Characterization of a novel VIIl sub-genotype of Newcastle disease virus circulating in Iran

ABSTRACT

Newcastle disease is an economically important and highly contagious disease affecting wild and domestic avian species. Despite extensive vaccination efforts within the poultry industry, Newcastle disease virus (NDV) outbreaks causing significant economic losses still occur. Rural chickens may act as a potential reservoir of NDVs for commercial poultry due to poor biosecurity and inadequate vaccination. The aim of this study was to investigate the phylogenetic relationship and molecular characterization of eight NDVs isolated from backyard poultry in Iran during 2011–2013. The complete coding sequence of fusion (F) and haemagglutinin-neuraminidase (HN) genes of eight NDVs were determined and compared with other published NDVs. Based on inter-population distances and phylogenetic topology between available NDV categories, Iranian isolates formed a novel VIIl sub-genotype distinct from previous groups designated in genotype VII. Furthermore, both F and HN genes of the Iranian isolates shared high nucleotide sequence similarity with viruses isolated in China. All viruses analysed contained a polybasic cleavage site motif (¹¹¹G/RRRQKR↓F¹¹⁷), indicating that all isolates could be categorized as a virulent pathotype. No mutation was observed in the neutralizing epitopes of the F protein. Analysis of amino acids associated with neutralizing antigenic sites within the HN protein revealed that all isolates exhibited a unique amino acid (Q) at position 347. These results emphasize the need for strengthening the biosecurity measures implemented on village flocks and practicing a mandatory vaccination programme for local poultry. Moreover, continuous monitoring of NDVs in different species of birds can help to gain more knowledge about the evolution of this virus and prevent future panzootics.

KEYWORDS:

• Newcastle disease virus
• F gene
• HN gene
• phylogenetic analysis
• genotyping
• chicken
• Iran

Introduction

Since the first outbreaks of Newcastle disease (ND) in the late 1920s, ND has become enzootic or a cause of regular epizootics for poultry producers in most parts of Asia, Africa and Central and South America. In less-developed countries, sporadic epizootics occur sometimes despite the use of prophylactic vaccination in the poultry industry (Alexander & Senne, 2008; Miller & Koch, 2013). ND is caused by the Newcastle disease virus (NDV), a member of the family Paramyxoviridae of the genus Avulavirus (ICTV, 2016). Based on pathogenicity, strains of NDV are categorized into four main pathotypes in decreasing order of virulence in chickens: velogenic, mesogenic, lentogenic and asymptomatic enteric strains. Subsequently, velogenic strains are subdivided into two forms of neurotropic and viscerotropic according to the main clinical signs (Miller & Koch, 2013).

The NDV genome consists of six genes, encoding six structural nucleocapsid (N), matrix (M), phosphoprotein (P), fusion (F), haemagglutinin-neuraminidase (HN) and large polymerase (L) proteins (Chambers et al., 1986). Moreover, two non-structural proteins of V and W, produced through RNA editing of the P mRNA, are only present in virally infected cells (Miller et al., 2010). Although pathogenicity of NDV is a polygenic trait, F and HN genes are the possible major contributing factors of virulence (Dortmans et al., 2011). A multi-basic cleavage site within the F protein is a crucial virulence determinant of NDV. Cleavage of the F precursor, F0 protein, into F1 and F2 by cellular proteases is required for the entry of NDV into the host cells (Glickman et al., 1988). Different viruses have different amino acid sequences in this site and consequently the number of these basic amino acids determines the type of proteases affecting them (Collins et al., 1993). HN is another key virulence factor. It acts as a multifunctional molecule and interacts with the F protein, required for promoting membrane fusion (Sergel et al., 1993). Additionally, these two proteins that form spike-like projections on the external surface of the viral envelope elicit protective neutralizing antibodies against NDV (Collins et al., 1993; Miller & Koch, 2013).

