Molecular characteristics and evolutionary analysis of a very virulent infectious bursal disease virus

Infectious bursal disease virus (IBDV) poses a significant threat to the poultry industry. Viral protein 2 (VP2), the major structural protein of IBDV, has been subjected to frequent mutations that have imparted tremendous genetic diversity to the virus. To determine how amino acid mutations may affect the virulence of IBDV, we built a structural model of VP2 of a very virulent strain of IBDV identified in China, vvIBDV Gx, and performed a molecular dynamics simulation of the interaction between virulence sites. The study showed that the amino acid substitutions that distinguish vvIBDV from attenuated IBDV (H253Q and T284A) favor a hydrophobic and flexible conformation of β-barrel loops in VP2, which could promote interactions between the virus and potential IBDV-specific receptors. Population sequence analysis revealed that the IBDV strains prevalent in East Asia show a significant signal of positive selection at virulence sites 253 and 284. In addition, a signal of co-evolution between sites 253 and 284 was identified. These results suggest that changes in the virulence of IBDV may result from both the interaction and the co-evolution of multiple amino acid substitutions at virulence sites.

Infectious bursal disease virus (IBDV) is an avian virus that specifically infects the bursa of Fabricius in young chickens. It induces the formation of acute lesions of lymphatic tissues and the destruction of lymphocytes in the bursa. Because of its immunosuppressive effects, IBDV severely impairs the immune response to many vaccines and causes a general failure of immunity [1][2][3]. The virus first emerged in the 1950s [1] and has caused egregious economic ruin to the modern poultry industry [4]. Two serotypes of IBDV, serotype 1 and serotype 2, have been characterized, but only serotype 1 is pathogenic to young chickens [5]. IBDV serotype 1 has undergone tremendous genetic diversification since it was first identified such as classical virulent strains (cvIBDV) [6], antigenic variant strains [7], and very virulent strains (vvIBDV) [8]. The vvIBDV strains, which have emerged and spread in the last 20 years, result in mortality rates as high as 60% in young chickens, thus posing severe challenges to avian epidemic control and prevention in China and throughout the world.
IBDV belongs to the genus Avibirnavirus, in the family Birnaviridae, and contains two double-stranded RNA segments, segment A (3.2 kb) and segment B (2.8 kb) [6]. Segment A encodes four viral proteins: VP2, VP3, VP4, and VP5 [9][10][11]. VP2, the major structural protein [11], assembles into trimers, and then 260 of these trimers aggregate to form the IBDV capsid. The folding of the amino acid chain of VP2 forms three domains, the protrusion domain (P domain), the shell domain (S domain), and the helical base domain (B domain). The P domain is located on the viral surface, perpendicular to the capsid [12,13]. According to the crystal structure of VP2 of an attenuated IBDV strain [12], neutralization escape epitopes and cell culture adaptation/virulence sites are located in the flexible loops of the P domain, which are highly exposed and away from the virion surface. Several viral infection studies have shown that the hypervariable region of VP2 is involved in virulence changes and that the amino acid residues at positions 253 and 284 of VP2 are the main determinants of the virulence and cell tropism of IBDV [14,15]. In other studies, residues 253, 279, and 284 were found to be important for cell culture adaptation/virulence [16,17]. Virulence declined rapidly when vvIBDV with amino acid mutations at positions 253 and 284 in the VP2 hypervariable region was passaged serially in chicken embryo fibroblast cells [18], while the very virulent pathogenicity was restored after a reverse mutation assay. Recent virulence and immunology studies have shown that the amino acid residue at position 253 is critical to the alteration of virulence of IBDV, while the residues at positions 279 and 284 are not directly related to viral virulence but are involved in cell adaptability and replication efficiency [19,20]. All these studies suggest that a few amino acid mutations at key positions in the VP2 hypervariable region can have significant impacts on the virulence of IBDV.
High-resolution structures of both VP2 and IBDV strains with different virulences are critical for a full understanding of how amino acid mutations may change the virulence of the virus. In 2005, the structure of IBDV Ct, an attenuated vaccine strain of IBDV, was resolved at 7 Å resolution by X-ray crystallography [12]. Later, the cryo-electron microscopy structure of Soroa IBDV, another attenuated strain of IBDV, was reported at 10 Å resolution [13]. Amino acid substitutions in a structural protein can markedly affect virulence, which is an important factor in viral adaptive evolution. Nevertheless, no structure of vvIBDV with substitutions at virulence sites in the hypervariable region of VP2 (e.g., Q253 and A284) is currently available. Gene sequence and evolutionary relationship analyses can provide a novel understanding of virulence proteins [21] and help to predict potential virulence and adaptive selection sites of variant viruses [22,23].
In this study, we built a pseudo-atomic structural model of VP2 from vvIBDV Gx, a very virulent strain of IBDV identified in China, and characterized the interaction of two virulence sites, Q253 and A284, by molecular dynamics simulation to reveal the molecular basis of IBDV virulence. In addition, we performed a population sequence analysis of the hypervariable region of VP2 to analyze the positive selection and co-evolution of virulence sites during the genetic diversification of IBDV. These results could shed new light on the molecular mechanisms of IBDV virulence and evolution.

