بررسی توالی کامل ژن های فسفوپروتئین و ماتریکس ویروس نیوکاسل جدا شده از مزارع مرغ شمال شرق ایران

Introduction

Newcastle disease (ND) is an avian viral infection that causes substantial economic losses in the poultry industry worldwide (Yune & Abdela, 2017). Newcastle disease virus (NDV), a.k.a. APMV-1, is a member of the Paramyxoviridae family, genus Orthoavulavirus (Yune & Abdela, 2017). The genome of NDV codes for six main viral proteins identified as NP, P, M, F, HN, and L (Flint et al., 2020).

F protein plays a pivotal role in the pathogenicity of NDV (OIE, 2012; Wang et al., 2017^a). This protein maintains various determinants such as neutralizing epitopes, hypervariable regions, and fusion peptides that have essential roles in the establishment of the viral infection (Sergel-Germano et al., 1994; Yusoff et al., 1989; Wang et al., 2016; Wang et al., 2017^b; Thi Huong et al., 2019).

The gene coding for the polycistronic phosphoprotein (P) is the most variable gene of NDV. Some genomic variations in this region help NDV escape from the host's immune response in previously vaccinated birds (Kattenbelt et al., 2006).

On the other hand, the gene coding for M protein endures minimum sequence variations, suggesting it is a valuable target for virus typing (Seal et al., 2000). However, evidence shows that the fundamental region in M protein plays characteristic roles in the pathogenicity of the virus (Coleman & Peeples, 1993; Kattenbelt et al., 2006; Duan et al., 2014^a; Duan et al., 2014^b;).

There is no effective treatment available for ND. Therefore, measures to prevent the disease's occurrence and outbreak, such as biosecurity approaches and vaccination, are the most effective approaches in keeping the disease under control and minimizing the economic losses (Yune & Abdela, 2017). However, the cross-protection from vaccination between the different genotypes of NDV is minimal (Miller et al., 2007). Accordingly, up-to-date knowledge of the NDVs prevalent and/or evolving in a region as well as maintaining efficient models to predict the possible introduction of non-native NDVs into the region is critical in this regard.

The Northeast region of Iran is a major producer and exporter of poultry products to other countries (Mohaddes, 2009). Common and rather long borders with countries with poor biosecurity measures necessitate consistent and dynamic surveillance of ND in this geographic region. The genomic analyses of NDV in Iran have so far been focused on regions harboring HN and cleavage site of F gene, and studies on other important genomic regions such as those coding for P and M proteins are limited. Here, we report the full genomic sequences of the P and M loci as well as that of the coding region (419 base pairs) of the F gene in a viral isolate obtained from an outbreak of ND in a broiler farm in the Northeast of Iran and analyze it. Further, phylogenetic and epidemiological relationships between the studied isolate and previously reported strains worldwide were surveyed. Our multiple alignments and phylogenetic analyses identified this virus of velogenic VII.1.1 genotype grouped with those reported from China.

Materials and Methods

History

An outbreak of Newcastle disease was reported in the Northeast of Iran in 2011. Sick chickens showed clinical manifestations of virulent ND. Briefly, the affected vaccinated broiler chickens displayed high mortality (between 40-80%) and severe symptoms of torticollis, lateral twisting of the head and neck, gasping and respiratory distress, green bile pigment, and white urates in faces. Vaccinated layer flocks also showed nervous symptoms, egg drop production, and deformed-shell eggs up to 70%. Gross pathology revealed cecal tonsils necrosis accompanied by proventriculus multifocal hemorrhages.

Virus Propagation and Primary Identification

Ten brain samples from a broiler flock with a high mortality rate (40% in 2 weeks), severe nervous clinical symptoms, and enteric lesions were collected for virus isolation. The samples were transferred on ice to the poultry clinic of the School of Veterinary Medicine, the Ferdowsi University of Mashhad, where they were processed for in vivo pathogenicity verification and virus propagation according to the standard guideline (OIE, 2012). The procedures outlined by OIE are used to verify the NDV virulence in vivo (OIE, 2012). For this, criteria such as mean death time (MDT) in 9-day-old embryonated specific-pathogen-free (SPF) eggs and intracerebral pathogenicity index (ICPI) in 1-day-old chicks (Gallus gallus) were considered. Accordingly, the sam-ples were homogenized, pooled, and used to inoculate allantoic fluids of 10-day-old Specific Pathogen Free (SPF) embryonated chicken eggs (Razi institute, Karaj) (OIE, 2012). Three to four days post-inoculation, embryos were found dead with severe and petechial hemorrhages. The allantoic fluids were harvested, and viral infection was confirmed by HA assay, as described previously (OIE, 2012).

