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

Background: Duck hepatitis A virus type 3 (DHAV-3) is one of the most harmful pathogens in the duck industry. However, the molecular mechanism underlying DHAV-3 infection in ducklings is remain poorly understood. To elucidate the genetic regulatory network for miRNA-mRNA and the signaling pathways involved in DHAV-3 infection in ducklings, we conducted global miRNA and mRNA expression profiling of duckling liver tissues infected with lethal DHAV-3 using high-throughput sequencing. Results: We found 156 differentially expressed miRNAs (DEMs) and 7717 differentially expressed miRNAs (DEGs) between mock-infected and DHAV-3-infected duckling livers. A total of 19,606 miRNA-mRNA pairs with negatively correlated expression patterns were identified in miRNA-mRNA networks constructed on the basis of these DEMs and DEGs. Moreover, immune-related pathways including the cytokine-cytokine receptor interaction, apoptosis, Toll-like receptor, Jak-STAT, and RIG-I-like receptor signaling pathway were significantly enriched through analyzing functions of mRNAs in the network in response to DHAV-3 infection. Besides, apl-miR-32-5p, apl-miR-125-5p, apl-miR-128-3p, apl-miR-460-5p, and novel-m0012-3p were identified as potential regulators in the immune-related signaling pathways during the DHAV-3 infection. Conclusions: To our knowledge, this is the first report on integrated analysis of miRNA-seq and mRNA-seq in DHAV-3-infected ducklings. The results indicated the important roles of miRNAs in regulating immune response genes and revealed the immune related miRNA-mRNA regulation network in the DHAV-3-infected duckling liver. Our findings may provide valuable information to further investigate the roles of miRNAs and their target genes in DHAV-3 replication and pathogenesis; additionally, they may offer clues for further understanding host-virus interactions. 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

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][21][22][23][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][26][27][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.

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.

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][27][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 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 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.
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][32][33][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 upregulated. 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][39][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.
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]. 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)

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

Virus and animals
Specific-pathogen-free (SPF) duck embryos were obtained from the Laboratory Animal lethal median dose (LD 50 )/mL.

RNA isolation and virus loads detection
Total RNA of all samples was isolated using TRIzol ® reagent (Invitrogen, Carlsbad, CA, 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]