Phylogenetic analysis of highly pathogenic avian influenza A ( H 5 N 8 ) virus outbreak strains provides evidence for four separate introductions and one between-poultry farm transmission in the Netherlands , November 2014

Citation style for this article: Bouwstra RJ, Koch G, Heutink R, Harders F, van der Spek A, Elbers AR, Bossers A. Phylogenetic analysis of highly pathogenic avian influenza A(H5N8) virus outbreak strains provides evidence for four separate introductions and one between-poultry farm transmission in the Netherlands, November 2014. Euro Surveill. 2015;20(26):pii=21174. Available online: http://www.eurosurveillance.org/ViewArticle.aspx?ArticleId=21174


Introduction
pathogenic avian Influenza (HPAI) is a viral disease of poultry causing high mortality [1].In 2003, HPAI A(H7N7) virus was detected on 241 Dutch farms [2].After the epidemic was stopped, syndromic surveillance and monitoring programmes were initiated to enable early detection of the introduction of any H5 and H7 avian influenza viruses on poultry farms.Since 2003, only low pathogenic avian Influenza viruses of diverse subtypes were detected.Up to 2014, HPAI H5 viruses had not been detected in the Netherlands.
In China, the H5N8 HP virus had been isolated in 2009-10 as part of a monitoring programme [3].From January 2014, the HPAI H5N8 virus spread very rapidly in South Korea.Genetic sequence analysis indicated that virus isolates from infected farm ducks in South Korea and dead Baikal teals in the surrounding area strongly resembled the earlier Chinese isolates [4].

Outbreak 1
Dead chickens from an indoor farm with 150,000 laying hens (Farm 1), in the municipality of Hekendorp in the province of Utrecht, were submitted for necropsy to the Dutch Animal Health Service in mid-November 2014, with an anamnesis of exponentially increasing mortality in the past days.Swab samples from these chickens were forwarded to the Central Veterinary Institute part of Wageningen UR, Lelystad, which is standard procedure when there is exponentially increasing mortality and inconclusive pathology.The farm was situated next to a river and in the middle of peat land with an abundant presence of wild waterfowl.Exponentially increasing mortality had been observed by the poultry farmer in one of six poultry houses in the days preceding the notification.

Outbreak 2
On 19 November, HPAI suspicion was raised on an indoor farm with 43,000 laying hens (Farm 2), in the municipality of Ter Aar in the province of Zuid-Holland.Three days previously, the poultry farmer saw exponentially increasing mortality, with birds showing conjunctivitis and ruffled feathers in one of three poultry houses.Egg production started to drop on 18 November and two days later, egg production decreased by 20%.Veterinarians of a specialist team from the Dutch Food and Safety Authority (NVWA) visited the farm on 19 November, took samples from clinically sick birds and the farm was quarantined.Like Farm 1, this poultry farm was also situated in an area with an abundant presence of wild waterfowl.Farm 2 was situated 21 km north-north-west of Farm 1.

Outbreak 3
One day later (20 November), HPAI suspicion was raised on a pullet-rearing farm with 11,100 birds (Farm 3), in the municipality of Kamperveen in the province of Gelderland, housing hens and cocks in two separate houses.This farm was located in a stretch of farmland boarded by a lake on one side and the IJssel river on the other.In the area, large quantities of wild waterfowl were present and it is known to be a congregating place for wild birds before they take off for their journey across the lake to the polder area of the province of Flevoland.Farm 3 was situated 90 km north-east of Farm 2 and 92 km north-east of Farm 1. Increased mortality was first observed on 18 November on Farm 3, as well as a temporary decrease in egg production.On 20 November, during a visit by the competent veterinary authorities after notification of the suspicion, increased mortality was seen in both poultry houses and the birds showed typical signs of HPAI [1].

Outbreak 4
In the evening of 19 November, a poultry farmer submitted several dead meat ducks to his veterinarian.The poultry farm, which housed 15,000 meat ducks each in two poultry houses indoors, was situated 550 m from Farm 3. The carcasses of the meat ducks were submitted the following day for necropsy to the Dutch Animal Health Service, who notified the Dutch Food and Safety Authority.Swab samples from the ducks at necropsy were forwarded, in the framework of the Dutch early-warning system for avian influenza, to the Central Veterinary Institute part of Wageningen-UR.

Outbreak 5
The fifth outbreak was reported on 29 November following exponentially increasing mortality in the four preceding days in an indoor farm housing 29,000 layer hens.This farm was located in Zoeterwoude about 30 km west of Farm 1 and 20 km south-west of Farm 2. The five outbreak locations are shown in Figure 1.
At the Central Veterinary Institute part of Wageningen UR, swab samples of sick chickens and ducks from the five outbreak farms were tested following our accredited procedures.All five farms were confirmed to be positive for HPAI H5N8.Backward and forward tracing of possible dangerous contacts in the framework of the standard epidemiological investigation by the Dutch Food and Consumer Product Safety Authority (i.e.transport, professional visitors such as advisors, veterinary practitioner, etc.; possibly contaminated materials delivered to the farm such as feed, bedding etc.; possibly contaminated transport vehicles, etc.) by the Dutch Food and Safety Authority revealed no indication for dangerous contacts of such connections between HPAI H5N8 outbreaks in Asia, Germany, the United Kingdom and the Netherlands occurring between October and December 2014 [5].Moreover, no links were discovered between the outbreaks in the Netherlands.In addition, recent analysis from our group showed that the Dutch virus (A/Ch/Netherlands/14015526) has high similarity to two South Korean and one Japanese strains [6].
The question arose as to whether the outbreaks on the five Dutch poultry farms were caused by separate virus introduction or by transmission between farms.In an attempt to answer this question, we sequenced the complete virus RNA genome obtained from several animals from each farm.The aim of our study was to assess possible routes of transmission of the virus by sequence and temporal phylogenetic analysis.

Genome sequence analysis
RNA was extracted from cloacal and oropharyngeal pools of five samples originating from clinically affected hens positive in a screening matrix-gene realtime PCR [7,8], which detects all avian influenza virus subtypes.The positive samples were then checked for the presence of notifiable subtypes (H5 and H7) by realtime PCR as recommended by the European Union reference laboratory in Weybridge, United Kingdom [9,10].Haemagglutinin (HA) and neuraminidase (NA) sequence analysis was based on PCR fragments that were generated according to the so-called KHA PCR [11] and PanHA [12] and PanNA [13] protocols.The sequence of the HA gene revealed polybasic amino acids -RRRKR*GLF -at the HA cleavage site, a motif typical for HPAI viruses.In addition, HA and NA sequence results showed that the virus was of the H5N8 subtype.
The sequencing results of the cleavage site also revealed high similarity to the German outbreak strain (November, 2014) [14].However, as complete genome sequencing is necessary to have more insight into the origin and emergence of this virus into Europe and specifically into the Netherlands, we amplified all eight RNA genome segments of the outbreak viruses using universal eight-segment primers and directly sequenced [15].Purified amplicons were sequenced at high coverage (average > 1,000 per nucleotide position) using the Nextera library preparation method and subsequently sequenced using Illumina MiSeq paired-end 150 base pairs sequencing (Illumina, San Diego, CA, United States).Quality control-passed sequence reads of high quality were iteratively mapped on resulting consensus sequences using Bowtie2 [16] starting against the South Korean H5N8 genome sequence (GenBank accession numbers KJ511809-KJ511816) to generate a majority (> 80% evidence) consensus sequence of all segments.The consensus sequences were compared with de novo assembled sequences reads using SPAdes-v3 [17] and no significant differences were detected.The majority consensus sequences were submitted to the Global Initiative on Sharing All Influenza Data (GISAID) (EPI_ISL_167905, EPI_ISL_174349, EPI_ISL_174350, EPI_ISL_174351, EPI_ISL_174352) and subsequently for all nucleic acid sequences, a basic molecular phylogenetic analysis was performed using the maximum likelihood method based on the Tamura-Nei model using a gamma distributed nucleotide substitution rate [18].Between two and four pools of samples per farm were sequenced.

