Introduction

Duck viral hepatitis (DVH) is caused by the duck hepatitis A virus (DHAV), which has fast-spreading and highly infectious disease in ducklings (Sui et al. 2022). There are DHAV-1, DHAV-2, and DHAV-3, three types in DHAV, which is an Avihepatovirus genus (Rohaim et al. 2021). Among three types of DHAV, DHAV-1 is the most widespread and can cause neurological symptoms and acute hepatitis in ducklings (Li et al. 2022a, b, c). The other two types were only found in Asia (Feher et al. 2021). In addition, DHAV-1 infection can also cause egg drop in laying ducks (Lan et al. 2019). In detail, DHAV-1 infection can cause opisthotonos and severe liver damage, leading to rapid death of ducklings and threaten duck production. It is necessary to study the mechanism of DHAV-1-infected host target cells to better control DHAV-1 infection and reduce the economic losses of the duck industry.

RNA-sequencing (RNA-seq) technology is widely used in clinical research, biological research, and related drug development, which can comprehensively analyze the transcription and regulation of genes (Seyama et al. 2022). MicroRNAs (miRNAs) are small-noncoding RNAs involved in gene regulation and cell signaling, which are found in almost all eukaryotes (Qi et al. 2022). It is reported that miRNAs regulate most transcriptional genes by binding to the 3′ untranslated regions (Aass et al. 2022). miRNAs regulate mRNA through a variety of biological pathways such as immunity and metabolism (Movassagh et al. 2022). miRNAs play a role in a variety of biological processes, including viral inhibition of the host immune response, host response to infection, and virus–host interaction (Chatterjee et al. 2022). For instance, miRNA-4776 can regulate the production of influenza A virus in infected cell via regulating NFKBIB expression (Othumpangat et al. 2017). In addition, miRNA-223 regulated the inflammatory response induced by SARS-CoV infection by regulating host mRNAs (Morales et al. 2022). Meanwhile, miR-122 inhibits virus replication by upregulating the expression of type I interferon (Zhang et al. 2022).

Recently, combined mRNA-miRNA analysis has helped advance our understanding of miRNA function and disease mechanisms. For example, after association analysis of mRNA and miRNA in H9N2 avian influenza A virus-infected chicken, miR-126-5p regulates antiviral innate immunity by targeting TRAF3 (Wang et al. 2022). Especially, mRNA and miRNA profiling of DHAV-1-infected primary duck embryo fibroblast cells revealed miR-222a as a key factor in anti-DHAV-1 (Sui et al. 2021). However, the hepatic is the main target organ for DHAV-1-infected ducklings. The mRNA and miRNA expression profiles of DHAV-1 directly infected with DEHs were limited. In this work, we analyzed the expression of mRNA and miRNA in DEHs infected with DHAV-1 by high-throughput sequencing. Additionally, the association analysis of mRNA and miRNA can further understand the potential metabolic and signal pathway in DHAV-1-infected DEHs.

Materials and methods

Cell culture

DEHs were isolated using 14-day-old SPF duck embryos and cultured as previously reported (Qiu et al. 2022). The DEHs were cultured in DMEM (BasaIMedia, China) at 37 °C with 5% CO2, supplemented with 10% fetal calf serum (TIANHANG, China), 100 U/mL penicillin, and 100 μg/mL streptomycin.

TCID50 of DHAV-1

DEHs were cultured in 96-well plates at 37 °C with 5% CO2. DHAV-1 stock solution was diluted and replaced according to tenfold series, and incubated for 12–96 h to observe the existence of cytopathic effect (CPE). Reed Muench method is used to calculate TCID50.

DHAV-1 infection

DEHs were infected using the DHAV-1 LQ2 strain (0.01 multiplicity of infection (MOI) was given by the Shandong Institute of Poultry in China) at 37 °C for 2 h. The virus solutions were discarded, and DEHs were washed by D-Hank’s. Next, the DEHs were cultivated with 1% DMEM at 37 °C with 5% CO2. The DEHs were collected and observed at 12, 24, 36, and 48 h post-infection (hpi) after virus infection.

Sample collection and RNA quality control (QC)

The DEHs were split up into cell control (CC) group and virus control (VC) group, and each group was performed in triplicates (CC1, CC2, and CC3; VC1, VC2, and VC3). The total RNA in all groups was isolated utilizing Trizol (Vazyme, China) at 36 hpi. In addition, RNA QC was evaluated by Standard Sensitivity RNA Analysis Kit and Fragment Analyzer to specifications.

RNA sequencing

The mRNAs were enriched and used as a template to synthesize double-stranded cDNA. Subsequently, the linked products were amplified by PCR. Then, the mRNA libraries were constructed using the quality-qualified RNA. Finally, the mRNA libraries were further sequenced and inspected by DNBSEQ-T7 Genome Analyzer (MGI, China) and Agilent 2100 Bioanalyzer.

