Comparative Transcriptomic Analysis of Primary Duck Hepatocytes Provides Insight into Differential Susceptibility to DHBV Infection

Primary duck hepatocytes (PDH) displays differential susceptibility to duck hepatitis B virus when maintained in the media supplemented with fetal bovine serum or dimethyl sulfoxide (DMSO) which has been widely used for the maintenance of hepatocytes, and prolonging susceptibility to hepadnavirus. However the mechanism underlying maintenance of susceptibility to hepadnavirus by DMSO treatment remains unclear. In this study, a global transcriptome analysis of PDHs under different culture conditions was conducted for investigating the effects of DMSO on maintenance of susceptibility of PDH to DHBV in vitro. The 384 differential expressed genes (DEGs) were identified by comparisons between each library pair (PDHs cultured with or without DMSO or fresh isolated PDH). We analyzed canonical pathways in which the DEGs were enriched in Hepatic Fibrosis / Hepatic Stellate Cell Activation, Bile Acid Biosynthesis and Tight Junction signaling. After re-annotation against human genome data, the 384 DEGs were pooled together with proteins belonging to hepatitis B pathway to construct a protein-protein interaction network. The combination of decreased expression of liver-specific genes (CYP3A4, CYP1E1, CFI, RELN and GSTA1 et al) with increased expression of hepatocyte-dedifferentiation-associated genes (PLA2G4A and PLCG1) suggested that in vitro culture conditions results in the fading of hepatocyte phenotype in PDHs. The expression of seven DEGs associated with tight junction formation (JAM3, PPP2R2B, PRKAR1B, PPP2R2C, MAGI2, ACTA2 and ACTG2) was up-regulated after short-term culture in vitro, which was attenuated in the presence of DMSO. Those results could shed light on DHBV infection associated molecular events affected by DMSO.


Introduction
Hepatitis induced by hepatitis B virus (HBV) infection remains a major public health problem with over 350 million chronically infected carriers worldwide [1]. The lack of an appropriate cell model supporting HBV infection has been a major limitation for the study of the pathogenesis of hepadnavirus infection. NTCP has been demonstrated to be the functional receptor of HBV on the cell membrane recently [2]. However, more detailed interactions between hepadnavirus and hepatocytes need to be further investigated. Because of their susceptibility to DHBV infection, primary duck hepatocytes (PDH) serve as a valuable cell model to elucidate the key features of hepadnaviral infection [3][4][5][6][7]. A reverse transcribed RNA intermediate in the replication of hepadnavirus was first described in the DHBV-PDH infection model [3]. However, cultured in vitro over a period of time, PDH lose the hepatocyte phenotype and their susceptibility to DHBV infection rapidly diminishes [4,8]. Only 10%-20% cells retain susceptibility to DHBV infection 6 d after plating under in vitro culture condition. The rapid loss of susceptibility to DHBV infection has limited the application of PDH in hepadnavirus study [8].
Several modified culture methods have previously been developed for the maintenance of the primary hepatocyte, including the addition of a low concentration of DMSO to the culture medium [9,10], using plates pre-coated with rat-tail collagen [11], or supplementation with cell factors (epidermal growth factor or hepatocyte growth factor) [9]. Among these methods, the addition of DMSO has been widely used for the maintenance of hepatocytes, and prolonging susceptibility to hepadnavirus [12]. It shows differential susceptibility to DHBV infection between PDH maintained in culture media supplemented with 5% fetal bovine serum and 1.5% dimethyl sulfoxide (DMSO) [8]. It has been revealed that maintenance of polarization phenotypes and prevention of tight junction formation in primary human hepatocytes or HepaRG cells by DMSO is essential for the entry of the HBV under in vitro culture conditions [13]. However the mechanism underlying maintenance of susceptibility to hepadnavirus by DMSO remains to be elucidated.
Next-generation sequencing technology provides a useful tool to understand the activation patterns of the cellular response to external stimuli. Here we conducted a global transcriptome analysis of PDH under different culture conditions to investigate the effects of DMSO on maintenance of susceptibility of PDH to DHBV in vitro. The differential expressed genes that were identified by comparison between PDHs cultured with or without DMSO, were pooled together with proteins belonging to hepatitis B pathway for the construction of a protein-protein interaction (PPI) network. In the network, the expression of a cluster DEGs associated with tight junction formation was up-regulated after short-term culture in vitro, while the upregulation was attenuated in the presence of DMSO. Supplementation with DMSO also relieved the declining of anti-virus associated DEGs expression in PDHs in vitro. Those results could shed light on DHBV infection associated molecular events affected by DMSO.

