HRV16 Infection Induces Changes in the Expression of Multiple piRNAs

Small RNAs can be categorized into three types based on their structures, origin, processing factors and the Argonaute proteins to which they bind, including microRNA (miRNA), small interfering RNA (siRNA), and Piwi-interacting RNA (piRNA) (Lipovich et al. 2010; Hombach and Kretz 2016). Piwi-interacting RNAs (piRNAs) are a class of small non-coding RNAs containing 25–31 nucleotides in length (Aravin et al. 2006; Girard et al. 2006; Sarkar et al. 2017). piRNAs are distinguished from microRNAs by being slightly longer, mainly expressed in the germline and binding to Piwi-class opposed to Agoclass Argonaute proteins (Aravin et al. 2006; Girard et al. 2006; Siomi et al. 2011). piRNAs are initially found to bind to Piwi family proteins to form Piwi/piRNA complexes. These complexes play important biological functions in the silencing of transposons, maintenance of genome stability and integrity in germ cells, regulation of germ cell development and differentiation, and formation of gametes (Rouget et al. 2010; Huang et al. 2013; Gou et al. 2015). Previous studies show that Piwi/piRNA complexes epigenetically and post-transcriptionally silence transposable elements in germline stem cells (Aravin and Bourc'his 2008; Aravin et al. 2008; Kuramochi-Miyagawa et al. 2008; Ku and Lin 2014). Furthermore, piRNAs are associated with the development and progression of several types of cancers (Cheng et al. 2011, 2012; Cui et al. 2011). A study has found that piRNA-823 in extracellular vesicles is associated with the production of the vascular endothelial growth factor (VEGF) and its proangiogenic activity in multiple myeloma cells (Yan et al. 2015; Li et al. 2019). Recent studies indicate that piRNAs play a vital role in physiological and pathological processes at the transcriptional or post-transcriptional level through piRNAs/piwi complex-mediated transcriptional gene silencing (TGS) and post-transcriptional gene silencing (PTGS) (Liu et al. 2019). The piRNAs/piwi complex directly binds with some proteins by piRNAs or the piwi protein PAZ domain to facilitate multi-protein interactions or to alter their subcellular localization (Liu et al. 2019).

Human rhinovirus (HRV) is one of the most important cold-causing pathogens in humans. HRV causes 50% or more upper respiratory tract infections and is the primary pathogen of acute lower respiratory tract infection in children. HRV is not only a key inducer of pediatric asthma (Stenberg-Hammar et al. 2017), but also associated with aggravation of acute exacerbations of chronic lung diseases, including asthma, bronchiectasis, and chronic obstructive pulmonary disease (COPD) (Gunawardana et al. 2014; Dijkema et al. 2016; Miller et al. 2016; Makris and Johnston 2018). However, the pathogenesis of HRV is unclear entirely. Whether piRNAs change during the HRV infection or whether piRNAs participates in the HRV pathogenesis remains unclear. Here we performed highthroughput sequencing on piRNAs of H1-HeLa cells infected with HRV16 at 12 h, 24 h and 36 h post-infection (p.i). Differential expression piRNAs were identified, their locations and functions were further predicted. Our data provide important insights into the interactions between piRNAs and HRV16 infection.

HRV16 strain (ATCCⓇ VR-283^TM) was propagated at MOI of 5 in H1-HeLa cells with 2% serum for maintenance and infection. H1-HeLa cells were purchased from American type culture collection (ATCC) and maintained in Dulbecco's Modified Eagle Medium (DMEM, Gibco) supplemented with 10% fetal calf serum (FCS) and antibiotics. Cells were harvested and analyzed at 12 h, 24 h, and 36 h p.i.

RNA was extracted from H1-HeLa cells infected with HRV16 at 12 h, 24 h, 36 h p.i and no-infected H1-HeLa cells through Trizol methods. For each sample, we subjected three micrograms (μg) of total RNA to PO treatment (Kirino and Mourelatos 2007). Total RNA (3 μg) and NaIO₄ (10 mmol/L) were mixed and kept at 0 ℃ for 40 min, then 5 μL of 1 mol/L rhamnose was added to quench the unreacted NaIO₄ at 0 C for 30 min. Fifty-five μL of 2 mol/L Lys-HCl (pH 8.5) was then added and the mixture was incubated at 45 C for 90 min for b-elimination. The treated RNAs were purified using a standard isopropanol precipitation method.

Small RNA (sRNA) library was constructed using the SOLiD Small RNA Expression Kit (Ambion). First, sRNAs were ligated to the 5' and 3' adaptors, and then reverse transcribed into cDNAs. A small amount of RNase was added to the cDNAs to remove any untranscribed RNAs. The cDNAs were amplified by PCR to obtain sufficient high-quality cDNAs for the construction of the sequencing library. After the library was constructed, the quality of the library was assessed using Bioanalyzer 2100, and the library was then sequenced using high-throughput sequencing. A substantial amount of sRNA sequences (16–31 nt) were obtained after removal of linkers and lowquality data. Most of these sRNAs were 22 nt in length, which is the typical length of sRNA generated from the Dicer pathway, and are hence well-conserved. After the removal of < 18 nt sequences, the final total purified sequences were used for piRNAs analysis.

Approximately 1.5 × 10⁴ cells were inoculated with HRV16 (MOI = 5) in a 96-well cell culture plate. Infected or uninfected cells of three repetitive wells were collected at 12 h, 24 h, and 36 h. Total cellular RNA of each well was isolated with Trizol Reagent. cDNAs of 21 piRNAs were gotten using stem-loop RT primers. The specific forward primer of each piRNA and a universal reverse primer were used to quantitatively analyze the expression of piRNAs on CFX96 (Bioread USA) using U6 small nuclear RNA (U6 snRNA) as a control. The relative amount of expressed piRNA was also calculated by the log₂ conversion of 2^-△△Ct.

The sequencing data of the four samples were analyzed in this study through TPM (Transcripts Per Kilobase Million). Comparable data were obtained by determining non-redundant sequences and copy number. piRNA prediction was performed on Piano using the SVM algorithm and transposon information (Wang et al. 2014). Genome distributions of piRNAs were produced by filtering reads with a length of 25–31 nts respectively and mapping the genome positions of the 5′ end of each read. Known piRNAs were verified by the Bowtie software with piRBase, and the secondary structure/folding of new piRNA sequences was predicted using UNAFold. The expression of piRNAs was normalized and compared using the t-test. The differences were considered statistically significant at P < 0.05. Differential gene expression analysis based on the negative binomial distribution was performed using DESeq2 as described (Love et al. 2014). We use adjusted P-value ≤ 0.1 and the absolute value of Log₂Ratio ≥ 1 as the default threshold to judge the significance of gene expression difference. Heatmap of piRNAs, which concurrently expressed among H1-HeLa cells at 12 h, 24 h and 36 h p.i, was generated by MultiExperiment Viewer.

To determine the overall differences in piRNA type and expression, we compared the common and unique piRNA sequences among H1-HeLa cells at 12 h, 24 h, and 36 h p.i. The results not only exhibited the typical size distribution of small ncRNAs (around 25–31 nts) (Fig. 1A), but also showed signatures of 5′1U and 10A biases as reported (Miesen et al. 2016a; Leggewie and Schnettler 2018), suggesting that both primary and secondary piRNAs produced when H1-HeLacells were infected with HRV16. Compared with normal control cells, there were 25, 523 of total 22, 151, 664, 36, 229 of total 24, 362, 486 and 25, 009 of total 22, 726, 546 upregulated piRNAs, and 29, 472 of total 22, 151, 664, 28, 906 of total 24, 362, 486 and 29, 541 of total 22, 726, 546 downregulated piRNAs in H1-HeLa cells at 12 h, 24 h and 36 h p.i, respectively (Fig. 1B, Supplementary Table S1). In particular, the infected H1-HeLa cells had 35, 264 upregulated and 25, 855 downregulated piRNAs at 24 h p.i compared with 12 h p.i; 25, 564upregulatedand26, 782downregulatedpiRNAsat36 h p.icomparedwith 12 h p.i; and 25, 268upregulated and 34, 172 downregulated piRNAs at 36 h p.i compared with 24 h p.i.

Figure 1.  Overall difference in piRNA expression among HRV16- infected H1-HeLa cells. A The overall nucleotide length distribution of RNA-sequencing results. B Upregulated and downregulated piRNAs in H1-HeLa cells at 12 h, 24 h and 36 h post-HRV16 infection. Black column indicates upregulated piRNAs. Grey column indicates downregulated piRNAs.

To screen for piRNAs induced from HRV16 infection, differentially expressed piRNAs at different time points post infection were analyzed. The expression of piRNAs was normalized by the filtering of known and highly-homologous sequences. Normalized piRNAs with a log2 ratio < 1 were removed. Compared with normal control cells, 120, 99 and 165 piRNAs were differentially expressed in H1-HeLa cells at 12 h, 24 h and 36 h p.i, respectively (Supplementary Table S2). In particular, 87 piRNAs were differentially expressed at 24 h p.i compared with 12 h p.i (Supplementary Table S3); 170 piRNAs were differentially expressed at 36 h p.i compared with 12 h p.i (Supplementary Table S3), and 174 piRNAs were differentially expressed at 36 h p.i compared with 24 h p.i (Supplementary Table S3). Our results showed that the expression of 21 piRNAs was concurrently changed among H1-HeLa cells at 12 h, 24 h and 36 h p.i (Fig. 2).

Figure 2.  Heatmap of 21 piRNAs concurrently expression among H1- HeLa cells at 12 h, 24 h and 36 h p.i by MultiExperiment Viewer.

Of these 21 piRNAs, seven piRNAs expression was significantly changed at the three infection time points, which have been reported in piRBase. The overall expression of hsa_piR_020391, hsa_piR_000805 and hsa_piR_000651 was downregulated relative to normal cells. Three of these piRNAs were higher than normal cells, which were respectively hsa_piR_017716, hsa_piR_004307 and hsa_piR_001311 (Fig. 3). Interestingly, the expression of has_piR_020326 was upregulated at 12 h p.i but downregulated from 1.49 to -2.38 and -1.18 compared with normal cells (Fig. 3).

Figure 3.  The RNA-sequencing result of seven known piRNAs expression levels in HRV16-infected H1-HeLa cells at the 12 h, 24 h and 36 h post infection.

Our results also found that 14 new piRNAs were differently expressed in infected-H1-HeLa cells relative to normal H1- HeLacellsat12 h, 24 hand36 hp.i.ThesenewpiRNAswere named as novel_pir97924, novel_pir78110, novel_pir78107, novel_pir78097, novel_pir78094, novel_pir76584, novel_pir73855, novel_pir70108, novel_pir70106, novel_pir 46604, novel_pir33182, novel_pir15900, novel_pir105705, and novel_pir105700, respectively. After the expression of these piRNA transcripts in infected H1-HeLa cells relative to normal cells was calculated, we found that the overall expression of these new piRNAs was downregulated, except novel_pir78107, novel_pir76584, novel_pir70106, which were slightly upregulated at 24 h p.i (Fig. 4).

Figure 4.  The RNA-sequencing result of 14 new piRNAs expression levels in HRV16-infected H1-HeLa cells at the 12 h, 24 h and 36 h postinfection.

To confirm the expression of above piRNAs identified by sequencing, we performed stem-loop RT-qPCR to examine the expression of those piRNAs in HRV16-infected H1-HeLa cells. Of these seven piRNAs, hsa_piR_001311, hsa_- piR_004307 and hsa_piR_017716, appeared trend of upregulation at 12 h, 24 h and 36 h p.i (Fig. 5). Expression of hsa_piR_004307 upregulated continuously over 1.5-fold as high as that of normal cells at 12 h, 24 h and 36 h p.i (Fig. 5). Expression of hsa_piR_001311 and hsa_piR_017716 reached over 1.5-fold of normal cells at 24 h and 36 h p.i. However, expression of hsa_piR_020391, hsa_piR_000651, hsa_- piR_000805, and hsa_piR_020326 were lower than normal cells and present a trend of downregulation at 12 h, 24 h and 36 h p.i (Fig. 5). The results showed that changes in piRNAs expression levels were consistent with RNA-seq data. Similarly, we also examined the expression of 14 new piRNAs using stem-loop RT-qPCR in HRV16-infected H1-HeLa cells. Of these novel piRNAs, similar to high-throughput results, the expression of all novel piRNAs was lower than normal cells at 12 h, 24 h and 36 h p.i (Fig. 6). The foldchange of expression level of new piRNAs that compared to normal cells were shown in Table 2.

Figure 5.  The expression of seven known piRNAs in HRV16- infected H1-HeLa cells at 12 h, 24 h and 36 h p.i. was validated by RT-qPCR. The infected and normal cells of three repetitive wells were collected at 12 h, 24 h, and 36 h post-infection. Total cellular RNA of each well was isolated with TRIzol Reagent. cDNAs of piRNAs were gotten using stem-loop RT primers. RT-qPCR was performed using specific forward primer of each piRNA and a universal reverse primer. U6 was used as an endogenous control for data normalization. The relative amount of expressed mRNA was calculated by the log2 conversion of 2^-△△Ct.

Figure 6.  The expression of new piRNAs in HRV16-infected H1-HeLa cells at 12 h, 24 h and 36 h p.i. was validated by RT-qPCR. The infected and normal cells of three repetitive wells were collected at 12 h, 24 h, and 36 h post-infection. Total cellular RNA of each well was isolated with TRIzol Reagent. cDNAs of piRNAs were gotten using stem-loop RT primers. RT-qPCR was performed using specific forward primer of each piRNA and a universal reverse primer. U6 was used as an endogenous control for data normalization. The relative amount of expressed mRNA was also calculated by the log2 conversion of 2^-△△Ct.

piRNA prediction was performed using Piano (Wang et al. 2014), and known piRNAs were verified by the Bowtie software. Secondary structure/folding of new piRNA sequences were predicted by UNAFold. The chromosome (chr) location of seven piRNAs was analyzed through piRBase. We found that hsa_piR_004307 and hsa_piR_020391 were located on chr 10 and 18, respectively. hsa_piR_001311 was located on chr 1 and 3, hsa_piR_020326 was located on chr 2, 8 and 9. Interesting, hsa_piR_020326recognized6sequences ofchr 9. hsa_piR_001311 recognized two sequences of chr 1. hsa_- piR_017716, hsa_piR_000805 and hsa_piR_000651 are located onchr 15. hsa_piR_000805was alsolocated onchr 13 (Table 1).

piRNA prediction found that novel_pir97924 is located on chr 8 and recognizes Long interspersed nuclear elements 1 (LINE-1). The novel_pir78110, novel_pir78107, novel_pir78097, novel_pir78094 and novel_pir76584, which are associated with the DNA/ hobo of Drosophila, Ac of maize and Tam3 of snapdragon (hAT)-Charlie transposon, are located on chromosome 3. The novel_pir73855, related to the long terminal repeat (LTR)/ Endogenous Retrovirus-K (ERVK) repetitive element, is located on chr 3. Both novel_pir70108 and novel_pir70106, which are associated with the LTR/ERVL-MaLR repetitive element, are located on chr 2. Novel_pir46604, which is related to the LTR/ERV1 repetitive element, are located on chr 18. The novel_pir33182 is located on chr 14 and is associated with the LTR/ERV1 repetitive element. The novel_pir15900 is located on chr 11 and is associated with the DNA/hAT-Tip100 repetitive element. Both novel_pir105705 and novel_pir105700 are located on the X chromosome and are associated with the LINE/L1 repetitive element (Table 2). Our results indicated that HRV16 infection-induced change of many new piRNA to recognize chromosomes.

SncRNAs are non-protein-coding RNAs that are shorter than 50 nt. These RNAs alone or in complexes with other proteins play certain biological functions. As a type of sncRNA, piRNAs along with miRNAs and siRNAs represent the three major mechanisms of small RNA interference in most animals (Malone and Hannon 2009). piRNAs are primarily found in intergenic regions that contain large numbers of transposons and repetitive sequences, which are otherwise known as piRNA clusters (Moyano and Stefani 2015). Compared with miRNAs and siRNAs, piRNAs are longer and are not derived from double-stranded hairpin RNAs, but from single-stranded RNAs transcribed from the piRNA gene cluster near the centromere. Though piRNAs are highly diverse in sequence and variable in length, characteristics of piRNA sequences preferentially harbor a Uridine(U) at their 50-most position (Brennecke et al. 2007; Stein et al. 2019). This 1U-bias has also been observed for other classes of small RNAs (Iwasaki et al. 2015). A report showed that there is a 60.9% bias of 1U in the sense reads for transposable elements (TEs), a 84.3% bias of 1U in antisense reads (George et al. 2015). However, in this study, we identified 28.2% bias of 1U from a total of 173, 084 novel piRNAs. Whether HRV infection disturbs on the formation of 1U bias needs further study.

piRNAs have diverse biological functions, including blockage of transcription and translation, silencing of transposons, regulation of development, and modulation of the transmission of genetic information. One of the major functions of piRNA is the regulation of retrotransposons during germ cell development. Methylation is one of the most effective means of retrotransposon regulation in normal cells. However, during the early stages of germ cell and embryonic development, epigenetic reprogramming provides versatility to the development of these cells through large-scale demethylation of the genome. This, in turn, activates the transcription of a large number of retrotransposons RNAs and enhances transposition activity. piRNA is the most effective means of transposon regulation at early stage, as it ensures that the integrity and stability of the genome are maintained when the genome is passed to the offspring (Ernst et al. 2017). Recent studies using C. elegans have shown that piRNAs can regulate almost all mRNAs in germ cells. The mode of binding between piRNA and its target gene is similar to miRNA, where a low degree of mismatch is permitted. Since the binding efficiency of piRNA is dependent on its binding affinity and not richness, piRNA can still exert its functions when it is at an extremely low concentration. Besides, it was found that transgenic alteration of the recognition and matching between piRNA and its target gene resulted in continuous, stable and high expression of the unregulated target gene, which demonstrates the broadness and effectiveness of piRNA-mediated regulation of gene expression in germ cells (Shen et al. 2018; Zhang et al. 2018). Furthermore, germ cell development can be severely impaired if the PIWI protein is mutated or the piRNA pathway is blocked (Watanabe et al. 2011; Gou et al. 2017). These studies demonstrate that piRNA regulates the transposition of transposons in multiple aspects and levels.

One of the major mechanisms by which piRNA inhibits the transposition of transposons is through histone modification at specific sites (Sienski et al. 2012; Le Thomas et al. 2013). A previous study showed that Piwi protein could directly interact with histone H1, which in turn increases H1 concentration at transposon gene sites and causes chromatin contraction (Watanabe and Lin 2014; Iwasaki et al. 2016). Moreover, piRNA can inhibit the transposition of LINE-1 DNA by inducing the latter's methylation (Aravin and Bourc'his 2008; Newkirk et al. 2017). In Mus musculus sperm cells, piRNA precursors transcribed from piRNA clusters bind to PIWI proteins to form piRNA-PIWI complexes. Upon entrance into the cell nucleus, recognition of the primary transcript of LINE-1 by these complexes leads to the recruitment of chromatin modifiers, which then results in the methylation of transposon-related histone H3K9 and demethylation of H3K4. The change in histone signaling subsequently recruits the DNA methyltransferase Dnmt3A/Dnmt3L complex, leading to the demethylation of LINE-1 DNA and, ultimately inhibiting LINE-1 transcription (Aravin et al. 2008). Our results showed that HRV16 infection induced the downregulation of several new piRNAs. We predicted that novel_pir97924, novel_pir105705 and novel_pir105700 could recognize LINE-1 since all of them were downregulated at the three infection time points. Therefore, we speculate that rhinovirus infection reduced these piRNAs to facilitate upregulation of LINE-1 transcription.