A recent unified nomenclature and classification system of NDVs has been designed based on the open reading frame (ORF) of the F gene, which is supported by phylogenetic topology and evolutionary distances between established genotypes or sub-genotypes (Diel et al., 2012). In this system, NDV is categorized into two classes, class I and class II. Class I consists of a single genotype, containing avirulent viruses usually isolated from waterfowl and shorebirds worldwide. Class II viruses are typically derived from wild birds and poultry species and are classified as either virulent or avirulent. The virus strains commonly used as vaccine seeds are grouped into genotype II. Viruses of genotypes I–XVIII form class II, which are recovered from both poultry and wild birds globally. Genotypes V–VII are to blame for the recent outbreaks in different continents (Diel et al., 2012; Miller et al., 2015). Some other genotypes have also been isolated but their occurrence is rare. Currently, genotype VII strains are the spillover of most outbreaks in the Middle East and Asia. These viruses are capable of causing significant economic losses, even in vaccinated poultry flocks, and expand their host range, sometimes resulting in disease in waterfowl. (Diel et al., 2012; Miller et al., 2015; Dimitrov et al., 2016).

Since the first detection of ND in 1951, this disease has been recovered from different species of birds in Iran (Sohrab, 1974). Despite extensive use of vaccination, outbreaks of NDV are observed in Iranian commercial poultry (Hosseini et al., 2014; Mehrabanpour et al., 2014). Non-commercial poultry and waterfowl seem to be the natural reservoirs for both virulent and avirulent NDVs (Schelling et al., 1999; Miller et al., 2015; Wajid et al., 2017). Due to insufficient biosecurity measures, NDVs are more likely exchanged between terrestrial poultry and wild birds in both directions, allowing new potential virulent strains to emerge. Subsequently, these terrestrial poultry play a crucial role in evolution of NDVs (Miller et al., 2015). Little is known of the molecular characteristics and genotype distribution of NDVs in the country, especially the NDV isolates from backyard poultry. Moreover, other studies were limited to partial or full sequencing of the HN or F genes, mostly in commercial poultry (Kianizadeh et al., 2002; Fathi et al., 2008; Ebrahimi et al., 2012; Hosseini et al., 2014; Samadi et al., 2014; Sabouri et al., 2016). Therefore, it is imperative to investigate the evolution and properties of these viruses. In the current study, the coding sequences of F and HN genes of eight Iranian NDVs, circulating in backyard poultry from 2011 to 2013, were sequenced and the sequences were analysed phylogenetically and molecularly.

Materials and methods

Virus isolation and identification

Eight NDVs were isolated from tracheal and intestinal samples of severely ill or dead backyard poultry, which were referred to the pet birds’ clinic of the Faculty of Veterinary Medicine, the University of Tehran, during the period 2011–2013. The involved backyard poultry were all domestic chickens. The chickens affected showed coughing, rattling when breathing, central nervous system signs and general signs of illness. All viruses were propagated in the allantoic cavity of 10-day-old pathogen-free embryonated chicken eggs, as described previously (Terregino & Capua, 2009). The allantoic fluids were harvested and confirmed to contain NDVs by reverse transcription-polymerase chain reaction (RT-PCR) method (Kant et al., 1997) and then frozen in aliquots at −70°C for further analysis.

Pathogenicity studies

According to OIE standards, pathogenicity of the studied viruses was determined using mean death time (MDT) and intracerebral pathogenicity index (ICPI) in 10-day-old pathogen-free embryonated chicken eggs and 1-day-old chicks, respectively (OIE, 2008).

RNA isolation, reverse transcription-PCR and sequence analysis of F and HN genes

Viral RNA was extracted from virus-infected allantoic fluid using RNX-Plus^™ solution (CinnaGen, Tehran, Iran) according to the manufacturer’s protocol. Reverse transcription was conducted with a commercial cDNA synthesis kit (RevertAid^™ first strand cDNA synthesis kit; Thermo Scientific, MA, USA) based on random priming. ORFs of both genes were amplified by Pfu DNA polymerase (Vivantis, Selangor, Malaysia) with two pairs of overlapping primers designed specifically for each gene (primer sequences are available upon request). After agarose gel electrophoresis, PCR products were purified with a gel extraction kit (GeneJET^™ gel extraction kit; Thermo Scientific, MA, USA) and sequenced by a commercial DNA sequencing company (Bioneer, Daejeon, South Korea) in both directions.

