Updated unified phylogenetic classification system and revised nomenclature for Newcastle disease virus

1The International Committee on Taxonomy of Viruses has recently created three genera, named Orthoavulavirus, Metaavulavirus, and Paraavulavirus, within a new subfamily Avulavirinae of the family Paramyxoviridae. Viruses of genus Avian orthoavulavirus 1 (AOAV-1) (formerly designated as Avian avulavirus 1 (AAvV-1)), commonly known as Avian paramyxoviruses 1 (APMV-1) or Newcastle disease viruses (NDV, used hereafter for the purposes of this paper), cause infections in a wide range of domestic and wild birds worldwide. Newcastle disease (ND), caused by virulent NDV, is highly contagious and can be devastating, particularly in immunologically naïve poultry. Virulent strains are defined by the World Organisation for Animal Health (OIE) as viruses that have an intracerebral pathogenicity index of 0.7 or higher (2.0 is maximum) or a fusion cleavage site with multiple basic amino acids and phenylalanine at position 117. Newcastle disease has a global impact and 109 of 200 member countries have reported the disease to the OIE in the last five years. Between 2006 and 2009, ND was ranked the 8th most important wildlife disease and the 3rd most significant poultry disease. The importance of NDV infections to avian health has long been recognized, and the viruses have been the subject of considerable scientific investigations over the past several decades. The need for NDV classification has resulted in a variety of methods for isolate identification and differentiation. These early techniques were mainly based on biological properties of the virus, such as pathogenicity, plaque formation, thermostability, analyses of structural polypeptides, and hemagglutination inhibition patterns using monoclonal antibodies.

2The broad circulation of NDV in poultry populations led to significant genetic diversity of the virus and constant emergence of NDV variants. Given the clinical and economical relevance of NDV to the poultry industry and the broad use of live ND vaccines worldwide, sequencing and phylogenetic analysis became the methods of choice for the characterization of NDV strains circulating in the field. With the incorporation of molecular methods into viral research, the collective knowledge of the NDV genetic makeup has greatly improved. To track NDV evolution and genetic diversity, several NDV molecular classification systems have been developed. In the late 1980s, two related studies based on analyses of the fusion (F) and hemagglutinin-neuraminidase (HN) genes' sequence diversity and HN gene sizes, proposed the existence of three NDV lineages, namely A, B and C. In 1996, Ballági-Pordani et al. suggested the classification of NDV isolates into six distinct genotypes (I to VI) based on restriction fragment length polymorphism analyses. This system was validated and further improved by phylogenetic analysis of partial F gene sequence data, and additional NDV genotypes have since been identified. A subsequent study based on complete genome sequences identified the existence of three different NDV genome sizes (15,186, 15,192, and 15,198 nucleotides) and separated NDV isolates into two classes – class I and class II. Each class was further divided into multiple genotypes. Another classification system, based on partial F gene analysis, was defined by Aldous and colleagues and separated NDV isolates into six distinct lineages (1 to 6) and 13 sub-lineages, with an additional lineage 7 proposed later. These systems pioneered the classification efforts for NDV and provided important information on NDV evolution and genetic diversity. However, each system was based on different approaches and lacked objective criteria for the classification of NDV isolates into different genetic groups. These deficiencies created inconsistencies in NDV classification, made discussions and comparisons among studies difficult, and generated confusion in the assignment of existing and new genetic groups.

3To overcome these challenges, a unified and objective NDV classification system was proposed in 2012. This system utilized the complete F gene coding sequences and incorporated a number of objective criteria for classification of NDV, including: i) phylogenetic topology; ii) inter-populational evolutionary nucleotide distances; iii) branch support; and iv) epidemiological independence of at least four isolates per sub/genotype. The use of system and criteria led to the classification of class I NDV isolates into a single genotype (genotype 1) containing three sub-genotypes, while class II encompassed viruses delineated into 15 genotypes (I to XV) and multiple sub-genotypes. This system was widely adopted, and its use resulted in the identification of three more genotypes (XVI, XVII, and XVIII). In addition, multiple sub-genotypes have been proposed within many of the existing genotypes.

4While the system provided objective criteria for classification of NDV, the ample circulation and constant evolution of NDV led to the almost simultaneous identification and uncoordinated naming of new genotypes. In other instances, new genotypes were created without applying all of the proposed criteria, or the classification was completed using a limited sequence dataset (often using partial F gene sequences) rather than using a curated dataset of all available complete NDV F gene sequences. For these reasons, there are inconsistencies in the naming and classification of some newly proposed genotypes and some sub-genotypes bear identical names while describing different and unrelated viruses. Furthermore, with the increased surveillance efforts and the improvement of sequencing technologies, the amount of available sequences has vastly increased, adding to the already complex challenge of studying the molecular evolution of NDV and the relationships among its isolates.

