Phylogenetic assessment reveals continuous evolution and circulation of pigeon-derived virulent avian avulaviruses 1 in Eastern Europe, Asia, and Africa

Sample collection and isolates background data

The samples and viruses studied here were collected from different species (pigeons, chickens and quail) in Egypt, Pakistan, South Korea, Ukraine and Bulgaria, representing three different continents (Africa, Asia and Europe). During April 2015, one hundred sixty-seven (n = 167) oropharyngeal and cloacal/fecal swabs were collected in Egypt from apparently healthy Rock Pigeons. Samples were collected from pigeon lofts in El Fayoum and Qena provinces, and a live bird market (LBM) in Cairo named Souq al-Goma’a (also named Souq Sayeda Aisha or Friday Market). The swabs were placed immediately into tubes with 3 ml of brain-heart-infusion broth (Difco, New Zealand) supplemented with penicillin G (10,000 IU /ml), amphotericin B (20 μg/ml), and gentamycin (1000 μg/ml). After collection, the samples were labelled, stored on ice, transported to the lab where kept frozen at −76 °C, and shipped on dry ice to the Southeast Poultry Research Laboratory (SEPRL) of the United States Department of Agriculture (USDA). Additionally, three more samples were collected from pigeons in Pakistan in 2015 and shipped to SEPRL, USDA. Six additional NDV from repositories in Bulgaria (n = 2), Ukraine (n = 2), and South Korea (n = 2), were also sent to SEPRL for further characterization. The Bulgarian and the Ukrainian viruses were passaged in eggs 3 times and 2 times, respectively. The Korean isolates were passaged 2 times at SEPRL and passage information was missing prior to their receiving. Nine more NDV, isolated between 2014 and 2016 from healthy and diseased pigeons (pet and zoo birds), were studied at the Quality Operations Laboratory (QOL) at the University of Veterinary and Animal Sciences (UVAS), Lahore, Pakistan.

Virus isolation, virus propagation and intracerebral pathogenicity index test

Initial screening of all Egyptian samples employing the NDV and avian influenza matrix gene rRT-PCR assays [24, 27], revealed that out of 167 oropharyngeal and cloacal swabs, 71 samples had cycle threshold (Ct) values ≤35 (40.3%) in the NDV assay and all samples were negative in the avian influenza assay. At SEPRL, thirty-one (selected to achieve representativeness of all locations) of the 71 rRT-PCR positive Egyptian samples and three Pakistani samples from 2015 (designated 21A, 22A and 25A) were selected for further studies and inoculated into 9-to-11-day-old SPF ECE, following standard procedures [28]. The SPF ECE and chickens used to characterize these viruses were from the SEPRL SPF White Leghorn flock. The allantoic fluids from both eggs with embryo mortality and embryos alive at the end of the incubation period were collected and tested by hemagglutination assay. All the hemagglutinating agents were confirmed as NDV using a HI assay with NDV specific antiserum [3]. Additionally, The NDV obtained from repositories in Ukraine, Bulgaria and South Korea were propagated into 9-to-11-day-old SPF ECE following the same procedures [28]. Assessment of the virulence in vivo of three selected viruses was done by the ICPI test using one-day-old SPF chickens following established procedures [3].

Complete fusion protein gene sequencing

RNA from the two South Korean isolates was extracted from infected allantoic fluid using TRIzol LS Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions at SEPRL. RT-PCR was performed to amplify the complete F gene using the Superscript™ III One-step RT-PCR kit with Platinum Taq DNA polymerase (Invitrogen, Carlsbad, CA, USA), per manufacturer’s instructions using a set of F-gene specific primers (4331F/5090R, MSF1/NDVR2, 4927F/5673R and 5491F/6341R) [29] (see Additional file 1: Table S1). The PCR amplicons were processed and sequenced as described previously [30]. RNA isolation and nucleotide sequencing of nine Pakistani viruses were performed in Quality operations laboratory at University of Veterinary and Animal Science in Pakistan as follows: RNA extraction and RT-PCR F-gene amplification were performed as described above using a set of F-gene specific primers (see Additional file 1: Table S1). The amplicons were electrophoresed using 1% agarose gel and purified using QIAquick® Gel Extraction Kit (Qiagen, Valencia, CA, USA). The purified products were sequenced using the ABI 3130 automated sequencer (Applied Biosystem Inc., Foster City, CA, USA), as described previously [30].