Molecular and phylogenetic analysis

All sequence data were compiled and edited using BioEdit Package, version 7.2.5 (Hall, 1999). Blast homology searches were used to retrieve sequences related to our queries from the GenBank databases. Multiple nucleotide and deduced amino acid sequence alignment were carried out by Clustal W. A preliminary phylogenetic analysis was conducted based on complete coding sequences of the F and HN genes together with 200 other GenBank sequences using the maximum likelihood method with the general time-reversible model and 500 bootstrapping replications in MEGA6 software (version 6.06) (Tamura et al., 2013).

Genotype definition

Classifications of genotypes and sub-genotypes were designated based on the topology of the phylogenetic tree and mean inter-population evolutionary distances between 200 complete F gene coding sequences using the criteria previously described (Diel et al., 2012; Miller et al., 2015). To assign different genotypes and sub-genotypes, cutoff values of mean distance >10% and 3–10% were used, respectively (Diel et al., 2012). To conclude whether the eight NDV strains were really classified as a novel sub-genotype, the F gene sequences were analysed with another 192 sequences from genotype VII. The mean evolutionary distances between sub-genotypes were computed using the maximum composite likelihood method with the gamma distribution model (shape parameter = 1) and 500 bootstrap replicates in MEGA6 software (version 6.06) (Tamura et al., 2013). The codon positions included in the analysis were the 1st + 2nd + 3rd + noncoding. All positions containing gaps and missing data were eliminated.

Nucleotide sequence accession numbers

All nucleotide sequences analysed in this study are available in GenBank under the following accession numbers: KU201408-KU201415 and KU201416-KU201423 for the F and HN genes, respectively.

Results

Pathogenicity of isolates

MDT and ICPI of the isolates were determined in accordance with instructions in the OIE manual (OIE, 2008). All eight NDV isolates had MDT of <60 h and ICPI of >1.5, and therefore were classified as velogenic strains. Moreover, multiple basic amino acids were identified in these viruses at the C-terminus of the F2 protein, which is considered as an essential factor for virulent strains.

Homology analysis

Comparison of the nucleotide and deduced amino acid sequences of F genes showed that all eight isolates had noticeably high identities with 98.8–99.9% and 98.6–99.8%, respectively. The HN genes of the Iranian isolates shared 98.8–99.9% nucleotide and 99.3–100% amino acid sequence homology. The nucleotide sequence homology observed between F genes of the tested isolates and other available complete ORFs of Iranian strains was 89.7–90.6% (data not shown). However, ORFs of the HN gene of other Iranian strains were not available for comparison.

The BLAST analysis of F and HN genes of the eight Iranian isolates are shown in Table 1. Both studied genes displayed the highest nucleotide sequence identities to Chinese strains isolated from some members of Anseriformes, indicating the possible origin of the current circulating viruses.

Table 1. Identities of F and HN genes between Iranian isolates and other viruses at the nucleotide level.

Virus	Gene	More similar to
Studied isolates	F KU201408-KU201415	Go/CN/JSG0210/02 JF340367 (96.9–97.2%)^a	MDk/CN/FP1/02 FJ872531 (96.9–97.2%)	MDk/CN/LGD/1/05 KM885167 (96.9–97.2%)
	HN KU201416-KU201423	Go/CN/JSG0210/02 JF340367 (96.0–96.4%)^b	Dk/CN/WF00D/00 FJ754272 (96.0–96.3%)	MDk/CN/FP1/02 FJ872531 (95.9–96.3%)

Note: Av: avian; Ck: chicken; Dk: duck; Fo: fowl; Go: goose; Ma: mallard; MDk: Muscovy duck; Pc: peacock; Pg: pigeon; Po: poultry; AU: Australia; CN: China; HN: Honduras; HU: Hungary; IN: India; ID: Indonesia; IR: Iran; IE: Ireland; IT: Italy; MG: Madagascar; MX: Mexico; NA: Namibia; PK: Pakistan; PE: Peru; US: United States.

