-
Human rhinoviruses (HRVs) are respiratory viruses that are commonly detected in patients with otitis media, croup, bronchiolitis, or pneumonia and are associated with exacerbation of lung diseases[1]. Though most infected patients are asymptomatic or only suffer from mild symptoms, HRVs have recently been linked to exacerbation of chronic obstructive pulmonary disease (COPD), asthma development, and severe bronchiolitis in infants and children and to fatal pneumonia in elderly and immunocompromised adults[2]. Several studies have found that wheezing episodes or lower respiratory tract infections are more common in patients with rhinovirus C than in patients with rhinovirus A or rhinovirus B[3]. HRV infection can cause a considerable economic burden in terms of medical visits and school and work absenteeism[2,4].
The miRNAs are a class of highly conserved small non-coding RNAs. The miRNAs are key post-transcriptional regulators of gene expression[5]. The miRNAs regulate the expression levels of about 60% of protein-coding genes in the human genome[6]. The miRNA-mediated control of infection with both DNA and RNA viruses has been described in a wide variety of hosts. Several types of miRNA-target interactions have been described to play crucial roles in the regulation of viral replication, maintenance of latency and/or reactivation, immune evasion, pathogenicity, and cell transformation.
It is therefore not surprising that HRV infection may regulate the expression of specific miRNAs for efficient replication as an ever-evolving pathogenicity and survival strategy. To investigate the miRNA expression levels in HRV-infected cells and the reaction network between miRNAs and their targets, miRNA expression profiles in HRV16-infected H1-HeLa cells were constructed through small RNA (sRNA) next-generation sequencing (NGS). Biological implications of 37 differentially expressed miRNAs were analyzed through GO and KEGG pathway analyses. Subsequently, interaction networks between 37 differentially expressed miRNAs and their targets were established by mirDIP and Navigator. The prediction results showed that QKI and NFAT5 interacted with 12 and 11 miRNAs, respectively. BNC2, CELF2, LCOR, MBNL2, MTMR3, NFIB, PPARGC1, RSBN1, TRPS1, WDR26, and ZNF148 all inter acted with nine miRNAs. Our results revealed how HRV affected miRNA expression during infection.
To investigate the changes in miRNA levels in HRV16-infected cells, HRV16 was propagated in H1-HeLa cells at a multiplicity of infection of 1. Infected cells were maintained in medium containing 2% FBS following 1 h of adsorption at 37 °C as in a previous study. After H1-HeLa cells were infected with HRV16 for 12, 24, or 36 h, total RNA was isolated from HRV16-infected and uninfected H1-HeLa cells at 12, 24, and 36 h post-infection (pi) using TRIzol Reagent (Thermo Fisher Scientific, USA). Using sRNA deep sequencing, a total of 22,151,664, 24,362,486, and 22,726,546 clean reads were obtained at the three time points for HRV16-infected H1-HeLa cells, respectively, and 22,090,435 clean reads were obtained for the control cells. All clean reads were analyzed based on sRNA databases such as miRBase, Rfam, siRNA, piRNA, snoRNA, etc., to determine their locus of origin, resulting in 2.2 × 107, 2.4 × 107, and 2.3 × 107 miRNA sequence tags in HRV16-infected cells’ libraries, respectively, and 2.2 × 107 miRNA sequence tags in the control cells’ library.
Subsequently, differentially expressed sRNAs between HRV16-infected and uninfected H1-HeLa cells were analyzed at each time point by the Exp Diff method. The results showed that 2,795, 2,797, and 2,779 miRNAs were differentially expressed at 12, 24, and 36 h, respectively. Using an appropriate filter, 327, 287, and 487 differentially expressed miRNAs were found at 12, 24, and 36 h, respectively. Among these, 41, 41, and 90 miRNAs were novel, respectively.
With Data Encryption Standard (DES), hierarchical clustering was performed. All differentially regulated miRNAs between HRV16-infected and uninfected H1-HeLa cells at 12, 24, and 36 h pi are shown in Figure 1. In total, 53 miRNAs were differentially expressed at all three timepoints. Of these, 37 met the fold change ≥ 1.5 criterion, with 17 miRNAs being upregulated and 20 miRNAs being downregulated. The expression of both hsa-miR-18b-5p and hsa-miR-100-5p was increased ≥ 4.9-fold compared to control cells (Table 1). The expression of hsa-miR-18b-5p, hsa-miR-100-5p, hsa-miR-181b-3p, hsa-miR-93-3p, and hsa-miR-3614-3p gradually decreased during infection. Particularly, the expression of hsa-miR-3614-3p in HRV16-infected cells at 36 h pi was 2.58-fold lower than that in uninfected cells (Table 1). In total, 16 miRNAs were found for the first time. Among these, 15 miRNAs were down regulated, and novel_mir334 showed a slight increase (Table 2). In conclusion, the expression of miRNAs in HRV16-infected H1-HeLa cells was comprehensively analyzed.
Figure 1. Hierarchical clustering. The x-axis represents samples. The y-axis represents DESs. Coloring indicates fold change values (upregulated: red; downregulated: blue).
No. miRNA 12 h 24 h 36 h 1 hsa-miR-18b-5p 9.51 8.64 8.49 2 hsa-miR-100-5p 7.4 6.9 4.91 3 hsa-miR-181b-3p 4.12 3.4 3.06 4 hsa-miR-4659a-3p 3.81 4.23 4.15 5 hsa-miR-5001-5p 3.43 4.89 3.32 6 hsa-miR-93-3p 3.35 2.96 1.66 7 hsa-miR-99b-3p 2.93 3.34 3.74 8 hsa-miR-196a-3p 2.83 2.04 2.49 9 hsa-miR-99a-5p 2.42 2.59 1.91 10 hsa-miR-423-3p 2.39 2.39 1.82 11 hsa-miR-3614-3p 2.37 1.55 −2.58 12 hsa-miR-1287-5p 2.28 2.52 2.79 13 hsa-miR-708-3p 2.16 2.33 2.63 14 hsa-miR-365a-5p 2.16 1.79 2.79 15 hsa-miR-767-3p 1.98 2.11 1.61 16 hsa-miR-548ab 1.58 3.65 3.37 17 hsa-miR-422a 1.49 3.72 4.87 18 hsa-miR-195-5p −1.51 −1.69 −3.64 19 hsa-miR-92a-1-5p −1.53 −4.75 −4.62 20 hsa-miR-708-5p −1.66 −2.75 −3.91 21 hsa-miR-324-5p −1.66 −1.71 −1.54 22 hsa-miR-200b-3p −1.74 −1.86 −2.49 23 hsa-miR-140-5p −1.88 −2.17 −2.06 24 hsa-miR-135a-3p −2.01 −1.56 −1.73 25 hsa-miR-374b-3p −2.3 −3.48 −2.4 26 hsa-miR-516a-5p −2.51 −6.07 −6.07 27 hsa-miR-3184-5p −2.59 −3.65 −2.01 28 hsa-miR-30b-3p −2.74 −2.08 −1.93 29 hsa-miR-151b −2.88 −4.83 −4.75 30 hsa-miR-3065-3p −3.03 −3.61 −3.73 31 hsa-miR-1283 −3.16 −2.67 −3.45 32 hsa-miR-23a-5p −3.53 −2.36 −2.18 33 hsa-miR-624-5p −3.56 −1.7 −1.71 34 hsa-miR-130a-5p −3.88 −5.58 −4.42 35 hsa-miR-522-3p −3.95 −3.28 −2.92 36 hsa-miR-196a-5p −4.46 −4.61 −2.42 37 hsa-miR-500b-5p −4.52 −4.52 −3.99 Table 1. Differentially expressed known miRNAs at 12, 24, and 36 h after HRV16 infection as determined by sRNA sequencing
No. miRNA 12 h 24 h 36 h 1 novel_mir803 −4.16 −3.74 −3.99 2 novel_mir801 −4.69 −2.38 −4.32 3 novel_mir715 −5.39 −2.87 −4.34 4 novel_mir69 −3.65 −5.15 −3.96 5 novel_mir65 −4.80 −2.51 −5.97 6 novel_mir64 −4.68 −2.74 −3.07 7 novel_mir61 −3.44 −3.25 −2.42 8 novel_mir59 −5.62 −5.25 −3.85 9 novel_mir56 −3.98 −3.44 −3.41 10 novel_mir471 −2.62 −3.16 −1.73 11 novel_mir468 −6.33 −3.46 −4.74 12 novel_mir334 1.47 1.21 1.34 13 novel_mir19 −1.50 −1.39 −2.03 14 novel_mir143 −1.20 −1.49 −1.25 15 novel_mir1076 −5.14 −3.50 −5.03 16 novel_mir10 −1.73 −1.14 −1.45 Table 2. Differentially expressed novel miRNAs at 12, 24, and 36 h after HRV16 infection as determined by sRNA sequencing
Since the 37 known miRNAs and 16 novel miRNAs have a similar change trend, the results indicate that these miRNAs simultaneously participate in biological regulation in virus-infected cells.
To verify the sRNA sequencing results, RNA was isolated from HRV16-infected cells and 37 miRNAs were validated by RT-qPCR using miRNA-specific TaqMan assays (CT Bioscience, China). Actin mRNA was used as an endogenous control for data normalization. The relative expression of mRNA was calculated by the 2−ΔΔCt method, and the experiment was performed three times. RT-qPCR results were consistent with the sequencing results, though fold change values of RT-qPCR were lower those of our sRNA sequencing analysis (Table 3). The high correlation between our sequencing and RT-qPCR results confirmed that HRV16 infection caused wide spread changes in miRNA expression levels.
No. miRNA 12 h 24 h 36 h 1 hsa-miR-18b-5p 1.85 1.82 1.55 2 hsa-miR-100-5p 1.53 2.46 1.87 3 hsa-miR-181b-3p 1.49 1.95 1.42 4 hsa-miR-4659a-3p 1.57 1.50 1.40 5 hsa-miR-5001-5p 1.34 1.66 1.55 6 hsa-miR-93-3p 1.68 1.79 1.50 7 hsa-miR-99b-3p 1.44 1.83 2.35 8 hsa-miR-196a-3p 1.52 2.19 1.80 9 hsa-miR-99a-5p 1.84 2.60 1.69 10 hsa-miR-423-3p 1.53 1.78 1.46 11 hsa-miR-3614-3p 1.49 1.50 1.67 12 hsa-miR-1287-5p 1.47 2.44 2.05 13 hsa-miR-708-3p 1.83 2.10 1.40 14 hsa-miR-365a-5p 1.57 1.49 1.52 15 hsa-miR-767-3p 1.75 1.54 1.53 16 hsa-miR-548ab 1.19 1.58 1.71 17 hsa-miR-422a 1.51 1.57 1.53 18 hsa-miR-195-5p −1.36 −1.25 −1.99 19 hsa-miR-92a-1-5p −1.38 −1.25 −2.12 20 hsa-miR-708-5p −1.48 −1.33 −2.28 21 hsa-miR-324-5p −1.49 −1.36 −2.38 22 hsa-miR-200b-3p −1.50 −1.36 −2.60 23 hsa-miR-140-5p −1.59 −1.39 −1.25 24 hsa-miR-135a-3p −1.61 −1.46 −1.26 25 hsa-miR-374b-3p −1.67 −1.57 −1.27 26 hsa-miR-516a-5p −1.69 −1.58 −1.28 27 hsa-miR-3184-5p −1.73 −1.61 −1.28 28 hsa-miR-30b-3p −1.76 −1.79 −1.37 29 hsa-miR-151b −1.79 −2.32 −1.40 30 hsa-miR-3065-3p −1.82 −1.47 −1.45 31 hsa-miR-1283 −2.74 −1.51 −1.52 32 hsa-miR-23a-5p −1.76 −1.50 −1.58 33 hsa-miR-624-5p −2.36 −1.35 −1.62 34 hsa-miR-130a-5p −1.95 −1.57 −1.62 35 hsa-miR-522-3p −2.08 −1.99 −1.72 36 hsa-miR-196a-5p −2.68 −1.18 −1.77 37 hsa-miR-500b-5p −2.89 −2.32 −1.84 Table 3. Differentially expressed miRNAs as confirmed by RT-qPCR at 12, 24, and 36 h after HRV16 infection
To further evaluate the biological implications of different miRNAs, Gene Ontology (GO) enrichment analysis of miRNA target genes was conducted. GO terms were classified into the biological process, cellular component, and molecular function categories. The top five enriched GO terms in the biological process category were cellular process, single-organism process, metabolic process, biological regulation, and regulation of biological process (Supplementary Figure S1, available in www.besjournal.com). GO analysis revealed that many upregulated miRNAs were involved in the regulation of immune processes, such as inhibition of antigen processing, presentation of exogenous peptide antigen via MHC class I/II, and the cytokine-mediated signaling pathway. Moreover, many miRNAs were found to be involved in the innate immune system, such as the Toll-like receptor (TLR) signaling pathway, the TRIF-dependent TLR signaling pathway, the MyD88-dependent/independent TLR signaling pathway, and type I interferon production.
Figure S1. GO functional classification of target genes of miRNAs induced by HRV16. The x-axis shows the means number of DEGs (the number is presented by its square root value). The y-axis shows GO terms. All GO terms are grouped into three categories: blue, biological process; green, cellular component; orange, molecular function.
After GO analysis, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of miRNA target genes. At the three infection time points miRNA target genes showed similar trends in metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases (Supplementary Figure S2, available in www.besjournal.com). The top three enriched KEGG pathways were pathways in cancer, the PI3K-Akt signaling pathway, and the MAPK signaling pathway. These results confirmed that the targets of miRNAs were involved in a wide variety of physiological processes.
Figure S2. KEGG classification of target genes of miRNAs induced by HRV16. The x-axis shows the number of DEGs. The y-axis shows the secondary KEGG pathway terms. All secondary pathway terms are grouped in top pathway terms indicated in different colors.
Next, 37 miRNAs whose expression was induced by HRV16 infection were inputted into the microRNA Data Integration Portal (miRDIP;
http://ophid.utoronto.ca/mirDIP ). Interactions between miRNAs and their predicted mRNA targets were then visualized as networks using NAViGaTOR v2.14 (http://ophid.utoronto.ca/navigator ). Interaction networks between 37 known differentially expressed miRNAs and significant gene targets identified by miRDIP were generated. Only the most highly connected miRNAs were used to build the networks.The target genes interacting with less than seven miRNAs were excluded, and the target genes regulated by at least seven miRNAs were retained. In concentric ellipses, the targets were regulated by 7, 8, 9, 11, and 12 miRNAs (Figure 2). The prediction results showed that the QKI gene was regulated by 12 miRNAs (Figure 2, purple circle) and NFAT5 was regulated by 11 miRNAs (Figure 2, blue circles). In total, 11 targets were recognized by nine miRNAs (Figure 2), 10 targets were regulated by 8 miRNAs (Figure 2), and 32 targets were regulated by 7 miRNAs. These genes played important roles in the regulation of cellular development, proliferation, differentiation, and cellular growth. Our results further demonstrated that HRV16 infection induced changes in the levels of cellular miRNAs to control multiple cellular processes, including cellular proliferation, differentiation, and development, as previously reported[7]. However, hsa-4481-5p, hsa-500b-3p, hsa-5088-3p, hsa-6825-5p, andhsa-6882-5p did not get high scores in the evaluation of gene recognition by mirDIP.
Figure 2. The miRNA-target gene regulatory interaction network induced by HRV16. Interaction networks between 37 known differentially expressed miRNAs and targets were established by mirDIP and Navigator. After exclusion of targets interacting with one to six miRNAs, targets recognized by at least seven miRNAs were used to construct miRNA–target gene networks. In concentric ellipses, the targets were regulated by 7, 8, 9, 10, and 11 miRNAs (red circles, orange circles, light green circles, green circles, blue circles, and purple circles, respectively). The miRNAs are represented by gray triangles on the periphery. The number of genes recognized by both hsa4659a-3p and hsa516a-5p and other miRNAs was less than six.
Subsequently, targets of maximum six upregulated miRNAs were predicted. TARDBP was recognized by five miRNAs, including hsa-miR-18b-5p, hsa-miR-100-5p, hsa-miR-99a-5p, hsa-miR-93-3p, and hsa-miR-3614-3p. BNC2, MBNL2, FKBP5, TRIB2, MTMR3, and MBNL1 were recognized by four miRNAs. In total, 35 targets were recognized by 3 miRNAs, and 236 targets were recognized by 2 miRNAs. We also found that hsa-miR-18b-5p and hsa-miR-767-3p recognized 174 and 105 target genes, respectively. The hsa-miR-422a, hsa-miR-100-5p, and hsa-miR-99a-5p recognized 92, 72, and 68 genes, respectively. The hsa-miR-1287-5p, hsa-miR-93-3p, hsa-miR-548ab, and hsa-miR-196a-3p recognized 23, 21, 16, and 13 target genes, respectively. The hsa-miR-423-3p, hsa-miR-3614-3p, hsa-miR-18b-6p, and hsa-miR-767-4p recognized eight, seven, four, and four targets, respectively. The hsa-miR-423a and hsa-miR-100-6p recognized two targets. The hsa-miR-99a-6p had only one target. The hsa-miR-1287-6p did not recognize any target.
In conclusion, a specific miRNA expression profile upon HRV16 infection was constructed by sRNA sequencing got three time points of infection. We found that HRV16 infection induced expression of 37 known miRNAs and 16 novel miRNAs in H1-HeLa cells. An interaction network between miRNAs and effector molecules during HRV16 infection has been constructed. Our results suggested tight coordination between transcriptional or post-transcriptional regulation of an miRNA and its target genes in HRV16-infected H1-HeLa cells. The miRNAs may serve as biomarkers to treat asthma and prevent asthma and COPD exacerbation. Our present results represented the first step in elucidating the roles of miRNAs during HRV infection, which will be our next research focus.
Conflict of Interest The authors declare that they have no competing interests.
Acknowledgments HAN Jun designed and conceived the experiments; SONG Qin Qin, WANG Yan Hai, and WANG Xin Ling performed the experiments, analyzed data, and wrote the original draft; SHI Bing Tian, WANG Rui Fang, SONG Juan, WANG Wen Jun, XIA Dong, XIA Zhi Qiang, WEI Qiang, and HAN Jun supervised the project. HAN Jun and WEI Qiang revised the manuscript.
HTML
22038+Supplementary Materials.pdf
|
|
Quick Links
DownLoad:
