Virus Detection and Subtyping

We extracted influenza virus RNA from tracheal or cloacal swabs from dead wild birds by using MagNa Pure 96 (Roche, Basel, Switzerland). For commercial poultry farms, we used pools of 5 samples from clinically affected animals. We tested samples by using a matrix-gene real-time PCR, which detects all avian influenza virus subtypes, as described previously (9). We subtyped positive samples by using an H5-specific real-time PCR as recommended by the European Union reference laboratory (11). We determined the sequence of the HA cleavage site and the N subtype by using Sanger sequencing (9).

Complete Genome Sequencing and Analysis

We purified influenza virus RNA by using the High Pure Viral RNA kit (Roche), amplified the RNA by using universal 8-segment primers, and then directly sequenced the RNA, as described previously (12). We sequenced purified amplicons at high coverage (average >1,000/nt position) using the Nextera XT DNA Library Preparation kit (Illumina, San Diego, CA, USA) and MiSeq paired-end 150-bp sequencing (Illumina). We mapped the reads by using the ViralProfiler-Workflow, an extension of the CLC Genomics Workbench (QIAGEN, Hilden, Germany). We generated consensus sequences by using a reference-based method. For generation of a defined set of influenza reference sequences, we clustered all complete genome sequences of Eurasia isolates in the GISAID EpiFlu database (https://www.gisaid.org) (years 2000–2016) per segment at 85% identity by using the CD-HIT-EST algorithm (13). We selected cluster representatives for each gene segment to generate a reference set. We first mapped reads to this reference set, and subsequently remapped them to the selected reference sequence. Finally, we extracted the consensus sequence of the complete virus genome. The ViralProfiler-Workflow also assigns the HA cleavage site and coding region of the proteins. We submitted sequences generated in this study to the GISAID database (Technical Appendix 1 Table 1).

Phylogenetic Analysis

We performed phylogenetic analysis of the complete genome sequences for each genome segment separately. We generated alignments with the closest related viruses obtained from the GISAID database (Technical Appendix 1 Table 2) by using Muscle in MEGA6 (14). We performed molecular phylogenetic analysis by using the neighbor-joining method. We modeled the rate variation among sites with gamma distribution (shape parameter = 1) and used 1,000 bootstrap replicates to estimate branch support.

We manually concatenated the 8 gene segment alignments to generate a single alignment, which we used to construct phylogenetic networks by applying the median-joining method implemented in NETWORK as described (15). This model-free method uses a parsimony approach, based on pairwise differences, to connect each sequence to its closest neighbor. It allows creation of internal nodes (median vectors), which could be interpreted as unsampled or extinct ancestral genotypes to link the existing genotypes in the most parsimonious way.

Molecular Dating

For each gene segment alignment, the simplest evolutionary model fitting the dataset was the Hasegaw–Kishino–Yano model with gammadistributed rates (16). We estimated nucleotide substitution rates by using Bayesian Markov chain Monte Carlo methods (17) until all parameters converged (chain lengths of 30 million delivered estimated sample sizes >500). We performed the analysis by using BEAST version 1.8.4 software (18) and strict or relaxed uncorrelated molecular clocks that we calibrated by using the sample isolation dates. We treated all gene segments separately to identify reassortants. We sampled the 30 million generations every 3,000 generations. We constructed a maximum clade credibility tree to summarize all 10,000 trees after removing the initial 10% burn-in by using TreeAnnotator (17). We visualized time-scaled phylogenetic trees by using FigTree version 1.4.2 (19) and further annotated the trees by using iTOL version 3.5.3 (20).

Detection of HPAI H5N8 in Wild Birds and Poultry

Figure 1

Figure 1. Geographic distribution of wild and captive birds and commercial poultry infected with highly pathogenic avian influenza A virus, the Netherlands, 2016. A) Location of dead wild birds infected with highly pathogenic...

The 2016 H5N8 epidemic in the Netherlands started with occurrence of hundreds of dead tufted ducks (Aythya fuligula) around a lake near Monnickendam, in the central part of the country. These tufted ducks tested positive for H5N8 on November 10, 2016. Sequencing of the HA gene showed the polybasic amino acid sequence PLREKRRKR*GLF in the cleavage site, confirming its high pathogenicity. As the epidemic continued, wild birds around several lakes and water-rich areas in the central part of the country, and later also in the northern part, were infected. We plotted the geographic locations where the HPAI H5N8 virus was detected (Figure 1, panel A). Tufted ducks were the most affected bird species in November and Eurasian wigeons (Anas penelope) in December (21). However, many species of wild birds became infected, including diving and dabbling ducks as well as mute swans, grebes, gulls, and buzzards (21). An HPAI virus of subtype H5N5 was detected in 2 birds, a tufted duck found dead near Monnickendam, and a mute swan near Groningen.

The first introduction of H5N8 into a commercial poultry holding with Peking ducks was detected on November 25, in the municipality of Biddinghuizen. The farm was located in the central part of the country, a few kilometers from the lake where dead H5N8-positive wild birds were detected earlier. Severe clinical signs, such as lethargy, neurologic symptoms, and sudden death, were observed at the farm. Samples from the ducks tested positive by PCR, and the HPAI cleavage site (PLREKRRKR*GLF) was confirmed by sequencing. Several days later, samples from 2 other duck farms located within a distance of 3 km tested positive for H5N8. After these initial outbreaks, 5 other poultry holdings were affected in the central and northern parts of the country: 4 farms with laying hens, in Abbega, Boven-Leeuwen, Hiaure, and Zoeterwoude; and 1 duck farm, in Kamperveen. In addition, 4 introductions into captive birds were detected (1 black swan, 2 ducks, and 1 chicken). We noted the geographic locations of the infected commercial farms and captive bird holdings (Figure 1, panel B). All farms, with the exception of Boven-Leeuwen, were located in water-rich areas where dead H5N8-postive wild birds were found.

Genetic Differences between Wild Bird and Poultry H5N8 Viruses

We performed deep sequencing to determine the complete genome sequence of 44 wild bird viruses, the viruses from 8 commercial poultry farms, and viruses from 4 captive birds in the Netherlands in 2016. We studied the genetic relationship between these viruses by median-joining network analysis, which showed that the viruses divide into 3 major clusters (Technical Appendix 1 Figure 1). Cluster A contained mainly Eurasian wigeons, whereas various wild bird species were found in clusters B and C (Technical Appendix 2 Figure 1). The viruses detected at the commercial poultry farms and in captive birds were located in clusters A and B, with the exception of those from Abbega. The network analysis showed no correlation between virus cluster and the geographic location of the virus. These results suggest that at least 3 separate introductions of the H5N8 virus occurred in the Netherlands, after which the virus likely spread to different parts of the country. In addition to the major clusters, we detected multiple separate viruses and small clusters, which suggests that more independent introductions of the H5N8 virus have occurred during the epidemic.

Figure 2

Figure 2. Median-joining network showing the genetic relationship between highly pathogenic avian influenza A viruses subtype H5N8 isolated from commercial poultry farms (red circles) and the most identical wild bird viruses (blue circles)...