^aHomology was calculated based on the F gene ORF (4550-6211).

^bHomology was calculated based on the HN gene ORF (6418-8133).

Phylogenetic analysis

Phylogenetic trees were constructed with eight sequences generated in the present study and other NDVs based on the ORFs of F and HN genes. The initial phylogenetic tree revealed that all F genes of the eight studied NDVs were grouped into a new sub-genotype of genotype VII (Figure 1). A further phylogenetic analysis using 192 strains from genotype VII proved that all recent Iranian NDVs fell into sub-genotype VIIl in a well-supported cluster (99% bootstrap value) (Supplementary Figure 1). According to the criteria proposed by Diel et al. (2012), the distances found between the eight new Iranian NDV isolates and other sub-genotypes of VII ranged between 0.04 and 0.1, with the highest and lowest distance values associated with the VIIk and VIId sub-genotypes, respectively (Table 2). However, other Iranian isolates, found in GenBank, grouped into genotype XIII together with other NDVs isolated from Pakistan and India (Diel et al., 2012). In the HN phylogenetic tree, viruses isolated from backyard poultry formed a new group close to the sub-genotype VIIb of genotype VII (Figure 2).

Figure 1. Phylogenetic tree of NDV F gene sequences from Iranian isolates based on the ORF (4550-6211). The phylogenetic tree was generated using the maximum likelihood method with the MEGA6 programme (version 6.06) (Tamura et al., 2013). The number near the nodes displays the percentage of bootstrap values of 500 replicates. Bootstrap values ≥60% are shown in this figure. The classification system used according to Diel et al. (2012) has been shown in this figure. Vertical lines indicate genotypes and sub-genotypes of class II. The viruses obtained in this study and previously published Iranian sequences are marked with squares and circles, respectively. Abbreviations can be found in the legend to Table 1.

Figure 2. Phylogenetic tree of NDV HN gene sequences from Iranian isolates based on the ORF (6418-8133). The phylogenetic tree was generated using the maximum likelihood method with the MEGA6 programme (version 6.06) (Tamura et al., 2013). The number near the nodes displays the percentage of bootstrap values of 500 replicates. Bootstrap values ≥60% are shown in this figure. The classification system used according to Diel et al. (2012) has been shown in this figure. Vertical lines indicate genotypes and sub-genotypes of class II. The viruses obtained in this study are marked with squares. Abbreviations can be found in the legend to Table 1.

Table 2. Estimates of mean interpopulational evolutionary diversity among new and previously designated sub-genotypes of genotype VII based on complete coding sequences of the F gene.

	VIIb	VIId	VIIe	VIIf	VIIg	VIIh	VIIi	VIIj	VIIk	VIIl
VIIb (n = 27)		0.004	0.005	0.006	0.007	0.010	0.010	0.003	0.010	0.006
VIId (n = 51)	0.041		0.003	0.005	0.006	0.008	0.009	0.005	0.010	0.004
VIIe (n = 16)	0.049	0.037		0.005	0.007	0.007	0.008	0.006	0.009	0.005
VIIf (n = 5)	0.57	0.047	0.042		0.008	0.007	0.008	0.007	0.008	0.006
VIIg (n = 4)	0.092	0.078	0.087	0.093		0.011	0.010	0.007	0.011	0.008
VIIh (n = 20)	0.104	0.094	0.088	0.079	0.135		0.009	0.009	0.010	0.009
VIIi (n = 37)	0.096	0.087	0.083	0.078	0.129	0.091		0.010	0.010	0.009
VIIj (n = 26)	0.030	0.048	0.055	0.066	0.096	0.109	0.103		0.010	0.006
VIIk (n = 6)	0.121	0.115	0.108	0.101	0.156	0.112	0.104	0.125		0.010
VIIl (n = 8)	0 . 057	0.041	0.049	0.059	0.094	0.102	0.099	0.062	0.127

