Differential miRNA expression profiling of highly pathogenic avian influenza virus H5N1 infected chicken lungs reveals critical microRNAs, biological pathways and genes involved in the molecular pathogenesis.

Dear Editor,

Avian influenza is a highly contagious viral infection affecting the respiratory system. MicroRNAs (miRNAs) are small, regulatory, endogenous, non-coding RNAs of ∼22 nt that regulate the gene expression of the target mRNAs by cleavage or translational repression. miRNAs are connected with the host response during avian influenza virus (AIV) infection (Wang et al., 2009, 2012). In this study, we used miRNomics approach to understand the complex host-pathogen interaction during the HPAIV H5N1 infection in chicken.

Briefly, six specific pathogen free (SPF) chickens, aged four weeks, were divided into two groups. One group was intranasally inoculated with 1 ​× ​106 EID50 of H5N1 virus (A/chicken/Navapur/7972/2006) and the other group was mock-infected with phosphate-buffered saline. Lung samples were collected post 18 h of the challenge and subjected to sequencing on ABI SOLiD platform. The detailed methodology is described in Supplementary file. On sequencing, the small RNA libraries on SOLiD generated 57,210,641 and 52,848,558 raw reads for mock-infected and AIV infected samples, respectively. The experiment identified 297 and 201 mature miRNAs in the mock-infected and the AIV-infected samples, respectively. This included 200 mature miRNAs to be common in both groups. In addition to this, 97 miRNAs were observed in mock-infected group and one miRNA (gga-miR-1670) was observed in the AIV-infected group (Supplementary Table S1). Out of 126 differentially expressed miRNAs observed between the AIV infected and mock groups, 36 miRNAs were found to be upregulated (UR) with a positive log fold change and 90 were downregulated (DR) with negative log fold change in the AIV infected group (Supplementary Table S2). The distribution of differentially expressed mature miRNAs were mapped onto individual chromosomes (Supplementary Fig. S1).

The top ten upregulated (gga-miR-1793, gga-miR-1728-3p, gga-miR-1689∗, gga-miR-7b, gga-miR-1614∗, gga-miR-1664, gga-miR-1718, gga-miR-1629, gga-miR-1618 and gga-miR-1397∗) and downregulated (gga-let-7b, gga-let-7c, gga-miR-181b, gga-miR-133c, gga-miR-30a-3p, gga-miR-1788-5p, gga-miR-455-3p, gga-miR-449b∗, gga-miR-24 and gga-miR-449c∗) miRNAs were validated by quantitative real time PCR (qPCR) using universal reverse primer and miRNA specific forward primer (Supplementary Table S3) using NCode™ miRNA First-Strand cDNA Synthesis and qRT-PCR Kits (Invitrogen, Carlsbad, USA). On qPCR validation, all log fold changes of the upregulated miRNAs were found to be positive and those of the downregulated were observed to be negative thereby, validating the findings of the next generation sequencing data (Fig. 1A). Though the trend of the fold change was found to be the same, the magnitude of fold change, as identified by NGS and by qPCR, were different.

Fig. 1

A The comparative fold change for the top ten upregulated and downregulated miRNAs. The trend of the fold changes was found to be similar direction but with a difference in magnitude. B The candidate novel miRNAs identified in this study. The three novel candidate miRNAs (candidate 1, candidate 2 and candidate 3) with MFEI of 0.96, 0.95, and 0.98 were identified in the mock-infected chicken sample. miRanalyzer pipeline was used for predicting secondary structure followed by validation by MIReNA algorithm.

The small RNA reads which were present five or more times in the data were considered for novel miRNA prediction. Three novel mature miRNAs (candidate 1, 2 and 3) were identified by miRanalyzer pipeline (Hackenberg et al., 2011) (Fig. 1B) followed by validation by the MIReNA algorithm. Minimum Folding Energy Index (MFEI) was calculated as MFEI ​= ​AMFE/(G ​+ ​C)%; Adjusted Minimum Folding Energy (AMFE) = (MFE/length of a pre-miRNA sequence) ​× ​100. Structures having an MFEI value greater than or equal to 0.85 were predicted as novel miRNA.