5To address these issues, a large group of international scientists and collaborators from 29 laboratories from all continents (except Antarctica) and representing all OIE reference laboratories for ND, was established in early 2014. Members of this consortium have convened during the OFFLU meeting in April 2014 in Paris, France, the 9th International Symposium on Avian Influenza in April 2015 in Athens, Georgia, USA, the 2nd International Avian Respiratory Disease Conference in May 2018 in Athens, Georgia, USA, and through multiple teleconferences. Three working groups were assigned, and the consortium was set out to revise the existing NDV classification systems and to propose updated classification and nomenclature criteria for NDV. The work described here was performed with the aims to: i) establish unified criteria for NDV sequence collection and curation; ii) generate and maintain updated sequence datasets; iii) perform comprehensive phylogenetic analyses using the generated datasets; iv) propose unified criteria for NDV classification; v) propose naming criteria for NDV genotypes; and vi) provide curated complete and partial datasets for future public use. Here, in-depth analyses of the genetic diversity of NDV utilizing all available complete F gene sequences are presented. An updated NDV classification system with new nomenclature criteria is proposed. In addition, curated reference datasets for phylogenetic inference, software recommendations, and guidelines for using high-performance computer clusters are provided.

6Three working groups (WG1 to 3) were formed within the consortium. WG1 was tasked with defining criteria for sequence data quality (e.g. size, coverage, recombination, and identical sequences) and generating the datasets that will be used in the analyses. The focus of WG2 was to make the datasets publically available for future application. The WG3 was tasked to evaluate methods for phylogenetic analyses utilizing the generated datasets and to propose unified NDV classification and nomenclature criteria. This WG was also tasked with selecting representative sequences for a smaller class II dataset for rapid genotype identification of new isolates (named the “pilot” dataset).

7Due to its variability and biological function (virulence determinant), the F gene is commonly targeted in sequencing efforts. Accordingly, the complete F gene coding nucleotide sequences of all available class I and class II NDV isolates were downloaded from GenBank of the National Center for Biotechnology Information. Only sequences that represented >99% of the full-length open reading frame of the F gene (≥1645 nucleotides, nt) were selected and aligned using Multiple Alignment with Fast Fourier Transformation (MAFFT v7.221.3) as implemented in the Galaxy platform. A total of 2966 sequences that had submission dates up to February 21st 2019 were obtained, resulting in two separate initial datasets – class I (n = 431) and class II (n = 2535). Similarly, all available complete genome sequences of class II NDV were downloaded from GenBank and aligned using the same tools (due to the scarcity of class I complete genomes, these were not analyzed). The leader and tail sequences and intergenic regions were trimmed, and the coding sequences for all six genes were concatenated. A total of 651 class II complete genome sequences submitted on or before Feb 21st 2019, were obtained. Metadata like host, year and location of isolation, and isolate name for all sequences were collected from GenBank annotations, where available, or from the respective peer-reviewed publications when the data was not provided with the sequence submission.

8Sequences identified as 100% identical and originating from the same species, country, laboratory, and/or outbreak or having identical names and dates were removed from the datasets (one representative from each duplication pair/group was kept). Each sequence was evaluated and all gaps and insertions that caused alignment shifts were deleted. Sequences that were found to originate from man-made clones, mutant viruses, vaccines, and chimeric viruses were removed from the datasets. Sequences obtained from viruses identified as spillovers or escapes (recent viruses identical or almost identical to vaccine and standard strains originally isolated decades before) were also excluded from further analyses. In addition, all sequences were subjected to recombination analysis using the RDP4 program as described previously, and sequences with recombination events were removed from the datasets.

9Three Nigerian and five Bulgarian historical viruses from underrepresented genotypes were obtained from the repository of the Southeast Poultry Research Laboratory of the USDA and were sequenced in this study: (chicken/Nigeria/Kano/N52/899/1973, chicken/Nigeria/Plateau/N53/900/1973, chicken/Nigeria/FLD/N54/901/1973, pullet/Bulgaria/Plovdiv/1153/1959, chicken/Bulgaria/ElovDol/1156/1981, pigeon/Bulgaria/Septemvriitsi/1157/1982, chicken/Bulgaria/Furen/1159/1988, pigeon/Bulgaria/NovoSelo/1161/1995). RNA extraction, next-generation sequencing, and genome assembly were performed as described previously.