Notes: The number of base substitutions per site from averaging the overall sequence pairs between groups is shown. Standard error estimates are shown above the diagonal and were obtained by a bootstrap procedure (500 replicates). Analyses were conducted using the maximum composite likelihood model. The rate variation among sites was modelled with a gamma distribution (shape parameter = 1). The analysis involved 200 nucleotide sequences. Codon positions included were 1st + 2nd + 3rd + Noncoding. All positions containing gaps and missing data were eliminated. There were a total of 1662 positions in the final data set. Evolutionary analyses were conducted in the MEGA6 programme (version 6.06) (Tamura et al., 2013). The average evolutionary distance of NDV isolates obtained in the present study with other sub-genotypes from genotype VII are shown in bold.

Characterization of the F gene

All the characterized NDVs, with the exception of Ck/IR/SMV-7/13, had the same ¹¹¹G-R-R-Q-K-R¹¹⁶ amino acid motif at the cleavage site of the F protein, which is characteristic of virulent NDVs in chickens (Table 3). The Ck/IR/SMV-7/13 virus had a G111R mutation, an additional basic amino acid, resulting in the motif of ¹¹¹R-R-R-Q-K-R¹¹⁶. Potential glycosylation sites with the N-X-T/S motif (in which X may be any amino acid except proline) were identified. All eight isolates possessed six conservative glycosylation sites in comparison to other NDVs, five of which were in the F1 portion and one was in the F2 part of the molecule (Table 3). Examination of amino acid sequences at neutralizing epitopes revealed that all isolates had conserved amino acids in comparison to other representative NDVs (Table 3).

Table 3. F gene comparison of the studied Iranian strains with other NDVs.

Virus strain	Cleavage site (111–117)	Glycosylation site	Neutralizing epitopes
			72	74	75	78	79	157–171	343
Studied isolates KU201408-KU201415 (VIIl)	G/RRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTNG	L
Ck/CN/FMW/06 GU564399 (VIIb)	RRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
MDk/CN/FP1/02 FJ872531 (VIId)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Go/CN/WF00G/00 FJ754273 (VIId)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Go/CN/NA-1/06 DQ659677 (VIIe)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	R	A	SIAATNEAVHEVTDG	L
Pg/CN/ND-03-018/03 GQ338309 (VIIf)	GKRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/CN/SRZ03/03 EU167540 (VIIg)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/ID/bali/020/10 HQ697261 (VIIh)	GRRRKR*F	85, 191, 366, 447, 471, 541	D	E	A	R	A	SIAATNEAVHEVTNG	L
Ck/ID/Kudus/017/10 HQ697259 (VIIi)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/IR/EMM-4/10 JQ267582 (XIIIa)	GRRRKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/IR/EMM-5/11 JQ267581 (XIIIa)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/IN/Bareilly/01/10 KJ577585 (XIIIa)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/BYP/PK/10 JN682210 (XIIIb)	GRRQKR*F	85, 191, 366, 447, 471, 541	D	Q	A	K	A	SIAATNEAVHEVTDG	L
Fo/US/B1/47 AF309418 (II)	GGRQGR*L	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L
Ck/US/Lasota-AF/46 AY845400 (II)	GGRQGR*L	85, 191, 366, 447, 471, 541	D	E	A	K	A	SIAATNEAVHEVTDG	L

Note: Abbreviations can be found in the legend to Table 1.

Characterization of the HN gene

Comparison of amino acids constituting the neutralizing epitopes of the HN protein revealed that all of the isolates tested had glutamine at residue 347, which is unique to these isolates (Table 4) in comparison to the other 200 selected HN sequences. Analysis of the HN sequence of the viruses showed the presence of at least five potential glycosylation sites (N-X-T/S) with a unique glycosylation site at residue 144, in which four of them were conserved in other examined strains. The two characterized viruses (Ck/IR/SMV-3/11 and Ck/IR/SMV-4/12) exhibited loss of the N-glycosylation site at amino acid 341, which was not observed in other viruses (Table 4). Residues E401, R416 and Y526 seem to be the key molecular determinants of receptor binding. All the isolates possessed the same amino acids in these sites, except for the Ck/IR/SMV-7/13 isolate, which had a mutation of E401K. The significance of this alteration is not known (Table 4).