Complete genome sequencing using next-generation sequencing (NGS)

Total viral RNA of six Egyptian (positive in virus isolation assay), two Ukrainian, three Pakistani, and two Bulgarian viruses were extracted from the infected allantoic fluids using QIAmp® Viral RNA Mini Kit (Qiagen, USA) according to manufacturer’s instructions. The recovered RNA was quantified using the Qubit® RNA HS Assay Kit (Life Technologies, USA) in the Qubit® Fluorometer instrument (Invitrogen, USA). Newcastle disease virus RNA was captured and enriched using Sera-Mag beads (GE Healthcare Life Sciences, USA), and three biotin-labeled oligonucleotide probes targeting three different positions in the NDV genome: 1) 8-AGA GAA TCT GTG AGG TAC GA/3Bio/, 2) 5905-TTC TCA AGT CAT CGT GAC AG/3Bio/, and 3) 12226-CCC TGC ATC TCT CTA CAG/3Bio/. Reverse transcription was performed using the Moloney Murine Leukemia Virus Reverse Transcriptase (M-MLV RT) kit (Thermo Scientific, USA) according to the manufacturer’s instructions. The cDNA products were recovered and purified using the Agencourt® RNAClean® XP beads (Beckman Coulter, USA) according to the manufacturer’s instructions and quantified using Qubit® ssDNA Assay Kit (Thermo Fisher Scientific, USA) in the Qubit instrument. The purified cDNA products were tagmented and amplified for NGS by using 1 ng/5 μl (0.25 ng/μl in water) of the cDNA product employing the Nextera XT DNA Library Preparation Kit (Illumina, USA) following manufacturer’s protocol. The two Bulgarian samples were processed using the KAPA Stranded RNA-Seq Library Preparation Kit for Illumina platforms (Kapa Biosystems, USA) according to the manufacturer’s instructions.

The distribution size and concentration of DNA in the prepared libraries were checked on a Bioanalyzer 2100 and Qubit instrument using Agilent High Sensitivity DNA Kit (Agilent Technologies, Germany) and Qubit® dsDNA HS Assay Kit (Life Technologies, USA), respectively. Paired-end sequencing (2 × 250 base pairs) of the generated libraries was performed on an Illumina MiSeq instrument using the 500 cycle MiSeq Reagent Kit version 2 (Illumina, USA). Raw sequence data were analyzed and assembled using MIRA version 3.4.1 [31] within a customized workflow on the Galaxy platform [32] as described previously [33]. Short internal gaps at the 3′ UTR of the nucleoprotein gene were closed using Sanger technology with primers designed using the sequences obtained from NGS (see Additional file 1: Table S1). The 5′ and 3′ ends of the genomes reported here were sequenced and confirmed as described previously [34].

Collection of sequences and phylogenetic analyses