10A total of 1956 complete F gene sequences remained after the curation of the datasets and were used for the phylogenetic analyses – class I (n = 284, Supplemental Table S1) and class II (n = 1672, Supplemental Table S2, “larger” class II dataset). Neighbor-joining (NJ), maximum-likelihood (ML), and Bayesian methods were used to infer phylogenetic relationships and construct phylogenetic trees of both classes. Neighbor-joining trees (Maximum Composite Likelihood model) were constructed using MEGA6 with 1000 bootstrap replicates. Codon positions included were 1st, 2nd, 3rd, and noncoding, and all positions containing gaps and missing data were eliminated. Maximum-likelihood trees based on general time-reversible (GTR) model (goodness-of-fit measured by the corrected Akaike information criterion) were constructed by using RaxML version 8.2.11 with 1000 bootstrap replicates and the following command: -s input_filename.phy -n output_tree_name -m GTRGAMMAI -f a -x 123 -N 1000 -p 456. A discrete Gamma distribution (Γ) was used to model evolutionary rate differences among sites and the rate variation model allowed for some sites to be evolutionarily invariable (I). Utilizing the same datasets and parameters, RaxML trees were also constructed through the CIPRES Science Gateway. A class II complete genome tree was constructed using the same ML method, utilizing the curated dataset of concatenated coding sequences (n = 491). Class I and class II complete F gene maximum-likelihood trees based on the GTR model were also constructed using MEGA6. Bayesian analyses were performed for comparison using the BEAGLE-enhanced parallel version of MrBayes v3.2. For both class I and II datasets, the GTR substitution model with a gamma-distributed variation of rates and a proportion of invariable sites were used. All other parameters were set to the MrBayes default settings. Independent runs of a 10 million generation length Markov Chain Monte Carlo (MCMC) simulation were performed sampling every 1000 states. Run convergence was assessed using Tracer v1.6 (http://tree.bio.ed.ac.uk/software/tracer) to ensure an effective sampling size larger than 200 for all parameters. After discarding non-convergent runs and removal of at least 10% burn-in, the runs were combined with MrBayes to produce a maximum clade credibility tree and Bayesian posterior probabilities (BPP) were estimated. Trees were visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree). Topological congruence between trees was compared through visual inspection for each sub/genotype.

11The final class I and class II datasets used for phylogenetic reconstruction were also used to estimate the mean inter-populational distances (nucleotide distances between sub/genotypes). The estimates of average evolutionary distances were inferred using MEGA6. Analyses were conducted using the Maximum Composite Likelihood model. The rate variation among sites was modeled with a gamma distribution (shape parameter = 1).

12After performing the analyses described above and identifying the main NDV genotypes, a smaller (pilot) dataset of class II sequences (n = 125, Supplemental Table S3, “smaller” class II dataset) was parsed from the complete dataset used in the phylogenetic analyses. This pilot dataset included representative sequences from all identified sub/genotypes and was phylogenetically analyzed with the ML method as described above to ensure that the topology inferred using the larger class II dataset would be maintained if fewer isolates were used. To explore the option to use sub-trees for genotype identification of new isolates, sub-trees were built separately, utilizing the ML method and using all sequences in each class II genotype. To ensure congruent topology with the biggest class II tree, these sub-trees were rooted to historical isolates ancestral to the respective genotypes.

13In all analyses, the criteria put forth by were used to initially assess the diversity of NDV and served as the foundation to update and propose a consensus classification system. The criteria accepted here, include the previously proposed criteria of clustering of sequences from one sub/genotype into a monophyletic branch, average nucleotide distance among genotypes (above 10%), and the presence in each group of at least four independent viruses from distinct outbreaks and without a direct epidemiologic link. To avoid excessive delineation of sub-genotypes, the cut-off value for nucleotide distance between these was increased to 5% (instead of 3%). In order to increase the reliability and stringency of the phylogenetic clustering, the bootstrap support value cut-off at the nodes that define sub/genotypes was raised to 70% or above. The updated consensus classification criteria are summarized in Table 1. Similar classification criteria have also been used by the World Health Organization/World Organisation for Animal Health/Food and Agriculture Organization H5N1 Evolution Working Group.

14When possible, names of genotypes were maintained with Arabic numerals used in class I and Roman numerals used in class II. The lowercase Latin letter system used to name sub-genotypes was replaced with a numerical-decimal system using Arabic numerals. Dichotomous splitting was used at every defining node (at which separation into sub-genotypes was done) using the numerals 1 and 2. It is assumed that each node has two immediate descendants (two higher order sub-genotypes). For example, within genotype VII (names of sub/genotypes used in this paragraph do not correspond to the classification proposed in this study and are used solely for demonstration purpose), the first two sub-genotypes will become VII.1 and VII.2 (see Supplemental Fig. S1 for examples). At the next node (closer to the tips or higher order) within sub-genotype VII.1, the sub-genotypes that are one order higher will become VII.1.1 and VII.1.2, and respectively at the next node within sub-genotype VII.2, the two further sub-genotypes will be named VII.2.1 and VII.2.2. The updated nomenclature criteria are presented in Table 2 and examples are provided in Supplemental Fig. S1, Supplemental Fig. S2, and Supplemental Fig. S3.

15The complete F-gene sequences (n = 8) of NDV obtained in this study were submitted to GenBank and are available under the accession numbers MH996897 to MH996904.

16A total of 1956 complete F gene sequences remained after the curation of the datasets, including 284 class I sequences and 1672 class II sequences. After utilizing the dataset curation criteria described above, 1010 sequences were removed from the initial datasets (class I: n = 147; class II: n = 863). In the class I dataset, 128 duplicate sequences, 17 recombinant sequences and 2 sequences that resulted from cloning were identified and removed. In the class II dataset, 419 duplicates, 170 recombinant forms, 56 man-made vaccine clones and 218 spillovers or escapes of vaccine/standard viruses were identified and removed. These numbers highlight the importance of carefully assessing the NDV datasets prior to phylogenetic inference. It has been reported, for example, that spillovers of vaccine and of older standard NDV strains can significantly impact analyses, resulting in skewed inferences. In addition, the increasing number of chimeric sequences as identified by the recombination analysis has implications in virus classification and impacts topology of phylogenetic trees. Although recombination has been suggested to play a role in NDV evolution and increasing numbers of chimeric sequences have been deposited in GenBank, the natural occurrence of recombination events remains disputable. Selection and the inherent error rate of the viral RNA polymerase are believed to be the main forces driving the evolution and diversity of NDV. The actual role and contribution of recombination to this process remains to be established. However, to provide a reliable phylogenetic inference based on clearly identified genotypes, we suggest that recombinant sequences always be identified through analysis with RDP (or alternative appropriate tool) and removed from the dataset before performing an analysis and a classification tree construction. A similar approach has already been adopted for porcine circovirus type 2 and infectious bronchitis virus.

17The phylogenetic analyses utilizing the NJ, ML, and Bayesian methods revealed high topology congruence within class I NDV (Supplemental Fig. S4, A-C). All class I trees were rooted to the oldest class I NDV isolate - EF564833/Canada Goose/USA(OH)/78/1987. The existence of a single genotype (i.e. genotype 1) was confirmed as previously determined by. Minor topological differences among a few unclassified sequences and lack of branch support (<70 bootstrap values) were observed when using the NJ method; however, the existence of a single class I genotype 1 was consistent between the used methods. Based on the ML method (selected due to easier and faster tree construction and because it outperformed NJ), there are three (n = 3) sub-genotypes within genotype 1 of class I (Supplemental Fig. S4B). Utilizing the new naming criteria, these are sub-genotypes 1.1.1 (former 1a), 1.1.2 (former 1b), and 1.2 (former 1c and 1d) (Supplemental Table S1, Table 3). For visualization purposes and to provide an overview of the topology of class I NDV, a condensed ML tree without taxa names is featured in Fig. 1A. The increase of the cut-off for differentiation of sub-genotypes to 5% resulted in the merging of former sub-genotypes 1c and 1d into one single sub-genotype 1.2 (distance between 1c and 1d - 3.68%) (Supplemental Table S4). Of note, there are five unclassified viruses within class I (Supplemental Fig. S4B, UNCL 9 to 13) that cluster together but are genetically distant from each other (9 to 11%) and have the potential to become members of multiple separate genotypes if sequences from additional epidemiologically independent viruses related to them become available. It is also possible that these sequences represent extinct viral lineages, which does not warrant their further classification into genotypes. These observations highlight the importance of having at least four epidemiologically independent sequences in order to classify NDV isolates into a new sub/genotype. Although there is evidence of ongoing evolution among the viruses of class I NDV, the evolutionary diversity within this class remains low. The majority of sequenced class I viruses are isolated from wild birds, however, poultry detections are not uncommon. The circulation of class I viruses mainly in wild birds, which are not vaccinated for ND, thus keeping the immune pressure on these viruses low, may explain the relatively lower genetic diversity of these viruses when compared to class II. While there has been one report of a virulent class I NDV, all remaining sequences in the class are from viruses of low virulence. It is likely that the low virulence nature of class I viruses in wild birds and poultry and their low incidence in poultry, has negatively influenced efforts to sequence and characterize them. This may explain the restricted spatial distribution of the majority of the available sequences (the U.S. and China, and a few sequences from Japan and Europe). In addition, not all RT-PCR assays are able to detect class I isolates and their circulation may have been missed as they do not commonly cause clinical disease. Additional sequencing of historical and prospective class I NDV isolates will provide a more detailed characterization of the genetic diversity of this group of viruses.

18The class II phylogenetic trees, constructed using the NJ, ML, and Bayesian methods, presented some significant differences, particularly between NJ (Supplemental Fig. S5A) and the other two methods (ML and Bayesian) (Supplemental Fig. S5B and S5C, respectively). Notably, lack of branch support was observed at many defining nodes of the NJ tree, and the identified monophyletic branches differed compared to the remaining two methods (e.g. in genotypes I, VI, and XIII) (Supplemental Fig. S5, A–C). With the exception of one small branch (former sub-genotype Vd) and low support between two genotype I sub-genotypes, the trees obtained using the ML and the Bayesian approach were generally consistent with well-defined monophyletic branches and strong branch support (≥70 bootstrap and ≥ 0.9 posterior probability values, respectively) at most defining nodes (Supplemental Fig. S5B and S5C, Supplemental Table S2). Since it statistically outperformed NJ, and due to its ease of use, greater speed and wide utilization and accessibility, the ML method was selected as a method of choice for class II NDV phylogenetic inference. Compared to class I, class II is more diverse, contains a range of non-virulent to virulent viruses, and the complete analyses identified at least 20 distinct genotypes (I to XXI, genotype XV that contains only recombinant sequences was excluded from the final analyses). Viruses previously classified as members of sub-genotype Va were separated into a new genotype XIX. Viruses previously classified as five of the sub-genotypes of genotype VI also formed two new genotypes, namely XX and XXI. All class II trees were rooted to the clade of genotype IV, consistent with molecular clock trees. The evolutionary distances between genotypes are presented in Table 4.

19Consistent with previous classifications, genotypes I, V, VI, VII, XII, XIII, XIV, and XVIII were confirmed to be further divided into sub-genotypes (Supplemental Fig. S5B, Table 3). The number of existing sub-genotypes in genotypes I, XIV, and XVIII, remained the same (four, two, and two, respectively), and their new names were assigned based on the updated nomenclature criteria (see Table 5, Supplemental Table S2, Supplemental Fig. S5B). No sub-genotypes were identified in genotypes II, III, IV, VIII, IX, XI, XVI, XIX, and XX. Sequences from eight viruses were not assigned to any of the identified genotypes (“UNCL” 1 to 8, see Supplemental Table S2). For an overview purposes and to provide an overall picture of the topology of class II NDV, a condensed phylogenetic tree without taxa names is presented in Fig. 1B. Utilizing the higher nucleotide distance cut-off of 5%, however, resulted in no sub-genotypes in genotypes X and XVII (distances of 4.83 and 4.59, respectively, Supplemental Table S4).

20The use of the updated classification and nomenclature criteria led to the identification of two sub-genotypes in genotype V, designated V.1 (former Vb) and V.2 (former Vc). Previously identified sub-genotype Vd lacked branch support and sufficient number of independent isolates and, following the newly proposed classification and naming criteria, was assigned to the lower order (i.e. genotype V). Eight more isolates from Europe and South and North America from the 1970s and 1980s could not be assigned to any of the genotypes (lack of branch support and branches not monophyletic) and were also designated as members of the lower order genotype V (Supplemental Fig. S5B).

21Undoubtedly, genotype VI is the most diverse among all NDV genotypes. While the lineage system separated this group into four sub-lineages with multiple further groups, 14 sub-genotypes were identified based on Diel's system. Although viruses classified into two new genotypes (XX and XXI) were taken out of genotype VI, after utilizing the updated criteria suggested here, this genotype was still found to contain seven sub-genotypes (Supplemental Fig. S5B). Former sub-genotypes VIa and VIn were merged into sub-genotype VI.2.1.1.1, as they did not pass the 5% criterion for nucleotide distance between sub-genotypes (3.46%, Supplemental Table S4). Previously identified sub-genotypes VIj and VIk were just above the 5% cut-off (Table 3) and were named VI.2.1.1.2.1 and VI.2.1.1.2.2, respectively. Former sub-genotypes VIb, VIe, VIf, VIh were confirmed and renamed accordingly (Table 5, Supplemental Table S2). Considering that many genotype VI NDV have been part of the panzootic in Columbiform birds that emerged in the Middle East more than four decades ago, and is still ongoing, it is not surprising that this is the most diverse group of NDV. Unlike any other NDV genotype, genotype VI viruses have been isolated in all continents, except Antarctica. Their potential to infect chickens has been demonstrated, but only rarely, and they appear to be highly adapted to some Columbiform birds. Viruses that belong to genotype VI are often referred to as pigeon paramyxoviruses 1 (PPMV-1). These viruses can be distinguished as unique antigenic variants of NDV by the patterns produced in a hemagglutination inhibition (HI) assay when tested with a panel of monoclonal antibodies. This panel, although described in the OIE manual, is not widely available and used by NDV reference laboratories. Not all genotype VI viruses have been confirmed to be PPMV-1.

22The use of the updated classification criteria significantly impacted genotype VII landscape. As many former sub-genotypes did not fulfill the distance (VIIb, VIId, VIIe, VIIj, VIIl), branch support (VIIh, VIIi, VIIk), and/or number of independent isolates (VIIk) criteria, these were merged into single sub-genotypes, resulting in a total of three genotype VII sub-genotypes (Table 3, Table 5, Supplemental Fig. S5B, Supplemental Table S4). The viruses responsible for the fourth NDV panzootic grouped together, and based on nucleotide distance, were classified into a single genotype (VII.1.1), combining former sub-genotypes VIIb, VIId, VIIe, VIIj, and VIIl. An exception is former sub-genotype VIIf that was classified as a separate sub-genotype, namely VII.1.2. The groups of viruses involved in the fifth NDV panzootic (VIIh and VIIi), that affected Indonesia, Asia, the Middle East, Europe, and Africa, lacked branch support, and were merged into a single sub-genotype VII.2. A group of five other sequences from viruses isolated in Namibia in 2016, previously identified as sub-genotype VIIk, were also assigned to VII.2 due to lack of branch support and insufficient number of independent isolates (Supplemental Fig. S5B). Interestingly, all these sub-genotype VII.2 groups met the nucleotide distance criterion with high genetic distances of 9.40% and 11.17% (Supplemental Table S4); however, all classification criteria need to be fulfilled before a clade is named as a separate sub-genotype. Of note, this branch (VII.2) including these three groups of viruses (former VIIh, VIIi, and VIIk) is almost 10% (9.83%, Table 3) distant from the remaining genotype VII viruses, and with the continuing evolution of NDV, will likely fulfill the criteria for a standalone genotype in the near future.

23Genotype XIII was found to consist of four sub-genotypes (Supplemental Fig. S5B). Two previously assigned sub-genotypes, were further split into two sub-genotypes each – XIIIa into XIII.1.1 and XIII.1.2, and XIIIb into XIII.2.1 and XIII.2.2, respectively (Supplemental Fig. S5B) because of the addition of new sequences from recent studies. Viruses isolated from different countries in Africa between 1995 and 2015 formed sub-genotype XIII.1.1, while viruses from Sweden, Russia, India, and Iran from 1997 to 2011 were classified in XIII.1.2. Viruses isolated in the last decade from Pakistan and India formed XIII.2.1 and XIII.2.2, respectively. Although very distant from the XIII.2.1 and XIII.2.2 sequences (12.63%), a branch of nine sequences from viruses from India lacked node support and number of independent isolates to be identified as separate sub-genotype and were assigned to the lower order branch (i.e. XIII.2). Similar to the viruses from genotype VII/lineage 5 (to which they were previously assigned), genotype XIII viruses are widely distributed on at least three continents. High genetic distances were observed between clades of genotype XIII (Table 3), and the continuous NDV evolution and addition of new genetic data will likely change the classification of this genotype.

24Genotype XIX is one of the three newly identified genotypes. This genotype contains all viruses previously classified as sub-genotype Va. Genotype XIX viruses are almost exclusively associated with outbreaks in double-crested cormorants in North America and Canada and additional isolates have been recovered from samples taken from pelicans and gulls in close proximity to cormorants during outbreaks.

25Genotype XX contains some of the oldest available NDV isolates, previously identified as members of genotype VI (former sub-genotype VIc). Notably, genotype XX outgroups genotype VI, and all genotype XX viruses have been isolated from chickens. This new genotype XX includes viruses from Asia and Europe isolated between 1985 and 2011.

26Genotype XXI contains viruses isolated from chickens and pigeons in different Asian, European, and African countries between 2005 and 2016. Similarly to genotype XX, genotype XXI viruses were previously assigned as members of genotype VI. At the root of genotype XXI is a clade of chicken viruses, isolated in Ethiopia between 2011 and 2012, that have previously been assumed to form a separate sub-genotype VIl. However, as the group lacked a sufficient number of independent isolates, these viruses were assigned to the lower order, namely XXI. Three sub-genotypes were identified within genotype XXI – XXI.2, XXI.1.1, and XXI.1.2 corresponding to former sub-genotypes VIi, VIg, and VIm, respectively (see Table 5, Supplemental Fig. S5B).

27Construction of the class II phylogenetic trees with any of the used methods and utilizing all available complete fusion gene sequences occasionally resulted in lack of branch support, presence of polytomies, and disruption of monophyletic branches (i.e. previously identified genotypes not grouping into a single branch) (Supplemental Fig. S5 A–C). For example, viruses of genotype IV did not fall into a monophyletic branch, there was lack of branch support between groups in genotype VII, and polytomies were present in different branches. These issues will likely be resolved when additional sequences from these groups become available. Indeed, the analysis using the concatenated complete genome coding sequences overcame the flaws observed when the complete fusion gene sequences were used for the analyses (Supplemental Fig. S6). Despite the fact that complete genomes provide more reliable phylogenetic inferences, the lack of complete genome sequences in all NDV sub/genotypes and the complexity of the analysis, the use of complete genome phylogenetic reconstruction to resolve NDV classification issues or for definitive classification, is not warranted at this time. Certainly, with the increasing use of high-throughput sequencing technologies and the development of more powerful computational resources, it is logical to suggest that the availability of complete genomes will increase in the near future. This will probably allow this classification to be revisited with more complete genomic data. In addition, genome size is no longer needed as a criterion for separation of NDV isolates into classes as there are sequences in genotype XIX of class II that have the same genome size as class I NDV (15,198 nucleotides). As the complete fusion gene mean nucleotide distance between the two classes is 59.25% and the individual distances between all class II genotypes and class I is above 40%, a cut-off of 40% nucleotide distance of the F gene between NDV classes is suggested here.

28Building a phylogenetic tree with a lower number of sequences, that still provides topology consistent to the larger tree, significantly decreases the time for the tree construction. To ensure that the topology of the class II genotypes would be maintained if fewer isolates were used, we performed a ML phylogenetic reconstruction with a smaller “pilot” dataset (n = 125, Supplemental Table S3). The generated tree (named “pilot” tree) is depicted in Fig. 2 and the congruence analysis showed complete consistency with the ML tree inferred using the larger class II fusion gene dataset (Supplementary Fig. S5B). Therefore, this pilot dataset can be utilized for rapid preliminary genotype identification of new isolates. A similar approach has been previously suggested by Aldous et al.. However, for naming new genotypes and sub-genotypes, building a tree using the complete curated F gene dataset is required.

29To further facilitate fast characterization of newly isolated NDVs, sub-trees utilizing all sequences within each genotype of class II that contains sub-genotypes, were constructed (Supplemental Fig. S7 A–I). To maintain consistent topology with the larger class II tree and to provide inference of branch support at the first defining node, each sub-tree was rooted with at least two sequences. A list of sequences to be used for rooting of each genotypes' sub-tree is provided in Supplemental Table S5. Proper sub-tree rooting appeared to be crucial for consistent topology between the larger tree and the sub-trees. Attempts to root the sub-trees with random sequences from different genotypes resulted in topologies largely incongruent to the full tree (data not shown). In addition, lower branch support was observed in some trees and for definitive classification, building a tree with the complete F gene dataset is recommended.

30Assuming that each node has two immediate descendants (partitions into two branches), a dichotomous naming system using the Arabic numerals 1 and 2 was used at every node at which sub-genotypes were separated. Arguably, strict following of the classification and nomenclature criteria proposed here, will prevent duplicate names of future sub-genotypes, as has happened before with the previous systems. In addition, when using this approach, the names of the sub-genotypes of higher order (closer to the tips), bear ancestral information for the sub/genotypes they split from. Of note, the use of the new naming criteria resulted in complex names for two of the genotype VI sub-genotypes – VI.2.1.1.2.1 and VI.2.1.1.2.2. However, the viruses from these two genotypes are not widely distributed, and their names would likely not have negative impact on the global use of the system. Similar decimal systems are used for classification of swine influenza and highly pathogenic H5N1 avian influenza viruses. A detailed comparison between the nomenclature used in the lineage system, the genotype system and the nomenclature suggested here, is provided in Table 5.

31The worldwide circulation and constant evolution of NDV will, expectedly, lead to emergence of new NDV genetic variants that will form clades fulfilling all criteria for separate genotypes. To ensure consistent naming, specific nomenclature criteria need to be utilized for newly identified genetic diversity. Updated guidelines for naming new genotypes are also outlined in Table 2. Newly identified virus diversity (a group of viruses undescribed before), that meets all classification criteria, will be classified as new genotype and will receive a subsequent Roman numeral (e.g. currently XXII is the next available) (example in Supplemental Fig. S2, red color). If a new sub-genotype of viruses is identified, but still falls into a higher order (not at the first defining node) within an existing genotype, this new sub-genotype receives the next consecutive sub-genotype numerical address for the respective level of order in the phylogeny. This is necessary to avoid re-numbering of all existing sub-genotypes that are of higher order (closer to the tip). For example, if a new sub-genotype that outgroups VII.1.1 and VII.1.2 is identified but it is still within sub-genotype VII.1, this new sub-genotype will be named VII.1.3 (see Supplemental Fig. S3 for an example). Existing sub-genotypes that fulfill the genotype criteria will not be designated with a new genotype name and will continue bearing the name assigned in the current classification (examples in Supplemental Fig. S2, green color).

32The increased number of available NDV sequences has led to a higher demand of computational resources. To this end, in addition to the widely used MEGA, and PhyML software, for example, tools like FastTree, Garli, and RaxML, have been extensively used for faster phylogenetic inferences. However, some of these utilize the Shimodaira-Hasegawa test to estimate nodal support (FastTree, PhyML), which requires multiple additional steps to calculate the traditional bootstrap value, commonly used in NDV phylogenetic analyses. Based on its increased productivity and reliable results, the RaxML tool (v.8.2.11) was utilized for ML inference in this study. Interestingly, when the invariable sites option is used with RaxML, the tool automatically sets the model optimization precision (in likelihood units) to 0.001 to avoid unfavorable effects caused by simultaneous optimization of gamma distribution and invariable sites. Although the combination GTR+ Γ + I has been widely used, it is suggested that distinct approaches to incorporate rate heterogeneity (e.g. Γ + I) should not be used at the same time and concerns have been raised that this might cause problems during the model parameter optimization process. Indeed, GTR+ Γ + I is often estimated as the best-fit model for NDV datasets based on the corrected Akaike information criterion. However, in most cases, trees constructed with this model using the MEGA software lack branch support at established defining nodes or contain multiple polytomies. Disabling the option for invariable sites in MEGA results in better branch support and resolving most polytomies (data not shown). Goodness-of-fit of models has to be established periodically in the future, as adding new isolates to the datasets may result in change of the best-fit model.

33The advances in molecular sequencing technologies, the increasing surveillance efforts, and the significantly larger amount of genetic data collected worldwide, have transformed the field of phylogenetic inference into a computational science. Therefore, faster tools, specifically ones that can be used in high-performance computer (HPC) environments, are needed. In this study, to test the utilization of HPC for NDV phylogenetic inference, the class II ML tree using the bigger complete F gene dataset was built again through the free CIPRES Science Gateway using RaxML-HPC (v.8.2.10) that allows multi-threading by message-passing interface. The use of the supercomputer shortened the construction of a fully congruent tree to less than two hours (compared to two days when using 30 cores, and two weeks when using a single computer). Step-by-step guidelines with all settings needed for ML tree construction through CIPRES (the same settings can be used for any of the freely available HPC) are provided in Supplemental Fig. S8 (provided as screenshots from http://www.phylo.org). Additional guidelines and video tutorials regarding the use of the CIPRES Science Gateway are available online (http://www.phylo.org/index.php/help/).

34The nucleotide distance between genotypes III and IX was estimated to be below 10% (8.62%) (Table 3) with topology not fully resolved. Some sequences, previously identified as members of genotype IV, did not form a monophyletic branch (Supplemental Fig. S5B). However, genotypes III and IV contain historical viruses that are believed to no longer naturally circulate. For consistency with the existent literature, the nomenclature and classification of these viruses were maintained. In addition, the genotype IV viruses from Nigeria could have evolved into a separate genotype as they formed a separate branch, but have probably become extinct as similar viruses have not been isolated since 1980.

35To facilitate future analyses, all complete fusion gene datasets, generated in this study, are provided as Supplemental materials – see Supplemental Table S1, Supplemental Table S2, Supplemental Table S3. The resulting alignments are also provided as Supplemental materials (Supplemental Table S1, Supplemental Table S2, Supplemental Table S3). Updated datasets will be deposited at the GitHub repository (https://github.com/NDVconsortium/NDV_Sequence_Datasets). The utilization of these alignments is highly recommended for NDV classification needs.

36The classification and nomenclature system proposed here is a result of the collaboration of 40 experts from 29 institutions on six continents, including all OIE NDV reference laboratories. The revised system includes objective criteria for identification of new NDV isolates and naming NDV sub/genotypes. Curated, up-to-date, complete F gene class I and class II NDV sequence datasets are provided for public use. In addition, a pilot dataset and rooting guidelines that allow rapid preliminary genotype identification of new isolates are included. To increase the speed of phylogenetic inference and to enable the use of consistent methods between laboratories, detailed guidelines for HPC use are also provided. In order to avoid super-delineation, more stringent criteria were proposed, and these resulted in the merging of multiple former sub-genotypes. However, the increased viral diversity due to constant NDV evolution and the significant increase of available sequences led to the naming of three new genotypes. The global adoption of the proposed consensus system will facilitate future studies on NDV evolution and epidemiology, and will make comparisons of results obtained across the world easier. Newcastle disease viruses, as any other RNA viruses, are rapidly evolving, and periodic (e.g. five years) revisiting of the consensus classification and nomenclature system proposed here, led by a working group from this consortium, is warranted.

37The following are the supplementary data related to this article.