Phylogeographic distribution of very virulent infectious bursal disease virus isolates in the Iberian Peninsula

Abstract

Viral population dynamics of very virulent infectious bursal disease virus (vvIBDV) field strains isolated in the Iberian Peninsula since the first outbreak in the 1990s have been analysed. Low levels of genetic variability and a global purification selection pattern were reported in 480 base pairs of the hypervariable region of the VP2 gene, indicating a lack of a selection-driven immune escape in the evolutive pathway of the virus. The viral population structure of vvIBDV strains in the Iberian Peninsula showed a strong relationship between geography and phylogeny, with two main groups observed. A global comparison among vvIBDV strains also showed an association with sequences from the same country. The low variability, the strong purifying selection and the geographical pattern observed point to a picture where the virus evolves slowly, occupying the same geographical niche for a long time. The scenario depicted fits well with the biological features of the virus: being able to remain viable for long periods of time due to a strong environmental resistance, and as an immunosuppressive agent, capable per se of annihilating temporally the immune system of the host.

Introduction

Infectious bursal disease (IBD) or “Gumboro disease” was reported for the first time nearly 50 years ago. However, nowadays it is considered to be endemic in most of the poultry-producing countries, and it is one of the major causes for economic losses in the poultry industry (Van den Berg et al., 2000). Infectious bursal disease virus (IBDV; genus Avibirnavirus, Fam. Birnaviridae) is the aetiological agent of this acute and highly contagious viral infection that only affects young chickens (Müller et al., 1979). Two serotypes (1 and 2) of IBDV have been recognized, but only serotype 1 has been demonstrated to be pathogenic (Lukert & Saif, 1997). IBDV particles have a non-enveloped, icosahedral capsid and present a bi-segmented (A and B) dsRNA genome. The larger segment A (around 3.4 kb) generates VP2, VP3, VP4 and VP5 proteins, and the smaller genome segment B (around 2.8 kb) encodes the viral RNA polymerase VP1 (Van den Berg, 2000). The external surface of the capsid is composed by trimeric sub-units formed by the VP2 protein that has been identified as the host-protective antigen (Brown et al., 1994; Cui et al., 2003; Müller et al., 2003; Van den Berg, 2000). Monoclonal neutralizing antibodies have been shown to bind to VP2, within a central hypervariable region comprising a highly hydrophobic region with two small hydrophilic peaks at each terminus, between amino acid residues 206 and 350 (Bayliss et al., 1990). Because it contains antigenic determinants and virulence markers, the VP2 protein has been commonly used in phylogenetic studies (Van den Berg et al., 2004).

In 1987, very virulent strains of infectious bursal disease virus (vvIBDV) were isolated in Europe, and quickly spread to Africa, Asia and South America (Lukert & Saif, 1997). Numerous reports of these acute IBD cases indicated that vvIBDV exhibited a high genetic and antigenic homogeneity (Van den Berg et al., 2004, and references cited therein). In the Iberian Peninsula, the first acute IBD outbreaks were reported in the early 1990s, and the pathologic and genetic studies of the isolated viruses demonstrated vvIBDV to be involved (Pagès-Mante et al., 1991; Majó et al., 2002). In 2002, a second outbreak of acute IBD started in the north-east region of Spain, and rapidly spread through the whole peninsula (Dolz et al., 2005). Genetically, the vvIBDV isolates of the 2002 outbreak were closely related with the ones isolated in the early 1990s. Since then, despite the intensive vaccination programmes and the strict biosecurity measures implemented, vvIBDV strains have become endemic in the Iberian Peninsula and acute IBD cases frequently occur.

Since the usage and availability of nucleic acid sequences have been extended in the veterinary field, numerous studies have characterized IBDV isolates all over the world (i.e., Brown et al., 1994; Etterradossi et al., 2000; Zierenberg et al., 2000; Liu et al., 2002; Rudd et al., 2002; Hernandez et al., 2006), generally following outbreak events. However, few studies have attempted to investigate population structure and viral dynamics both at a longitudinal and geographical scale seeking for clarification of IBDV evolution. Since 2002, an ongoing epidemiological surveillance has been carried out at the Centre de Recerca en Sanitat Animal (Bellaterra, Spain) by means of clinical samples submitted from suspicious IBD cases. The present work analyses the viral population dynamics of vvIBDV field strains isolated in the Iberian Peninsula since the first outbreak in the 1990s. A global comparison among vvIBDV strains reported worldwide is also assessed.

