Cytokine-cytokine receptor interactions in the highly pathogenic avian influenza H5N1 virus-infected lungs of genetically disparate Ri chicken lines

Objective The highly pathogenic avian influenza virus (HPAIV) is a threat to the poultry industry as well as the economy and remains a potential source of pandemic infection in humans. Antiviral genes are considered a potential factor for HPAIV resistance. Therefore, in this study, we investigated gene expression related to cytokine-cytokine receptor interactions by comparing resistant and susceptible Ri chicken lines for avian influenza virus infection. Methods Ri chickens of resistant (Mx/A; BF2/B21) and susceptible (Mx/G; BF2/B13) lines were selected by genotyping the Mx dynamin like GTPase (Mx) and major histocompatibility complex class I antigen BF2 genes. These chickens were then infected with influenza A virus subtype H5N1, and their lung tissues were collected for RNA sequencing. Results In total, 972 differentially expressed genes (DEGs) were observed between resistant and susceptible Ri chickens, according to the gene ontology and Kyoto encyclopedia of genes and genomes pathways. In particular, DEGs associated with cytokine-cytokine receptor interactions were most abundant. The expression levels of cytokines (interleukin-1β [IL-1β], IL-6, IL-8, and IL-18), chemokines (C-C Motif chemokine ligand 4 [CCL4] and CCL17), interferons (IFN-γ), and IFN-stimulated genes (Mx1, CCL19, 2′-5′-oligoadenylate synthase-like, and protein kinase R) were higher in H5N1-resistant chickens than in H5N1-susceptible chickens. Conclusion Resistant chickens show stronger immune responses and antiviral activity (cytokines, chemokines, and IFN-stimulated genes) than those of susceptible chickens against HPAIV infection.


INTRODUCTION
Avian influenza (AI) is a virus of the genus influenzavirus A in the Orthomyxoviridae family of viruses [1]. Based on genetics and disease severity, AI viruses are classified either as highly pathogenic avian influenza (HPAI) or lowly pathogenic AI [2]. In particular, a high number of deaths were observed in birds caused by HPAI [3]. Moreover, the influenza A virus subtype H5N1, type of HPAI virus (HPAIV), is a threat to the poultry industry and the economy and remains a potential source of pandemic infection in humans [4].
Mx dynamin like GTPase (Mx) proteins are members of the dynamin family of large GTPases. They prevent viral RNA replication by inhibiting the trafficking or activity of viral polymerases [5]. Several previous reports showed that only the asparagine (Asn-AAT) polymorphism at the 631st position triggered antiviral activity, while Mx proteins carrying a serine (Ser-AGT) at that position did not suppress viral growth [6]. In addition, the major histocompatibility complex (MHC) haplotype can affect the antiviral activity of the host [7]. Previous research has shown a significant association of the major histocompatibility complex class I antigen BF2 (BF2) haplotype in chicken MHC class I with resistance or susceptibility to a number of pathogens, including Marek virus [8], and AI virus [9]. Furthermore, chickens with the BF2-B13 haplotype had higher mortality than those that have the BF2-B21 haplotype [9].
One of the key determinants of the severity and outcome of AI virus infection is the regulation of the host innate immune response [10]. Cytokines and chemokines play crucial roles in the balance of the immune system. Previous studies in various animals have shown that tissue damage and host death are caused by cytokine and chemokine dysregulation [11]. In chickens, HPAIV H5N1 induced the excessive expression of cytokines and chemokines in lung tissues [12]. Furthermore, cytokine-cytokine receptor interactions were significantly increased in Fayoumi and Leghorn chicken lines after infection with HPAIV [13]. The regulation of cytokines and chemokines has been shown to be important in host defense against HPAIV.
In this study, we used Ri chickens, a local chicken breed in Vietnam, as an experimental animal [14]. Chickens resistant or susceptible to HPAIV were distinguished by genotyping their Mx and BF2 genes. We infected Ri chickens with HPAIV H5N1 and analyzed gene expression patterns in the lung tissue using high-throughput RNA sequencing. In particular, we analyzed the expression of genes related to cytokine-cytokine receptor interactions between resistant and susceptible Ri chickens.