To analyze the genetic relationship between H5N8 viruses detected at commercial poultry farms and in wild birds, we identified the most identical wild bird virus for each farm (Technical Appendix 1 Table 2) and performed median-joining network analysis (Figure 2). The geographic locations of the poultry farms and the dead wild birds infected with the most related virus are shown (Figure 1, panel C). Epidemiologic investigation of the affected farms revealed no dangerous contacts. However, the 3 infected farms in Biddinghuizen were detected within a period of 7 days and are located within 3 km of one another. The median-joining network (Figure 2) showed that the complete genome sequence of the viruses from the farms in Biddinghuizen are highly similar, containing at most 5 nt changes. This fact indicates that the farms were either infected by direct transmission between farms or by separate introduction of virus from the same outside source. The viruses detected at the other 5 farms differed between 14 and 59 nt and are likely the result of separate introductions from wild birds (Figure 2). Closely related wild bird viruses with 5–9 nt changes were identified near the farms in Biddinghuizen, Hiaure, Kamperveen, and Zoeterwoude at 3–20 km distance. For the farm in Boven-Leeuwen, the closest related wildbird virus was found at a distance of 48 km and contained 16 nt differences. For the farm in Abbega, the most similar wild-bird virus contained 56 nt changes and was found at a distance of 52 km. We did not identify a more closely related virus, suggesting that the variation of wild bird viruses introduced into the Netherlands was larger than revealed in this study.

The comparison of wild bird and poultry sequences identified no specific mutations related to adaptation of the virus to poultry. Finally, the consensus sequence of the 2016 H5N8 virus was compared with the consensus sequence of the virus detected in the Netherlands in 2014, when 5 commercial poultry holdings were infected (results not shown). This analysis identified 97 nt differences, resulting in 18 aa changes in the HA gene segment between the consensus sequences of these viruses. In addition, we found 56 nt (14 aa) differences in the NA segment. In the 6 internal segments, we found a total of 142 aa differences in the open reading frames. The numerous genetic differences between the viruses indicate that the H5N8 virus found in the Netherlands in 2016 results from a new introduction and therefore is not from continuous circulation of the 2014 virus.

Phylogenetic Analysis of the H5N8 Viruses

Figure 3

Figure 3. Schematic representation of the reassortant highly pathogenic avian influenza A virus subtypes H5N8 and H5N5 detected in the Netherlands, 2016. The Russia-Mongolia ancestor is shown in green. NL clusters 1 (blue)...

To study the origin of the H5N8 viruses found in the Netherlands in 2016, we performed a detailed phylogenetic analysis for all gene segments individually (Technical Appendix 2 Figure 2). This analysis confirms the differences between the 2014 and 2016 H5N8 viruses, which are found separate clusters. The closest relatives of the 2016 viruses in the Netherlands are the Russia–Mongolia H5N8 clade 2.3.4.4 viruses that were first detected in May 2016 in dead wild birds near Uvs-Nuur Lake. However, we identified several reassortments, which we have represented schematically (Figure 3). Five of the gene segments (PB 1, HA, NA, MP, and NS) are highly similar to the Russia–Mongolia viruses. However, the closest relatives of PB2 in the GISAID database are LPAI viruses found in Bangladesh and Russia (Technical Appendix 1 Table 1). The closest relatives of the NP gene segment are found in LPAI viruses in the Netherlands and the Republic of Georgia. An LPAI H7N9 virus detected in a commercial chicken holding in the Netherlands in June 2016 also carried a closely related NP segment.

Phylogenetic analysis of the PA gene segment showed that the 2016 H5N8 viruses in the Netherlands divert into 2 separate clusters, designated PA I and PA II. Different subtypes of LPAI viruses found in Russia and Asia were identified as closest relatives for the 2 PA clusters (Technical Appendix 1 Table 3). Most of the viruses found in wild birds and all viruses detected in commercial poultry and captive birds are in cluster PA I. Cluster PA II contains 12 viruses isolated from various wild bird species and is represented by cluster C in the network analysis (Technical Appendix 2 Figure 1). These findings show that no correlation exists between PA cluster and bird species or geographic location (Technical Appendix 2 Figure 1 panels A and B; Figure 1, panel D). The viruses in the Netherlands evolved from the Russia–Mongolia H5N8 viruses by obtaining novel PB2 and NP segments and 2 variants of the PA segment. These results demonstrate that the H5N8 HPAI epidemic in the Netherlands is caused by 2 distinct reassorted viruses.

Phylogenetic Analysis of Two H5N5 Viruses