Our prediction also found that novel_pir78110, novel_pir78107, novel_pir78097, novel_pir78094 and novel_pir76584 are associated with the DNA/hATCharlie transposon. Novel_pir73855, novel_pir70108, novel_pir70106, novel_pir33182, novel_pir46604 and novel_pir15900 are related to the repetitive element LTR/ ERVK, LTR/ERVL-MaLR, LTR/ERV1 and DNA/hATTip100. Given that these piRNAs were downregulated, it suggests that rhinovirus infection may upregulate the expression of retrotransposons. And whether the upregulation of retrotransposons transcripts promotes pathogenesis of HRV infection or virus replication required further studies to confirm.

In addition to transposons, sufficient evidence supports the antiviral effect of piRNA pathway in host strategies for protection against invading viruses. It is worth noting that the genome analysis of Aag2 mosquito cells shows that a valuable 'genome record' of host virus interaction is closely related to the production of piRNAs, which leads to a genetic immune system based on piRNA (Whitfield et al. 2017). The system is believed to be a defense mechanism against viruses that have infected Anopheles aegypti in the past and is proposed to provide long-term virus specific protection for future infections. The piRNA pathway has been shown to be a contributor to antiviral immunity in mosquitoes or derived cells (Morazzani et al. 2012; Schnettler et al. 2013). The piRNA has also been found to origin from a number of mosquito-borne viruses, including Dengue virus (DENV), Chikungunya virus, Semliki Forest virus (SFV) and Rift Valley fever virus (Morazzani et al. 2012; Schnettler et al. 2013). Further, Piwi4 has been demonstrated to inhibit replication of SFV (Schnettler et al. 2013), and both Ago3 and Piwi5 benefit the generation of specific piRNAs to prohibit the proliferation of DENV and SINV in Aag2 cells (Miesen et al. 2015, 2016a, b ); Recent findings showed that the human piwil 2 (hili) bound to tRNAs to inhibit HIV replication in activated CD4⁺ T cells (Peterlin et al. 2017). Our previous study found that Coxsackievirus B3 infection-induced changes in many piRNAs expression in HeLa cells (Yao et al. 2020). However, we don't know whether piRNAs effect on gene expression of H1-HeLa cells to benefit viral replication.