Temporal phylogenetic analysis
The five fully sequenced Dutch H5N8 sequences were each aligned with 22 H5N8 sequences obtained from GISAID using Muscle in MEGA6 [19].For each of the eight segment alignments, the simplest evolutionary model fitting the dataset was the Hasegaw-Kishino-Yano model with gamma distributed rates [20].Nucleotide substitution rates were estimated using Bayesian Markov Chain Monte Carlo methods [21].Analysis was performed using the programme BEAST v1.8.1 [22] using strict or relaxed uncorrelated molecular clocks that were calibrated using the sample isolation dates.All genome segments were treated as virus partitions with individual substitution and clock models and analysed for 30 million generations, sampling every 3,000 generations.Effective sample sizes were checked using Tracer 1.6 [23]: the values were far above the minimum threshold of 200.A maximum clade credibility tree was constructed to summarise all 10,000 trees after 10% burn-in using TreeAnnotator [21].The final time-scaled phylogenetic tree was visualised and annotated using FigTree 1.4.2[24].Three or eight gene segment alignments were manually concatenated to generate a single alignment that was used to construct phylogenetic networks using the median-joining method implemented in the programme NETWORK as described by Bataille et al. [25].This model-free method uses a parsimony approach, based on pairwise differences, to connect each sequence to its closest neighbour, and allows the 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.

Results
The number of nucleotide differences between sequences of all eight segments of viruses detected at the five Dutch outbreak farms varied from four for the viruses from the two farms in Kamperveen to 51 nucleotides resulting in 12 amino acid differences between the viruses from Ter Aar and Zoeterwoude (Figure 2).
The median-joining network shows that all five sequences were derived from one or more calculated ancestors and were not descendants of each other.To determine the relationship with other European isolates, we aligned the five Dutch sequences with 22 sequences obtained from the EpiFlu database.Isolate ID numbers and providers are listed in the Table.

Figure 2
Median-joining network analyses of five highly pathogenic avian influenza A(H5N8) sequences, the Netherlands, November 2014 The median-joining network was constructed from the combined sequence of eight gene segments data.This network included all the most parsimonious trees linking the sequences.Each unique sequence genotype is represented by a yellow circle sized relative to its frequency in the dataset.Numbers refer to the position of the mutation within the combined sequences.Red circles represent median vectors.The sequence detected in samples from the duck farm in Kamperveen (Outbreak 4) was heterogeneous at position 12,486 (A or G).Analyses of the 27 H5N8 isolates showed diversification into three lineages that diverged between 2009 and 2010 (Figure 3).

Calculation of the time of most recent common ancestor reveals that the origin of the H5N8 viruses occurred between 9 September 2007 and 20 June 2008 (Supplementary Table 1 [26]
). Calculations were also performed using a strict clock or uncorrelated relaxed clocks, with essentially similar outcomes (Supplementary Tables 2 and 3 [26]).The viruses isolated in Europe are located within lineage I, forming a separate sublineage.Lineage I diversified between the end of June and end of December 2013 into three sublineages of viruses that caused outbreaks in poultry and wild birds in South Korea and China in the beginning of 2014.Viruses isolated in Europe and Japan (Chiba) in the autumn of 2014 seem to have evolved from viruses that circulated in wild birds, including Baikal teals, in South Korea in beginning of 2014 The common ancestor of the isolates from Europe and Chiba was estimated to have emerged somewhere between 15 July and 8 August 2014 (mean: 28 July).
Four of the five Dutch strains were differently located in lineage I: A/Ch/Nl/Ter_Aar and A/Ty /Germ-MV are closely related and diverted around the end of August 2014.In contrast, A/Ch/NL-Hekendorp, A/Dk/Eng, A/ Ch/NL-Kamperveen and A/Ch/NLZoeterwoude diverged slightly earlier, around 8 August.Both viruses isolated from Kamperveen diverged only in the first week of November, consistent with the possible transmission between the two farms (Figure 3).
Trees of each eight gene segments of the 27 H5N8 viruses were constructed using the minimal likelihood method using a Tamura-Nei substitution model, gamma-distributed substitution rate and 1,000 bootstraps.The five Dutch sequences, (A/Ch/ NL-Hekendorp, A/Ch/NL-Ter Aar, A/Ch/NL-Kamperveen, A/Dk/NL-Kamperveen, A/Ch/NL-Zoeterwoude), A/Dk/ Eng, A/Germ/-MV, A/Ty/Engl and A/Dk/Chiba, were located in the same subgroup of lineage I and clustered in a similar fashion to each other (Figure 4).