Two HPAI H5N5 viruses were identified in dead wild birds in the Netherlands. For the tufted duck isolate (t_dk/NL-Werkendam/16014159–001/2016) PB1, PB2, HA, NP, MP, and NS are closely related to the H5N8 viruses detected in the Netherlands in 2016. The closest relatives of NA in the GISAID database are LPAI N5 viruses found in the Netherlands, Europe, and Asia in previous years (Technical Appendix 1 Table 3; Technical Appendix 2 Figure 1). For the PA segment, this H5N5 virus does not fall in either cluster PA I or PA II. The closest relatives of PA are LPAI viruses found in the Netherlands and Asia during 2012–2014. For the mute swan isolate (m-swan/NL-Groningen/16015825–001/2016), all genome segments are highly similar to an H5N5 virus identified in feces of a wild bird in Kamchatka, Russia, in October 2016. Phylogenetic analysis of the HA, PA, NP, MP, and NS segments shows that this H5N5 virus clusters with the Russia–Mongolia H5N8 viruses and not with the H5N8 2016 viruses in the Netherlands. The closest relatives of NA in the GISAID database are LPAI N5 viruses found in Asia in previous years. A similar N5 segment was also found in H5N5 viruses detected in Croatia, Italy, and Poland during the H5N8 epidemic in 2016–2017. The PB2 and PB1 segments cluster with various LPAI viruses detected earlier in Russia, Mongolia, and Europe. The mute swan virus thus is a direct descendant of a reassorted H5N5 virus from Russia, given that it was also previously detected in Kamchatka. We have depicted the genetic constitution of the 2 H5N5 viruses (Figure 3). These results indicate that 2 genetically distinct H5N5 viruses have been found in the Netherlands.

Molecular Dating of H5N8 Reassortment Events

Figure 4

Figure 4. Schematic representation of the molecular dating analysis of highly pathogenic avian influenza A virus subtypes H5N8 and H5N5 detected in the Netherlands, 2016. Time-scaled phylogenetic trees are shown for 4 gene...

To estimate the timing of the reassortment events that led to the emergence of the H5N8 viruses in the Netherlands, we performed molecular dating (Technical Appendix 2 Figure 3). We estimated the time of most recent common ancestor for all gene segments by using a relaxed molecular clock (Table). A schematic representation of the time-scaled phylogenetic tree for the HA, PB2, NP, and PA segments is shown (Figure 4). The H5N8 viruses in 2016 diverted from the 2014 viruses in April 2012 (Figure 4, panel A [node 1]), from the Russia–Mongolia H5N8 viruses in March 2016 (node 2), and we estimated the common ancestor of the viruses in the Netherlands to have occurred in August 2016 (node 3). Similar results were obtained for the PB 1, NA, MP, and NS gene segments (Table; Technical Appendix 2 Figure 3). The PB2 segment was transferred from LPAI viruses in June 2016 (Figure 4, panel B [node 3]) and the NP segment in August 2016 (Figure 4, panel C [node 3]). Two different PA segments were detected in H5N8 viruses in the Netherlands. The viruses in cluster PA I obtained the PA segment from LPAI viruses in April 2016 (Figure 4, panel D [node 4]) and in cluster PA II in August 2016 (Figure 4, panel D [node 5]). We also performed calculations by using a strict molecular clock, with similar outcomes (results not shown). Thus, the H5N8 viruses in the Netherlands were generated by multiple reassortment events with LPAI viruses, a process that was likely completed in August 2016. Because the H5N8 virus was introduced in the Netherlands in November 2016, this analysis suggests that the reassortment events took place in wild birds in Russia-Mongolia. Molecular dating for H5N5 shows that the common ancestor of the mute swan isolate and the Kamchatka isolate already circulated in January 2015. The N5 gene segments of the tufted duck and mute swan isolate diverted in 1989 (node A). Multiple reassorted viruses thus were introduced in the Netherlands by the autumn migration of wild birds from their breeding grounds in Russia-Mongolia.