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 105.78 copies (1μg cDNA)-1 at 12 hpi, and 108.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).