Table 4. HN gene comparison of the studied Iranian strains with other NDVs.

Virus strain	Glycosylation site	Receptor binding site			Neutralizing epitopes
		401	416	526	193–201	263	287	321	332–333	345–353	356	494	513–521	569
Studied isolates KU201416-KU201423 (VIIl)	119, 144, 341^a, 433, 481	E/K	R	Y	LSGCRDHSH	K	D	K	GK	PDQQDYQIR	K	D	RVTRVSSSS	D
Ck/CN/XJ-2/97 JN618348 (VIIb)	119, 341, 433, 481, 508	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	D
MDk/CN/FP1/02 FJ872531 (VIId)	119, 341, 433, 481, 508	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	D
Go/CN/ZJ1/00 AF431744 (VIId)	119, 341, 433, 481, 508, 538	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	D
Go/CN/NA-1/06 DQ659677 (VIIe)	119, 341, 433, 481	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	D
Pg/CN/ND-03-018/03 GQ338309 (VIIf)	119, 341, 433, 481, 508	E	R	Y	LSGCRDHSH	K	D	K	GK	PDGQDYQIR	K	D	RVTRVSSSS	G
Ck/CN/Hebei/01/06 KC542895 (VIIg)	119, 341, 433, 481, 508	E	R	Y	LSGCRDHSH	K	D	K	GK	PDKQDYQIR	K	D	RVTRVSSSS	D
Ck/ID/bali/020/10 HQ697261 (VIIh)	119, 341, 433, 481, 508, 538	E	R	Y	LSGCRDHSH	R	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	D
Ck/ID/Kudus/017/10 HQ697259 (VIIi)	119, 341, 433, 481, 538	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	N
Ck/IN/Bareilly/01/10 KJ577585 (XIIIa)	119, 341, 433, 481	E	R	Y	LSGCRDHSH	K	D	K	GK	PDEQDYQIR	K	D	RVTRVSSSS	V
Ck/BYP/PK/10 JN682210 (XIIIb)	119, 341, 433, 481	E	R	Q	LSGCKDHSH	K	D	K	GK	PDDQDYQIQ	K	D	RVTRVSSSS	V
Fo/US/B1/47 AF309418 (II)	119, 341, 433, 481, 538	E	R	Y	LSGCRDHSH	N	D	K	GK	PDEQDYQIR	K	G	RITRVSSSS	D
Ck/US/Lasota-AF/46 AY845400 (II)	119, 341, 433, 481, 538	E	R	Y	LSGCRDHSH	N	D	K	GK	PDEQDYQIR	R	G	RITRVSSSS	D

Note: Abbreviations can be found in the legend to Table 1.

^aTwo isolates (SMV-3 and SMV-4) displayed loss of the glycosylation site at this position.

Discussion

Considering the capability of NDV to infect all domestic and wild bird species, investigation of backyard poultry as a potential reservoir of this virus is crucial. Backyard poultry are reared for either companionship or small-scale meat and egg production to provide a source of income in most developing countries (Alexander, 2001). Due to insufficient biosecurity measures and poor or absent disease control strategies, they might play a prominent role in the transmission of NDV to commercial poultry and wild birds (Schelling et al., 1999; Shekaili et al., 2015; Wajid et al., 2017). Simultaneously with the poultry improvement plan introduced in the 1950s, NDV caused severe disease in chickens with high mortality and soon after spread across the country (Sohrab, 1974). From the beginning, vaccination of commercial poultry in Iran has been used to bring ND under control; however, sporadic cases occurred despite the implementation of these prophylactic programmes even with induction of high haemagglutinin-inhibition antibody levels (Sohrab, 1974; Hosseini et al., 2014). Hence, the study of molecular and phylogenetic analysis of the HN and F genes encoding two major NDV immunogenic proteins is worthwhile to identify the common genotypes circulating in such a reservoir and investigate the future ND outbreaks.