MicroRNA sequencing

The miRNA libraries were built and sequenced by DNBSEQ-T7 Genome Analyzer. Finally, the raw data were used to identify and predict miRNAs.

DEG and DEM analysis

We use DESeq2 software for DEG and DEM analysis. A smaller false discovery rate (FDR) value indicates a larger FoldChange, indicating a more significant difference in expression. Analysis of DEGs and DEMs with biological duplication is defined as | FoldChange |≥ 2 and p-adj ≤ 0.05.

Function and pathway analysis

To understand the potential function and pathways of DEGs and DEMs in DEHs, GO and KEGG were analyzed using clusterProfiler software. In addition, the screening criteria for significant results are p-adj < 0.05.

Combined mRNA-miRNA analysis

The mRNA-miRNA networks were created and visualized by Cytoscape 3.9.1, which based on mRNA and miRNA correlation.

RT-qPCR validation

The sequencing results of some screened mRNAs and miRNAs were verified by RT-qPCR. The relative PCR primers were designed based on the NCBI (Table 1). Foremost, total RNAs were extracted with Trizol and quantified with Nanodrop (Thermo, USA). Subsequently, the cDNA was formed using the HiScript® II RT supermix for qPCR Kit (Vazyme, China) and PCR Amplifier. For the miRNA, the cDNA was extracted by miRNA first-strand cDNA synthesis kit. In addition, the mRNA and miRNA analyses were achieved using qPCR mix kit, miRNA qPCR mix kit (Vazyme, China), and Light Cycler 96 instrument (Roche) and quantified using the 2−ΔΔCt method (Xie et al. 2018).

Table 1 Primer sequence
Full size table

Statistical analysis

All data were expressed as means ± standard deviations (SD), and the statistical tests were performed by IBM SPSS Statistics 26.0 and GraphPad Prism 9.0. The p value < 0.05 was considered statistically significant.

Results

TCID50 of DHAV-1 in DEHs

The TCID50 of DHAV-1 was 10−4.48/0.1 mL (Table 2).

Table 2 Observation results of DHAV-1 infection with DEHs
Full size table

Replication of DHAV-1 in DEHs

CPE and the relative expression of DHAV-1 in DEHs were detected at different time points after infection with DHAV-1. The observation results revealed that CPEs in infected DEHs were visible at 24 hpi and became more pronounced at 36 hpi, and cell fragmentation and detachment were severe at 48 hpi (Fig. 1a). Additionally, the relative expression of DHAV-1 increased rapidly between 24 and 36 hpi (Fig. 1b). However, the infected DEHs were severely shed at 48 hpi. Consequently, in order to ensure definite infection rates and eliminate ultra-CPE, we chose to collect DEHs at 36 hpi for further mRNA and miRNA sequencing.

Fig. 1
figure 1

Characteristics of DEHs following DHAV-1 infection. a Cell morphology of DEHs following DHAV-1 infection at 12, 24, 36, and 48 hpi. b The relative quantification of DHAV-1 on DEHs at 12, 24, 36, and 48 hpi. Values within a column without the same superscripts (a–d) are significant different (p < 0.05)

Full size image

RNA-seq data analysis

During DHAV-1 infection, the expression of mRNA in DEHs was applied to analyze by the Illumina high-throughput sequencing platform. More than 12.3 billion raw data were obtained of each sample, and clean reads were obtained after filtering and QC. The Q20 and Q30 were > 92.00%, and the GC was > 47.00%, showing the accuracy of the RNA-seq data (Table 3).

Table 3 Overview of RNA-seq data of each group
Full size table

miRNA sequencing data analysis

The miRNA expression profile of DEHs infected with DHAV-1 was obtained by the high-throughput sequencing platform. After filtering, clean reads surpassed 10 million and comprised > 93% of the mapped reads in reference genome (Table 4). Most miRNA lengths were 23 nt, which is in line with the typical miRNA length (Fig. 2a) (Jia et al. 2022). Furthermore, the first nucleotide bias of miRNAs indicated that a large proportion of first nucleotide biases were uridine (U), which is consistent with the miRNA typical characteristic (Fig. 2b) (Xie et al. 2015).

Table 4 Overview of miRNA sequencing data of each group
Full size table
Fig. 2
figure 2

Analysis of the miRNA data. a miRNA length distribution. b The nucleotide bias of known miRNAs. c The nucleotide bias of novel miRNAs

Full size image

DEG and DEM analysis

To understand the changes of DEHs infected by DHAV-1, the DEGs and DEMs in two group were compared. There were 2298 upregulated and 1112 downregulated DEGs in CC vs. VC comparison (Fig. 3a). Meanwhile, 97 upregulated and 45 downregulated DEMs were expressed (Fig. 3b). Notably, the hierarchical clustering of DEGs and DEMs revealed a strong correlation between the intra-group samples (Fig. 3c, d).