Primer Design

Three and two sets of primers were designed to amplify P and M genes, respectively (Table 1) using Primer 3 Web (version 4.1.0) (Koressaar et al., 2018). They were designed to generate overlapping amplicons to have the sequence information covering full-length genes. Primers reported by Aldous et al. (2003) were used to amplify the cleavage site of the F gene.

Table 1. Specifications of the primers used to amplify P and M genes.

RNA Extraction and RT-PCR

RNA isolation kit-III (Dena zist, Iran) was used for RNA extraction from Allantoic fluids. The RNA was then reverse transcribed into cDNA using AcuuPower®CycleScript RT PreMix kit (Bioneer Corp-oration, South Korea). The cDNAs were used to amplify the genes of interest (P, M, or F genes) using AccuPower®PCR PreMix kit (Bioneer Corporation, South Korea). The PCR conditions were as follows: 95°C/ 4min, 25x (95 °C/30 s, 58/40 s for M gene fragment 1, 55/45 s for F gene, and 60°C/40 s for all the other amplifications, 72°C/1 min), 72°C/10 min.

Electrophoresis and Sequencing

The PCR products were gel separated using agarose 1.2%, and the amplified products were excised and extracted using S-1050-1 kit (Dena zist, Iran). The purified products were sequenced in both directions using Sanger sequencing by Macrogen Com-pany (South Korea).

Bioinformatics Analyses

The sequenced data were evaluated, and the corresponding full gene sequences were reconstructed using CLC Main Workbench (version 5.5) (CLC Genomics Workbench). The sequences were confirmed to be that of the corresponding genes, i.e., P, M, or F genes, using the BLAST function of the National Center for Biotechnology Information (NCBI) database.

Nucleotide sequences were deposited in the GeneBank database with the accession numbers MG017-442, MG017443, and JQ344320 for M, P, and partial sequence of F genes, respectively (The accession number of complete cds of F gene is available on request).

MEGA 7 (version 7.0.26) program (Kumar et al., 2016), set at Maximum Likelihood, K2+G model, and 1000 bootstrapping, was used to construct the phylogenetic trees. Initially, a 419-bp nucleotide sequence of the F gene obtained for the studied isolate and the related sequences for other 47 ND strains previously reported in Iran was used to generate the phylogenetic tree. Subsequently, other phylogenetic trees were generated based on our isolate's P, M, and F gene sequences in conjunction with those of more than 100 other NDV strains classified by Dimitrov et al. (2019). Three class I NDVs with large distances were used as outgroups. It is noteworthy that the complete sequence of F was analyzed. However, the related data are not shown in this study.

The nucleotides and their translated amino acid sequences aligned by Clustal X (version 2.1) (Higgins & Sharp, 1988) were analyzed using the bioinformatics programs BioEdit (version 7.2.6) (Hall, 1999) and MEGA7 (version 7.0.26) (Hall, 1999) and MEGA7 (version 7.0.26) (Kumar et al., 2016). Applications were used to conduct similarity and distance evaluation, estimate nucleotide substitution rate, and compare amino acid composition.

Results

Primary Virus Identification and Pathogenicity Evaluation

Based on clinical findings indicative of ND, brain samples belonging to 10 chickens from one broiler farm were pooled and used for further molecular characterization (see Material and Methods). We subjected the isolate to MDT and ICPI assays, which were determined to be <60 h and >1.5, respectively. This, plus the F protein amino acid composition (Table 2), showed that our isolate was a velogenic strain. The studied isolate was used to inoculate the allantoic fluid of 10-day-old SPF chicken eggs. Embryos in the inoculated eggs died with severe petechial hemorrhages 3-4 days postinoculation. HA assay on the allantoic fluids of the inoculated embryos showed an antigen titer of 6 HA units per 25 μL.

Table 2. The amino acid substitutions in fusion peptide and cleavage site of F protein in the studied isolate and some other similar sequences submitted in GenBank. Isolates are represented by the nucleotide and protein GenBank accession numbers host name, country of isolation, and year of isolation. The isolate subject of the current study is shown in bold.

^a The consensus sequence was obtained from 100 sequences represented in GenBank 

^b Identical amino acid as the consensus sequence

Homology Analyses

BLAST function of the NCBI database was used to run pairwise alignment analyses. For this purpose, we first used the coding sequence of the F gene (419 bps), and this analysis revealed high homology (98.33-99.76%) of our isolate to other VII.1.1 genotypes isolated between 2011 and 2019 in Iran (Table 3). Further, the VII.1.1 genotypes of Chinese origin were the closest foreign strains to our isolate (Table 3). However, two VII.2 isolates (MG871466 and MT254060) recently reported in 2017 and 2019 from Iran were closer to variants reported from Pakistan (Table 3) than our isolate.

