Integrated analysis of miRNA-seq and mRNA-seq reveals immune related miRNA-mRNA regulation network in the DHAV-3-infected duckling liver

Abstract

Background

Duck viral hepatitis is a highly fatal, contagious, and rapidly spreading viral infection in young ducklings [1]. At least five viruses can cause this disease. The most widely distributed causative agent is duck hepatitis A virus (DHAV), which is classified as genus Avihepatovirus in the family Picornaviridae [2]. Based on phylogenetic analysis and virus neutralization tests, DHAV has been classified into DHAV-1, the most widespread serotype [3]; DHAV-3, a serotype isolated in Korea and China [4, 5]; and DHAV-2, a serotype isolated only in Taiwan [6]. With the widespread use of the officially approved DHAV-1 live vaccine in 2013, a higher rate of DHAV-3 causing DHAV infections has been observed in China since that time [7].

DHAV-3 causes acute hepatitis and induces typical duckling liver lesions, which is characterized by petechial hemorrhages of the liver surface [8]. The results of detecting DHAV-3 loads in the liver exhibited that the virus could replicate rapidly in the liver [9]. In addition, the highest viral titer was detected in the liver among DHAV-3 infected organs in vivo [10]. More importantly, the liver is a major site for the regulation of immune and inflammatory responses, plays a critical role in defending against invasive pathogens [11, 12]. However, there is limited information available about miRNA combined with gene expression profiles in duckling liver in response to DHAV-3 infection.

MicroRNAs (miRNAs) are a class of endogenous, small RNAs (21–24 nt) that possess a regulatory role in mRNA translation [13, 14]. They direct the regulation of gene expression by binding partially complementary target sites in the 3’ untranslated region (3’UTR) of mRNAs, cause post-transcriptional inhibition and facilitate target gene knockout and degradation [15]. It has been proved that miRNAs have significant effect on the proliferation of viruses and viral-host interactions [16-18]. For instance, previous reports have shown that several RNA viruses, including dengue virus (DENV), Japanese encephalitis virus (JEV) and West Nile virus (WNV), can be regulated in their replication by host miRNAs [18, 19]. On the other hand, it has been proved that miRNAs can regulate expression of host immune-related genes, and thereby inhibit or activate the downstream signaling pathway and mediate the anti-viral immune response [20-24].

Nowadays, high-throughput sequencing is an efficient means to obtain the sum of all miRNAs and mRNA products of a specific tissue or cell in a specific state. After the advent of genome-wide miRNA expression profile using deep sequencing, several studies in duck have been performed on miRNAs to examine their roles in growth, reproduction and virus infection [25-28]. However, there is currently no report about the integrated analysis of miRNA and mRNA expression responding to DHAV-3 infection in duckling liver. In this study, miRNA-seq and mRNA-seq were conducted to provide a global picture of the different mechanisms activated as response to DHAV-3 infection in duckling liver. The results may pave the way to understand the pathogenesis of DHAV-3 and the mechanisms of host-virus interactions.

Methods

Virus and animals

Specific-pathogen-free (SPF) duck embryos were obtained from the Laboratory Animal Center in the Harbin Veterinary Research Institute at the Chinese Academy of Agricultural Sciences (HVRI; Harbin, China). All SPF ducklings used in this study were hatched from SPF duck embryos and housed in isolators. All animal experiments followed the National Institutes of Health guidelines for the performance of animal experiments.

The highly virulent DHAV-3 SD strain used in this study was isolated in Shandong Province by our own laboratory. The virus was propagated in 1- to 3-day-old SPF ducklings. The virus titer was calculated by the Reed and Muench method and determined to be a 10^6.3 lethal median dose (LD₅₀)/mL.

DHAV-3 infection and sample collection

Twenty 1-day-old SPF ducklings were randomly chosen and intramuscularly (IM) inoculated with 0.2 mL of lethal DHAV-3 SD strain under 10⁵ dilutions. At 3, 6, 12, 24 and 36 hour post-infection (hpi), three infected ducklings were euthanized by decapitation, respectively. Their liver tissues were collected, washed with ice-cold PBS and frozen in liquid nitrogen immediately, then stored at –80˚C until the need of total RNA extraction for qRT-PCR, transcriptome and miRNA sequencing analysis.

RNA isolation and virus loads detection

Total RNA of all samples was isolated using TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instruction. The quantity and quality of the extracted total RNAs were detected and assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA).

Viral copies in total RNA were measured using TaqMan-based RT-qPCR established in our laboratory. The qPCR was carried out on an ABI Step One Plus thermocycler (Applied Biosystems, Foster City, CA, USA) with a final volume of 20 μL using the One Step PrimeScript RT-PCR Kit (TaKaRa, Dalian, China) according to the manufacturer’s instructions. The probe and primer were listed in Additional file 13, their final concentrations were 0.4 μM and 0.2 μM respectively. In addition, the amplicon of the forward and reverse primer pair was cloned into the pMD18-T vector (TaKaRa, Dalian, China) using as the DHAV-3 standard plasmid for constructing the standard curve of virus copy. Thermal cycling conditions included 42 °C for 5 min, 95 °C for 10 s, 40 cycles of 95 °C for 5 s and 60 °C for 34 s. Negative controls contained PCR-grade water. All reactions were conducted in triplicates.

miRNA-Seq and data analysis