Fig. 3
figure 3

Analysis of the DEGs and DEMs. Volcano plot of DEGs in CC vs. VC (a) and DEMs in CC vs. VC (b). Hierarchical clustering of DEGs in CC vs. VC (c) and DEMs in CC vs. VC (d)

Full size image

Functional and pathway enrichment analysis

The functions and pathway of DEGs and DEMs were found using GO and KEGG analyses. In short, DEGs were related primarily to transmembrane transport, ion transport, immune response, and so on (Fig. 4a, b). Simultaneously, the KEGG pathway analysis found that the DEG pathways were mostly in connection with organismal systems, environmental information processing, and metabolism, such as RIG-I like receptor (RLR) pathway, arachidonic acid metabolism, cytokine and cytokine receptor, and so on (Fig. 4c, d).

Fig. 4
figure 4

GO and KEGG enrichment analysis based on DEGs of CC vs. VC group. a GO enrichment classification histogram. b Scatter diagram of GO enrichment. c KEGG enrichment classification histogram. d Scatter diagram of KEGG enrichment

Full size image

Moreover, DEM target genes were forecasted to characterize the biological functions of miRNAs. GO analysis revealed that the target genes were mainly related to molecular functions such as oxidoreductase activity, ATPase-coupled intramembrane, and voltage-gated calcium channel activity (Fig. 5a, b). Meanwhile, cGMP-PKG, cAMP, and calcium signaling pathway were the top KEGG pathways related to the target genes, indicating that the energy metabolism was significantly changes (Fig. 5c, d).

Fig. 5
figure 5

GO and KEGG enrichment analysis based on DEM target genes of CC vs. VC group. a GO enrichment classification histogram. b Scatter diagram of GO enrichment. c KEGG enrichment classification histogram. d Scatter diagram of KEGG enrichment

Full size image

mRNA–miRNA correlation network construction

To clarify the co-expression of mRNA and miRNA, the PPI network was constructed (Fig. 6a). In addition, a miRNA can concurrently target multiple genes, and the expression of miRNAs and mRNAs was multi-directional. Later, the differential expression of mRNA-miRNA pairs in the main pathways was visual (Fig. 6b, c). For instance, the downregulated novel-mir163 regulates TRAF3, and the upregulated gga-miR-132c-3p regulates MFN2.

Fig. 6
figure 6

mRNA–miRNA correlation networks in DHAV-1-infected DEHs. a All differential expression of mRNA-miRNA pairs. b Differential miRNA-mRNA interaction networks in RIG-I-like receptor signaling pathway. c Differential miRNA-mRNA interaction networks in mitophagy. Red: upregulated miRNA. Yellow: downregulated miRNA. Blue: upregulated mRNA. Green: downregulated mRNA.

Full size image

DEG and DEM validation

To test the accuracy of the RNA-seq and miRNA-seq results, the pivotal antiviral-related and mitophagy-related genes and miRNAs were examined. Briefly, the level of selected mRNAs was consistent with the RNA-seq results (Fig. 7a). Simultaneously, the tendency of selected miRNAs was in conformity with the miRNA-seq results (Fig. 7b). Moreover, the validation results of miRNA target genes were consistent with RNA-seq results (Fig. 7c).

Fig. 7
figure 7

Validation the mRNAs and miRNAs by RT-qPCR. a The expression of DEGs between RNA-seq and RT-qPCR. b The expression of DEMs between miRNA-seq and RT-qPCR. c The expression of DEM target genes between RNA-seq and RT-qPCR. *p < 0.05, **p < 0.01

Full size image

Discussion

As a viral pathogen, DHAV-1 has emerged as the leading cause of ducklings’ death (Xie et al. 2019). DHAV-1 is the most ordinary serotype and gives rise to serious acute hepatitis in ducklings (Liu et al. 2021). The typical symptom is punctate or ecchymotic bleeding on liver surface (Zhang et al. 2020). The current research revealed that m6A modification motifs and unique motifs in attenuated DHAV-1 infected duckling livers by m6A-Seq (Wu et al. 2022). Nevertheless, the DHAV-1-infected target cell molecular mechanism is unclarified. Some researchers found that miRNA plays a necessary role in the mutual effect between virus and host (Chirayil et al. 2018). Herein, we analyzed the mRNA and miRNA in DHAV-1-infected DEHs to illuminate duckling-DHAV-1 interaction mechanisms.