The sequence analysis of F genes revealed that all the studied isolates were closely related to each other (sequence homology of 98.8–99.9%). Moreover, sequence comparison of the HN gene also showed the highest percentage of homology between the analysed strains ranging from 98.8% to 99.9%. They also shared a lower homology (about 90%) with those of Iranian sequences previously published. All these results suggest that the studied isolates might have originated from a common ancestor.

Unlike previous studies, all eight Iranian strains showed a high percentage identity of over 97% in both genes from Chinese viruses recovered from species of the order Anseriformes (Muscovy duck, goose and duck). In the early studies on partial sequences of NDV isolates from different regions of Iran, it was reported that Iranian isolates were closely related to VOL95, a Russian NDV isolate from 1995 which was a possible origin of the outbreaks of ND in different parts of Iran in 1995–1998 (Kianizadeh et al., 2002; Fathi et al., 2008). Since 1997, virulent strains of NDV have been frequently reported in waterfowl in China, and it was previously assumed that this virus was not capable of inducing overt clinical diseases in these species (Liu et al., 2003; Shi et al., 2011; Ababneh et al., 2012). The majority of the viruses isolated from these hosts consisted of lentogenic viruses while virulent ones were rarely identified. Different studies demonstrated that waterfowl could act as a natural reservoir of NDVs (Chong et al., 2010; Miller et al., 2015). Due to the flyways of migratory birds over Iran and lack of biosecurity measures for non-commercial poultry, interaction between wild aquatic birds and rural chickens might have occurred. Consequently, new recombinant viruses emerged and infected commercial and backyard poultry (Spalatin & Hanson, 1975; Kelleher et al., 1985). Moreover, agricultural trade between Iran and other countries could be considered another risk factor, causing the introduction of these new strains.

Cleavability of the precursor F0 into F1 and F2 by host cell proteases is considered to be the main molecular basic of pathogenicity in NDVs. Lentogenic viruses possess a monobasic amino acid motif at the F cleavage site, ¹¹²G-R/K-Q-G-R↓L¹¹⁷, which are cleaved extracellularly by trypsin-like enzymes restricted to enteric or respiratory epithelia. Cleavage of the mesogenic and velogenic strains with the multi-basic amino acid motif ¹¹²R/G/K-R-Q/K-K/R-R↓F¹¹⁷ occurs intracellularly by ubiquitous furin-like enzymes (Glickman et al., 1988; Collins et al., 1993). In accordance with pathogenicity assays (MDT and ICPI tests), the F protein of all eight isolates harboured multiple basic amino acids (¹¹¹G/R-R-R-Q-K-R↓F¹¹⁷), a feature of virulent strains, which is similar to the Iranian sequences previously published (Fathi et al., 2008; Hosseini et al., 2014; Samadi et al., 2014).

Changes in glycosylation sites can have dramatic effects on biological activity, virulence and infectivity of viruses (Vigerust & Shepherd, 2007). Analysis of potential N-glycosylation sites (N-X-S\T) in the F protein revealed that all the analysed isolates have a similar pattern of glycosylation sites that are conserved in all strains (Yusoff & Tan, 2001; Samal et al., 2012). Mutation in these sites can increase replication and virulence of NDV by altering the fusogenic nature of the F protein (Samal et al., 2012). The potential glycosylation sites of HN were located at residues 119, 144, 341, 433 and 481; however, two strains (SMV-3 and SMV-4) have lost a glycosylation site at position 341. Moreover, the studied isolates exhibited a unique potential glycosylation site at residue 144 in comparison to other NDVs. It has been reported that mutation in glycosylation sites exerts an influence on pathogenicity of NDVs by changing the biological function of the HN protein and modulating the virulence of the virus (Panda et al., 2004).