Virus propagation and histopathology
Groups of 3-week-old specific-pathogen-free (SPF) chickens, 10 chickens per group, were inoculated with vvIBDV Gx (1,000 ELD 50 (median embryo lethal dose) per 0.1 mL) in the nasal mucosa and conjunctiva. Bursae of Fabricius were collected three to five days later, when chickens showed clinical lesions or were dying. Lesions were fixed and embedded in paraffin blocks. Cross sections were cut with a microtome, stained with hematoxylin and eosin, and examined under a light microscope. DMEM (Dulbecco's Modified Eagle's Medium) and an attenuated strain, IBDV Gt, were used to infect SPF chickens as negative controls. All animal experiments were approved by the Harbin Veterinary Research Institute (HVRI) of the Chinese Academy of Agricultural Sciences (CAAS) and were performed in accordance with animal ethics guidelines and approved protocols.

Virus purification
Fresh bursae were ground in 10 volumes (w/v) of phosphate-buffered saline (pH 7.2). The cell lysate was centrifuged at 4,800×g for 30 min at 4°C to remove cell debris. The supernatants were pelleted through a 20% glycerol cushion at 186,000×g for 2 h at 4°C. The virus was then purified by density gradient centrifugation on 30%-45% sucrose at 125,000×g for 12 h at 4°C and 27% cesium chloride at 125,000×g for 24 h at 4°C. Fractions containing virions were collected and dialyzed against phosphatebuffered saline (pH 7.2).

vvIBDV Gx VP2 homology modeling
The structure of vvIBDV Gx VP2 was built by the I-TASSER (Iterative Threading ASSEmbly Refinement) server [24]. The full-length vvIBDV Gx VP2 model was constructed by multiple-threading alignment and iterative template fragment assembly simulations.

Molecular dynamics simulation
Molecular dynamics (MD) simulations of the effects of amino acid mutations at positions 253 and 284 on the IBDV VP2 structure were performed with nanoscale molecular dynamics (NAMD) [25], based on the crystal structure of cvIBDV VP2 (253H, 284T) [26]. Four systems, i.e., the wild type, two single mutations (H253Q, T284A), and a double mutation (H253Q/T284A), were subjected to structural dynamics simulations over a 100 ns simulation period after a 14 ns pre-equilibration period.

