The DHAV-1 HA5 strain utilized in this study was isolated from a one-week-old duckling in China, which exhibited characteristic clinical symptoms and pathological changes of DVH. The HA5 strain was propagated in eleven-day-old embryonating specific-pathogen-free (SPF) duck eggs and subsequently stored at −80 °C for further experiments. The SPF duck eggs were sourced from Shandong Haotai Experimental Animal Breeding Co., Ltd. Additionally, one-day-old SPF shelducks were procured from the same supplier and maintained in isolators until use.

Thirty two-day-old SPF shelducks were randomly selected and intramuscularly inoculated with a lethal dose of the DHAV-1 HA5 strain, each receiving 3 × 10⁵ ELD₅₀. At 6, 12, 24 and 48 hpi, three ducklings from the infected group were euthanized, and fresh liver tissue samples were collected. A control group, treated with an equivalent dose of isotonic sodium chloride. For each time point, there were three independent biological replicates, with one duckling per replicate, and no pooling of samples was performed. All liver tissues were immediately washed with ice-cold PBS, frozen in liquid nitrogen, and stored at −80 °C until processed for total RNA extraction for RT-qPCR and transcriptome sequencing analysis.

Liver tissue samples collected at various time points were subjected to ten-fold serial dilutions in sterile phosphate-buffered saline containing antibiotics. Following incubation, allantoic fluid was collected from both dead embryos and those surviving for 5 days post-inoculation. Embryos were examined for the presence of characteristic lesions. The 50% embryo infectious dose (EID₅₀) was determined according to the method of Reed and Muench.

mRNA sequencing was performed by Biomarker Technologies Co., Ltd. (Beijing, China). Total RNA from three ducks per group was extracted using TRIZOL reagent (Takara, China) according to the manufacturer’s protocol. The quantity and quality of the RNA were assessed, and only high-quality RNA was selected for library construction. Strand-specific RNA-Seq libraries were prepared using the Illumina Stranded mRNA Prep kit, ensuring the preservation of transcriptional directionality. Each library contained an average 300 bp cDNA insert flanked by adapter sequences. Sequencing was then carried out on an Illumina HiSeq 2,500 platform (Illumina, Inc.; San Diego, CA, USA) in 125 bp paired-end (PE125) mode. A typical coverage of at least 20–30 million clean reads per sample was achieved to ensure robust quantification of gene expression.

To ensure the reliability of the analyzed data, raw sequencing reads underwent rigorous filtering to eliminate low-quality and contaminated sequences, yielding a set of clean reads. Quality assessment of the filtered data was conducted using FastQC (v0.11.9) (17). The high-quality sequences were aligned to the Anas platyrhynchos reference genome (GenBank assembly accession: GCF_015476345.1, ZJU1.0) using the STAR (v2.7.9a) (18). We utilized the NCBI Anas platyrhynchos Annotation Release 104. However, we noted that potential gaps in microchromosome assembly and GC-rich promoter regions remains a known limitation in the duck genome which may impact the detection of low-abundance transcripts. Transcript reconstruction was carried out with Cufflinks software (19). The expression of each gene was calculated according to the reads per kilobase per million reads (RPKM). DEGs were screened using DESeq2 (v1.26.0), with selection criteria set at a false discovery rate (FDR) < 0.01 and |log2FoldChange| > 1 (20). The log2FoldChange threshold was applied globally across all time points, and DEGs were then compared and tracked across different time points for downstream analyses to observe temporal expression patterns. Functional enrichment analyses (COG, GO, and KEGG) were executed using clusterProfiler. In order to control the errors in these large-scale comparisons, multiple-testing correction was applied using the Benjamini-Hochberg (BH) procedure. GO terms and KEGG pathways were considered significantly enriched only if the adjusted p-value was < 0.05.