All available NDV complete genome and complete fusion protein gene coding sequences were downloaded from GenBank and curated, resulting in two large datasets (n = 331 and n = 1406, respectively). Together with the sequences obtained in the current study, all sequences were aligned using Multiple Alignment with Fast Fourier Transformation (MAFFT v7.017) [35] as implemented in the Geneious software v8.1.4 [36]. The datasets were used for two preliminary maximum likelihood phylogenetic analyses using MEGA6 [37]. Based on evolutionary relatedness to the viruses sequenced here, two smaller datasets (including newly obtained sequences) of complete fusion protein gene coding sequences (n = 82) (see Additional file 1: Table S2) and complete genome sequences (n = 83) (see Additional file 1: Table S3) were parsed from the larger datasets. Sequences of selected representative isolates from other NDV genotypes were also included in each of the datasets. The coding regions of the complete genome and complete fusion gene were used to construct final phylogenetic trees using MEGA6. To select best-fit substitution model, the Bayesian Information Criterion (BIC) and corrected Akaike Information Criterion (AICc) values were estimated using MEGA6. The General Time Reversible (GTR) model as implemented in MEGA6 with a discrete gamma distribution (4 categories [+G, parameter = 0.6519 for the full fusion gene tree and 0.5984 for the complete genome tree]) with 1000 bootstrap replicates was used in all data analysis. The phylogenetic trees were visualized and edited using tree explorer implemented in MEGA6 and branch lengths are proportional to the differences between the isolates. The evolutionary distances were inferred using MEGA6 and showed as the average number of base substitutions per site. Analyses were conducted using the maximum composite likelihood model [38] with a gamma distribution (shape parameter = 1) of rate variation among sites. For all analyses, the codon positions included were the 1st, 2nd, 3rd, and noncoding and positions containing gaps and missing data were eliminated from the datasets. The classification criteria proposed by Diel et al. for naming sub-genotypes and genotypes were followed in the current study [12].

Real-time reverse transcription polymerase chain reaction (rRT-PCR)