DHAV-1 infection causes the dysregulation of mRNAs expressions on DEHs. A total of 3410 DEGs were found in DHAV-1-infected DEHs. Subsequently, DEGs were forecasted for further functional enrichment analysis. GO analysis specified that the DEGs were centralized in multiple biotic processes, including transmembrane transport, ion transport, and immune response. Moreover, KEGG analysis revealed that DEGs involved in PI3K-Akt pathway, Ras pathway, RLR pathway, and mitophagy. The RLRs are necessary for discovering viral RNA and starting the innate immune response (Kell and Gale 2015). Briefly, in RNA viruses’ infections, the viral RNA binds to RLRs after entering the host cells and stimulates the activation of RIG-I or MDA5 (Zhang et al. 2021). The sensitized RIG-I or MDA5 recruits the MAVS at the mitochondrial membrane to promote the constitution of MAVS filaments (Li et al. 2021). Whereafter, MAVS filaments recruit TBK1, leading to the assembly of a complex crucial to the excitation of IRF3 and IRF7, thereby activating type I interferon responses and forming an antiviral state (Wang et al. 2021a, b, c). A previous study implicates RLR pathway as a main antiviral pathway against Kaposi’s sarcoma-associated herpesvirus infection (Zhao et al. 2018). In our study, several genes of RLR pathway, including DDX58 (RIG-I), DHX58 (LGP2), IFIH1 (MDA5), STING1, TRAF3, and IRF7, were significantly increased whether RNA-seq or RT-qPCR. The recent study demonstrated that long noncoding RNAs play essential roles in the RLR pathway by regulating related genes, defending DHAV-1 infection (Sui et al. 2022). For instance, DDX58 is a crucial factor in the cytoplasmic pattern recognition receptor family that initiates innate immune responses by detecting viral RNA (Frietze et al. 2022). Generally, IFIH1 (MDA5) can detect RNA from particular viral species, such as Picornaviridae (Stok et al. 2022). Hence, we speculate that RLR pathway is the key pathway for DEHs against DHAV-1 infection.

DHAV-1 infected motivates the dysregulation of miRNAs. Specifically, in DHAV-1-infected DEHs, a majority of the miRNA lengths ranged from 21 to 24 nt and most first nucleotide bases were U, indicating that the DEH miRNA characteristics were consistent with previous research (McGeary et al. 2022). In CC vs. VC comparison, we found a total of 142 DEMs. Then, the DEM target genes were predicted and performed functional annotation. In a word, DEM target genes have enriched energy metabolism. Moreover, mitochondria are the main organelles responsible for energy metabolism and ATP production in cells (Yan et al. 2021). Notably, not only DEGs but also DEM target genes has enriched mitophagy. In this study, the upregulated CALCOCO2 (NDP52) and OPTN and the downregulated PINK1 and MFN2 in mitophagy were found from RNA/miRNA-seq and PT-qPCR. Previous studies have reported that mitophagy is induced during virus infection, thereby destroying the function of mitochondria in managing inflammation and immune response (Li et al. 2022a, b, c; Wang et al. 2021a, b, c). For example, some investigators elaborated that syndrome coronavirus 2 ORF10 inhibits antiviral innate immune response by reducing MAVS through mitophagy (Li et al. 2022a, b, c). Additionally, another study indicated that influenza A virus induced mitophagy via interactions with TUFM and LC3B, thereby suppressing the antiviral innate immune response (Wang et al. 2021a, b, c). Therefore, the DHAV-1 infection may affect the antiviral pathway (RLR pathway) through mitochondria and mitophagy.

Next, the relevance between mRNA and miRNA was used for constituting the networks and implementing a comprehensively analysis. The correlation analysis identified 115 mRNA-miRNA pairs in CC vs. VC comparisons. Universally, majority of miRNAs could target more than one gene. Furthermore, the antiviral-related and mitophagy-related mRNA-miRNAs were established. The gga-miR-132c-3p (targeted by MFN2), gga-miR-6542-3p (targeted by PINK1), and novel-mir163 (targeted by TRAF3) were selected to verify, and the trend was consistent with miRNA-seq. However, these miRNAs have not been reported in antiviral studies. During RNA virus infection, MFN2 was responsible for IL-1β generation and interacted with MAVS to inhibit antiviral immunity (Silwal et al. 2021). It is reported that classical swine fever virus replication is accomplished by regulating PINK1 to induce mitophagy (Chengcheng et al. 2022). In addition, several studies elucidated that antiviral immunity is activated by catalyzing TRAF3 (Gao et al. 2021; Sun et al. 2020). Thus, these miRNAs may be prospective candidates for antiviral immunity.

In summary, this work is the first to identify mRNA and miRNA expression patterns in DHAV-1-infected DEHs for a further exploration of DHAV-1-target cell interactions. Moreover, we screen several dysregulated genes and miRNAs during DHAV-1 infection, involved in antiviral immunity and mitophagy. These findings provide serviceable mRNA and miRNA sequencing to better understand the changes of target cells after DHAV-1 infection. Further research will focus on the mechanisms of antiviral immunity and mitophagy during DHAV-1 infection, and find the potential therapeutic targets against DHAV-1.