Materials and Methods

Isolates

vvIBDV strains detected from clinical samples (bursa of Fabricius) submitted to the Centre de Recerca en Sanitat Animal between 2002 and 2009 were included (n=77). In addition, sequences from vvIBDV strains isolated in the 1990s (n=3) and 2002 (n=22) outbreaks, published in previous studies (Majó et al., 2002; Dolz et al., 2005), were included in the comparison. Names, origins and accession numbers of the sequences included in the study are reported in Table 1 (n=102). In addition, 65 sequences available at GenBank (Supplementary Table 1, online only) from different countries around the world were used to study global geographic patterns in vvIBDV.

Table 1.  Origin, year of isolation and accession number for the 102 partial VP2 Iberian sequences used in the study.

Origin/strain	Year of isolation	Accession number	Origin/strain	Year of isolation	Accession number
Catalunya			Andalucía
CATVG248	1990	AY083925	AND13	2002	AY770583
CAT03	2002	AY770581	AND31	2002	AY770593
CAT33	2002	AY770591	AND37	2002	AY770583
CAT35	2002	AY770581	AND38	2002	AY770581
CAT40	2002	AY770583	AND39	2002	AY770583
CAT41	2002	AY770583	AND84	2008	JF682281
CAT43	2002	AY770581	AND121	2009	JF682293
CAT47	2002	AY770581	AND144	2009	JF682298
CAT48	2002	AY770581	AND147	2009	JF682299
CAT49	2002	AY770590	AND1472	2009	JF682315
CAT60	2004	JF682245	Navarra
CAT728	2004	JF682247	NAV781	2008	JF682278
CAT34	2006	JF682260	NAV783	2008	JF682279
CAT2650	2006	JF682257	NAV118	2009	JF682291
CAT54	2008	JF682274	NAV177	2009	JF682304
CAT73	2008	JF682276	La Rioja
CAT74	2008	JF682277	RIO119	2009	JF682292
CAT106	2008	JF682288	C. La Mancha
CAT113	2009	JF682290	CM6145	1990	AY083930
CAT1004	2008	JF682286	CM244	2009	JF682318
CAT1005	2008	JF682287	C. León
CAT1010	2008	JF682283	CLe1296	2004	JF682244
CAT1013	2008	JF682285	CLe20	2005	JF682253
CAT1014	2008	JF682284	CLe1836	2005	JF682249
CAT124	2009	JF682294	CLe271	2006	JF682258
CAT168	2009	JF682300	CLe272	2006	JF682259
CAT186	2009	JF682308	CLe584	2007	JF682261
CAT197	2009	JF682310	CLe61	2008	JF682271
CAT204	2009	JF682312	CLe06	2008	JF682270
CAT213	2009	JF682313	Galicia
CAT234	2009	JF682316	GAL50	2002	AY770592
CAT1741	2009	JF682301	GAL222	2005	JF682254
CAT2142	2009	JF682314	GAL224	2005	JF682255
Illes Balears			GAL284	2005	JF682248
IB5939	1990	AY083926	GAL59	2007	JF682262
IB21	2005	JF682251	GAL61	2007	JF682265
IB24	2005	JF682256	GAL27	2008	JF682273
València			GAL83	2008	JF682280
VAL46	2002	AY770582	GAL136	2009	JF682295
VAL14	2005	JF682250	GAL141	2009	JF682297
VAL96	2008	JF682282	GAL2511	2009	JF682319
VAL109	2009	JF682289	GAL2512	2009	JF682320
VAL2141	2009	JF682314	Portugal
Murcia			POR241	2009	JF682317
MU1124	2004	JF682246	Undetermined
MU11	2008	JF682272	SP01	2002	AY770581
Aragón			SP02	2002	AY770581
AR42	2002	AY770581	SP27	2002	AY770581
AR140	2009	JF682296	SP32	2002	AY770581
AR192	2009	JF682309	SP44	2002	AY770582
AR202	2009	JF682311	SP15	2005	JF682252
AR1742	2009	JF682302	SP80	2007	JF682263
AR1744	2009	JF682303	SP81	2007	JF682264
AR1781	2009	JF682305	SP82	2007	JF682266
AR1782	2009	JF682306	SP83	2007	JF682267
AR1783	2009	JF682307	SP84	2007	JF682268
			SP85	2007	JF682269