Experimental animals and HPAIV infection
Fourteen-day-old post-hatching Ri chickens were kept in specific-pathogen-free conditions and observed daily for signs of disease and mortality. The care and experimental use of the chickens were approved by the Ministry of Agriculture and Rural Development of Vietnam (TCVN 8402: 2010/TCVN 8400-26:2014).
For the Mx gene, Ri chickens that have the non-synonymous adenine single nucleotide polymorphism (SNP) at residue 631 were genotyped as resistant while guanine was genotyped as susceptible (Supplementary Figure S1). In addition, for the BF2 gene, those that had the B21 genotype were determined to be resistant and those that had the B13 genotyped were susceptible. Taking these together, HPAIVresistant Ri chickens had the genotype Mx(A)/BF2 (B21) and HPAIV-susceptible Ri chickens had the genotype Mx(G)/ BF2(B13). For HPAIV infection, a total of 10 Ri chickens (five resistant and five susceptible) were inoculated intranasally with the collected allantoic fluid including 10 4 egg infectious dose (EID50) of A/duck/Vietnam/QB1207/2012 (H5N1) based on the OIE guidelines [15]. All experiments, including chicken management, HPAIV infection, and sample collection, were conducted at the National Institute of Veterinary Research, Vietnam.

Sample collection and total RNA extraction
At three days post-infection, the lung tissues were collected from five chickens per group, according to the WHO Manual on Animal Influenza Diagnosis and Surveillance. Total RNA from lung tissue was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer' s instructions.