The sequencing results of selected mRNAs were validated using RT-qPCR. Primers for RT-qPCR were designed based on sequences obtained from the NCBI database (Table 1). Total RNA was extracted from the collected tissue samples using MolPure® TRIeasy Plus Total RNA Kit (Cat No. 19211; Yeasen, Shanghai, China) and quantified with a Nanodrop spectrophotometer (Thermo Fisher Scientific, USA). cDNA synthesis was performed using the HiScript® II RT Supermix for qPCR Kit (Vazyme, China) following the manufacturer’s protocol. RT-qPCR reactions were carried out in a 20 μL volume using the LightCycler (Roche Diagnostics, Mannheim, Germany) and the SYBR Green PCR Kit (Takara, Dalian, China) according to the manufacturer’s instructions. The PCR cycling conditions were set as follows: an initial denaturation at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s, and extension at 60 °C for 30 s. Each sample was analyzed in triplicate.

For statistical analysis, the relative expression levels of target genes in the infected and control groups were calculated using the 2^−∆∆Ct method. Expression levels were normalized to the housekeeping gene GAPDH, which served as an endogenous control, and results were expressed as fold changes in gene expression.

EID₅₀ of DHAV-1 in duck liver was measured at various time points following infection with the DHAV-1 HA5 strain. The viral titer increased progressively from 6 hpi, peaking at 24 hpi, before beginning to decline (Figure 1).

We conducted transcriptional profiling of livers infected with the HA5 strain. RNA-seq was employed to analyze the total RNA from the livers of SPF shelducks infected with the HA5 strain at various time points, as well as from mock-infected controls. The DEGs identified (FDR < 0.01 and |log2FoldChange| > 1) are listed in Supplementary Table S1.

At 6 hpi, 12 DEGs were identified, with an equal distribution of 6 upregulated and 6 downregulated genes. By 12 hpi, the number of DEGs increased to 73, comprising 67 upregulated and 6 downregulated genes. The most substantial transcriptional response was observed at 24 hpi, with 4,067 DEGs, including 1981 upregulated and 2086 downregulated genes. At 48 hpi, the number of DEGs decreased to 1,054, with 543 upregulated and 511 downregulated genes (Table 2). The progressive increase in DEGs during the course of DHAV-1 infection, peaking at 24 hpi, suggests that the infection intensifies between 12 to 24 hpi and then begins to subside after 24 hpi. The upregulated and downregulated genes are depicted in the volcano plots (Figures 2A1–A4). To distinguish between transient and sustained host responses, a cross-time-point overlap analysis was conducted (Figure 2C). We identified 35 genes that were consistently differentially expressed across the 12, 24, and 48 hpi time points. This persistent core signature was predominantly composed of genes involved in innate immune signaling and antiviral defense (e.g., IRF3, MX, CCL4 and IFIT5), indicating a sustained attempt by the host to suppress DHAV-1 replication.

DEGs were annotated using the COG database to determine the biological processes involved at each time point. At 6 hpi, the COG analysis revealed that most DEGs were involved in amino acid transport and metabolism, translation, ribosomal structure, and biogenesis (Figure 2B1). These findings suggest early alterations in protein synthesis and metabolic adaptation, likely as an immediate host response to viral entry.

By 12 hpi, the range of affected functional categories expanded significantly. Key processes included amino acid transport and metabolism, transcription, replication, recombination and repair, inorganic ion transport and metabolism, general function prediction only, and signal transduction mechanisms (Figure 2B2). The involvement of transcriptional and replication-related processes indicates the host’s attempt to modulate gene expression and potentially manage viral replication.

At 24 hpi, the host response became more complex. DEGs were associated with energy production and conversion, amino acid transport and metabolism, lipid transport and metabolism, translation, ribosomal structure, and biogenesis, transcription, replication, recombination and repair, secondary metabolites biosynthesis, transport and catabolism, general function prediction only, and signal transduction mechanisms (Figure 2B3). This broad range of functional categories suggests a multifaceted host response, involving metabolic shifts, cellular stress responses, and signal transduction pathways.