Polymerase chain reaction and sequencing

For all suspicious cases, viral RNA was extracted with the NucleoSpin RNA Virus kit (Macherey-Nagel, Düren, Germany) and amplified using a reverse transcriptase-polymerase chain reaction assay for the VP2 gene, as previously described by Dolz et al. (2005). Sequences of 480 base pairs from each of the hypervariable region in the VP2 gene were obtained. The analysed region included 159 amino acid residues, from positions 202 to 360, containing the hydrophilic peaks responsible for the IBDV antigenicity.

Phylogenetic and population studies

The phylogenetic relationships among the 102 Iberian IBDV sequences were estimated with a maximum parsimony network of the nucleotide sequences by the statistical parsimony method of Templeton et al. (1992) using the TCS computer program (Clement et al., 2000). Patterns of nucleotide diversity distribution among regions were estimated by a hierarchical nested analysis of molecular variance (AMOVA; Excoffier et al., 1992) of the frequency distribution of sequences and their pairwise divergence at three hierarchical levels: within populations (ϕst), among populations within groups (ϕsc) and among groups (ϕct). Two analyses were performed using Arlequin software (version 3.0: Excoffier et al., 2005): one considering a temporal design and a second considering a geographical design. The temporal design separated sequences among outbreaks: during the 1990s, in 2002, and after 2002. In the geographic design, two large areas that contained different regions of the Iberian Peninsula were considered: the Mediterranean region (Catalunya, València, lles Balears, Murcia, Andalucía) together with the north-eastern region (Aragón, Navarra, Rioja; M/E region); and the Atlantic and central-western regions (Portugal, Galicia, Castilla-la-Mancha, Castilla-León; A/W region). A matrix of nucleotide divergence between regions was calculated using MEGA 4 (Tamura et al., 2007) and the gamma Tamura–Nei model. This distance matrix of nucleotide divergence was used to perform principal coordinate (PCO) analysis using the PCO program (Anderson, 2003). The PCO analysis determines a small set of synthetic variables (the coordinates) that explain the variability of the original dataset. The first principal coordinate accounts for as much variability in the data as possible, and each following coordinate accounts for as much of the remaining variability as possible. Finally, a neighbour-joining tree including Iberian and world vvIBDV nucleotide sequences was constructed using MEGA 4 (Tamura et al., 2007) and the gamma Tamura–Nei model. The confidence of internal branches was estimated by means of 1000 bootstrap replicates.

Selection pressure