Small RNA (sRNA) libraries have been generated from DHAV-3-infected and mock-infected duckling livers at 24 hpi. The sequencing of miRNAs was performed by Gene Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq^TM 2500.The raw reads were first trimmed, then mapped to the GenBank database (Release 209.0), Rfam database (11.0) and the reference duck genome (CAU_duck1.0). After the sRNAs screened with the above processes were excluded, the remaining sRNAs were searched against miRbase database (Release 21.0) to identify known miRNAs. In addition, all of the unannotated tags were aligned with reference genome. According to their genome positions and hairpin structures predicted by software Mireap_v0.2, the novel miRNA candidates were identified. In order to profile the differentially expressed miRNAs in DHAV-3-infected libraries versus mock-infected libraries, the miRNA expression level was calculated and normalized to transcripts per million (TPM). The p-value was adjusted using the Benjamini-Hochberg method [54, 55]. A corrected p-value < 0.05 and | log2 (fold change) | ≥ 1 were set as the threshold parameters for the significantly DEMs. Miranda (v3.3a), TargetScan (Version 7.0) and RNAhybrid (v2.1.2) & svm_light (v6.01) were used to predict the target genes. Then the intersection set of the three software’s results were chosen as predicted miRNA target genes.

mRNA-Seq and data analysis

The RNA samples used for the sRNA sequencing were the same as those used for the transcriptome sequencing. After total RNA extraction and DNase I treatment, mRNA was enriched by Oligo(dT) beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcripted into cDNA with random primers. Second-strand cDNA were synthesized. The cDNA fragments were purified, end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq^TM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China). Raw sequence data were assessed and filtered. Clean reads were aligned to the reference duck genome (CAU_duck1.0). The reconstruction of transcripts was carried out with software Cufflinks [56]. The expression of each gene was calculated according to the reads per kilo bases per million reads (RPKM). To identify DEGs, the edgeR package (http://www.rproject.org/) was used. In addition, | log2 (fold change) | ≥ 1 and the false discovery rate (FDR) < 0.05 were used as thresholds to define significant difference in gene expression

Construction of miRNA-mRNA regulatory networks

To elucidate the interaction network of miRNA-mRNA with positive and negative correlations, the DEMs and DEGs were constructed a miRNA-mRNA regulatory network. Expression correlation between miRNA-mRNA was evaluated using the Pearson correlation coefficient (PCC). Pairs with PCC < - 0.7 and p < 0.05 were selected as co-expressed negatively miRNA-mRNA pairs. And the visualization of miRNA-mRNA network was conducted through Cytoscape software (v3.6.0).

Functional analysis

To assess functional enrichment, Gene ontology (GO) functional analysis and a Kyoto Encyclopedia of Genes Genomes (KEGG) pathways analysis against the mRNAs in the network were implemented [57, 58]. The significantly enriched GO or KEGG terms were analyzed using p-value < 0.05.

Quantitative RT-PCR verification

The RNA samples used for the high-throughput sequencing assays were also used for the qPCR assay. The DEMs were validated by stem-loop qRT-PCR [59, 60]. Briefly, 1 μg of total RNA was reverse-transcribed to cDNA. cDNA was synthesized with stem-loop reverse-transcribed primers for DEMs validation, using PrimeScript^TM RT Master Mix (TaKaRa, Dalian, China). cDNA was synthesized with oligo(dT) primer for DEGs validation, using ReverTra Ace^® qPCR RT Kit (TOYOBO, Osaka, Japan). Next 2 μL cDNA was used for qPCR according to PowerUp^TM SYBR^® Green Master Mix kit (Applied Biosystems, Foster City, CA, USA) instructions. All primers are showing in Additional file 13. The small nuclear RNA (snRNA) U6 was used as endogenous internal control gene for miRNA and GAPDH was used for mRNA. ABI Step One Plus thermocycler (Applied Biosystems, CA, USA) was used to conduct qRT-PCR. In addition, reactions were performed in triplicate, and the relative expression levels were quantified by 2^-ΔΔCt method [61].

Results

Characterization of DHAV-3 infection in duckling liver

One-day-old ducklings infected with virulent DHAV-3 exhibited visibly ecchymosis hemorrhages throughout the liver at 24 hpi (Figure 1C), and generally showed typical clinical signs, such as metal depression and anorexia. Mortality occurred within 24-36 hpi. We detected DHAV-3 loads in the livers using TaqMan qRT-PCR. DHAV-3 replicated rapidly in the liver and reached 10^5.78 copies (1μg cDNA)^-1 at 12 hpi, and 10^8.62 copies (1μg cDNA)^-1 at 24 hpi (Figure 1E). On the basis of these results, we chose duckling livers at 24 hpi to perform the RNA-Seq analysis.

Analysis of small RNA libraries

To determine the miRNA expression pattern in response to DHAV-3 in the duckling liver, two sRNA libraries from mock-infected and DHAV-3-infected groups were constructed. Following high-throughput sequencing, the total number of raw reads collected from uninfected and infected duckling liver were 11,826,502 and 14,689,775, respectively (Table 1). After removing low-quality sequences, adapter sequences, and short reads smaller than 18 nt, clean reads were obtained (Figure 2A and Table 1). All of the clean reads were annotated and classified as rRNA, snRNA, snoRNA, snoRNA, tRNA, exon-sense, exon-antisense, intron-sense, intron-antisense, miRNA and repeat (Table 1). The length distribution of the clean reads were mainly between 21 and 24 nt (Figure 2B), which is consistent with previous reports [26-28]. These results indicate that miRNAs have been enriched successively from the two libraries.

Identification of known and novel miRNAs

To identify known miRNAs, we aligned sRNA from our libraries to the known mature miRNAs and their precursors in miRBase 21.0 database. A total of 349 and 291 known miRNAs were identified in the mock- and DHAV-3-infected liver libraries, respectively (Additional file 1). A number of unannotated sRNAs were present in each library (Table 1), and these sRNAs were matched to duck genome (CAU_duck1.0) for predicting novel miRNA candidates. A total of 109 and 34 novel miRNAs were predicted in the mock-infected and DHAV-3-infected duckling liver libraries using MIREAP_v0.2 software. These novel miRNAs were shown in Additional file 2.

Different expression analysis of miRNAs

For expression comparison between DHAV-3-infected and mock-infected libraries, miRNA sequences were analyzed through log2 (ratio) test and Chi-square 2×2 tests based on their normalized reads. Following significant differences standard (p ˂ 0.05 and |log2 (fold change)| ≥ 1), 156 differentially expressed miRNAs (DEMs) were detected in the two libraries (Figure 3A, Additional file 3). Comparing with the mock-infected library, 102 miRNAs were upregulated and 54 miRNAs were downregulated in the DHAV-3-infected library (Figure 3A, Additional file 3). Venn analysis was utilized to precisely examine DEMs (Figure 4). There were 59 miRNAs representing the numbers of miRNAs only expressed in DHAV-3 infected ducklings when compared with those in the control (Figure 4).

Target genes prediction for miRNAs

In order to understand the molecular function and biological processes of miRNAs during DHAV-3 infection in duckling liver, three independent algorithms, RNAhybrid，Miranda and TargetScan were used to predict the mRNA targets. A total of 26,886 genes for 398 known miRNAs and 119 novel miRNAs were predicted as potential miRNA targets. GO analysis revealed these predicted target genes were involved in biological process, cellular component and molecular function (Additional file 4). To explore the roles that DEMs may play in regulatory networks, we also assigned the putative miRNA targets to the KEGG pathways using the KEGG Orthology Based Annotation System (KOBAS). The results showed that 6,746 candidate genes were annotated for 241 biological processes. Many immune-related pathways, such as apoptosis (ko04210), ubiquitin mediated proteolysis (ko04120), FoxO signaling pathway (ko04068), NOD-like receptor signaling pathway (ko04621), p53 signaling pathway (ko04115), and RIG-I-like receptor signaling pathway (ko04622) were significantly enriched (Additional file 5). The results indicated that DEMs might play an important role in the virus-host interactions during the DHAV-3 infection in duck.

Global transcriptome profiles

In parallel with miRNA profile, we also explored the global changes in gene expression associated with DHAV-3 infection in order to assess the effect of DEMs on their predicted targets. Six cDNA libraries representing the livers of duckling in the mock-infected group (C1, C2, C3) and those in the DHAV-3-infected group (SD1, SD2, SD3) were constructed and subjected to Illumina sequencing. Overviews of the sequencing data are shown in Additional file 6. Q20 and Q30 values were all > 95%, and the percentage of GC content was similar, indicating that the data was accurate and reliable. Approximately 70% of the clean reads mapped to the duck reference genome (CAU_duck1.0) (Table 2). 

Analysis of DEGs

Genes with at least a twofold change (| log2 (fold change) | ≥ 1) and the FDR < 0.05 were considered DEGs. As a result, a total of 7717 DEGs including 6358 up-expressed and 1359 down-expressed were generated in the DHAV-3-infected duckling liver compared with those in the mock-infected ducklings at 24 hpi (Figure 5A and Additional file 7). It is notable that there are more up-regulated genes than down-regulated genes at 24 hpi. Overall, DHAV-3 infection had a significant impact on the global gene expression profile.

To ensure DEGs function, we performed GO enrichment analysis (Additional file 8). The DEGs were mainly enriched in 244 biological processes, mainly including the regulation of immune system process (GO:0002682), regulation of response to stimulus (GO:0048583), and lymphocyte activation (GO:0046649), of which the immune-related GO term accounted for the most (Additional file 9).

KEGG pathway enrichment analysis were employed to further define DEGs function in duckling liver after DHAV-3 infection (Additional file 10). The top 20 significant enrichment KEGG pathways are listed in Additional file 11 according to q-value < 0.05. Five functional categories were identified to potentially play important roles associated with DHAV-3 infection, including cytokine-cytokine receptor interaction (ko04060), Jak-STAT signaling pathway (ko04630), Toll-like receptor signaling pathway (ko04620), Influenza A (ko05164), and apoptosis (ko04210). These results suggested that genes in these pathways may play an important role in response to DHAV-3 infection in ducklings.

Integrated analysis of DEMs and DEGs

The correlation between miRNAs and mRNAs was performed on the basis of the miRNA and transcriptome sequencing. To construct miRNA-mRNA networks, the genes identified as putative targets of DEMs which were also differentially expressed in transcriptome were selected as the candidate target genes. According to the negative correlation principle of miRNA and mRNA, we anticipate an inverse relationship between the miRNAs and their target genes. On the basis of these criteria, 19,606 miRNA-mRNA interactions with the involvement of 155 DEMs and 4,484 DEGs were identified (Additional file 12).

Functional enrichment analysis of genes in these negatively correlated miRNA-mRNA pairs provided us with an overall clue of their functional roles during DHAV-3 infection in ducklings. Among the top 20 significantly enriched GO terms in biological process, immune- and signal-related terms were also the majority (Figure 6A). Pathway analysis helps to gain a better understanding of the biological function of genes. The KEGG pathway analysis demonstrated that the target genes of 155 DEMs were significantly associated with         cytokine-cytokine receptor interaction (ko04060), apoptosis (ko04210), Toll-like receptor signaling pathway (ko04620), FoxO signaling pathway (ko04068), and Jak-STAT signaling pathway (ko04630) (Figure 6B). Partial miRNA-mRNA interaction network of these pathways which are associated with DHAV-3 and host interactions were generated in Figure 7, using the Cytoscape network construction software.

Quantitative RT-PCR validation of significant DEMs and DEGs

Quantitative RT-PCR was performed to further investigate the differentially expressed miRNAs and mRNAs from sequencing data. Ten DEMs (including nine known miRNAs and a novel candidate miRNA) and ten DEGs were randomly selected for validation. Overall, the results analyzed by qPCR were in accordance with the high-throughput sequencing data, confirming the reliability of the sequencing technique used in this work (Figure 3B and Figure 5B).

Discussion

High-throughput sequencing was used as a powerful tool to reveal the expression profiles of miRNAs and genes in animals and plants, especially to identify differentially expressed miRNAs and genes under physiological perturbation. The combination of the miRNA profile with transcriptome has great significance to the study of miRNA regulatory mechanisms during the virus infecting procedure. Although some researchers had investigated changes of miRNA expression profile during the process of some viruses infected duck cells [27, 28], to our knowledge, the miRNA expression profile in ducklings exposed to DHAV-3 has not been identified. To this end, in the present study high-throughput sequencing was conducted to identify miRNAs and genes involved in the response to DHAV-3 in duckling liver. Moreover, we attempt to construct the miRNA-mRNA regulatory network according to the DEMs and DEGs datasets. Subsequently, differential expression of 156 miRNAs and 7717 mRNAs with 19,606 negative miRNA-mRNA pairs in DHAV-3-infected duckling liver were identified successfully. As shown in the miRNA-mRNA regulatory network (Figure 7), a single miRNA can target multiple genes. For example, apl-miR-125-5p can synchronously target RIG-I (retinoic acid inducible gene-I), DHX58 (DEXH box polypeptide 58), and TIRAP (toll/interleukin-1 receptor domain-containing adapter protein). On the contrary, one gene can be targeted by several miRNAs, such as SOCS1 (suppressor of cytokine signaling 1) in Figure 7 which can be simultaneously targeted by apl-miR-101-3p, apl-miR-128-3p, and apl-miR-148-3p. However, the network showed in Figure 7 is only a small part of the entire network. Apparently, the miRNA-mRNA regulatory network is extraordinarily complex.

MicroRNA-mediated changes in gene expression have been demonstrated to modulate viral replication, antiviral immune responses, viral latency, and pathogenesis [18, 29, 30]. In order to highlight the miRNAs related to interactions between DHAV-3 and duckling, we employed pathway analysis of mRNAs in the negative miRNA-mRNA network to understand their biological function and molecular mechanism. In performing KEGG pathway analysis for 4,484 mRNAs involved in the negative miRNA-mRNA network, immune-related pathways, including cytokine-cytokine receptor interaction, apoptosis, Toll-like receptor, FoxO, Jak-STAT, and RIG-I-like receptor signaling pathways were found to be significantly enriched. Similarly, pathway enrichment analysis for 7,717 DEGs of transcriptome showed almost the same results. In addition, this pathway enrichment result is consistent with another study that investigated transcriptome of the DHAV-3 infected duckling at 48 hpi [9]. Previous reports have shown that these pathways are involved in the host anti-virus processes [31-34]. Briefly, miRNAs have an effect on the regulation of immune responses and may be related to DHAV-3 replication at early infection in ducklings.

The innate immune system is the first line of defense against invading viruses [35, 36]. To date, accumulating data have proven a particular role for miRNAs in modulating many levels of the innate immune response [35, 37]. Pattern-recognition receptors (PRRs), such as Toll-like receptors (TLRs) and RIG-I-like receptors, recognize pathogens and initiate a serious of signaling programs that execute the first line of host innate immune responses. In our study, the transcriptional expressions of TLR2 and TLR7 were significantly up-regulated. Meanwhile, the negative miRNA-mRNA network revealed that there were 12 significantly down-regulated miRNAs targeting these TLRs. It has been reported that some miRNAs, like miR-140-5p and well-described miR-155, could target TLRs to regulate the innate immune response during the virus infection that had subsequent impact on virus replication [38-40]. Upon recognition of pathogens, TLRs recruit a specific set of adaptor molecules, such as MyD88 (myeloid differentiation factor 88) and TIRAP, then initiate downstream signaling events that leads to the secretion of inflammatory cytokines, type I IFN, chemokines, and antimicrobial peptides [41]. In the present study, we found that apl-miR-129-3p, apl-miR-148-5p, and apl-miR-3529-5p were down regulated upon DHAV-3 infection, and they all targeted MyD88 (Figure 7). In addition, significantly up-regulated TIRAP after DHAV-3 infection was targeted by 7 miRNAs in our network in Figure 7. Previous reports have shown that MyD88 was targeted by miR-7 and miR-21 during several viruses infection that led these viruses to evade the immune surveillance system to enhance them proliferation [42-44]. It is important to note that, apl-miR-32-5p had multiple target genes in the Toll-like receptor signaling pathway, including TLR2, TOLLP (toll-interacting protein), TRAF3 (TNF receptor-associated factor 3), and TBK1 (TANK-binding kinase 1). So did apl-miR-301-3p (targeting TLR7, TIRAP, TBK1, and CCL4) and novel-m0012-3p (targeting TLR7, TOLLP, TRAF3, and TBK1). This indicates that these three miRNAs may play a regulative role in the Toll-like receptor signaling pathway.

In the RIG-I-like receptor signaling pathway, the downregulated apl-miR-125-5p and apl-miR-148-5p were observed to bring about the upregulation of RIG-I in the DHAV-3 infected liver of ducklings. RIG-I recognizes viral RNA in the cytosol of most cell types that executes a critical role in the detection of RNA viruses [32]. Similarly, MDA5, melanoma differentiation-associated gene 5, which is another virus RNA detector was targeted by apl-miR-128-3p, apl-miR-27-3p, and apl-miR-499-3p. There are researches demonstrated that gga-miR-142-5p and miR-34b-5p could target MDA5 to promote IBDV (infectious bursal disease virus) and ALV-J (avian leukosis virus subgroup J) replication, respectively [45, 46]. In our constructed miRNA-mRNA network, four miRNAs (e.g. apl-miR-32-5p, novel-m0012-3p) were found to target RIG-I signaling-related molecule TIRM25 (tripartite motif containing 25). Moreover, MAVS (mitochondrial antiviral signaling protein), which functions as a platform for innate antiviral signal transduction and positively regulates type I interferon (IFN) production, were targeted by apl-miR-20-3p and apl-miR-460-5p (Figure 7). Detailed analysis revealed that miRNAs suppressed the expression of MAVS, thereby inhibiting MAVS-mediated NF-kappa B and IRF3 signaling, decreasing type I IFN and antiviral gene production, thus facilitating viral replication [47]. Briefly, DEGs with their target miRNAs in the RIG-I-like receptor signaling pathway may have important roles in regulating the innate immune system during DHAV-3 infection.

IFNs provide a robust line of host innate immune defense against viral infection [48]. Previous report has shown that IFN-α and IFN-β expressions were induced in the DHAV-3 infected duckling liver [10]. However, little is known about the mechanisms for miRNA regulation of IFN induction in DHAV-3 infection in ducklings. At the type I IFN receptor level, multiple miRNAs have been reported to target interferon alpha/beta receptor (IFNAR) to regulate type I IFN signaling [35, 49, 50]. In the current study, the expression of IFNAR1 and IFNAR2 had a 2.15-fold and a 4.88-fold increase, respectively. Apl-miR-455-5p and apl-miR-499-3p which targeted IFNAR1, plus novel-m0213-5p which targeted IFNAR2 were significantly down-regulated after DHAV-3 infection. Besides, in our negative miRNA-mRNA network, five miRNAs (e.g. apl-miR-128-3p, apl-miR-455-5p) targeted the signal transducer and activator of transcription 1 (STAT1) which is a member of JAK-STAT pathway. As a transcription factor, STAT1 plays an important role in antiviral responses by inducing expression of antiviral IFN-stimulated genes (ISGs) [51]. Hou et al. found that miR146a targeted STAT1 and attenuated the production of type I IFN-induced antiviral factors [52]. Moreover, it has been previously proven that miRNAs can also indirectly modulate type I IFN signal transduction by targeting SOCS1, a negative regulator of the JAK-STAT pathway [53]. Our study revealed that the negative regulator SOCS1 was targeted by apl-miR-101-3p, apl-miR-128-3p and apl-miR-148-3p in DHAV-3-infected duckling liver library. These results indicate that miRNAs regulate at many stages of the IFN response, such as through regulation of their cognate receptor components IFNAR and downstream signal transduction pathways including STATs. This may give an insight into how the post-transcriptional gene regulation mechanisms of IFN signaling in DHAV-3 infection.

It is important to note that several miRNAs discussed above, such as apl-miR-32-5p, apl-miR-125-5p, apl-miR-128-3p, apl-miR-460-5p, and novel-m0012-3p, can target vital genes that not only involved in the Toll-like receptor and the RIG-I-like receptor signaling pathway, but also involved in their downstream signaling pathways (Table 3). These results indicate the possibility that these miRNAs are vital regulators participating in the immune-related signaling pathways in DHAV-3-infected ducklings. And further investigations are needed to clarify the underlying mechanisms on how these miRNAs regulate their target genes and their roles in promoting or inhibiting DHAV-3 replication.

Conclusions

To our knowledge, this study is the first exploration to simultaneously characterize miRNAs and mRNAs of the duckling in response to DHAV-3 infection. Numerous significantly differentially expressed miRNAs and mRNAs were screened and identified. According to our data of integrated analysis, the negative miRNA-mRNA network were constructed. The functional enrichment analysis of genes involved in the network provided informative results which could help us articulate the host-virus interactions during DHAV-3 infection. For resistance to viral infection, the duckling mobilized many immune-related miRNAs and genes that took part in several pathways, such as Toll-like receptor signaling pathway, RIG-I-like receptor signaling pathway, and Jak-STAT signaling pathway. Besides, we found some miRNAs (e.g. apl-miR-32-5p, apl-miR-125-5p, apl-miR-128-3p, apl-miR-460-5p, and novel-m0012-3p) may play an important role in regulating the host defense response and having effect on DHAV-3 replication. We believe that these data will contribute to further studying the pathogenic mechanism of DHAV-3.

Abbreviations

DHAV-3: Duck hepatitis A virus type 3; DEMs: Differentially expressed miRNAs; DEGs: Differentially expressed miRNAs; miRNAs: MicroRNAs; 3’UTR: 3’ untranslated region; DENV: Dengue virus; JEV: Japanese encephalitis virus; WNV: West Nile virus; RIG-I: Retinoic acid inducible gene-I; DHX58: DEXH box polypeptide 58; TIRAP: Toll/interleukin-1 receptor domain-containing adapter protein; SOCS1: Suppressor of cytokine signaling 1; PRRs: Pattern-recognition receptors; TLRs: Toll-like receptors; MyD88: Myeloid differentiation factor 88; TOLLP: Toll-interacting protein; TRAF3: TNF receptor-associated factor 3; TBK1: TANK-binding kinase 1; IBDV: Infectious bursal disease virus; ALV-J: Avian leukosis virus subgroup J; TIRM25: Tripartite motif containing 25; MAVS: Mitochondrial antiviral signaling protein; IFN: Interferon; IFNAR: Interferon alpha/beta receptor; STAT1: Signal transducer and activator of transcription 1; ISGs: IFN-stimulated genes; SPF: Specific-pathogen-free; hpi: Hour post-infection; sRNA: Small RNA; TPM: Transcripts per million; PCC: Pearson correlation coefficient; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes Genomes; snRNA: small nuclear RNA

Declarations

Ethics approval and consent to participate

The study was approved by the Animal Ethics Committee of Jiangsu Academy of Agricultural Sciences. The experimental procedures were performed strictly according to the approved guidelines.

Consent for publication

Not applicable.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was financially supported by the National Key R&D Program of China (2017YFD0500804), China Postdoctoral Science Foundation (2018M642195), and Jiangsu Planned Projects for Postdoctoral Research Funds (2018K069B). The funders had no role in the experimental design and samples collection, as well as data analysis and interpretation, and in writing of this manuscript.

Authors' contributions

FW and XZ conceived and designed the experiments; FW, FL and XF performed the majority of experiments; JC, CL, QP, and HS participated part of the experiments; FW prepared and revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thanks Gene Denovo Biotechnology Company for performing the library preparation and small RNA sequencing.

References

1. World organization for animal health. Duck Virus Hepatitis. In: Manual of Diagnostic Tests and Vaccines for Terrestrial Animals. 2019. http://www.oie.int/fileadmin/Home/eng/Health_standards/tahm/3.03.08_DVH.pdf.
2. Kang M, Roh JH, Jang HK. Protective efficacy of a bivalent live attenuated vaccine against duck hepatitis a virus types 1 and 3 in ducklings. Vet Microbiol. 2018;214:108-122. doi: 10.1016/j.vetmic.2017.12.018.
3. Swayne DE, Glisson JR, Mcdougald LR, Nolan LK, Suarez DL, Nair VL, editors. Diseases of poultry, 13th edition. New Jersey: Wiley-Blackwell Press; 2013. p. 417-463.
4. Kim MC, Kwon YK, Joh SJ, Kim SJ, Tolf C, Kim JH, et al. Recent Korean isolates of duck hepatitis virus reveal the presence of a new geno- and serotype when compared to duck hepatitis virus type 1 type strains. Arch Vrol. 2007;152:2059-2072. doi: 10.1007/s00705-007-1023-0.
5. Li J, Bi Y, Chen C, Yang L, Ding C, Liu W. Genetic characterization of duck hepatitis A viruses isolated in China. Virus Res. 2013;178:211-216. doi: 10.1016/j.virusres.2013.10.007.
6. Tseng CH, Knowles NJ, Tsai HJ. Molecular analysis of duck hepatitis virus type 1 indicates that it should be assigned to a new genus. Virus Res. 2007;123:90-203. doi: 10.1016/j.virusres.2006.09.007.
7. Wen X, Zhu D, Cheng A, Wang M, Chen S, Jia R, et al. Molecular epidemiology of duck hepatitis a virus types 1 and 3 in china, 2010-2015. Transbound Emerg Dis. 2018;65:10-15. doi: 10.1111/tbed.12741.
8. Yugo DM, Hauck R, Shivaprasad HL, Meng XJ. Hepatitis Virus Infections in Poultry. Avian Dis. 2016;60:576-588. doi: 10.1637/11229-070515-Review.1.
9. Zhang XL, Cao C, Liu Y, Qi HH, Zhang WJ, Hao CX, et al. Comparative liver transcriptome analysis in ducklings infected with duck hepatitis a virus 3 (DHAV-3) at 12 and 48 hours post-infection through RNA-seq. Vet Res. 2018;49:52. doi: 10.1186/s13567-018-0545-7.
10. Zhang XL, Cao C, Qu ZH, Zhang WJ, Liu Y, Qi HH, et al. Pathogenicity of duck hepatitis a virus type 3 and innate immune responses of the ducklings to virulent DHAV-3. Mol Immunol. 2018;95:30-38. doi: 10.1016/j.molimm.2018.01.007.
11. Doherty DG. Immunity, tolerance and autoimmunity in the liver: A comprehensive review. J Autoimmun. 2016;66:60-75. doi: 10.1016/j.jaut.2015.08.020.
12. Strnad P, Tacke F, Koch A, Trautwein C. Liver - guardian, modifer and target of sepsis. Nat Rev Gastroenterol Hepatol. 2017;14:55-66. doi: 10.1038/nrgastro.2016.168.
13. Ambros V. MicroRNAs: tiny regulators with great potential. Cell. 2001;107,:823-826.
14. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281-297.
15. Kloosterman WP, Plasterk RH. The diverse functions of microRNAs in animal development and disease. Dev Cell. 2006;11:441-450. doi: 10.1016/j.devcel.2006.09.009.
16. Yao YX, Nair V. Role of virus-encoded microRNAs in avian viral diseases. Viruses. 2014;6:1379-1394. doi: 10.3390/v6031379.
17. Diogo P, Azevedo-Pereira JM. The role of microRNAs in the pathogenesis of herpesvirus infection. Viruses. 2016;8:156. doi: 10.3390/v8060156.
18. Trobaugh DW, Klimstra WB. MicroRNA regulation of RNA virus replication and pathogenesis. Trends Mol Med. 2017;23:80-93. doi: 10.1016/j.molmed.2016.11.003.
19. Slonchak A, Shannon RP, Pali G, Khromykh AA. Human MicroRNA miR-532-5p exhibits antiviral activity against West Nile virus via suppression of host genes SESTD1 and TAB3 required for virus replication. J Virol. 2016;90:2388-2402. doi: 10.1128/JVI.02608-15.
20. Dang L, Teng M, Li HZ, Ma SM, Lu QX, Hao HF, et al. Marek's disease virus type 1 encoded analog of miR-155 promotes proliferation of chicken embryo fibroblast and DF-1 cells by targeting hnRNPAB. Vet Microbiol. 2017;207:210-218. doi: 10.1016/j.vetmic.2017.06.015.
21. Fu M, Wang B, Chen X, He Z, Wang Y, Li X, et al. MicroRNA gga-miR-130b suppresses infectious bursal disease virus replication via targeting of the viral genome and cellular suppressors of cytokine signaling 5. J Virol. 2018;92:e01646-17. doi: 10.1128/JVI.01646-17.
22. Hazra B, Kumawat KL, Basu A. The host microRNA miR-301a blocks the IRF1 mediated neuronal innate immune response to Japanese encephalitis virus infection. Sci Signal. 2017;10:eaaf5185. doi: 10.1126/scisignal.aaf5185.
23. Smith JL, Jeng S, McWeeney SK, Hirsch AJ. A MicroRNA screen identifies the Wnt signaling pathway as a regulator of the interferon response during Flavivirus infection. J Virol. 2017;91:e02388-16. doi: 10.1128/JVI.02388-16.
24. Wang B, Fu M, Liu Y, Wang Y, Li X, Cao H, et al. Gga-miR-155 enhances type I interferon expression and suppresses infectious burse disease virus replication via targeting SOCS1 and TANK. Front Cell Infect Microbiol. 2018;8:55. doi: 10.3389/fcimb.2018.00055.
25. Yu DB, Jiang BC, Gong J, Dong FL, Lu YL, Yue HJ, et al. Identification of novel and differentially expressed MicroRNAs in the ovaries of laying and non-laying ducks. J Integr Agric. 2013;12:136-146. doi: 10.1016/S2095-3119(13)60214-2.
26. Chen XY, Ge K, Wang M, Zhang C, Geng ZY. Integrative analysis of the Pekin duck (Anas anas) microRNAome during feather follicle development. BMC Dev Biol. 2017;17:12. doi: 10.1186/s12861-017-0153-1.
27. Cui M, Jia RY, Huang J, Wu XL, Hu ZQ, Zhang XC, et al. Analysis of the microRNA expression profiles in DEF cells infected with duck Tembusu virus. Infect Genet Evol. 2018;63:126-134. doi: 10.1016/j.meegid.2018.05.020.
28. Wu XL, Jia RY, Zhou JK, Wang MS, Chen S, Liu MF, et al. Virulent duck enteritis virus infected DEF cells generate a unique pattern of viral microRNAs and a novel set of host microRNAs. BMC Vet Res. 2018;14:144. doi: 10.1186/s12917-018-1468-2.
29. Samir M, Vaas LAI, Pessler F. MicroRNAs in the host response to viral infections of veterinary importance. Front Vet Sci. 2016;3:86. doi: 10.3389/fvets.2016.00086.
30. Bernier A, Sagan SM. The diverse roles of microRNAs at the host-virus interface. Viruses. 2018;10:440. doi: 10.3390/v10080440.
31. Takaoka A, Hayakawa S, Yanai H, Stoiber D, Negishi H, Kikuchi H, et al. Integration of interferon-α/β signalling to p53 responses in tumour suppression and antiviral defence. Nature. 2003;424:516-523. doi: 10.1038/nature01850.
32. Chan YK, Gack MU. RIG-I-like receptor regulation in virus infection and immunity. Curr Opin Virol. 2015;12:7-14. doi: 10.1016/j.coviro.2015.01.004.
33. Santhakunar D, Rubbenstroth D, Martinez-Sobrido L, Munir M. Avian interferons and their antiviral effectors. Front Immunol. 2017;8:49. doi: 10.3389/fimmu.2017.00049.
34. Zhang LT, Liu J, Qian LP, Feng Q, Wang XF, Yuan YK, et al. Induction of OTUD1 by RNA viruses potently inhibits innate immune responses by promoting degradation of the MAVS/TRAF3/TRAF6 signalosome. PLoS Pathog. 2018;14:e1007067. doi: 10.1371/journal.ppat.1007067.
35. Forster SC, Tate MD, Hertzog PJ. MicroRNA as type I interferon-regulated transcripts and modulators of the innate immune response. Front Immunol. 2015;6:334. doi: 10.3389/fimmu.2015.00334.
36. Kawai T, Akira S. Toll-like receptors and their crosstalk with other innate receptors in infection and immunity. Immunity. 2011;34:637-650. doi: 10.1016/j.immuni.2011.05.006.
37. Momen-Heravi F, Bala S. miRNA regulation of innate immunity. J Leukoc Biol. 2018;103:1205-1217. doi: 10.1002/JLB.3MIR1117-459R.
38. Zhang Y, Shao LY. Decreased microRNA-140-5p contributes to respiratory syncytial virus disease through targeting Toll-like receptor 4. Exp Ther Med. 2018;16:993-999. doi: 10.3892/etm.2018.6272.
39. Oliveira BRSM, Vieira FV, Vieira DD, da Silva SEL, Gameiro R, Flores EF, et al. Expression of miR-155 associated with Toll-like receptors 3, 7, and 9 transcription in the olfactory bulbs of cattle naturally infected with BHV5. J Neurovirol. 2017;23:772-778. doi: 10.1007/s13365-017-0564-6.
40. Sarkar N, Panigrahi R, Pal A, Biswas A, Singh SP, Kar SK, et al. Expression of microRNA-155 correlates positively with the expression of Toll-like receptor 7 and modulates hepatitis B virus via C/EBP-beta in hepatocytes. J Viral Hepatitis. 2015;22:817-827. doi: 10.1111/jvh.12390.
41. Kawai T, Akira S. The role of pattern-recognition receptors in innate immunity: Update on Toll-like receptors. Nat Immunol. 2010;11:373-384. doi: 10.1038/ni.1863.
42. Chen YN, Chen JB, Wang H, Shi JJ, Wu KL, Liu S, et al. HCV-Induced miR-21 Contributes to Evasion of Host Immune System by Targeting MyD88 and IRAK1. PLoS Pathog. 2013;9:e1003248. doi: 10.1371/journal.ppat.1003248.
43. Feng N, Zhou ZZ, Li YX, Zhao LF, Xue ZF, Lu R, et al. Enterovirus 71-induced has-miR-21 contributes to evasion of host immune system by targeting MyD88 and IRAK1. Virus Res. 2017;37:27-36. doi: 10.1016/j.virusres.2017.05.008.
44. Huang Y, Wang W, Xu ZQ, Pan JL, Zhao Z, Ren Q, et al. Eriocheir sinensis microRNA-7 targets crab Myd88 to enhance white spot syndrome virus replication. Fish Shellfish Immunol. 2018;79:274-283. doi: 10.1016/j.fsi.2018.05.028.
45. Ouyang W, Qian J, Pan QX, Wang JY, Xia XX, Wang XL, et al. gga-miR-142-5p attenuates IRF7 signaling and promotes replication of IBDV by directly targeting the chMDA5 ' s 3 ' untranslated region. Vet Microbiol. 2018;221:74-80. doi: 10.1016/j.vetmic.2018.05.018.
46. Li ZH, Luo QB, Xu HP, Zheng M, Abdalla BA, Feng M, et al. MiR-34b-5p Suppresses Melanoma Differentiation-Associated Gene 5 (MDA5) Signaling Pathway to Promote Avian Leukosis Virus Subgroup J (ALV-J)-Infected Cells Proliferaction and ALV-J Replication. Front Cell Infect Microbiol. 2017;7:17. doi: 10.3389/fcimb.2017.00017.
47. Xu TJ, Chu Q, Cui JX, Bi DK. Inducible MicroRNA-3570 Feedback Inhibits the RIG-I-Dependent Innate Immune Response to Rhabdovirus in Teleost Fish by Targeting MAVS/IPS-1. J Virol. 2018;92:e01594-17. doi: 10.1128/JVI.01594-17.
48. Chen K, Liu J, Cao XT. Regulation of type I interferon signaling in immunity and inflammation: A comprehensive review. J Autoimmun. 2017;83:1-11. doi: 10.1016/j.jaut.2017.03.008.
49. Papadopoulou AS, Dooley J, Linterman MA, Pierson W, Ucar O, Kyewski B, et al. (2012). The thymic epithelial microRNA network elevates the threshold for infection-associated thymic involution via miR-29a mediated suppression of the IFN-alpha receptor. Nat Immunol. 2012;13:181-187. doi: 10.1038/ni.2193.
50. Zhou CX, Yu Q, Chen L, Wang J, Zheng SD, Zhang JK. A miR-1231 binding site polymorphism in the 3’UTR of IFNAR1 is associated with hepatocellular carcinoma susceptibility. Gene. 2012;507:95-98. doi: 10.1016/j.gene.2012.06.073.
51. Cheon H, Yang J, Stark GR. The functions of signal transducers and activators of transcriptions 1 and 3 as cytokine-inducible proteins. J Interferon Cytokine Res. 2011;31:33-40. doi: 10.1089/jir.2010.0100.
52. Hou ZH, Han QJ, Zhang C, Tian ZG, Zhang J. miR146a impairs the IFN-induced anti-HBV immune response by downregulating STAT1 in hepatocytes. Liver Int. 2014;34:58-68. doi: 10.1111/liv.12244.
53. Piganis RAR, De Weerd NA, Gould JA, Schindler CW, Mansell A, Nicholson SE, et al. Suppressor of cytokine signaling (SOCS) 1 inhibits type I interferon (IFN) signaling via the interferon alpha receptor (IFNAR1)-associated tyrosine kinase Tyk2. J Biol Chem. 2011;286:33811-33818. doi: 10.1074/jbc.M111.270207.
54. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to Multiple Testing. J R Statist Soc. 1995;57:289-300.
55. Storey JD, Tibshirani R. Statistical signifcance for genomewide studies. Proc Nati Acad Sci U S A. 2003;100:9440-9445. doi: 10.1073/pnas.1530509100.
56. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562-578. doi: 10.1038/nprot.2012.016.
57. Young MD, Wakeeld MJ, Smyth GK, Oshlack A. goseq: Gene Ontology testing for RNA-seq datasets Reading data. Genome Biol. 2012;1-21.
58. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2007;36(Database):D480-D484. doi: 10.1093/nar/gkm882.
59. Mestdagh P, Feys T, Bernard N, Guenther S, Chen C, Speleman F, et al. High-throughput stem-loop RT-qPCR miRNA expression profiling using minute amounts of input RNA. Nucleic Acids Res. 2008;36:e143. doi: 10.1093/nar/gkn725.
60. Kramer MF. Stem-loop RT-qPCR for miRNAs. Curr Protoc Mol Biol. 2011;95:15.10.1-15.10.15. doi: 10.1002/0471142727.mb1510s95.
61. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt method. Methods. 2001;25:402-408. doi: 10.1006/meth.2001.1262.

Tables

Table1. Overview of miRNA-Seq in the mock-infected and DHAV-3-infected libraries.

Table 2. Major characteristics of mRNA libraries and reads mapping to the duck reference genome.

Table 3. Candidate immune-related miRNAs in the DHAV-3-infected duckling.