Differential susceptibility to DHBV infection of primary duck hepatocytes maintained under different conditions
In order to confirm susceptibility to DHBV infection, PDH isolated from Cherry Valley ducklings were maintained in Leibovitz's L-15 Medium supplemented with 5% FBS for 1d (Con-PDH), and then followed by culture in media with 5% FBS (FBS-8d PDH) or 1.5% DMSO (DMSO-8d PDH) for another 7 days (Fig 1A). The PDH under different culture conditions were infected with DHBV at a multiplicity of genome equivalents of 100, and maintained in media with 5% FBS for another 5 days. DHBV in the supernatant of Con-PDH was 3.29×10 5 copies/μl, 7.89×10 4 copies/μl in that of the DMSO-8d PDH, and only 6.72×10 3 copies/μl in that of FBS-8d PDH, determined by qPCR (Fig 1B). DHBV covalently closed circular DNA (cccDNA) and DHBV DNA in the core particles (core DNA) were extracted from the cells and analyzed by Southern blot assay. The highest level of cccDNA and core DNA were detected in Con-PDH, followed by DMSO-8d PDH. The lowest level was detected in FBS-8d PDH (Fig 1C  and 1D). It confirmed that treatment with 1.5% DMSO could maintain the susceptibility of the PDH to DHBV infection [8].