The existence of selective pressures along partial VP2 sequences was assessed by calculating the difference between non-synonymous (dN) and synonymous (dS) rates (dN – dS) both globally and per codon using the modified Nei–Gojobori method by MEGA 4 software. These differences were calculated with the SNAP web utility (http://hivweb.lanl.gov/content/hivdb/SNAP/WEBSNAP/SNAP.html).

Results

Nucleotide and amino acid variability

Among the 102 sequences of 480 base pairs each, 68 variable nucleotide positions were reported, most of which (n=55) were synonymous and only 13 were non-synonymous. The majority of the sequence changes were point mutations; more than one-half of the nucleotide (53%, 36 out of 68) and amino acid (61%, eight out of 13) changes were detected in only one or two strains. Only five nucleotide mutations were present in a large number of sequences (35C → T, 71G → A, 116C → T, 295G → A, 362G → A); four were silent, and only the nucleotide change in position 295 (VP2 position 1027) caused an amino acid substitution in residue 98 (VP2 position 299) (S → N). The percentage of nucleotide identity ranged from 94.7 to 100%.

Most of the amino acid codons showed a neutral selection (Figure 1). Codons identified under negative selection (55 residues, 34.6%) were mainly grouped in two regions of the VP2 protein (residues 212 to 272, and residues 311 to 351). Few positions (5.7%) showed a positive selection pressure; nearly all presented a value close to a zero, except position 299 where a frequent amino acid substitution was reported.

Figure 1.  Accumulated (line, values referred to the left axis) and individual (points, values referred to the right axis) differences between non-synonymous and synonymous rates (dN – dS) for the VP2 gene. Numbers on the horizontal axis represent the amino acid position.

Temporal and geographical analyses of Iberian vvIBDV sequences

The statistical parsimony network of the IBDV (Figure 2) reflected the low number of mutations previously described in the nucleotide and amino acid variability: a maximum of 14 steps separated the most distant Iberian strains. Isolates from the first outbreak in the 1990s were placed closely in the network, whereas sequences from the second outbreak in 2002, and more recent ones, were distributed all over the network, indicating the lack of a temporal structure. In contrast, nearly all of the strains from the Atlantic and north-western Iberian Peninsula (Galicia, Castilla-León, Portugal) clustered together and separated from the Mediterranean and north-eastern ones (Figure 3), which correlates with the nucleotide positions described above. This point was also supported by the AMOVA results (Table 2), where the temporal structure of the Iberian strains was not supported (Φ_ct=0.0775, P=0.4909), whereas a significant percentage of variation (24.45%) was allocated between the M/E and A/W groups (Φ_ct=0.2445, P=0.0091). In addition, the neighbour-joining tree among regions clearly separated the A/W regions from the M/E ones (data not shown). Finally, the first axis of the PCO analysis placed the A/W regions in the positive scores and the M/E regions in the negative scores (Figure 4).

Figure 2.  Statistical parsimony network based on the nucleotide sequences. Every branch represents a nucleotide mutation: black circles, mutations with no frequency; ovals, strains isolated during the 1990s; dashed squares, strains isolated in 2002; solid squares, strains isolated after 2002.

Figure 3.  Iberian Peninsula map, showing the relative frequencies of viral strain group Atlantic/West (A/W, black) and Mediterranean/East (M/E, grey) in every geographical area.

Figure 4.  PCO analysis. Numbers between brackets indicate the percentage of variation explained by every coordinate.

Table 2.  AMOVA results.

Source of variation	Degrees of freedom	Sum of squares	Variance components	% of variation	Related Φ statistics
Temporal design
Within populations	95	320.10	3.679	81.64	Φst=0.1066 (P=0.0000)
Among populations within groups	5	58.34	0.758	14.88	Φsc=0.1709 (P=0.0000)
Among groups	2	13.64	0.319	3.48	Φct=0.0775 (P=0.4909)
Total	102	392.07	4.118	100
Geographic design
Within populations	80	248.58	3.069	65.09	Φst=0.3491 (P=0.0000)
Among populations within groups	10	61.51	0.493	10.46	Φsc=0.1384 (P=0.0000)
Among groups	1	50.38	1.153	24.45	Φct=0.2445 (P=0.0091)
Total	91	360.47	4.715	100

Global geographic distribution patterns of vvIBDV sequences

When Iberian sequences were compared with vvIBDV from other parts of the world, a geographic clustering tendency was observed (Figure 5) despite the Iberian sequences being split into distinct branches. Nearly all sequences from the A/W Iberian group clustered together in a single branch, except three sequences that, similarly to the statistical parsimony network, were placed in the M/E group. Sequences from the M/E group clustered into several branches along the tree, but one-half (36 out of 77) clustered into a single homogeneous branch. Several terminal branches showed a strong association with sequences from the same country. In several countries all of the sequences clustered into a single clade (Indonesia, Uruguay, Venezuela, Malaysia, Nigeria, Brazil), and in other cases were distributed in two clades (Iran, China).

Figure 5.  Circular neighbour-joining tree based on 102 Iberian and 65 global vvIBDV nucleotide sequences. Colours indicate different geographic origins.

Discussion

In contrast with other dsRNA viruses, high nucleotide identity has characterized vvIBDV strains worldwide (>94%; Etterradossi et al., 2000). Indeed, the Iberian vvIBDV field strains showed low levels of genetic variability, alongside a global purification selection pattern in the studied fragment of the VP2 gene. Negative or purifying selection is reported when synonymous mutations (a nucleotide mutation that does not produce an amino acid change) outnumber non-synonymous mutations (a nucleotide mutation that produces an amino acid change) for a certain codon (Nei & Kumar, 2000). The negatively selected codons were mainly located in two regions roughly coincident with the major hydrophilic peaks (210 to 225 and 312 to 324) reported for the VP2 gene (Bayliss et al., 1990). Despite the regions under purifying selection being wider than the hydrophilic peaks, this is most probably a consequence of a hitchhiking effect commonly observed in regions under strong functional constraint (Barton, 2000). Considering that any change within these hydrophilic regions may potentially cause a change in the antigenic features of an isolate (Letzel et al., 2007; Jackwood & Sommer-Wagner, 2011), a strong purifying selection pattern as observed in this study could indicate a lack of a selection-driven immune escape in the evolutive pathway of the virus.

Partial VP2 sequences have been widely used in routine diagnosis and epidemiological studies worldwide (Jackwood & Sommer-Wagner, 2011). Nonetheless, in order to accurately infer an IBDV phenotype, genome fragments including segments A and B need to be analysed. Alternatively, in vivo assays should be carried out to confirm the very virulent phenotype. Keeping this limitation in mind, field isolates in the present study were classified as vvIBDV when clinical signs of acute IBD were observed in the birds and partial VP2 sequences clustered in the vvIBDV branch (data not shown). Many molecular studies in Gumboro have been restricted to the regional level: Spain (Pagès-Mante et al., 1991; Majó et al., 2002; Dolz et al., 2005), Uruguay (Hernandez et al., 2006), Croatia (Lojkic et al., 2008), People's Republic of China (Liu et al., 2002) or Brazil (Fernandes et al., 2009). In contrast, few studies seek for evolutive relationships on a global scale among vvIBDV strains (Eterradossi et al., 2000; Zierenberg et al., 2000). In order to determine the population structure present among vvIBDV strains, this study adopted two approaches: a general comparison among most (n=65) of the vvIBDV strains described worldwide, and a regional longitudinal analysis including 102 vvIBDV Iberian strains isolated during a 20-year period. According to all of the performed analyses—statistical parsimony network, AMOVA and PCO analysis—the viral population structure of vvIBDV strains in the Iberian Peninsula nowadays follows a geographic pattern, with two main groups observed: A/W and M/E. Nearly all of the strains isolated after 2002 clustered according to their geographical origin. When compared on a global scale, the vvIBDV strains clustered with a similar geographical-based pattern, despite some minor incongruence being shown, indicating a strong association between geography and viral diversity.

Very virulent IBDV exhibit a so-called vvIBDV genetic fingerprint in the VP2 hypervariable region at amino acid residues 222 (Ala), 242 (Ile), 256 (Ile), 294 (Ile) and 299 (Ser) (Van den Berg et al., 2004). However, the presence of a serine in the residue 299 was not a constant trait, although it is believed to be a common feature among European vvIBDV field strains (Domanska et al., 2004). Comparing the European and African vvIBDV strains, Etterradossi et al. (2000) found an asparagine (N) for an African isolate at position 299 while the remaining French, British and Dutch isolates presented a serine (S) in this position. Interestingly, in our study, isolates from the A/W group presented the serine residue, while most isolates from the M/E group showed asparagine. Nevertheless, two of the M/E isolates that do not present the asparagine residue were isolated during the first vvIBDV reports in the 1990s. Whether a new strain arrived from Africa and spread only in the Mediterranean area after the 1990s, or a de novo mutation restricted to the Mediterranean region occurred between the first (1990) and the second (2002) outbreaks, remains an open question.

Taken together, the low variability, the strong purifying selection and the geographical pattern observed point to a picture where the virus evolves slowly, occupying the same geographical niche for a long time. The depicted scenario fits in well with the biological features of the virus: being able to remain viable for long periods of time due to a strong environmental resistance (Van den Berg, 2000); and as an immunosuppressive agent, being capable per se of annihilating temporally the immune system of the host. Due to the strong nature of IBDV and the widespread distribution of the virus, prevention and control of this disease may be achieved by a combination of biosecurity measures and the correct vaccination programme.

Acknowledgements

K. Bertran held an FPU Pre-Doctoral grant from the Spanish Ministry of Education. M. Cortey received a grant from the Government of Catalonia.