Although previous research suggests that piRNAs are germline specific, the high-throughput small RNA sequencing data show widespread presence of piRNA in multiple somatic tissues in more than 130 fruit fly, mouse and rhesus macaque, and piRNA abundance of mouse pancreas and macaque epididymis was compatible with piRNA abundance in the germline (Yan et al. 2011). These findings indicate that piRNA-like molecules might play important roles outside of the germline. The role of somatic piRNAs in cancer biology is increasingly being recognized. Growing evidence shows that piRNAs strongly correlate with tumor cell malignant phenotype and clinical stage through maintaining genomic stability (Yan et al. 2011). Not only that, some evidence showed that expression of many piRNAs are wildely found in somatic tissues of human, flies, mice and macaques (Yan et al. 2011), but also participate in cancer development (Martinez et al. 2015). Though both level of expression and the number of piRNAs are lesser in somatic cells than in germline, piRNAs contribution to tumorigenesis is significant. Research results from Martinez et al. showed that 273 known piRNAs are found in somatic tissues (n = 508), and 522 cancer-type specific piRNAs appeared in tumor samples (n = 5, 752) from 12 tumor types (Gou et al. 2017). The functional role of piRNAs in nongonadal cells is a rapidly emerging field of research. Numerous piRNAs are dysregulated in tumor tissues, playing tumor-promoting or tumor-suppressor roles (Liu et al. 2019). 369 down-regulated piRNAs and 235 up-regulated piRNAs appear in 106 renal clear cell carcinoma tissues through piRNA microarray's study (Busch et al. 2015). Some invertebrates affect antiviral machinery, which piRNA pathway processes viral RNA into piRNAs (Miesen et al. 2016b). However, some virus infections, for example, CVB3, can also induce abnormal expression of piRNA in host cells (Yao et al. 2020). Our results showed that many dysregulated piRNAs appeared in HRV16-infected H1-HeLa cells again.