Transcriptome analysis of primary duck hepatocytes under three culture conditions
Transcriptomes of PDH with differential susceptibility to DHBV infection were analyzed using RNA-seq. The sequencing libraries were prepared in triplicates for PDH under three different culture conditions (Con-PDH, DMSO-8d PDH and FBS-8d PDH), and a total of 220.25 million raw reads were generated. After the removal of ambiguous and low-quality reads (Q20<20), 195.74 high-quality reads were obtained and 173.89 million were further mapped onto the Anas platyrhynchos genome (GenBank accession: PRJNA46621) ( Table 1).
Principal component analysis (PCA) was applied to analyze the sample data from PDHs under the different culture conditions (3 samples for each condition). It showed that transcriptome profiles of the samples from the same culture condition were clustered together (Fig 2A). We further compared the differentially expressed gene of PDH under three different culture condition (Con-PDH, DMSO-8d and FBS-8d). A Gene with an FDR (False discovery rate) adjusted P value less than 0.05 (t-tests) and at least 2-fold change in transcript level were deemed to be differentially expressed. By three comparisons between each library pair (FBS-8d vs Con-PDH, DMSO-8d vs Con-PDH and FBS-8d vs DMSO-8d PDH), 4190, 2249 and 2194 genes were identified as DEGs, respectively. A total of 384 DEGs were overlapped, as illustrated by a Venn-diagram ( Fig 2B). Of those overlapping DEGs, expression of 298 genes were up-regulated and that of 86 genes down-regulated (S1 File). The overlapping DEGs were further analyzed using K-mean clustering. The 273 up-regulated and 77 down-regulated DEGs (were assigned to two significant expression patterns as Con-PDH < DMSO-8d < FBS-8d (P value: 5.1×10 −135 , Fisher's exact test) and Con-PDH > DMSO-8d >FBS-8d (P value: 6.1×10 −3 , Fisher's exact test), respectively ( Fig 2C).
Expression of five genes in PDH under different culture conditions were selected for validation by quantitative reverse transcription PCR (RT-qPCR). Among them, two liver-specific expressed genes (alb encoding albumin and ambp for alpha-1-microglobulin/bikunin precursor) were identified as DEGs by RNA-seq, and three other genes (cpd, gldc and furin encoding carboxypeptidase D, glycine dehydrogenase and Furin, respectively) were not detected as DEGs but reported as proteins associated with DHBV infection [14][15][16][17]. Compared with that in FBS-8d PDH, PDH cultured under the culture condition supplemented with DMSO resulted in a 19-fold increase in expression of ALB and a 4-fold increase in AMBP determined by RNA- seq analysis (FPKM value), which was confirmed by RT-qPCR (Table 2). However, in vitro culture of PDH led to dramatic reduction of ALB and AMBP expression, whether or not DMSO was applied. For the DHBV-infection-associated-genes (cpd, gldc and furin), the presence of DMSO did not provide a distinct expression advantage over culture with FBS.

Key DEGs in DHBV infection revealed by a protein-protein interaction network
Hepadnavirus infection in hepatocytes is a complex, multi-step process mediated by diverse proteins. To determine links between the DEGs and hepadnavirus infection, a protein-protein interaction (PPI) network was constructed based on the KEGG database. Due to the limited interaction information for duck proteins, a human protein interaction database was applied to construct the network. The proteins encoded by DEGs were pooled with 113 proteins in hepatitis B (http://www.genome.jp/dbget-bin/get_linkdb?-t+orthology+path:map05161) as candidates of PPI network. The protein-protein interaction network illustrates 65 out of 384 proteins encoded by the DEGs and 48 out of 113 proteins encoded by genes belonging to hepatitis B pathway. In the network, 49 proteins encoded by DEGs showed a direct or indirect interaction with Hepatitis B pathway-associated proteins, and the others interact with each other among the DEGs. Most of the direct interactions were observed between proteins with similar alteration in expression under different culture conditions. Those up-regulated proteins interacted preferentially with proteins belonging to the Hepatitis B pathway. In the PPI network, 47 proteins encoded by DEGs were up-regulated (indicated as red nodes in Fig 4) and 18 ones were down-regulated (as green nodes). The combination of decreased expression of liver-specific genes (CYP3A4, CYP1E1, CFI, RELN and GSTA1 et al, Table 3) with increased expression of hepatocyte-dedifferentiationassociated genes (PLA2G4A and PLCG1) indicates that in vitro culture conditions results in the fading of hepatocyte phenotype in PDHs. The presence of DMSO alleviated the hepatocyte phenotype fading and maintained the hepatocyte characters. The discrepancy of susceptibility of PDH to DHBV under different conditions was consistent with expression pattern of ferritin (FTH1) that was used as a predictor of host response to hepatitis B virus infection [18]. DEGs involving in virus infection were observed in the PPI network including CREB3L1, FDPS and IRS2 (Table 3), in which CREB3L1 was reported to be associated with iCell hepatocyte (induced pluripotent stem cell derived hepatocytes from Cellular Dynamics International) susceptibility to HBV infection [19]. It was unexpected that three genes (RSAD2, CMPK2, and RIPK2) playing roles in anti-viral immune responses (Table 3) showed down-regulated expression in PDHs under in vitro culture.
Expression of seven DEGs associated with tight junction formation, including JAM3, PPP2R2B, PRKAR1B, PPP2R2C, MAGI2, ACTA2 and ACTG2, showed up-regulation pattern, which is in accordance with GO analysis of tight junction signaling pathway (P value: 1.25×10 −3 , Fisher's exact test, Fig 3). Expression of those DEGs in DMSO-8d PDH was lower   than that in FBS-8d PDHs, suggesting that DMSO-8d PDH had less tight junction formation that serves a blockade of HBV infection [13]. DHBV attachment to PDH with or without DMSO was carried out. DHBV were inoculated into PDH maintained under the three different culture conditions, and incubated at 4°C for 2 h. Anti-Pres/s of DHBV was used as control, which could block DHBV attaching to PDH. After washing of free virus, the attached DHBV were detected by qPCR. A decrease of DHBV attachment was observed in the presence of anti-Dpres/s for Con PDH. Attached DHBV in FBS-8d PDH (2.94-fold decrease) were, lower than that in DMSO-8d PDH (1.37-fold decrease), compared with Con PDH (Fig 5).

Discussion
Hepadnaviruses exhibit a strong liver tropism and preferentially replicate in hepatocytes, the predominant liver cell type. Primary cultured hepatocytes remain the most important cell model used to elucidate the molecular pathogenesis underlying HBV infection. However, it remains unclear why primary hepatocytes gradually lose susceptibility to DHBV infection under in vitro culture conditions. Supplementation with DMSO can prolong PDH susceptibility to hepadnavirus [8]. In this study, DEGs were identified in PDH cultured under conditions with or without DMSO with RNA-seq analysis. Our results indicated that multiple factors were affected by DMSO treatment, and differential cell-cell junction may be one of the reasons leading to differential susceptibility in PDH in vitro.
For screening the genes of interest, we compared transcriptome data from PDHs cultured with 1.5% DMSO or 5% FBS, and fresh isolated PDH was used as a control. To identify related DEGs associated with PDH susceptibility to DHBV infection, a Venn diagram was used to visualize overlapping DEGs across the three comparisons (FBS-8d/Con-PDH, FBS-8d/DMSO-8d, DMSO-8d/Con-PDH). Among the DEGs (fold change>2 and FDR adjusted p-value<0.05) observed in the three comparisons, a total of 384 overlapping genes were obtained. Meanwhile, significant changes in 350 overlapping genes (273 up-regulated and 77 down-regulated genes) were enriched in two expression patterns (Con-PDH<DMSO-8d<FBS-8d and Con-PDH>DMSO-8d>FBS-8d), which provided further confirmation of a relationship between the DEGs and PDH susceptibility to DHBV infection. In another point of view, FPKM values of the down-regulated gene, which may be helpful in DHBV infection, are supposed to be higher than that of up-regulated ones. Therefore, FPKM values of the overlapped DEGs in Con PDH were compared by t-test. It indicated FPKM values of down-regulated gene were higher than that of up-regulated ones (p<0.0001, S1 Fig). Differential DHBV attachment to PDH maintained under different conditions. PDH with differential susceptibility was pretreated with nocodazole and infected with DHBV at a MGE of 500. The DHBV attached to PDH were quantified by qPCR. Anti-Dpres/s or IgG were used as control. The asterisk represents a significant difference (t test, P value<0.05). To obtain more protein information, DEGs were matched against the human genome using a BLAST search and 247 uncharacterized proteins were further annotated. Then Proteins encoded by overlapped DEGs and those from hepatitis B pathway were pooled together to construct a protein-protein interaction network. In the PPI network, a cluster of tight junction associated proteins (JAM3, PPP2R2B, PRKAR1B, PPP2R2C, MAGI2, ACTA2 and ACTG2) were up-regulated, suggesting that the formation of tight cell-cell junction between hepatocytes cultured in vitro over time. Since it has been assumed that tight junction formation creates a physical barrier preventing virus access to basolateral localized receptor(s), cell-cell junction alteration may contribute to differential susceptibility to DHBV between PDHs under different conditions. Under the culture condition with DMSO, tight junction formation between PDHs seemed to be retarded, evidenced by lower expression of tight junction genes than that in the absence of DMSO. The result will drive us to focus on molecular pathway of tight junction formation and virus entry in further investigation of influence of DMSO on hepadnavirus infection.
It is admitted that impacts of DMSO involve in more events of virus life cycle. Therefore, we analyzed the DHBV attachment, replication and release in PDH under different culture conditions. We have studied the effect of DMSO on DHBV attachment in PDHs ( Fig 5)  Differential susceptibility of PDH to DHBV also could be due to differential virus defense response to DHBV infection in PDHs under different culture condition. We observed downregulated expression of anti-virus genes (RSAD2, CMPK2, and RIPK2), suggesting the collapse of cellular anti-virus immune concomitant with the loss of hepatic characters in PDH in vitro. It could rule out the possibility that weaken susceptibility of PDH to DHBV is the consequence of enhancement of cell defense against virus infection. In addition, expression of FTH1, a predictor of host response to hepatitis B virus infection in vivo [18] was down-regulated in FBS-8d PDH but partially recovered in DMSO-8d PDH. Expression of FTH1 could be used to predict availability of cell model for DHBV/HBV infection or to screen anti-virus infection chemicals.
We also analyzed the expression of genes which play important roles in the DHBV infection, such as cpd, gldc and furin. CPD is the docking receptor and binds to DHBV particles attached to heparin sulfate and co-factors on the hepatocyte membrane [14]. DHBV particles that bind to CPD are internalized into the cytoplasm, forming endosomes where the large protein on the DHBV membrane is cleaved by the Furin enzyme [15]. The cleaved large DHBV protein then binds to GLDC [16,17], followed by the membrane fusion of the DHBV particle and the endosome. NTCP is considered to be a functional receptor on the hepatocyte membrane and leads to the susceptibility of HepaG2 cell to HBV infection [2]. However, low FPKM value of NTCP was observed in PDH under different conditions (0.068, 0.36 and 0.28 in PDH-Con, DMSO-8d and FBS-8d, respectively) and cpd, gldc, furin and NTCP were not indicated as DEGs both by RNA-seq analysis and RT-qPCR.
Despite the candidate genes revealed in this study, further studies are still needed to investigate the association between these candidate genes and PDH susceptibility to DHBV infection. The mechanisms underlying the contribution of these genes to DHBV infection will be helpful for establishing new strategy against hepadnavirus infection.

Conclusion
In this study, a global comparative transcriptomic analysis has been conducted to investigate differential susceptibility of PDH which were maintained under three conditions. A PPI network was constructed by pooling the DEGs with proteins belonging to hepatitis B pathway. The combination of decreased expression of liver-specific genes with increased expression of hepatocyte-dedifferentiation-associated genes suggested that the fading of hepatocyte phenotype of PDHs were resulted by in vitro culture conditions. At the same time, the expression of seven DEGs associated with tight junction formation was up-regulated, which was attenuated by DMSO. The study paves the way to fully understand mechanism underlying DHBV life cycle and to establish new strategy against hepadnavirus infection.

Ethics statement
All animal experiments with Cherry Valley ducklings that we performed in this study were approved by the Institutional Animal Care and Use Committee (IACUC) of Fudan University (NO.20150305).

Experimental animals and cell lines
The 1-day-old Cherry Valley ducklings (Anas domesticus) were purchased from Breeding Center of Shanghai Institute of Veterinary Medical Sciences, Shanghai, China. The ducklings without of DHBV infection were used for PDH preparation, DHBV in the serum of ducklings were detected by qPCR with the specific primers (based on DHBV genome NC_001344.1), listed in Table 4.
LMH-D2 cells (a generous gift of Dr. William Mason, Institute for Cancer Research, Philadelphia, Pa. USA), were maintained in DMEM F12 with 10% FBS and constitutively secreted DHBV virions into the supernatant [20].

PDH preparation and culture
The PDH were isolated from the livers of 1-day-old cherry valley ducklings as previously described [21]. Briefly, DHBV-negative ducklings were anesthetized with 0.75% Table 4. Primers for detection of DHBV DNA and differential expressed genes. pelltobarbitalum natricum (10 mg/kg) and the portal vein was exposed. The liver was perfused with 80 ml of hepatocyte perfusion medium (Invitrogen, Carlsbad Table 4. The qPCR amplification was performed with the Takara Probe PCR Kit on an Eppendorf RealPlex 4, using 50 cycles of a 2-step PCR program.

DHBV infection and detection of viral cccDNA and core DNA
The PDH were maintained in Leibovitz's L-15 Medium supplemented with 5% FBS for 1 d (Con-PDH), and then followed by culture in media with 5% FBS (FBS-8d PDH) or 1.5% DMSO (DMSO-8d PDH) for 7 days. DHBV infection was carried out at a MGE of 100 by incubating overnight at 37°C. Then the inoculum was removed and L-15 medium supplemented with 5% FBS was added to the cells. On the 5 th day post-DHBV infection, PDH and culture supernatants were collected for isolation of DHBV cccDNA, and core DNA.
To isolate DHBV cccDNA, PDH were lysed in 1 ml of cccDNA extraction buffer containing 50 mM Tris-HCl, pH 8.0, 150 mM NaCl, 10 mM EDTA and 1% SDS. Then 250 μl of 2.5 M KCl was added to the lysate, rotated overnight at 4°C and centrifuged at 18,000 × g for 30 min at 4°C. The DHBV cccDNA was extracted from the supernatant 3 times with phenol and once with chloroform. After adding 2.2 volumes of absolute ethanol and incubating at -80°C overnight, precipitated DHBV cccDNA was centrifuged at 18,000 × g for 20 min, washed 3 times with 75% Ethanol and dissolved in 35 μl of TE buffer.
The extracted DHBV core or cccDNA was separated in a 1.2% agarose gel overnight at 30 v. The DNA transferred onto nitrocellulose membrane via upward capillary was detected by p32-labeled probes which were synthesized with a random-primed DNA Labeling Kit (Roche Diagnostics, Basel, Switzerland); the signals were exposed to a Phosphor screen for 24 h, and visualized by a PhosphorImager (Fujifilm, Kanagawa, Japan).

Detection of DHBV attachment to PDH
The detection of DHBV attachment to PDH was described by Funk A et al [22]. In brief, PDH in 6-well plate were pretreated with 16.5 μM nocodazole for 1 h. The DHBV was inoculated at a MGE of 500 and incubated at 4°C for 2 h. DHBV premixed with anti-Dpres/s monoclonal antibody (1:1000) (a generous gift of Dr. William Mason in Fox Chase Cancer Center) or IgG (1:1000) were used as controls. After incubation with DHBV, the cells were washed with cold PBS for 5 times and lyzed with 200 μl lysis buffer (50 mM KCl, 10 mM Tris-HCl [pH 8.3], 0.45% Tween 20, 0.45% Nonidet P-40, 0.1 mg/ml proteinase K). The lysates were centrifuged at 16100 g for 10 min followed by DHBV DNA extraction with TIANamp Virus DNA/RNA Kit. Then DHBV DNA was quantified by qPCR.

Preparation of the PDH cDNA library, and sequencing
For transcriptome analysis, cDNA libraries were established for the PDH samples. The total RNA was extracted from PDH under three culture conditions (in Leibovitz's L-15 Medium supplemented with 5% FBS for 1 d, followed by culturing in media with 5% FBS or 1.5% DMSO for 7 d) using Trizol reagent following the manufacturer's protocol (Invitrogen, Carlsbad, CA, USA). Then, poly-A mRNA was isolated from 10μg total RNA using a TruSeq 1 RNA Sample Preparation v2 kit (Illumina Inc., San Diego, CA, USA) and was fragmented at 94°C for 8 min. The first and second cDNA strands were synthesized successively, followed by the purification of double-stranded cDNA using AMPure XP beads (Beckman Coulter, Brea, CA, USA). After end repair, cDNA was adenylated at the 3 0 end, ligated to the sequencing adapters and amplified by PCR. The quality and quantity of the cDNA library was determined with an Agilent DNA 1000 Kit and an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). PDH cDNA libraries were sequenced on the Illumina Solexa platform using the pairend approach and the Illumina data processing pipeline.

Processing and mapping of Illumina reads
In order to rule out sequencing errors, all raw reads obtained from the Illumina platform were pre-processed prior to the next step. Clean reads were obtained by removing the adaptor sequences, reads with >5% ambiguous bases (noted as N) and low-quality reads containing more than 20 percent of bases with qualities of <20. High-quality reads were mapped to the Anas platyrhynchos genome and transcripts were downloaded from ENSEMBL (http://www. ensembl.org) with TopHat v2.0.12 software (http://ccb.jhu.edu/software/tophat/index.shtml). TopHat allows multiple alignments per read and the default parameters were used. Cufflinks v2.2.1 software (http://cufflinks.cbcb.umd.edu/) was later used for analyses that included transcript assembly and FPKM value calculations to quantify gene expression; this program was also run with default parameters. Genes with an FDR adjusted P value less than 0.05 and a change in levels of at least 2-fold were considered to be differentially expressed.

Validation of DEGs by real-time PCR
For each RNA sample, 4 μg of total RNA was reverse-transcribed using the TIANScript II RT Kit (Tiangen, Beijing, China). The cDNA was diluted 1:10, and 2 μl of the cDNA was used as template in 20 μl of amplification reaction, using SYBR 1 Premix Ex Taq™ II according to the manufacturer's protocols (Takara, Dalian, China) on a Mastercycler ep realplex (Eppendorf, Hamburg, Germany). The specific primers for each gene were designed using Beacon designer 7.5 software (Palo Alto, CA, USA), and are listed in Table 4.

Canonical pathway analysis by IPA
The differentially expressed genes with corresponding expression values were uploaded to the IPA Software, and analyzed by canonical pathway. We used Fisher's exact test to select the significant pathways and the significance threshold was defined by a P value of 0.05. The significance of the pathway was indicated by the ratio of the number of genes form the dataset that mapped to the pathway to the total number of genes present in the canonical pathway.

Construction of protein-protein interaction network
A protein-protein interaction network was constructed using Cytoscape software. The candidate proteins are encoded by DEGs in the present study and from proteins involved in hepatitis B pathway in KEGG (http://www.genome.jp/dbget-bin/get_linkdb?-t+orthology+path:map05161). A total number of 113 genes which were reported involved in hepatitis B pathway were extracted. Due to the limited interaction information for duck proteins, DEGs identified by RNA-seq analysis were matched against the human genome using BLAST method, and matched ones were pooled together as candidate genes. The protein-protein interaction network was constructed according to the functional relationships annotated in the KEGG database. In the network, lines were used to represent protein-protein interactions including binding/association, phosphorylation, activation, and inhibition. Proteins encoded by up-regulated or down-regulated DEGs were indicated in red or green, respectively, while purple nodes were used to present proteins involved in the hepatitis B pathway.

Statistical analysis
T test was used to filter the differentially expressed genes, after the significant analysis and FDR analysis under the following criteria: i) Fold Change>2 or <0.5; ii) FDR<0.05. Fisher's exact test was used in K-means clustering and Canonical pathway analysis. A P value<0.05 was considered to be statistically significant.