High-throughput RNA sequencing and data analysis
Potentially existing sequencing adapters and low-quality bases in the raw reads were trimmed using Skewer version 0.2.2. The cleaned high-quality reads after trimming the lowquality bases and sequencing adapters were mapped to the reference genome using STAR software version 2.5. Since the sequencing libraries were prepared in a strand-specific manner using the Illumina strand-specific library preparation kit (Illumina Way, San Diego, CA, USA), the strandspecific library option was applied in the mapping process. The gene annotation of the chicken reference genome gg6 from the UCSC genome (https://genome.ucsc.edu) and the corresponding expression values were calculated as fragments per kilobase of transcript per million fragments mapped (FPKM). The differentially expressed genes (DEGs) between the two selected biological conditions were analyzed using the Cuffdiff software in the Cufflinks package version 2.2.1. Assignment to significance was determined using a threshold. The default threshold was set as |fold-change| ≥2 and p-val-ue≤0.05. To gain insight into the biological roles of the genes differentially expressed between the susceptible and resistant lines, a gene set overlapping test was performed between the analyzed DEGs and functionally categorized genes, including the biological processes determined from the gene ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG) pathways, and other functional gene sets using g: Profiler version 0.6.7.

Cytokine-cytokine receptor interaction enrichment analysis
KEGG pathway analysis of the DEGs revealed that most DEGs were related to cytokine-cytokine receptor interactions. Therefore, cytokine-cytokine receptor interactions were analyzed in this study. We used the Search and Color Pathway tool of KEGG Mapper (https://www.genome.jp/ kegg/tool/map_pathway2.html) to identify the upregulated or downregulated genes between resistant and susceptible Ri chicken lines. The mutual interactions of the cytokine sig-naling pathway-related genes were analyzed using the Pathway Studio program (Ariadne Genomics, Rockville, MD, USA).

Quantitative real-time polymerase chain reaction validation for cytokine-cytokine receptor interaction genes
Quantitative real-time polymerase chain reaction (qRT-PCR) was performed to verify the DEGs obtained from RNA sequencing. Before cDNA synthesis, 2 μg of total RNA (R1D3, R2D3, R3D3, and R5D3 for resistant line and S1D3, S3D3, and S5D3 for susceptible line) was treated with 1.0 U DNase I (Thermo Fisher Scientific, Waltham, MA, USA) to remove potentially contaminating genomic DNA. cDNA synthesis was performed using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, USA) according to the manufacturer's instructions. The primer sequences of 15 genes and the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were designed using Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) ( Table 1). The qRT-PCR was performed in a LightCycler 96 system (Roche, Indianapolis, IN, USA) using the AMPI-GENE qPCR Green Mix Hi-Rox kit (Enzo Life Sciences, Inc., Farmingdale, NY, USA) according to the manufacturer's recommendations. Relative gene expression was calculated using the 2 -ΔΔCt method after normalization with chicken GAPDH [16]. All experiments were performed independently in triplicates.

Statistical analysis
Statistical analysis was carried out using SPSS software (IBM, SPSS 26.0 for Windows, Chicago, IL, USA), and p<0.05 was considered as a statistically significant expression level using Student's t-test. All qRT-PCR experiments were replicated independently three times, and the mean±standard error of the mean values for each group were reported.

Analysis of RNA sequencing
After H5N1 infection, we observed that chickens had ruffled hair; pulmonary edema, pneumorrhagia, and congestive lung. Among the lung tissue samples obtained from resistant and susceptible Ri chickens three days post H5N1 infection, six samples passed the quality control check and were subject high-throughput RNA sequencing (Table 2). Table 2 shows the statistics of the raw and clean reads of individual sample transcriptomes after sequence processing and analysis. The six libraries produced approximately 6 GB worth of cDNA sequences per chicken. After data filtering, approximately 21 million clean reads were obtained for each sample from the resistant and susceptible lines.
The number of mapped reads, percentages, and transcripts is shown in Table 3. More than 88.5% of the filtered reads from each library were mapped to the reference genome.

Gene ontology functional enrichment and KEGG pathway analysis
The calculated FPKM values were used to compare the mRNA expression levels between resistant and susceptible Ri chickens ( Figure 1). Among the 972 DEGs, 309 genes (32%) were significantly upregulated and 663 genes (68%) were downregulated in susceptible chickens compared to the resistant chicken lines (Supplementary Table S1).
To explore the function of these DEGs, we performed GO and KEGG pathway analyses. Among the 972 DEGs, 114 were associated with biological processes, seven were associated with molecular functions, and three were associated with cellular components (Figure 2; Supplementary Table S2). KEGG pathway analysis revealed six categories of immunerelated pathways, including cytokine-cytokine receptor interactions (33 DEGs), nucleotide-binding oligomerization domain like receptor (NOD-like receptor) signaling pathway (27 DEGs), influenza A (26 DEGs), Toll-like receptor signaling pathway (23 DEGs), herpes simplex infection (20 DEGs), and RIG-I-like receptor signaling pathway (11 DEGs) (Figure 3). Most of the DEGs were related to cytokine-cytokine receptor interactions.

Analysis of cytokine-cytokine receptor interactions related to the DEGs
Based on the RNA sequencing results and subsequent mapping to the KEGG pathway database, 33 DEGs were involved in cytokine-cytokine receptor interactions (Supplementary  Table S3). Among them, 32 genes were significantly downregulated, and one gene was significantly upregulated in susceptible chickens compared to the resistant chicken lines (Figure 4).
For the interaction analysis, 32 DEGs out of 33 genes were identified using the Pathway Studio program ( Figure 5). Of these, 28 genes showed strong interactions with each other, but four genes showed no interaction.

Quantitative real-time polymerase chain reaction analysis of genes associated with cytokine-cytokine receptor interactions
To validate the RNA sequencing results, the expression levels of 14 cytokines, chemokines, and receptors of the 33 DEGs between the two chicken lines were quantitatively determined via qRT-PCR ( Figure 6). After H5N1 infection, the expression levels of interleukin-1β (IL-1β), IL-6, IL-18, C-C       infected with HPAIV H5N1. RNA sequencing was conducted after infection and 972 DEGs were identified after comparing the transcriptome profiles of lung tissue obtained between resistant and susceptible chickens. KEGG analysis revealed that most of the DEGs were related to cytokine-cytokine receptor interactions.
Avian influenza viral pathogen-associated molecular patterns are recognized by host pattern recognition receptors (PRRs). Viral double-stranded RNA (dsRNA) and cytosineguanosine oligodeoxynucleotides, which form during the replication of AIV, are recognized by toll-like receptor 3 (TLR3) [17] and TLR21 [18], respectively, through adaptor proteins such as TIR-domain-containing adapter-inducing interferon and myeloid differentiation primary response 88 (MyD88), respectively [19]. These adaptors activate the transcription factor interferon regulatory factor 7 (IRF7) and nuclear factor kappa B (NF-κB) to induce cytokines, chemokines, and IFN-stimulated genes [20]. High expression levels of PRRs, MyD88, IRF7, and NF-κB were previously observed in H5N1-infected chickens [12]. Our results showed that the expression levels of TLR3, TLR21, IRF7, and MyD88 were higher in resistant chickens compared to susceptible chickens (Supplementary Table S1); therefore, we suggest that resistant Ri chickens respond more sensitively to HPAIV H5N1 infection than susceptible Ri chicken lines through increased expression of TLR3, TLR21, IRF7, and MyD88.
Pro-inflammatory cytokines cause inflammation and recruit other immune cells to the site of infection [21]. High expression levels of pro-inflammatory cytokines and chemokines have been found in humans and duck infected with H5N1 [10,11]. Similarly, high expression levels of pro-inflammatory cytokines and chemokines, including IL-1β, IL-6, IL-8, IL-18, CCL4, CCL17, and IFN-γ, were reported in chickens infected with H5N1 [12]. A significant increase in the levels of cytokines and chemokines exacerbates the inflammatory response, leading to apoptosis, multi-organ failure, and host death [22]. However, HPAIV infection was lethal in mice lacking tumor necrosis factor and IL-1 receptors and inhibition of the cytokine response does not protect mice against H5N1 influenza infection [23]. In this study, expression levels of pro-inflammatory cytokines and chemokines such as IL-1β, IL-6, IL-8, IL-18, CCL4, CCL17, IFN-γ, and receptors were upregulated in resistant chickens, compared to susceptible chickens, after HPAIV H5N1 infection ( Figure 6). However, the induction of cytokines and chemokines did not cause Ri chickens to die during experimentation. Therefore, we suggest that cytokine and chemokine induction is necessary to protect the host against HPAIV infection.
Type I interferons (IFN-α and IFN-β) trigger the expression of IFN-stimulated genes. IFNs and IFN-stimulated genes can inhibit viral replication by blocking virus entry into the host cells, binding to viral RNA to stop translation, and regulating host antiviral responses [24]. IFN-γ directly or indirectly inhibits viral replication by strictly regulating the production of nitric oxide or interfering with the onset of the RNase L pathway [25]. Moreover, several studies have shown that IFN-stimulated genes have antiviral activity [26,27]. The Mx gene inhibits the trafficking and activity of viral polymerases [5]. Viperin (RSAD2) inhibits newly synthesized influenza virions [28]. The CCL19, CCL21, and CCR7 chemokine axes are homeostatic [29]. Chicken interferon-inducible 2'-5'-oligoadenylate synthase-like (OASL) and RNase L restrict both viral and cellular RNA, preventing viral genome replication [30]. Wild-type duck OASL inhibits the replication of a variety of RNA viruses in vitro, including the influenza virus [25]. Protein kinase R (EIF2AK2) inhibits the translation of viral mRNAs, including those from influenza A viruses [26]. Interferon-induced proteins of the tetratricopeptide repeats (IFIT) protein family sequester viral nucleic acids [31]. Furthermore, the clinical results were ameliorated in chIFIT5-transgenic chickens after treatment with HPAIV and Newcastle disease virus [27]. As our results showed, although IFN-α and IFN-β were not found in the DEGs, we observed a higher expression of IFN-γ and IFNstimulated genes (Mx, CCL19, OASL, RSAD2, EIF2AK2, and IFIT5) in the resistant Ri chickens compared to the susceptible ones (Supplementary Table S1). Therefore, we suggest that the antiviral response of resistant chickens is higher than that of susceptible chickens.
In summary, by using RNA sequencing and qRT-PCR, we evaluated the differential expression of genes associated with cytokine-cytokine receptor interactions from the lung tissue of resistant and susceptible H5N1-infected Ri chickens. Interestingly, the expression of PRRs, cytokines, chemokines, and IFN-stimulated genes in resistant chickens were higher than those in susceptible chickens. These results suggest that resistant Ri chickens show a higher antiviral activity compared to susceptible Ri chickens, and this antiviral activity may be related to the expression of antiviral genes. Therefore, we propose that more studies on Ri chickens be done to compare their resistance and susceptibility to HPAIV infections. With this, there will be more studies to support advancing the development of chicken lines with disease resistance genes.