Allantoic fluids from all of the samples studied at SEPRL were subjected to the USDA-validated matrix gene rRT-PCR (M-gene assay), as well as to the USDA-validated fusion gene rRT-PCR assay (F-gene assay) described previously by Wise et al. [24]. An additional test was performed using pigeon-specific fusion protein gene assay as described previously by Kim et al. [39]. AgPath one-Step rRT-PCR Reagents (Thermo Scientific, USA) was used and the reactions were carried out in the Cepheid Real-Time Thermal Cycler (Life Science, USA). The samples that had a Ct value ≤35 were considered positive in all assays. Due to the failure of the validated F-gene rRT-PCR assay [24] and pigeon-specific fusion protein gene assay [39] to detect some of the genotype VI viruses studied here, the probe and the forward primer described by Kim et al., [39] were optimized as follows: F-4876 5′–[6-FAM] AAG CGY TTC TGT CTC YTT CCT CCT [BHQ_1]–3′ and (F + 4837) 5′- TGA TTC CAT CCG CAG GAT ACA AG -3′. Additionally, the reverse primer for the pigeon-specific F-gene assay (F-4837) was replaced with a new primer (F-4943) 5′- GCT GCT GTT ATC TGT GCC GA-3′. The optimized probe and newly designed primers were analyzed by OligoAnalyzer 3.1. (Integrated DNA Technologies, USA, https://eu.idtdna.com/calc/analyzer ) tool and checked for self-annealing, hairpin loops and heterodimers. The sequences of the optimized primers and probes along with the tested fifteen NDV and selected representatives from other genotypes were aligned and compared to the M-gene and F-gene primers and probes and the pigeon-specific fusion probe to identify variable sites that determined the different outcomes of the F-gene rRT-PCR assay [24, 39] (see Additional file 2: Figs. S1, S2, and S3).

Virus isolation and biological properties

Complete fusion protein gene and complete genome sequencing

The complete genome sequences of thirteen of the studied viruses were obtained at SEPRL, and the genome characteristics of these complete genomes are typical for NDV (See Additional file 1: Table S4). Additional eleven complete fusion gene coding sequences were also obtained and utilized for prediction of virulence and phylogenetic analysis. The fusion protein cleavage site of NDV is a major virulence determinant marker [40], and the deduced amino acid sequences of the fusion protein cleavage site revealed that all 24 isolates had multiple basic amino acids at positions 113–116, and a phenylalanine residue at position 117 (Table 1). The cleavage site motif of the Egyptian and Ukrainian viruses was ¹¹³RQKR↓F¹¹⁷ (n = 8). The Pakistani viruses (n = 12) had the following three motifs, ¹¹³RKKRF¹¹⁷, ¹¹³RQRR↓F¹¹⁷, and ¹¹³RQKR↓F¹¹⁷, sharing the last motif with the Bulgarian viruses (n = 2). The South Korean (n = 2) viruses had a ¹³RRKR↓F¹¹⁷ motif (Table 1). All of these motifs are characteristic of virulent NDV [3, 41].

Distance and phylogenetic analyses

In order to determine the phylogenetic relationship between the studied viruses and other NDV isolated worldwide, the complete fusion gene coding sequences obtained in the current work, along with sequences of highly related viruses, were used to construct a phylogenetic tree (n = 82) (Fig. 1 and Additional file 1: Table S2). The full fusion coding region of all available class II NDV and of genotype VI (n = 1430 and n = 281, respectively, both including the 24 sequences from the current study) was also used to estimate evolutionary distances (Tables 2 and 3). A second phylogenetic analysis was performed using 83 complete genome-concatenated coding sequences of viruses pertaining to class II (see Additional file 1: Table S3 and Additional file 2: Fig. S4). Both, the full fusion gene and the complete genome phylogenetic analyses displayed similar topology confirming the phylogenetic classification of the viruses studied here into different sub-genotypes related to viruses previously isolated in Asia, Africa and Eastern Europe. The phylogenetic analyses (Fig. 1 and Additional file 2: Fig. S4) demonstrated that these 24 NDV studied here the clustered with viruses of genotype VI. The isolates from pigeons grouped with previously described viruses from columbine birds, while the chicken isolates grouped with other viruses from chickens. The topology of the full fusion phylogenetic tree (Fig. 1) indicated the viruses from Egypt (n = 6) (0%–3.1% genetic distance among themselves) and Ukraine (n = 2) (1.9% genetic distance between themselves) isolated from pigeons grouped together and clustered within sub-genotype VIg. These viruses were closely related to viruses isolated from pigeons in Russia, Ukraine, Kazakhstan, and Nigeria during 2005–2014 (Fig. 1) [42–44]. These results were consistent with the results of the complete genome analysis (see Fig. S4). Four Pakistani viruses isolated during 2015–2016 clustered together in a separate monophyletic branch within sub-genotype VIg (Fig. 1). The genetic distance of these Pakistani viruses (4.7%) to the rest of the VIg viruses shows higher nucleotide diversity within the sub-genotype.

Sub-genotype (number of analyzed sequences)	VI a	VI b	VI c	VI e	VI f	VI g	VI h	VI i	VI Ethiopia
VI a (n = 165)
VI b (n = 10)	0.071
VI c (n = 19)	0.091	0.069
VI e (n = 16)	0.080	0.063	0.083
VI f (n = 16)	0.075	0.058	0.082	0.059
VI g (n = 25)	0.096	0.080	0.088	0.097	0.091
VI h (n = 9)	0.071	0.071	0.098	0.084	0.081	0.103
VI i (n = 4)	0.113	0.100	0.104	0.108	0.106	0.093	0.120
VI Ethiopia (n = 8)	0.125	0.099	0.106	0.119	0.119	0.107	0.123	0.125
VI m (n = 9)	0.111	0.091	0.099	0.112	0.101	0.085	0.112	0.107	0.120

^a The numbers of base substitutions per site from averaging over all sequence pairs between groups within genotype VI are shown. The analysis involved 281 nucleotide sequences. There were a total of 1647 positions in the final dataset

Genotype (number of analyzed sequences)	I	II	III	IV	V	VI	VII	VIII	IX	X	XI	XII	XIII	XIV	XV	XVI	XVII	XVIII
I (n = 139)
II (n = 171)	0.125
III (n = 10)	0.117	0.143
IV (n = 5)	0.104	0.130	0.084
V (n = 90)	0.191	0.206	0.178	0.148
VI (n = 272)	0.185	0.204	0.177	0.136	0.159
VII (n = 476)	0.180	0.213	0.168	0.139	0.159	0.132
VIII (n = 5)	0.144	0.165	0.133	0.100	0.129	0.119	0.122
IX (n = 35)	0.107	0.127	0.093	0.079	0.169	0.170	0.164	0.125
X (n = 11)	0.115	0.114	0.141	0.129	0.207	0.201	0.199	0.164	0.125
XI (n = 14)	0.201	0.224	0.191	0.130	0.229	0.234	0.239	0.197	0.172	0.224
XII (n = 9)	0.192	0.224	0.179	0.152	0.169	0.129	0.118	0.126	0.177	0.205	0.248
XIII (n = 44)	0.186	0.216	0.178	0.146	0.165	0.140	0.117	0.126	0.167	0.204	0.233	0.112
XIV (n = 56)	0.221	0.260	0.221	0.184	0.191	0.167	0.143	0.155	0.216	0.235	0.283	0.134	0.139
XV (n = 6)	0.145	0.132	0.135	0.109	0.163	0.145	0.112	0.127	0.108	0.159	0.210	0.149	0.144	0.178
XVI (n = 4)	0.173	0.199	0.169	0.129	0.168	0.165	0.169	0.129	0.161	0.189	0.233	0.170	0.166	0.199	0.167
XVII (n = 56)	0.182	0.222	0.187	0.152	0.171	0.151	0.132	0.138	0.174	0.212	0.235	0.124	0.120	0.132	0.155	0.184
XVIII (n = 18)	0.192	0.217	0.185	0.153	0.172	0.139	0.125	0.136	0.176	0.209	0.236	0.119	0.115	0.134	0.148	0.179	0.109
VI m (n = 9)	0.193	0.218	0.185	0.158	0.176	0.112	0.143	0.134	0.182	0.222	0.242	0.148	0.155	0.182	0.155	0.184	0.173	0.157

^a The numbers of base substitutions per site from averaging over all sequence pairs between genotypes are shown. The analysis involved 1430 nucleotide sequences. There were a total of 1596 positions in the final dataset

The remaining eight Pakistani viruses (0%–3.6% genetic distance among themselves, see Additional file 1: Table S5) isolated from pigeons during 2014–2016 did not cluster within any of the previously known sub-genotypes within genotype VI. These Pakistani isolates grouped with a Pakistani virus isolated from pigeon (KU885949/pigeon/Pakistan/MZS-UVAS/2014) creating a separate branch in the phylogenetic tree (Fig. 1) and were 8.5% to 12% distant from the rest of the sub-genotypes in genotype VI (Table 2). This new group of viruses fulfills all classification criteria set by Diel et al. [12] and was named as a novel, hitherto undescribed, sub-genotype of class II genotype VI, namely sub-genotype VIm. The viruses from the newly designated sub-genotype VIm were also more than 10% (11.2–24.2%) distant from all known NDV class II genotypes (Table 3).

The topology of the phylogenetic tree based on the full fusion gene sequences revealed that the Bulgarian virus (Chicken/Bulgaria/Mokresh/1982) grouped together with isolates from Ethiopia (during 2011–2012) [13] (Fig.1). However, the Bulgarian isolate was quite divergent from other members in this sub-genotype showing high genetic distance (10–11%). When this Bulgarian isolate was phylogenetically analyzed based on complete genome, it shared common ancestor with sub-genotypes VIg and VIm described above (see Additional file 2: Fig. S4). The other Bulgarian virus (chicken/Bulgaria/Dolno Linevo/1992), as well as both South Korean viruses obtained from a chicken and a quail, grouped within sub-genotype VIc with other chicken isolates from Sweden, Japan and China indicating that they belong to the older genotype VI sub-genotype (Fig. 1).

Rapid diagnostic and development of new primers and probe