Characterization of DHAV-3 infection in duckling liver
One-day-old ducklings infected with virulent DHAV-3 exhibited visible ecchymosis hemorrhages throughout the liver at 24 hpi (Figure 2C), and showed typical clinical signs, such as mental 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 2E). 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 was 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 3A 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 repeats (Table 1). The length distribution of the clean reads was mainly between 21 and 24 nt (Figure 3B), which is consistent with previous reports [26–28]. These results indicate that miRNAs had 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 the miRBase 22.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 the 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 are 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 4A, Additional file 3). Compared to the mock-infected library, 102 miRNAs were upregulated and 54 miRNAs were downregulated in the DHAV-3-infected library (Figure 4A, Additional file 3).
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 that 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.
Furthermore, 23 DEMs, including 18 known miRNAs and 5 novel miRNAs, were predicted to target the DHAV-3 genome (Additional file 6). Most of these miRNAs were predicted to target genes of DHAV-3 structural proteins. For example, the VP1 gene of DHAV-3 was predicted to be targeted by 7 miRNAs, the VP0 gene was targeted by 5 miRNAs, and the VP3 gene was targeted by 3 miRNAs. Moreover, 4 non-structural protein genes of DHAV-3, including the 2C, 3B, 3C and 3D genes, were predicted targets of 7 DEMs. In addition, the apl-miR-340-5p was predicted to target the 5′-untranslated region (5′-UTR) of the DHAV-3 genome.
Global transcriptome profiles
In parallel with the miRNA profile, we also explored the global changes in gene expression associated with DHAV-3 infection to assess the effect of DEMs on their predicted targets. Six cDNA libraries representing the livers of ducklings 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. An overview of the sequencing data is shown in Additional file 7. Q20 and Q30 values were all > 95%, and 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 DHAV-3-infected duckling livers compared with those in the mock-infected ducklings at 24 hpi (Figure 5A and Additional file 8). 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 determine functionality of the DEGs, we performed GO enrichment analysis (Additional file 9). 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 enriched terms (Additional file 10).
KEGG pathway enrichment analysis was used to further define DEGs function in duckling liver after DHAV-3 infection (Additional file 11). The top 20 significantly enriched KEGG pathways are listed in Additional file 12 based on a q-value < 0.05. Five functional categories were identified to potentially play a role in 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 be involved in response to DHAV-3 infection in ducklings.
Integrated analysis of DEMs and DEGs
The correlation between miRNAs and mRNAs was based on the miRNA and transcriptome sequencing. To construct miRNA-mRNA networks, the genes identified as putative targets of DEMs, which were also differentially expressed in the transcriptome, were selected as the candidate target genes. According to the negative correlation principle of miRNA and mRNA, we anticipated an inverse relationship between the miRNAs and their target genes. Based on these criteria, 19,606 miRNA-mRNA interactions with the involvement of 155 DEMs and 4,484 DEGs were identified (Additional file 13).
Functional enrichment analysis of genes in these negatively correlated miRNA-mRNA pairs provided us with an integrated picture 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 in the majority (Figure 6A). Pathway analysis help us to obtain 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). A partial miRNA-mRNA interaction network of these pathways associated with DHAV-3 and host interactions were generated in Figure 7, using Cytoscape network construction software.
Quantitative RT-PCR validation of significant DEMs and DEGs
Quantitative RT-PCR was performed to 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. This confirmed the reliability of the sequencing techniques used in this work (Figure 4B and Figure 5B).