The coding and non-coding regions of the P, M, and F genes were used in the second run. As shown in Table 4, this analysis produced similar data to the previous one. Furthermore, the analyses of the full sequence of F showed similar patterns (data are not presented in this study).

Distance Analyses

Distance analyses using the coding sequence (419 bps) of the F gene showed the same trends revealed by the homology analyses. The nucleotide distances between our isolate and other VII.1.1 sequences reported from Iran between 2011 and 2019 were less than 0.05 (Table 3). This value was higher when VII.1.1 isolates from other countries such as Russia, Japan, South Korea, Vietnam, and Taiwan were considered. However, it remained less than 0.05 for Chinese strains (Table 3). As expected, sequences derived from other VII sub-genotypes (VII.2 and VII.1.2) belonging to Iran and other countries were significantly distant from our isolate (Table 3).

Table 3. Identities (%) and distances of the 419-bp coding sequences of the F gene of the studied isolate with some other similar VII genotype isolates. The codon positions included were 1st+2nd+3rd+Noncoding. The sequences with the highest level of identity to the studied isolate are shown in bold.

As is shown in Table 4, distance analyses based on the coding and non-coding region of the P, M, and F genes were consistent with the data generated using the coding region of the F gene alone.

Phylogenetic Analyses

Phylogenetic analysis using sequences of F (Figures 1 and 2), P (Figure 3), and M (Figure 4) genes showed that our isolate is clustering with sub-genotype VII.1.1 of NDVs reported from various countries from 2000-2019; and especially with those reported in China and Iran. Our isolate was far from vaccine strains (KU665482 and KU886038) and other NDV genotypes (VII.2, XIII, XXI) previously reported from Iran.

Figure 1. Phylogenetic tree, based on the nucleotide sequences of 419-bp coding region of the F gene of the studied isolate, compared to the related sequences from some other isolates reported from Iran available until 2019. The maximum likelihood method and K2+G model with 1000 Bootstrap replicates were used. The taxa names include the GenBank accession number (Dimitrov et al. classification/old classification), hostname, country of isolation, strain designation, and year of isolation. The isolated subject of the current study is marked by a circle. Class I is assigned as an outgroup.

Table 4. Identities (%) and distances based on the genomic region of P and M and the 419-bp coding sequences of the F gene. The isolate reported in the present work is compared to some other VII genotypes. The codon positions included were 1st+2nd+3rd+Noncoding. The sequences with the highest level of identity to the studied isolate are shown in bold.

Figure 2. Phylogenetic tree, based on the nucleotide sequences of 419-bp coding region of the F gene of the studied isolate, compared to more than 100 NDV isolates with different genotypes submitted in GenBank. The maximum likelihood method and K2+G model with 1000 Bootstrap replicates were used. The taxa names include the GenBank accession number (Dimitrov et al. classification/old classification), hostname, country of isolation, strain designation, and year of isolation. The isolated subject of the current study is marked by a circle. Class I is assigned as an outgroup.

Figure 3. Phylogenetic tree, based on complete nucleotide sequences of P genes of the studied isolate, compared to more than 100 NDV isolates with different genotypes submitted in GenBank. The maximum likelihood method and K2+G model with 1000 Bootstrap replicates were used. The taxa names include the GenBank accession number (Dimitrov et al. classification/old classification), hostname, country of isolation, strain designation, and year of isolation. The isolated subject of the current study is marked by a circle. Class I is assigned as an outgroup.

Figure 4. Phylogenetic tree, based on complete nucleotide sequences of M genes of the studied isolate, compared to more than 100 NDV isolates with different genotypes submitted in GenBank. The maximum likelihood method and K2+G model with 1000 Bootstrap replicates were used. The taxa names include the GenBank accession number (Dimitrov et al. classification/old classification), hostname, country of isolation, strain designation, and year of isolation. The isolated subject of the current study is marked by a circle. Class I is assigned as an outgroup.

Table 5. Nucleotide substitutions pattern in the NDV isolate presented in the current work. Transition substitutions shown in bold and transversion substitutions are shown in italics.

Multiple Alignments of Phosphoprotein and Matrix Proteins

Multiple alignments of P and M translated sequences from our isolate with those reported from other countries showed, as expected, more transition than transversion mutations (Table 5), plus higher frequencies of third-place substitutions at each codon (data not shown). Further, this analysis identified Valine (V) to Isoleucine (I) and Alanine (A) to Valine (V) as the most common amino acid changes in both proteins (data not shown).

