Analyzing Molecular Traits of H9N2 Avian Influenza Virus Isolated from a Same Poultry Farm in West Java Province, Indonesia, in 2017 and 2023

Sample collection

Sampling was conducted at a commercial layer chicken farm in West Java Province, the original location of the AIV subtype H9N2 strain A/Layer/Indonesia/WestJava-04/2017 (H9N2) (GenBank Accession No. MG957203). After five years of vaccination with a homologous commercial vaccine at the farm, a decrease in production and bird mortality occurred in February 2023. Samples consisting of brain, trachea, and oviduct organs from three deceased chickens with symptomatic of H9N2 infection. The samples were from flock of 100,000 chickens at a commercial egg production farm. All organs from all chickens combined in one sample and then was tested by virus isolation and PCR methods.

Virus isolation

The sample was multiplied using standard laboratory procedures. Tissue of pooled organs of three chickens was combined with sterile phosphate-buffered saline (pH 7.4) from ThermoFisher (28348) containing antibiotics (200 g/ml penicillin and 100 g/ml streptomycin) and then centrifuged at 1,000 g and 4 °C for fifteen minutes. The remaining supernatant was filtered through a 0.45 μm membrane filter and injected intra-allantoically into 9-day-old SPF eggs. The SPF eggs were incubated at 37 °C for 48–72 hours, and daily monitoring of embryo mortality was performed.¹⁷ After 72 hours of incubation, all ECEs, regardless of their viability, were subsequently preserved at a refrigerated temperature of 8 °C overnight and then harvested. The collected specimens consisted of allantoic fluid, which was subsequently utilized for conducting a rapid HA assay to promptly detect AIV in early stages.¹⁶

Viral RNA was identified using a reverse transcription-polymerase chain reaction (RT–PCR) assay

The RNA was further extracted from the allantoic fluid that was collected using a total RNA mini kit reagent from Geneaid^®, Taiwan (RPD050). The RNA obtained was then analyzed for the presence of the H9N2 virus using the SensiFASTTM SYBR Lo-ROX Kit from Bioline^®, Taunton (BIO-94005). For the RT–qPCR assay, the primers H9-F: 5′-ATCGGCTGTTAATGGAATGTGTT-3′, H9-R: 5′- TGGGCGTCTTGAATAGGGTAA-3′¹⁸ and N2-F: 5′-CTCCAATAGACCCGTACTAT-3′, N2-R: 5′-CCTGAAGTCCCACAAAATAC-3′¹⁹ were utilized with the LongGene Q2000C (China). The thermal profile for gene amplification included two minutes of polymerase activation at 95 °C, which was followed by a total of 45 cycles of denaturation at 95 °C for 5 seconds, annealing at 60 °C for 10 seconds, and extension at 72 °C for 15 seconds. When the cycle threshold (Ct) value falls below 40, a positive result is indicated. In addition, RT–PCR assays were conducted to detect avian influenza subtype H9N2, Newcastle disease virus (NDV), and infectious bronchitis virus (IBV). RT–PCR was conducted according to the protocol used in a previous study.²⁰

Sequencing