The gene expression profile at 48 hpi (Figure 2B4) closely mirrored that of 24 hpi, indicating that the host had settled into a sustained response pattern. The persistence of DEGs in categories such as energy production, amino acid and lipid metabolism, and signal transduction mechanisms suggests ongoing efforts to maintain cellular homeostasis and manage viral load.

To elucidate the functional roles of the DEGs, GO analysis was conducted. This analysis categorized and annotated the DEGs into three primary domains: biological processes, cellular components, and molecular functions.

At 6 hpi, the overall number of DEGs was relatively small (Figure 3A). Nevertheless, these DEGs were enriched in biological processes such as cellular process, single-organism process, and metabolic process. For cellular components, the DEGs were mainly associated with cellular parts, cells, organelles, and membranes. In terms of molecular functions, the DEGs were predominantly involved in binding, catalytic activity, and transporter activity.

By 12 hpi, there was an increase in the number of DEGs (Figure 3B). In addition to the categories identified at 6 hpi, there was enrichment in biological regulation, response to stimulus, immune system process, membrane, organelle part, receptor activity, and nucleic acid binding transcription factor activity. This increase indicates a more complex host response, including regulatory and immune-related processes.

At 24 hpi, DEGs were enriched across a broad range of categories in all three GO domains (Figure 3C). This expansion reflects the host’s intensified response to the viral infection, activating various biological processes, cellular components, and molecular functions.

The number of DEGs decreased and stabilized at 48 hpi (Figure 3D). Despite the reduced number, the categories of enriched GO terms remained consistent, indicating a sustained but regulated response. The stabilization suggests the host’s adaptation to the infection, maintaining essential cellular and molecular functions while possibly minimizing damage.

At 6 hpi, the DEGs were enriched in pathways related to RIG-I-like receptor signaling, ECM-receptor interaction, cytokine-cytokine receptor interaction, and cytosolic DNA-sensing (Figures 4A1,B1). These pathways are critical for initiating antiviral responses and modulating the immune response. Additionally, significant enrichment in the focal adhesion pathway suggests that early cellular changes occur in both structure and intercellular communication.

By 12 hpi, the immune profile expanded significantly. In addition to the pathways identified at 6 hpi, DEGs showed further enrichment in several critical signaling cascades: MAPK signaling, Toll-like receptor signaling, NOD-like receptor signaling, and cytosolic DNA-sensing pathways (Figures 4A2, B2). The continued involvement of these pathways, underscores a robust and systemic activation of the innate immune response and inflammatory signaling.

At 24 and 48 hpi, the complexity and diversity of DEG involvement increased (Figures 4A3–A4). DEGs were further enriched in pathways associated with Jak–STAT signaling pathway, oxidative phosphorylation, peroxisome, valine, leucine, and isoleucine degradation, fatty acid metabolism, glycine, serine, and threonine metabolism and PPAR signaling pathway. Among these, oxidative phosphorylation and lipid metabolism are particularly relevant to liver pathology. Oxidative phosphorylation plays a critical role in energy production and mitochondrial function, which is often disrupted during viral infections, contributing to liver damage. Lipid metabolism pathways are also crucial, as disturbances in lipid homeostasis can lead to liver steatosis and inflammation, common features of viral-induced liver injury. Furthermore, the Jak–STAT signaling pathway is pivotal in the regulation of immune responses and inflammation. Aberrant activation or dysregulation of this pathway has been linked to chronic inflammation and fibrosis in liver disease. These findings suggest that the dysregulation of these pathways may significantly contribute to DHAV-1-induced liver damage.

We randomly selected six genes (MDA5, IL-12, IFN-α, TLR-7, IL-8, and IL-10), which are implicated in RIG-I-like receptor and Toll-like receptor signaling pathways, for validation via RT-qPCR within 24 h post-stimulation. This approach aimed to assess the consistency and reproducibility of the DEGs identified through transcriptome sequencing. The RT-qPCR results confirmed significant changes in the expression of these genes, with upregulation patterns aligning closely with those observed in RNA-Seq data (Figure 5). These findings support the reliability of transcriptome sequencing in identifying DEGs.