P nucleotide sequence in our isolate showed more variability than the M gene, and most of the variations in P had occurred within the N-terminal site of the protein (data not shown). Our isolate's nuclear localization sequence (NLS) of the M protein (aa.s 247-263) was identical to that of the other VII genotypes reported from Iran, China, and Jordan. Still, it showed variations compared to that of VII.2 isolates identified in Pakistan, Indonesia, and Malaysia (Table 6). Two clusters of basic amino acids at positions 250-251 and 262-263 are reported to have a main role in the nuclear localization of M protein (Mayo, 2002). However, our isolate showed no substitution at these positions. Further, the FPIV-like domain (residues 23-26) and residue 42 of the M protein of our isolate also showed no variations compared to those sequences from the velogenic strains we tested (Table 6). But our isolate showed variations at positions 257 (L→I) and 259 (R→E) compared to Colemen and Peeples sequences (Coleman & Peeples, 1993). Notably, characteristic changes were observed in our isolate compared to LaSota and B1 vaccine strains, with four substitutions at 247, 257, 259, and 263 residues of the M protein (Table 6).

Table 6. Sequences of the important domains of the M protein of NDVs reported around the world including those of the NDV from the current work. Isolates are represented by the nucleotide and protein GenBank accession numbers host name, country of isolation, and year of isolation. The isolate subject of the current study is shown in bold.

Multiple Alignments of the Fusion Protein

At the fusion peptide region of the fusion protein, the two live attenuated vaccinal strains, i.e., ANW-12438 and ANQ45239 commonly used in Iran, showed the greatest variations compared to our isolate (Table 2). At this region, the V→I substitution at either position of 121 or 125 was the most common variation observed between our isolate and the other strains reported from Iran.

At the cleavage site, similar to the other velogenic strains tested, our isolate maintained the 112RRQ-KRF117 sequence (Table2). This sequence is different from that of the lentogenic vaccinal strains of ANW12438 and ANQ45239, which is 112GRQ-GRL117.

Consistent with the comparisons at the fusion peptide region and cleavage site, the hypervariable region of the F protein in our isolate showed the highest variation to the two vaccinal strains (ANW12438 and ANQ45239) commonly used in Iran (Table 7). The amino acids at positions 11, 20, and 27 were the most variable residues in this region. The C →R substitution at position 27 was the most frequent variation observed in our isolate compared to VII.2 sequences reported from Iran and other countries (Table 7).

The neutralizing epitopes at the F protein were highly conserved among the tested isolates (Table 7), andXX1.2 strain (QDE10569) was the only Iranian strain with amino acid substitution at position 78 (K→R) compared to our isolate.

Discussion

ND is enzootic to many countries, and NDV is constantly evolving, generating new genotypes (Miller et al., 2010). Vaccination is the gold standard measure in keeping ND under control. However, the occurrence of the disease in hosts previously infected, and hence deemed resistant, as well as in the vaccinated flocks, raises concerns on the possible emergence of strains that may cause another panzootic (Miller et al., 2010; Miller et al., 2015). Further, it is known that the attenuated vaccine strains can gnomically recombine with the wild strain(s) to produce new genotypes previously absent in the area (Miller et al., 2015). Therefore, knowing the genomic structure of the strain(s) prevalent in a region is not only helpful in keeping the disease under control by choosing effective vaccine strain(s), but it also helps in keeping the virus reservoir low by minimizing the introduction of new strains into the region.

Genomic sequencing provides the highest genotyping resolution and helps identify the NDV strains prevalent in a region. Here, we present a 419 base long nucleotide sequence of the F gene and full genomic sequences of the loci coding for P and M proteins in an NDV isolated from the Northeast of Iran. Our sequencing data and the associated bioinformatics analyses confirmed the presence of the VII.1.1 strain in the region.