Figure 4a
Phylogenetic trees of nucleotide sequences of influenza A(H5N8) viruses Evolutionary history was inferred using the maximum likelihood method based on the Tamura-Nei model [18].The trees of polymerase 1 (panel A), polymerase B2 (panel B), polymerase A (panel C), haemagglutinin (panel D), nucleoprotein (panel E), neuraminidase (panel F), matrix protein (panel G) and non-structural protein (panel H), with the highest log likelihood shown.The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) is shown next to the branches.A discrete gamma distribution was used to model evolutionary rate differences among sites.The trees are drawn to scale, with branch lengths measured in the number of substitutions per site.Evolutionary analyses were conducted in MEGA6 [19].Viruses detected at the Dutch farms are shown in bold.

Figure 4b
Phylogenetic trees of nucleotide sequences of influenza A(H5N8) viruses Evolutionary history was inferred using the maximum likelihood method based on the Tamura-Nei model [18].The trees of polymerase 1 (panel A), polymerase B2 (panel B), polymerase A (panel C), haemagglutinin (panel D), nucleoprotein (panel E), neuraminidase (panel F), matrix protein (panel G) and non-structural protein (panel H), with the highest log likelihood shown.The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) is shown next to the branches.A discrete gamma distribution was used to model evolutionary rate differences among sites.The trees are drawn to scale, with branch lengths measured in the number of substitutions per site.Evolutionary analyses were conducted in MEGA6 [19].Viruses detected at the Dutch farms are shown in bold.

Figure 4c
Phylogenetic trees of nucleotide sequences of influenza A(H5N8) viruses Evolutionary history was inferred using the maximum likelihood method based on the Tamura-Nei model [18].The trees of polymerase 1 (panel A), polymerase B2 (panel B), polymerase A (panel C), haemagglutinin (panel D), nucleoprotein (panel E), neuraminidase (panel F), matrix protein (panel G) and non-structural protein (panel H), with the highest log likelihood shown.The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) is shown next to the branches.A discrete gamma distribution was used to model evolutionary rate differences among sites.The trees are drawn to scale, with branch lengths measured in the number of substitutions per site.Evolutionary analyses were conducted in MEGA6 [19].Viruses detected at the Dutch farms are shown in bold.

Figure 4d
Phylogenetic trees of nucleotide sequences of influenza A(H5N8) viruses Evolutionary history was inferred using the maximum likelihood method based on the Tamura-Nei model [18].The trees of polymerase 1 (panel A), polymerase B2 (panel B), polymerase A (panel C), haemagglutinin (panel D), nucleoprotein (panel E), neuraminidase (panel F), matrix protein (panel G) and non-structural protein (panel H), with the highest log likelihood shown.The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) is shown next to the branches.A discrete gamma distribution was used to model evolutionary rate differences among sites.The trees are drawn to scale, with branch lengths measured in the number of substitutions per site.Evolutionary analyses were conducted in MEGA6 [19].Viruses detected at the Dutch farms are shown in bold.