The bioinformatics prediction of target genes of differentially expressed miRNAs using miRDB (http://mirdb.org/) identified 7,317 and 10,089 unique targets for UR and DR miRNAs, respectively. The predicted target genes were used for gene ontology (GO) analysis using DAVID bioinformatics resources 6.7 (https://david.ncifcrf.gov/) with the default parameters. The GO analysis using DAVID 6.7 showed many significant GO terms to be enriched. Of the total GO terms, 281 GO terms were common to both groups, 67 and 171 GO terms were found to be exclusively in TUR and TDR miRNAs. The GO terms, which were common to both the groups, included immune and apoptosis (Supplementary Fig. S2), signaling (Supplementary Fig. S3), and gene expression (Supplementary Fig. S4). The GO terms observed for the targets of upregulated (TUR) miRNAs included negative regulation of transcription from RNA polymerase II promoter, B cell differentiation, histone acetylation, positive regulation of myeloid cell differentiation, lipopolysaccharide-mediated signaling pathway, positive regulation of protein ubiquitination, etc (Fig. 2A). The important GO terms obtained for targets of downregulated (TDR) miRNAs were grouped into immune-related (Fig. 2B), apoptosis (Fig. 2C) and signaling (Fig. 2D).

Fig. 2

The GO terms observed for the targets of downregulated miRNAs. A The exclusive GO terms observed for targets of upregulated miRNAs. The GO terms including positive regulation of myeloid cell differentiation, histone acetylation, B cell differentiation, transmembrane receptor protein serine/threonine kinase signaling pathway, and negative regulation of RNA polymerase II promoter were found to be enriched. B The exclusive immune related GO terms observed for targets of downregulated miRNAs. GO terms like regulation of myeloid leukocyte differentiation, B cell activation, leukocyte activation, lymphocyte activation, T cell activation regulation of immune system process are found to be enriched. C The exclusive apoptosis related GO terms observed for targets of downregulated miRNAs. The GO terms including regulation of neuron apoptosis, regulation of apoptosis, cell death were enriched. D The exclusive signaling related GO terms observed for targets of downregulated miRNAs. These include GO terms like steroid hormone receptor signaling pathway, positive regulation of MAP kinase activity, Ras protein signal transduction, regulation of MAPKKK cascade and cell-cell signaling.

Six KEGG pathways were affected exclusively by TUR and included adipocytokine signaling pathway, notch signaling pathway, SNARE interactions in vesicular transport, vascular smooth muscle contraction, etc. Additionally, eleven pathways were found common for both the groups including calcium signaling pathway, MAPK signaling pathway, TGF-beta signaling pathway, insulin signaling pathway, regulation of actin cytoskeleton, endocytosis, Wnt signaling pathway, adherens junction, ubiquitin-mediated proteolysis, ErbB signaling pathway, etc. Moreover, seven pathways were found to be altered by TDR miRNAs (Supplementary Fig. S5). The functional protein association network of most upregulated and downregulated miRNA (gga-miR-1793 and gga-let-7b) was analyzed using STRING v.11 (https://string-db.org/) and is depicted in Fig. 3. Analysis of targets of gga-miR-1793 identified INTERPRO protein domain death-like domain superfamily (IPR011029), death domain (IPR000488) and SMART protein domain death domain, found in proteins involved in cell death (SM00005) with genes nuclear factor kappa b subunit 1 (NFKB1), NFKB inhibitor alpha (NFKBIA), inhibitor of nuclear factor kappa B kinase subunit beta (IKBKB), interleukin 1 beta (IL1B), component of inhibitor of nuclear factor kappa B kinase complex (CHUK) and REL proto-oncogene, NF-KB subunit (REL) to be involved (Fig. 3A). The analysis of targets of gga-let-7b showed mTOR signaling pathway (KEGG:gga04150) and p53 signaling pathway (KEGG:gga04115) with genes insulin like growth factor (IGF1), insulin like growth factor binding protein 1 (IGFBP1), insulin like growth factor binding protein 3 (IGFBP3), insulin like growth factor binding protein 4 (IGFBP4), insulin like growth factor binding protein 5 (IGFBP5), and insulin like growth factor 1 receptor (IGF1R) to be involved (Fig. 3B).

Fig. 3

The functional protein association network analysis using STRING v.11 for gga-miR-1793 and gga-let-7b. A The protein-protein interaction analysis for targets of gga-miR-1793. INTERPRO protein domain death like domain superfamily (IPR011029), death domain (IPR000488) and SMART protein domain death domain, found in proteins involved in cell death (SM00005) with genes NFKB1, NFKBIA, IKBKB, IL1B, CHUK and REL were observed to be involved. B The protein-protein interaction analysis for targets pf gga-let-7b showed mTOR signaling pathway (KEGG: gga04150) and p53 signaling pathway (KEGG: gga04115) with genes IGF1, IGFBP1, IGFBP3, IGFBP4, IGFBP5, and IGF1R to be involved.

The miRNAs viz., gga-miR-155, gga-miR-455-5p, gga-miR-214, gga-miR-140, gga-let-7c, gga-let-7j were also reported in some other AIV infection (Wang et al., 2009; Peng et al., 2015). miR-155, let-7b, miR-27b, miR-9 and miR-223 were reported to be associated with an acute respiratory viral infection and targeting innate immunity related pathways and genes, while miR-30a, let-7c, miR-221, miR-9, miR-98 and miR-214 were regulating apoptosis pathways (Leon-Icaza et al., 2019). Likewise, miR-let-7c is in the list of top most DR miRNAs and reported to repress the hepatitis-C virus multiplication by activating antiviral interferon and blocking of hepatitis-C virus protease activity (Chen et al., 2019).

The immune GO terms included leukocyte homeostasis, myeloid cell differentiation, hemopoiesis, and immune system development were found to be enriched for both the groups. The first cells in the respiratory tract to respond to influenza A virus infection are the immune cells. The viral entry has altered leukocyte homeostasis and leukocytes are critical in helping the host to eliminate the virus (Lamichhane and Samarasingh 2019). B cells and influenza-specific T lymphocytes are critical in clearing the virus and the activities of these T lymphocytes along with antibodies determine the capacity of the host to eliminate AIV (Hufford et al., 2015). Myeloid cells, representing a wide range of innate leukocytes, are critical for the host control over the virus. They identify pathogen-associated molecular patterns initiating a signaling cascade that marks the production of cytokines (Stegelmeier et al., 2019). Functional analysis also revealed GO terms related to gene expression and one of the interesting GO terms to be observed was the regulation of transcription from RNA polymerase II promoter indicating that the virus gained control over the host as Influenza virus RNA-dependent RNA polymerase associated with cellular RNA Polymerase II via cap snatching to generate primers for viral transcription which ultimately, resulted in host shut-off (Walker and Fodor, 2019).

Mitogen activated protein kinase (MAPK) signaling was found to be enriched in both groups as well as observed common for KEGG pathway analysis. This activation of MAPK pathways might have sent some signal to the nucleus for the production of favourable proteins. The p38 MAPK was listed as critical in mediating viral sensor signaling cascade and involved in the expression of antiviral chemokines and cytokines. Targeting the kinases was proposed to have a vital role in the treatment of influenza infections (Meineke et al., 2019). Another interesting pathway was the mTOR pathway which was found to be enriched for the TDR miRNAs. This pathway was also found to be significant during the protein-association network analysis of gga-let-7b using STRING v.11. The viruses have evolved a variety of mechanisms to opt with the mTOR pathway to make the host cell work their way. The analysis of gga-let-7b also pointed towards the p53 signaling pathway and p53 activity involved in influenza virus-induced cell death (Turpin et al., 2005). Further, some previous studies summarized the key role of nuclear factor kappa B (NF-κB) and insulin like growth factor (IGF) in AIV infection. Interaction analysis of gga-miR-1793 also revealed the complex role of NF-κB in H5N1 infection. NF-κB has been reported to be linked with host innate immune defense but influenza viruses have attained the capacity to turn this to their benefit (Ludwig and Planz, 2008). In another interaction analysis of miR-let-7b miRNA, IGF was found to be a key regulator of complex host-pathogen interaction. This plays a critical role in the vaccine-triggered immune response (Yoon et al., 2017) and is proposed to be a therapeutic target for humans in response to an influenza outbreak (Li et al., 2019).

This study identified the miRNAs which were differentially expressed during the H5N1 infection in chicken and by employing bioinformatics approaches the molecular pathways were predicted to gain insight into the host response after the viral entry. However, studies should be undertaken to find out the expressed miRNAs at different time intervals post viral entry and validation of the miRNA targeted genes would be more helpful to gain in-depth knowledge about the molecular pathogenesis in chicken during H5N1 infection.

Footnotes

This work was funded by the (grant number BT/IN/Indo-UK/FADH/48/AM/2013). We acknowledge Director, ICAR-National Institute of High Security Animal Diseases, Bhopal, India and Director, ICAR-Indian Veterinary Research Institute, Izatnagar, India for the help and support to carry out this work. The authors declare that they have no conflict of interest. This article was approved by the Ethics Committee of ICAR-NIHSAD vide approval no. 48/IAEC/HSADL/10. All institutional and national guidelines for the care and use of laboratory animals were followed/the experiments were conducted under the guidelines of Committee for the Purpose of Control and Supervision of Experiments on Animals (CPCSEA), Ministry of Environment and Forests, Govt. of India.

All the data generated during the current study are included in the manuscript. Supplementary data to this article can be found online at https://doi.org/10.1016/j.virs.2022.03.004.