Genomic studies have already shown NDVs under genotype VII as the most important genotypes circulating in Asia (Yang et al., 1999; Mase et al., 2002; Tsai et al., 2004; Lee et al., 2004; Miller et al., 2010; Diel et al., 2012; Ebrahimi et al., 2012; Miller et al., 2015; Shohaimi et al., 2015; Dimitrov et al., 2016; Liu et al., 2019; Almubarak, A. I., 2019; Xiang et al., 2020). Earlier studies based on the genomic variations at the cleavage site of the F gene indicated VII.1.1 as the prevalent sub-genotype in Iran (Ebrahimi et al., 2012; Hosseini et al., 2014; Kiani et al., 2016; Boroomand et al., 2016; Esmaelizad et al., 2017; Alborz, I., 2018; Sabouri et al., 2018; Goudarzi et al., 2019; Molouki et al., 2019; Beheshtian et al., 2020; Allahyari et al., 2020). Our analyses based on comparing the sequences of our isolate with those previously reported in Iran also verified the persistent presence of genotype VII.1.1 from 2011-2019 in Iran. However, recent studies show that VII.1.1 is not a unique genotype in the region (Molouki et al., 2019), but genotypes/sub-genotypes VI, VII.2, XIII.1.1, XIII.1.2, XX1.1, and XXI.2 are also reported from different parts of the country (Ebrahimi et al., 2012; Mayahi & Esmaelizad, 2017; Jabbarifakhar et al., 2018; Dizaji et al., 2020). Additionally, a study suggested that some evolutionary events in NDVs happen in Iran and Indian subcontinent countries (Ebrahimi et al., 2012).

Phylogenetic analyses using partial sequencing of P, M, and F genes had previously suggested a close relationship between some Iranian NDV isolates and those reported from Russia (AY865652) and Pakistan (JN682210) and India (Kianizadeh et al., 2002; Langeroudi et al., 2014; Ahmadi et al., 2016; Boroomand et al., 2016; Rezaei Far et al., 2017). However, our analyses based on the hypervariable region of F and full P and M gene sequences, as well as some similar recent works (Ebrahimi et al., 2012; Langeroudi et al., 2014; Sabouri et al., 2018; Molouki et al., 2019), suggests that Iranian NDV isolates are closest to Chinese VII.1.1 strains (Figures 2, 3 and 4). In addition to the deepness of the genomic interrogations concerning the number of genes and the length of the genomic sequences used, the fact that these viruses have been isolated from different geographic locations within Iran may be a reason for differences in the genotypes reported from Iran.

Residues 257 and 259, known as hypervariable residues of the M protein, were the same sequence in our isolate and other VII.I.I strains previously reported (Table 6). However, various changes, especi-ally substitutions of basic amino acids (R and K), were observed in VII.2 and vaccinal strains compared to our isolate. Further, our isolate showed variations at positions 257 and 259 compared to Col-emen and Peeples sequences (Coleman & Peeples, 1993); but, changes to the basic amino acids at these positions must not have an important effect on the nuclear localization efficiency of the virus (Coleman & Peeples, 1993).

Two other critical clusters of basic amino acids at positions 250-251 and 262-263 of M protein and FPIV-like domain had completely conserved sequences in all tested genotypes. This propounds the view that M can also be considered a target for virus typing along with F gene analysis (Seal et al., 2000).

LaSota and B1, the two vaccinal strains commonly used in Iran, are among the most deviated strains compared to our isolate (Tables 2 and 7). This may explain the reduced protections from the vaccination programs.

Considering that exotic and backyard birds can act as a source of virus infection for industrial poultry flocks (Alexander, 2001; Derksen et al., 2018), isolation of and studying NDVs (virulent or avirulent) present in exotic and backyard birds can also be beneficial in rooting the source of the NDVs found in commercial flocks. In this context, it is essential to note that just two amino acid changes in the NDV cleavage site of F protein can convert an avirulent strain to a virulent one (Alexander, 2001; Derksen et al., 2018).

Table 7. The amino acid substitutions in hypervariable regions and neutralizing epitopes of F protein in the studied isolate and some other similar sequences submitted in GenBank. Isolates are represented by the nucleotide and protein GenBank accession numbers host name, country of isolation, and year of isolation. The isolate subject of the current study is shown in bold.

Conclusion

In conclusion, our work shows the presence of NDV VII.I.I genotype already prevalent in the Far East in Iran. Although the path for disseminating these viruses into the region is currently not obvious, this finding highlights the importance of biosecurity management in controlling ND. Further, our isolate and the other VII.I.I viruses reported in Iran are genetically distant from the vaccine strains currently used in Iran. Hence, it is also important to have dynamic surveillance to catalog and continuously investigate the NDVs circulating in each region. This will enable authorities to administer effective and appropriate vaccination programs that will be of maximum effectiveness in preventing ND and minimize the possibility of introducing a new strain(s) into the region.

Acknowledgments

We would like to thank Dr. Saeed Yaghfuri for his help with bioinformatics analysis and Mr. Ali Kargar for his technical assistance. This study was supported by a research grant from the School of Veterinary Medicine, Ferdowsi University of Mashhad, Iran.

Conflict of Interest

The authors declared no conflict of interest.