Discussion
Comparative analyses of all European and Asian H5N8 viruses using Bayesian Markov Chain Monte Carlo and median-joining network analyses suggest that four of the five outbreaks of HPAI H5N8 virus in the Netherlands were caused by separate introduction and not by farm-to-farm spread.In addition, the analyses suggest between-farm transmission between the third and the fourth outbreaks, both located in Kamperveen, at a distance of 550 m from each other, although we cannot entirely exclude that both outbreaks resulted from two separate introductions from the same source.If the virus had spread between farms after a single introduction, it is expected that the Dutch viruses would have diversified from a single node in the tree [25].Moreover, the five Dutch viruses had a maximum of 20 nucleotide substitutions in the PB2, HA and NA gene segments (Supplementary Figure 1 [26]) that must have been generated during circulation in poultry during nine days, if we assume that between-farm spread caused all outbreaks.However, during the HPAI H7N7 outbreak in the Netherlands, a maximum of 25 substitutions in the same three genes were generated in nine weeks and in Italy, 66 substitutions in nine months [27].On the basis of these numbers, we would have expected between three and six mutations between the five H5N8 outbreak virus sequences.In addition, the fact that H5-specific antibodies were not detected in animals of the outbreak farms (data not shown) excluded the possibility that the virus had circulated for some time unnoticed.Moreover, tracing of contacts by the Dutch Food and Safety Authority revealed no indication of epidemiological connections between the Dutch outbreaks and the farms in outbreaks 1, 2, 3 and 5 were located far from each other (16 to 112 km).
Remarkably, H5N8 viruses isolated from non-specified species of ducks (Anatidae) in Chiba, Japan, in November 2014 [28] most likely derived from the same precursor as that of viruses isolated in Europe in the same period.Rates of nucleotide substitution and time of most recent common ancestor analyses showed that the origin of the H5N8 Japanese and European H5N8 viruses dated back to late summer of 2014 based on the combined PB1, PB2, PA, HA, NP, NA, MP and N genome components in the dataset (Figure 3).Values for the time of most recent common ancestor and highest posterior density were robust whether a strict or relaxed clock (Supplementary Tables 1, 2 and 3 [26]) or different viruses (Supplementary Figures 2 and 3 [26]) were used in the analyses.The results are consistent with the hypothesis that the precursor virus arose in Siberia on breeding grounds where migratory birds from the East Atlantic Flyway and Asia Australia Flyways [29] may have mingled during the breeding season in 2014.Recently published HA and NA sequences of virus RNA detected in a wigeon that was shot in September 2014 in north-east Russia show that the HA and NA sequences of this virus are phylogenetically located near the node of the Eurasian and Chiba viruses [6] (Figure 4 panels D and F).
It is known that viruses that pass species barriers will show adaptation mutations [30].These adaptations become visible in phylogenetic analysis.The fact that European viruses diverged from the same ancestor as that of two Japanese viruses might indicate that the virus travelled from Asia to Europe in just a single or a few bird species.If a lot of bird species were involved in the transport of H5N8 from Asia to Europe, you would expect no Asian viruses in lineage I.In addition, none of the amino acid mutations in the five Dutch viruses (Figure 2) are known to affect virulence of the virus in mammals.
Pathogenicity studies showed that the virus was highly virulent for chickens, but mildly to moderately virulent for wild ducks [31], suggesting the potential for transport of the virus over large distances.From January 2014, HPAI H5N8 virus spread very rapidly in South Korea, initially mainly among farm ducks.At the time of the first outbreaks among farm ducks, a large number of dead Baikal teals (Anas formosa) -a species of migratory ducks -were found near the affected farms, leading to the hypothesis that the infection may have been carried by these migratory ducks].Genetic sequence analysis indicated that isolates of infected farm ducks in South Korea and dead Baikal teals in the surrounding area strongly resembled Chinese isolates from 2010 to 2013 [4], while it was also noted that isolates of HPAI H5N8 virus in South Korea was a product of reassortment of A/duck/Jiangsu/k1203/2010 (H5N8) and other avian influenza virus subtypes that co-circulated in birds in east Asia from 2009 to 2012 [32].Kang et al. very recently demonstrated in a pathogenicity study of various HPAI H5N8 virus isolates that the virus was moderately virulent in experimental infection trials with wild ducks (Anas platyrhynchos) and Baikal teals and did not result in serious disease and/or mortality [31].
These findings emphasise the clear need for the utmost attention concerning hygienic measures and biosecurity by poultry farmers to prevent introduction of disease agents into poultry houses.Somehow, the virus was brought into indoor poultry farms.For example, via persons with contaminated clothing/boots/materials/feed, or by vermin and flies.In addition, it would be wise at this moment in the epidemic to house poultry that normally use outdoor facilities in order to prevent exposure to possibly infected wild waterfowl and their excrement [33].
In conclusion, phylogenetic analysis of the Dutch HPAI H5N8 outbreak strains helped to unravel possible routes of transmission.These kinds of analyses, combined with other epidemiological and laboratory data, might provide tools to support specific preventive measures.

Figure 1
Figure 1Location of five highly pathogenic avian influenza A(H5N8) outbreaks in the Netherlands, November 2014

Figure 3
Figure 3Phylogenetic trees derived from complete genome sequences of highly pathogenetic avian influenza A(H5N8) viruses