key: cord-0867093-k31xtghv
authors: Wu, Wenzhe; Choi, Eun-Jin; Wang, Binbin; Zhang, Ke; Adam, Awadalkareem; Huang, Gengming; Tunkle, Leo; Huang, Philip; Goru, Rohit; Imirowicz, Isabella; Henry, Leanne; Lee, Inhan; Dong, Jianli; Wang, Tian; Bao, Xiaoyong
title: Changes of Small Non-coding RNAs by Severe Acute Respiratory Syndrome Coronavirus 2 Infection
date: 2022-02-23
journal: Front Mol Biosci
DOI: 10.3389/fmolb.2022.821137
sha: 9924c1b28b8d687c53eb6c54d68ebffbe6907134
doc_id: 867093
cord_uid: k31xtghv

The ongoing pandemic of coronavirus disease 2019 (COVID-19), which results from the rapid spread of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is a significant global public health threat, with molecular mechanisms underlying its pathogenesis largely unknown. In the context of viral infections, small non-coding RNAs (sncRNAs) are known to play important roles in regulating the host responses, viral replication, and host-virus interaction. Compared with other subfamilies of sncRNAs, including microRNAs (miRNAs) and Piwi-interacting RNAs (piRNAs), tRNA-derived RNA fragments (tRFs) are relatively new and emerge as a significant regulator of host-virus interactions. Using T4 PNK‐RNA‐seq, a modified next-generation sequencing (NGS), we found that sncRNA profiles in human nasopharyngeal swabs (NPS) samples are significantly impacted by SARS-CoV-2. Among impacted sncRNAs, tRFs are the most significantly affected and most of them are derived from the 5′-end of tRNAs (tRF5). Such a change was also observed in SARS-CoV-2-infected airway epithelial cells. In addition to host-derived ncRNAs, we also identified several small virus-derived ncRNAs (svRNAs), among which a svRNA derived from CoV2 genomic site 346 to 382 (sv-CoV2-346) has the highest expression. The induction of both tRFs and sv-CoV2-346 has not been reported previously, as the lack of the 3′-OH ends of these sncRNAs prevents them to be detected by routine NGS. In summary, our studies demonstrated the involvement of tRFs in COVID-19 and revealed new CoV2 svRNAs.

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 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-2positive 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.

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.

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 , followed by qRT-PCR.

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 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.

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 4 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) .

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 postinfection.

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 2 , 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′ CCTACTAAATTAAATGATCTCTGCTTTACT; reverse primer, 5′' CAAGCTATAACGCAGCCTGTA.

To evaluate sncRNAs expression, qRT-PCR was performed, as described previously 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 ) Wu et al., 2021) . The primers and 3′-RNA linker sequences are listed in Table 1 .

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 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 32 P-labeled DNA probe in ULTRAhyb-Oligo solution (Life Technologies, Grand Island, NY, United States), followed by 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.

Frontiers in Molecular Biosciences | www.frontiersin.org February 2022 | Volume 9 | Article 821137 5 membrane washing and image development. The 32 P-labeled DNA probe for tRF5-GluCTC was 5′-CGCCGAATCCTAACC ACTAGACCACCA-3′.

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 twotailed 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.

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) .

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.

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, topranked 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 .

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 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.

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 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. Figure 4F) , which was the most abundant tRF5 among the tested four tRF5s, confirming the liability of our newly developed qRT-PCR. 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 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.

Viral-derived sncRNAs are also an important family of sncRNAs (Tycowski et al., 2015) . To investigate whether SARS-CoV-2encoded 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-2encoded svRNAs are listed in Table 6 . 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.

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 ). 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 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.

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-ligationbased 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 , 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 virusspecific. Previously, we and others have shown that RSV, hepatitis B virus (HBV), and hepatitis C virus (HCV) infections lead to (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 miRNAslike 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 . 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.

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.

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.

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. 

Non-coding RNAs in Development and Disease: Background, Mechanisms, and Therapeutic Approaches

Cells and Culture Systems Used to Model the Small Airway Epithelium

Aberrant Methylation of T RNA S Links Cellular Stress to Neurodevelopmental Disorders

Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19

Syncytia Formation by SARS-CoV-2-Infected Cells

Fast-track Development of an In Vitro 3D Lung/immune Cell Model to

Verification of SARS-CoV-2-Encoded Small RNAs and Contribution to Infection-Associated Lung Inflammation

ELAC2, an Enzyme for tRNA Maturation, Plays a Role in the Cleavage of a Mature tRNA to Produce a tRNA-Derived RNA Fragment during Respiratory Syncytial Virus Infection

Small Noncoding RNAs: Biogenesis, Function, and Emerging Significance in Toxicology

Systematic Classification of Non-coding RNAs by Epigenomic Similarity

The Mechanics of miRNA-Mediated Gene Silencing: a Look under the Hood of miRISC

Stress Induces tRNA Cleavage by Angiogenin in Mammalian Cells

A Virus-Derived microRNA-like Small RNA Serves as a Serum Biomarker to Prioritize the COVID-19 Patients at High Risk of Developing Severe Disease

Phospho-RNA-seq: a Modified Small RNA-Seq Method that Reveals Circulating mRNA and lncRNA Fragments as Potential Biomarkers in Human Plasma

Virus-encoded microRNAs

Identification of microRNAs and Other Small Regulatory RNAs Using cDNA Library Sequencing

A Computational Approach for Predicting Role of Human MicroRNAs in MERS-CoV Genome

Human Lung Epithelial Cell Cultures for Analysis of Inhaled Toxicants: Lessons Learned and Future Directions

Vienna RNA Secondary Structure Server

Sex Hormone-dependent tRNA Halves Enhance Cell Proliferation in Breast and Prostate Cancers

Potential Prognostic Biomarkers of Lung Adenocarcinoma Based on Bioinformatic Analysis

West Nile Virus Encodes a microRNA-like Small RNA in the 3' Untranslated Region Which Up-Regulates GATA4 mRNA and Facilitates Virus Replication in Mosquito Cells

G-quadruplex Structures Contribute to the Neuroprotective Effects of Angiogenin-Induced tRNA Fragments

The UCSC Table Browser Data Retrieval Tool

Epigenetic Regulator miRNA Pattern Differences Among SARS-CoV, SARS-CoV-2, and SARS-CoV-2 World-wide Isolates Delineated the Mystery behind the Epic Pathogenicity and Distinct Clinical Characteristics of Pandemic COVID-19

A Novel Class of Small RNAs: tRNA-Derived RNA Fragments (tRFs)

H5N1 Influenza Virus-specific miRNA-like Small RNA Increases Cytokine Production and Mouse Mortality via Targeting poly(rC)-Binding Protein 2

A tRNA-Derived RNA Fragment Plays an Important Role in the Mechanism of Arsenite -induced Cellular Responses

MicroRNome Analysis Unravels the Molecular Basis of SARS Infection in Bronchoalveolar Stem Cells

Frontiers in Molecular Biosciences | www.frontiersin.org

Human Plasma and Serum Extracellular Small RNA Reference Profiles and Their Clinical Utility

SARS-CoV-Encoded Small RNAs Contribute to Infection-Associated Lung Pathology

Exogenous ACE2 Expression Allows Refractory Cell Lines to Support Severe Acute Respiratory Syndrome Coronavirus Replication

Influenza A Virus-Generated Small RNAs Regulate the Switch from Transcription to Replication

Methods and Matrices: Approaches to Identifying miRNAs for Nasopharyngeal Carcinoma

Automated 3D Structure Composition for Large RNAs

Pathological Significance of tRNA-Derived Small RNAs in Neurological Disorders

Regulation of Host Innate Immunity by Non-coding RNAs during Dengue Virus Infection

tRNA Wobble Modifications and Protein Homeostasis. Translation

Small Noncoding RNA and Cancer

The RNA-Induced Silencing Complex Is a Mg2+-dependent Endonuclease

Influenza Virus Infection Increases ACE2 Expression and Shedding in Human Small Airway Epithelial Cells

Small tRNA-Derived RNAs Are Increased and More Abundant Than microRNAs in Chronic Hepatitis B and C

Transfer RNA-Derived Fragments and tRNA Halves: Biogenesis, Biological Functions and Their Roles in Diseases

Viral Noncoding RNAs: More Surprises

Analysis of microRNA Transcriptome by Deep Sequencing of Small RNA Libraries of Peripheral Blood

Identification and Functional Characterization of tRNA-Derived RNA Fragments (tRFs) in Respiratory Syncytial Virus Infection

A Cytoplasmic RNA Virus Generates Functional Viral Small RNAs and Regulates Viral IRES Activity in Mammalian Cells

Broad Anti-coronavirus Activity of Food and Drug Administration-Approved Drugs against SARS-CoV-2 In Vitro and SARS-CoV In Vivo

World Health Organization

Non-Coding RNAs and Their Role in Respiratory Syncytial Virus (RSV) and Human Metapneumovirus

tRNA-Derived Fragments in Alzheimer's Disease: Implications for New Disease Biomarkers and Neuropathological Mechanisms

Expression of the SARS-CoV-2 ACE2 Receptor in the Human Airway Epithelium

Identification of Two Novel Functional tRNA-Derived Fragments Induced in Response to Respiratory Syncytial Virus Infection

tRNA-derived Small Non-coding RNAs in Human Disease

A Novel Coronavirus from Patients with Pneumonia in China

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).

We thank Darby L. Buck for the manuscript editing.

The 

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.