IBDV VP2 sequence analysis and phylogenetic tree construction
The cDNA sequences for VP2 of variant strains of IBDV (i.e., cvIBDV, vvIBDV) deposited in GenBank as of October 2014 were collected and aligned using ClustalW (http://www.genome.jp/tools/clustalw/). Duplicate sequences were detected with DAMBE [27] and discarded prior to analyses. Strains were grouped by geographic origin: Europe, Africa, West Asia, South Asia, East Asia, North America, or South America. The consensus sequences in each group were obtained separately and the pairwise distances by proportion of nucleotide substitutions were calculated using MEGA [28]. The nucleotide substitution matrix was used to construct the phylogenetic relationships among IBDVs from sevenregions. Analysis of amino acid variability in the hypervariable region of IBDV VP2 was performed using DAMBE software.

Positive selection and co-evolution analyses
The maximum likelihood method was used to estimate the ratio of the non-synonymous substitution rate (dN) to the synonymous substitution rate (dS) codon by codon (= dN/dS) [29,30]. An  value greater than one indicates that the codon is under positive selection. Two different algorithms, FEL (fixed effects likelihood) and SLAC (single likelihood ancestor counting), were used to calculate the likelihood [31]. All analyses were performed using the software HyPhy (http://www.hyphy.org) and the online tool Datamonkey (http://www.datamonkey.org/). Spidermonkey [32], a component of HyPhy, was used to identify pairs of sites that co-evolved during IBDV evolution. Summing the probabilities of all branches gives the overall probability Q of co-occurrence of a pair of amino acid substitutions. A threshold probability of 0.9 was set, which means a pair of amino acid substitutions was regarded as co-evolving if Q was larger than 0.9.

Virulence of vvIBDV Gx
Young SPF chickens were infected with purified vvIBDV Gx virions and the bursae of infected chickens were subjected to histopathology after 3 to 5 d [18,33]. In comparison with infection by attenuated IBDV Gt, infection by vvIBDV Gx induced severe bursal lesions and resulted in ~100% morbidity and ~60% mortality (Figure 1). Pathological ex- amination of infected tissues showed severe lymphocyte depletion, bursal atrophy, connective tissue hyperplasia, and regional fibrosis. In contrast, infection with IBDV Gt induced only subclinical pathological changes, including slight hemorrhage and ecchymosis. No bursal lesions were observed in the negative control (DMEM) group ( Figure  1B). These results show that the purified vvIBDV Gx virions were highly infectious to bursae of young chickens.
VP2, the major capsid protein of IBDV, is a key determinant of virulence, cell tropism, and pathogenic phenoltype [16]. We compared the VP2 sequences of vvIBDV Gx and attenuated IBDV Gt and found that the homology between them is greater than 97%, while two amino acid substitutions, H253Q and T284A, exist in the hypervariable region of vvIBDV Gx VP2 ( Figure 1A). In a previous study in which vvIBDV Gx was attenuated to IBDV Gt by blind passage, we identified residues 253 and 284 as possible key virulence sites of vvIBDV Gx [18]. These two amino acids were selected for the molecular dynamics simulation and structural analysis in this study.

Interactions between virulence sites
We built a pseudo-atomic structural model of vvIBDV Gx VP2 by homology modeling using I-TASSER, based on the crystal structure of cvIBDV VP2 (PDB: 2DF7) (Figure 2). The overall architectures of VP2 in these two strains of IBDV are quite similar (Figure 2A). Most of the conformational differences are found in four flexible loops (P BC , P HI , P DE , and P FG ) in the jellyroll  barrels of the P domain (Figure 2B). Both virulence sites of vvIBDV Gx VP2, Q253 and A284, are exposed at the tips of loops P DE and P FG .
To determine how the amino acid mutations at positions 253 and 284 in VP2 may contribute to virulence change, we performed a simulation of the dynamics of a VP2 trimer based on the crystal structure of cvIBDV VP2. Four sys-tems, including the wild type (253H, 284T), two single mutations (H253Q, T284A), and a double mutation (H253Q/ T284A), were subjected to long-timescale MD simulations. A total simulation time of 114 ns was used for each system, allowing 14 ns for equilibration and 100 ns for the subsequent isothermal-isobaric simulation.
According to the homology model, the distance between virulence sites 253 and 284 is slightly greater in vvIBDV Gx VP2 than in cvIBDV VP2 ( Figure 2B). In the MD simulation, residues H253 and T284 of cvIBDV VP2 were found to interact with each other via a hydrogen bond (Figure 3A), which could stabilize the connection between loops P DE and P FG . In contrast, this hydrogen bond is remarkably weakened when VP2 carries the mutation H253Q ( Figure  3B and 3C) and is completely abolished by the mutations T284A and H253Q/T284A (data not shown). These results suggest that the mutations at positions 253 and 284 in vvIBDV VP2 could disrupt the hydrogen bond between these two sites and lead to the relaxation of loops P DE and P FG , which may contribute to the receptor binding and infection process of vvIBDV.

Positively selected and co-evolving sites in IBDV VP2
We collected 726 VP2 coding sequences of variant IBDV strains reported in GenBank (including 93 strains of vvIBDV) and performed a positive selection analysis to identify positively selected sites. All sequences were initially pooled together and non-synonymous/synonymous substitution rate ratios (=dN/dS) were calculated, but no significant signal of positive selection was found. We then classified the strains into seven geographic regions, i.e., Europe (78 strains), East Asia (China, Korea, and Japan; 114 strains), West Asia (33 strains), South Asia (38 stains),   (Table 1). Analysis of the phylogenetic relationships among IBDV strains by geographic region indicates that the strains in Asia are more closely related to those in Europe than to those in other areas such as the Americas ( Figure 4A).
All of the positively selected sites identified by this analysis are located in the P domain of VP2. Most of them are present in the -barrel loops, i.e., P BC , P HI , P DE , and P FG , which are found on the exterior surface of IBDV and thus  Figure 4B). Interestingly, although the amino acid variability of the hypervariable region of VP2 is similar among vvIBDV in East Asia, cvIBDV in East Asia, and vvIBDV in other regions ( Figure 5), two sites important for the virulence of IBDV, Q253 and A284, were found to be positively selected in East Asia (Table 1). We then performed a co-evolution analysis of the VP2 sequences of vvIBDV and cvIBDV strains in East Asia. The analysis identified three significant pairs of co-evolved sites: 253 and 284 (Q=0.97), 284 and 294 (Q=0.95), and 242 and 279 (Q=0.93) ( Figure 4C). Interestingly, the pair of amino acid sites with the largest Q value, 253 and 284, is the same as that identified in the positive selection analysis. These results suggest a direct relationship or contact between these two virulence sites.

Discussion
In this study, we purified vvIBDV Gx virions from the bursae of infected young SPF chickens and analyzed the virulence of this strain of IBDV. Compared with the attenuated IBDV Gt, the purified vvIBDV Gx virions are highly infectious to bursae; they have mutations at amino acids 253 and 284 in VP2, which were previously reported to be important for the virulence of IBDV [15]. To explore the molecular basis of the contribution of these two mutations to virulence change, we built a pseudo-atomic structural model of vvIBDV Gx VP2 by homology modeling and compared it with the structure of cvIBDV VP2. The overall architectures of VP2 from vvIBDV and cvIBDV are quite similar. Nevertheless, structural differences, mostly in the P domain of VP2, were found. Although the amino acid sequences of variant strains of IBDV are generally well conserved, mutations are frequently observed, especially in the hypervariable region of VP2.
To further explore the effects of the H253Q and T284A mutations on the virulence of IBDV, we performed a mo-  lecular dynamics simulation of the structural changes that occur in VP2 with different mutations. One amino acid mutation, H253Q, replaces a relatively bulky side chain (histidine) with a more linear side chain (glutamine), which may increase the flexibility of loop P DE in the P domain of VP2. The other amino acid mutation, T284A, involves a change from a hydrophilic residue (threonine) to a hydrophobic residue (alanine), which could increase the hydrophobicity of the surface of loop P FG . A more flexible conformation of loops and a more hydrophobic environment could facilitate the interaction between virus and target cell and promote stronger binding to IBDV-specific receptors.
Our MD simulation showed that there is a hydrogen bond between H253 and T284 in cvIBDV. This hydrogen bond is weakened by the mutation H253Q and completely abolished in vvIBDV VP2 with both mutations (H253Q and T284A). This disruption of the hydrogen bond between residues 253 and 284 by amino acid substitutions could affect the conformational flexibility of the loops in the P domain of VP2 and likely alter the virulence of IBDV.
The degree of virulence could be important for the evolution and adaptation of IBDV in different hosts or environments. To determine whether the amino acid substitutions at the virulence sites were positively selected during the molecular evolution of IBDV, we classified the VP2 coding sequences of variant IBDV strains into different groups by geographic origin and performed a positive selection analysis. Amino acid sites with significant signals of positive selection were found in individual geographic regions. Most of these sites are in the exposed -barrel loops of VP2 and thus could contribute to the interaction between IBDV and its specific receptors. Interestingly, the two virulence sites, 253 and 284, were found to be positively selected in East Asia. These results indicate that the virulence sites could be critical in the adaptive evolution of the virus.
Our positive selection analysis suggested that the virulence sites of IBDV may be subjected to different pressures in different geographic regions. We identified positively selected sites when we performed the analysis for six of the seven individual geographic regions, but no positively selected sites with obvious signals were identified when we performed a similar analysis using pooled vvIBDV VP2 sequences from all six geographic areas. Although the different regions exhibit similar levels of VP2 amino acid variability ( Figure 5), positively selected sites could be detected only when the IBDV strains were grouped by geographic region (Table 1). This could be due to diverse adaptations to the specific environments found in the different geographic areas.
In our positive selection analysis, amino acid 253 was identified as a target of positive selection in most geographic regions (four of six), while site 284 was identified only in East Asia; moreover, site 284 had a lower dN/dS ratio (lower positive selection pressure) than site 253 (Table 1). These results suggest that the mutation at site 253 plays a key role in determining the virulence of IBDV, while the mutation at site 284 could be population specific and may play a minor role in virulence change. We calculated the frequencies of amino acid residues at positions 253 and 284 in variant strains of IBDV. Interestingly, most vvIBDV strains (96.6%) have Q and A at sites 253 and 284, respectively, while most IBDV strains with intermediate virulence (82.4%) have Q and T at these two sites, respectively. Furthermore, all avirulent strains have H and T at these two sites ( Table 2). These results also suggest that site 253 plays a more critical role in the virulence of IBDV.
A single amino acid substitution can have a great influence on virulence; however, virulence can be affected by multiple factors as well as multi-site interactions [34]. An example of a multi-site interaction that may play a role in determining virulence is the co-evolution of different sites. Three pairs of amino acids, including pair 253/284, were identified in our co-evolution analysis as co-evolving sites that may promote changes in the virulence of IBDV.
In summary, in this study we analyzed the structure of VP2 and characterized the virulence sites 253 and 284 in vvIBDV Gx, a strain identified in China, and explored the molecular mechanism of shifts in the virulence of IBDV caused by amino acid substitutions. Amino acids 253 and 284 in the hypervariable region of VP2 were found to be both positively selected and co-evolving, which suggests a diversity of virulence determination factors. Two mutations in VP2, H253Q and T284A, increase the virulence of vvIBDV Gx in comparison with avirulent strain, but amino acid 253 plays a more critical role. A future study combining structure and sequence analysis of different strains of IBDV will provide more insights into the mechanism of virulence and molecular evolution of this virus.