The confirmed HA of the H9 gene from RT–PCR was amplified for sequencing using the MyTaq One-Step RT–PCR kit from Bioline^®, Taunton (BIO-65049), and the primers HAp1-F: 5′-TCCACGGAAACTGTAGACACA-3′, HAp1-R: 5′-TTCTGTGGCTCTCTCCTGAAA-3′ and Hap2-F: 5′-AGGCCTCTTGTCAACGGTTT-3′, Hap2-R: 5′-CCAACGCCCTCTTCACTTTA-3′ were used.²¹ The Sanger sequencing method was performed by First BASE Laboratories, Malaysia, which separated the PCR products using electrophoresis and purified the desired band for sequencing. Using Bioedit v.7 (https://bioedit.software.informer.com/7.0/), the nucleotide and amino acid sequences of a recently isolated strain’s HA gene were determined, and ClustalW was used for alignment. A modern virus phylogenetic tree was constructed using MEGA v 7.0 (https://www.megasoftware.net/), the neighbor-joining method, and 1,000 alignment repetitions. The phylogenetic tree was compared the clade of the recent study isolate, a previous isolate A/Layer/Indonesia/WestJava-04/2017 (H9N2) (GenBank Accession No. MG957203) and other prototype H9N2 isolate strains of the all H9N2 clades from GenBank. The genetic distance between isolates and the topology of the phylogenetic tree were used to compared strains. The genetic distance between isolates and the topology of the phylogenetic tree were used to compared strains. The genetic similarity percentage of all study strains was calculated by using the maximum composite likelihood model implemented in MEGA v 7.0 (https://www.megasoftware.net/).²²

Several key regions were investigated in this investigation, including the receptor-binding sites (RBS) on the left and right edges of the binding pocket, as well as the cleavage site. In addition, the correlation between particular amino acid residues (54, 80, 106, 109, 113, 123, 125, 129, 130, 135, 137, 146, 147, 149, 150, 152, 164, 165, 178, 179, 182, 183, 188, 189, 194, and 216) and H9N2 virus antigenicity was analysed. The amino acid sequence of the HA gene from the most recent H9N2 virus was uploaded to http://www.cbs.dtu.dk/services/NetNGlyc/to evaluate potential N-glycosylation sites.²¹

Molecular Docking of Protein HA1 (2017 and 2023 isolates) with Ligands Neu5Ac2-3Gal and Neu5Ac2-6Gal

The simulation was conducted on a computer hardware system with the following specifications: 8.00 GB RAM and an Intel^® Core i5-2520 M CPU with a clock speed of 2.50 GHz. The software used in the study included Google Chrome v114.0.5735.199 for accessing the SwissDock website, which served as the protein–ligand docking server and ran automatically (http://swissdock.ch/docking). Additionally, BIOVIA Discovery Studio Visualizer v21.1.0.20298 (https://discover.3ds.com), Chimera v1.17.3 (https://www.cgl.ucsf.edu/chimera/download.html), AutoDock Tools v1.5.7 (https://autodock.scripps.edu/), PyMol v2.5.2 (https://pymol.org/), and LigPlot+ v2.2.8 (https://www.ebi.ac.uk/thornton-srv/software/LigPlus/download.html) were used for data analysis and visualization.

The simulation begins with the creation of a 3D model of the HA1 protein from the AIV, specifically the recent study isolate (HA1-2023) and the A/Layer/Indonesia/WestJava-04/2017 strain (HA1-2017). This is achieved using the PHYRE2 Protein Fold Recognition Server (http://www.sbg.bio.ic.ac.uk/). The ligands used in the simulation were obtained from https://pubchem.ncbi.nlm.nih.gov/. They are Neu5Ac2-3Gal (compound CID: 13832708) and Neu5Ac2-6Gal (compound CID: 53262334).

The docking was performed using the automated server facility on SwissDock, where the uploaded documents were in PDB format (.pdb) for the protein acting as the receptor and SYBYL MOL2 format (.mol2) for the molecules acting as ligands. Docking was carried out for the protein HA1-2017 with the Neu5Ac2-3Gal ligand, HA1-2017 with the Neu5Ac2-6Gal ligand, HA1-2023 with the Neu5Ac2-3Gal ligand, and HA1-2023 with the Neu5Ac2-6Gal ligand. Only models with ligands attached to the Leu216 residue region are observed and selected for further analysis.

The selection of the residue is based on the fact that the amino acid at position 216 determines the virus’s tendency to infect mammalian cells dominated by Neu5Ac2-6Gal or avian cells dominated by Neu5Ac2-3Gal.²⁰ The docking results were downloaded in the Chimera Web Data format (.xml) and run in Chimera software v1.17.3 to select the best docking model. The best model was chosen based on the Gibbs free energy (ΔG) value, where ΔG represents the stability parameter of the bond formed between the protein and the ligand. A more negative value indicates a more stable bond. The protein–ligand interactions are then evaluated, observing hydrogen bonding, hydrogen bond distances (Å), interacting residues and functional groups, and residues that interact with the ligand noncovalently.²³

Nucleotide sequence accession numbers

The partial CDS of the HA gene of A/chicken/Indonesia/MSL0123/2023 (H9N2) in this study was deposited in GenBank and received the accession number OR243721.

H9N2 subtype AI virus detection

The presence of viral RNA was identified by RT–qPCR specifically targeting H9 and N2. The results shown in Supplementary Figure 1a and Supplementary Figure 1b indicate that the recent isolate of this study is an AIV subtype H9N2. The H9 gene in a recent isolate was amplified with a Ct value of 28.50, while the N2 gene was amplified at a Ct value of 34.69.

Amplification of the Hemagglutinin (HA) gene

The results of HA gene amplification using primers HApar1 and HApar2 on a recent sample from this investigation are shown in Supplementary Figure 1c. The electrophoresis results revealed the presence of DNA bands at the 736 bp (for HApar1) and 712 bp (for HApar2) sites. Marker 100-2000 bp is the size of the molecular marker displayed in the illustration.

Amino acid analysis

The partial CDS of the HA gene from this study was deposited in GenBank, with the accession number OR243721 as A/chicken/Indonesia/MSL0123/2023 (H9N2). The amino acid sequence at the cleavage site of the recent isolate is PSRSSR↓GLF, while the receptor binding site has the motif PWTNTLY ( Table 1). Additionally, at position 217 on the left side of the receptor binding site (RBS), the recent isolate contained the amino acid M (methionine). It should be noted that this amino acid sequence is identical to A/Layer/Indonesia/WestJava-04/2017, which originated from the same farm.

In a recent study, the antigenic site I of the HA gene in the recent sample displayed the same motif as A/Layer/Indonesia/WestJava-04/2017 (H9N2). At site I, amino acids S (serine), K (lysine), and P (proline) were present. For site II, the motif consisted of amino acids D (aspartate), N (asparagine), and L (leucine) at positions 135, 183, and 216. In contrast, A/Layer/Indonesia/WestJava-04/2017 (H9N2) exhibited D (aspartate) at position 183. The analysis of the HA gene in the recent sample included an examination of potential glycosylation sites (PGS). The HA1 gene segment revealed the NXT/S motif (where X represents any amino acid except proline) at positions 11-13 (NST), 123-125 (NVS), 200-202 (NRT), 280-282 (NTT), 287-289 (NVS), and 295-297 (NCS). There were also several other amino acid differences observed at positions 23, 34, 53, 69, 72, 74, 114, 120, 163, and 179 ( Table 1).

Homology comparison

The results of the homology comparison are shown in Table 2. The recent isolate has a similarity of 95.84% with A/Layer/Indonesia/WestJava-04/2017 (H9N2), which was isolated from the same farm but in a different year. When analyzed against several representatives of the H9N2 subtype clades, the recent isolate has a similarity of approximately 72.24% with clade h9.1, 72.56% with clade h9.2, 74.72% with clade h9.3, and 81.80% with h9.4.1. These results are lower compared to the similarity of the recent isolate with viruses classified under clade h9.4.2. Based on homology comparison, a recent isolate, A/chicken/Indonesia/MSL0123/2023 (H9N2), belonged to subclade h9.4.2.5.

Phylogenetic tree

The phylogenetic tree showed that the samples were together with other isolates from Indonesia such as A/Layer/Indonesia/WestJava-04/2017 (H9N2), A/chicken/East_Java/M92_10/2017 (H9N2), A/chicken/East_Java/M92_24/2017 (H9N2), Vietnam A/muscovy_duck/Vietnam/LBM719/2014 (H9N2) and China such as A/chicken/Zhejiang/HE6/2009 (H9N2), A/chicken/Henan/LY-36/2013 (H9N2), A/chicken/Guangdong/LGQ02/2014 (H9N2) so that these were included in China, Vietnam, and Indonesia (CVI) clades ( Figure 1). Recent viruses belonged to subclade h9.4.2.5, according to a phylogenetic study.

Figure 1. Phylogenetic tree of the hemagglutinin (HA) gene of the AIV H9N2 virus was constructed.

Docking analysis

The target HA1 proteins of AIV subtype H9N2 isolated in 2017 (HA12017) and 2023 (HA12023) were chosen for the study because of their importance in viral attachment for the docking analysis, as there were no 3D structures of the protein, and its modeling was carried out. After submission of their amino acid sequences on Phyre2, the modeled proteins obtained were downloaded and visualized in Discovery Studio Visualizer. The Neu5Ac2-3Gal (Sia3) and Neu5Ac2-6Gal (Sia6) ligands used in this study are shown in Figure 2.

Figure 2. Modeled protein HA1 of H9N2 predicted by the Robetta server, visualized in Discovery Studio Visualizer.

Based on the docking results, the data obtained are presented in Table 3. The obtained results represent the best model with the most negative ∆G and the ligand positioned and bound near the Leu216 residue. The interactions between each protein and ligand are shown in Figure 3, with hydrogen bonding data and residues involved in hydrogen bonding interactions presented in Table 4. The interactions that occur can be either hydrogen bonding or nonbonding interactions that strengthen the ligand’s affinity with the protein.

Figure 3. Interactions occurring in HA12017 and HA12023 with ligand Sia3 and Sia6.

Ethical statement

The present research adheres to the guidelines outlined in the Indonesian Law on Animal Health Research (UU/18/2009, article 80). Due to the absence of live animals in this study, ethical approval was not necessary.