Most piRNA pathway proteins localize to specific cytoplasmic compartments, including nuage in animal germ cells, Yb bodies in the somatic ovarian follicle cells of flies, and the mitochondrial outer membrane in all phased piRNA-producing cells (Trenkmann 2019). The enrichment of the piRNA machinery in these subcellular structures may serve to increase the local concentration of specific proteins or protect piRNA precursors from housekeeping nucleases. Virus infection may prevent compartmentalization from entering the piRNA pathway. It is possible to explain that many piRNAs decreased in the H1-HeLa cells infected with HRV16. In the meantime, virus infection-induced increase of some piRNA through different piRNA biosynthesis mechanisms. Production of some piRNA precursor transcripts may be shunted to other sites of piRNA production (Zhang et al. 2012).

In this study, we detected and analyzed the expression of piRNAs in HRV16 infected H1-HeLa cells, and found HRV16 infection altered the expression of piRNAs compared to normal cells. Further analysis showed that rhinovirus infection could reduce some piRNAs, which may facilitate upregulation of LINE-1 transcription or upregulate retrotransposons' expression. Though a variety of human cells, including human embryonic kidney, diploid fibroblast lines from human embryonic lung, and the HEP- 2, KB, L132, and H1-HeLa cell lines, are available for growth of rhinovirus (Fields et al. 2006), only H1-HeLa cells were used to evaluate the expression of piRNAs in this study. If the expression of piRNAs from HRV16 infected different cells are compared, a comprehensive map of affected piRNA will be drawn.

This work was supported by the China MegaProject for Infectious Disease (2018ZX10102001, 2018ZX10711001, 2018ZX10734401 and 2018ZX10734404), and the SKLID Development Grant (2011SKLID104).

JH designed the experiments. JL, XW, YW, JS, QS carried out the experiments. JL, YW and JH analyzed and interpreted the data. JL, XW and JH wrote the manuscript. All authors read and approved the final manuscript.

The authors declare that they have no conflict of interest.

This article does not contain any studies with human or animal subjects performed by any of the authors.