Mutation in the neutralizing epitopes of the F and HN proteins of NDVs may produce antigenic variants and escape antibody neutralization (Chambers et al., 1988; Iorio et al., 1989; Neyt et al., 1989; Qin et al., 2008; Snoeck et al., 2013). Among the seven neutralizing epitopes of the F protein, all isolates showed no mutation compared to other strains of NDV. However, examination of amino acids in neutralizing epitopes of the HN protein revealed that all the studied isolates had one unique amino acid substitution at position 347 (Q) compared to other strains. A mutation of E347K in the neutralizing epitope of the HN protein resulted in a neutralization escape variant (Gotoh et al., 1988; Iorio et al., 1991). Whether this amino acid can have the same effect in the enhanced virulence of these isolates needs further investigation.

Three amino acid residues at positions 401 (E), 416 (R) and 526 (Y) that are considered to have a pivotal role in receptor binding of the HN protein (Connaris et al., 2002), were conserved in all the viruses analysed except for the E401K mutation in Ck/IR/SMV-7/13. Viruses with an Y526Q mutation in the HN protein of NDV had reduced functional activities and attenuated viral replication and virulence (Khattar et al., 2009). Whether mutation at residue 401 has similar effects is not known.

There used to be two conventional NDV genotyping systems based on the partial or entire sequence of the F gene until 2012 (Aldous et al., 2003; Czegledi et al., 2006; Kim et al., 2007). In 2012, Diel et al. proposed a new unified system for classification of NDVs attributable to considerable genetic variety within the NDV sub-genotypes (Diel et al., 2012). This unified system grouped NDVs into two distinct classes; class I included just a single genotype whereas class II was composed of 18 genotypes with subsequent dividing of genotypes I, V, VI, VII, XIII, XIV, XVII and XVIII into sub-genotypes (Diel et al., 2012; Miller et al., 2015; Dimitrov et al., 2016). Most diversity was observed in genotype VII viruses, which are associated with ND outbreaks in the Middle East and Asia. Genotype VII viruses likely appeared in the Far East in the 1990s and eventually spread to Asia, Europe, South Africa and South America. Currently, genotype VII viruses have been classified into sub-genotype VIIb and VIId–VIIk (Dimitrov et al., 2016; Molini et al., 2017; Xue et al., 2017). The sub-genotype VIId strains spread rapidly throughout the world and are considered as one of the major circulating genotypes (Zhang et al., 2014). All the eight characterized NDVs clustered into a novel VIIl sub-genotype in contradiction to earlier Iranian findings. Ebrahimi et al. categorized Iranian isolates as VIIb sub-genotype based on the entire ORF sequence of the F gene, whereas their isolates fell into sub-genotype XIIIa according to the recently developed genotyping (Ebrahimi et al., 2012). In another study, genetic analysis of partial F sequences of the Iranian isolates displayed that Iranian isolates belonged to genotype VII (Samadi et al., 2014). This discrepancy between genotyping of current and previous studies may have resulted from different genotyping systems used in various studies. In the phylogenetic tree of HN, the topology of the tree showed a slightly different pattern to the F gene phylogenetic tree in which the analysed isolates formed a new sub-genotype close to sub-genotype VIIb.

To our knowledge, this study is the first report to characterize NDVs isolated form rural chickens in Iran. All the eight analysed viruses belonged to a novel sub-genotype (VIIl) in genotype VII. Furthermore, these NDVs displayed the highest nucleotide sequence identities to Chinese isolates from waterfowl, suggesting that they might be considered as a source of NDVs in backyard poultry; however, the whole genome sequencing of different Iranian NDV strains will be crucial to determine genetic relatedness among NDVs from different geographical regions of the world. With regard to the substitutions observed in key positions of the studied HN genes, further studies are needed to evaluate these alterations in order to obtain more information on the virulence and antigenic change of these mutations. Since local chickens serve as a reservoir of NDV and pose a threat to commercial poultry, implementation of biosecurity measures and routine ND vaccination are required for this group of birds (Conan et al., 2012). Taken together, continuous surveillance in different species of birds would help to gain an insight into the evolution of NDVs and control of future panzootic viruses.

Acknowledgements

The authors would like to thank Dr Claudio L. Afonso and Dr Kiril M. Dimitrov for their technical assistance. We are also grateful to Dr Aidin Molouki for his useful comments on the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by a grant from the Research Council of the University of Tehran [No. 7508001/6/15].