Phylodynamic analysis and molecular diversity of the avian infectious bronchitis virus of chickens in Brazil.

Avian infectious bronchitis virus (IBV) is the etiological agent of a highly contagious disease, which results in severe economic losses to the poultry industry. The spike protein (S1 subunit) is responsible for the molecular diversity of the virus and many sero/genotypes are described around the world. Recently a new standardized classification of the IBV molecular diversity was conducted, based on phylogenetic analysis of the S1 gene sequences sampled worldwide. Brazil is one of the biggest poultry producers in the world and the present study aimed to review the molecular diversity and reconstruct the evolutionary history of IBV in the country. All IBV S1 gene sequences, with local and year of collection information available on GenBank, were retrieved. Phylogenetic analyses were carried out based on a maximum likelihood method for the classification of genotypes occurring in Brazil, according to the new classification. Bayesian phylogenetic analyses were performed with the Brazilian clade and related international sequences to determine the evolutionary history of IBV in Brazil. A total of 143 Brazilian sequences were classified as GI-11 and 46 as GI-1 (Mass). Within the GI-11 clade, we have identified a potential recombinant strain circulating in Brazil. Phylodynamic analysis demonstrated that IBV GI-11 lineage was introduced in Brazil in the 1950s (1951, 1917-1975 95% HPD) and population dynamics was mostly constant throughout the time. Despite the national vaccination protocols, our results show the widespread dissemination and maintenance of the IBV GI-11 lineage in Brazil and highlight the importance of continuous surveillance to evaluate the impact of currently used vaccine strains on the observed viral diversity of the country.

Introduction

Infectious bronchitis (IB) is an acute and highly contagious viral disease that affects domestic fowl (Gallus gallus) worldwide (Cavanagh, 2007). The etiological agent is the avian infectious bronchitis virus (IBV), a Gammacoronavirus from the family Coronaviridae (ICTV, 2016). IBV genome is a single positive sense RNA strand with approximately 27.6 Kb in length that encodes four structural proteins - nucleocapsid (N), membrane (M), envelope (E), and spike (S) - in addition to an RNA-dependent RNA polymerase and numerous accessory proteins (Jackwood, 2012). Among all the structural proteins, S is the most important for antigenic and immunogenic reasons. It is cleaved into the subunits S1 and S2 with approximately 535 and 625 amino acids, respectively. S1 glycoprotein is important in adsorption to the cellular receptor and virus entry into the host cell, besides inducing neutralizing antibodies. S1 gene is highly variable among the different viral strains and is directly related to the diversity of IBV antigenic and genetic groups (Cavanagh, 2007; Toro et al., 2012).

IBV genetic diversity, mainly in the S1 gene, was demonstrated in different poultry-producing regions of the world. Historically, Massachusetts (Mass) and Connecticut (Conn) serotypes were the first isolates in the 1940s and 1950s, respectively (Schalk and Hawn, 1931; Jungherr et al., 1956). Since then, several new IBV sero/genotypes have been identified and associated with the disease around the world (Cook et al., 2012; De Wit et al., 2011). However, for decades, the identification and classification of IBV genetic types were mostly performed without clear criteria regarding the nomenclature, methods to compare viral molecular data and exact genetic region to be analyzed. Consequently, subclades of the closely related virus, sometimes circulating in small geographical regions, have been assigned as new genotypes and more than 50 genetic groups were already reported (De Wit et al., 2011). Recently, a more definitive phylogeny-based classification system was proposed for the assignment of IBV strains. This study classified IBV in six genotypes, which are further divided into 32 genetic lineages, and provides reliable reference sequences (lineage prototypes) to guide the viral classification (Valastro et al., 2016).

Currently, Brazil is the main exporter and second-biggest producer of chicken meat in the world (USDA, 2016). In 2016, more than 13 million tons of chicken was produced and over 4 million tons of chicken products were exported to all continents. This high production is obtained in intensive raising systems, favoring the dissemination of respiratory infections (Bermudez, 2008) such as IBV, which causes great economic losses in Brazilian poultry flocks (broilers, breeders and layers) (Balestrin et al., 2014; Carranza et al., 2017; Colvero et al., 2015).

IB has been described in Brazil since the 1950s. However, a better genetic characterization of the field IBV isolates has started to be performed only in the last decade. It is known that a Brazilian variant (BR-I) is widely disseminated in the main poultry-producing regions of the country. Similar strains were also observed to circulate in Argentina and Uruguay and the whole genetic cluster was also identified as South America I (SA-I) (Marandino et al., 2015) and recently renamed as GI-11 (Valastro et al., 2016). Moreover, subclades of GI-11 were reported by distinct studies (Fraga et al., 2013; Villarreal et al., 2010) but the origin of these lineages and the role of recombination are still to be investigated.

Moreover, IBV molecular diversity in Brazil has been investigated by using diverse methods on different regions of S1 gene, making it difficult to compare results from different studies and impairing epidemiological surveillance efforts to track down IB outbreaks. In the view of the recently proposed system of IBV classification, the present study investigated the lineages circulating in Brazil and the role of recombination for the current observed genetic diversity. This study also applied phylodynamic methods to estimate the time of the most common recent ancestor and demographic history of the IBV field variants in Brazil and related sequences from South America.

Material and methods

Sequence dataset compilation and maximum likelihood (ML) analysis

All available IBV S1 gene sequences from Brazil were downloaded from GenBank. Alignment was performed with Mafft (Katoh and Standley, 2013) and visually inspected in AliView (Larsson, 2014). A reference sequence dataset for the genotypic classification of IBV was used as provided by Valastro et al. (2016). Due to different sizes in length and sequences that cover distinct and not overlapping regions of S1 gene, the genotyping of Brazilian sequences was performed in separate datasets when necessary.

To further analyze the global circulation of IBV strains isolated in Brazil, all international IBV S1 sequences with information for country and date of sampling were downloaded from GenBank. Alignment was performed with Mafft and trimmed aiming to keep the highest number of sequences from Brazil. RAxML (Stamatakis, 2014) was used to remove identical sequences and construct maximum likelihood (ML) trees. The general time reversible model (GTR) with gamma-distributed rate heterogeneity plus a proportion of invariable sites (GTR+G+I) was used as the optimal nucleotide substitution model as identified in the jModelTest program (Posada, 2008).

Recombination analyses

Analyses of recombination were performed for sequences grouped within the GI-11 lineage. The S1 gene fragment analyzed corresponds to nucleotide positions 8 to 550, according to the H120 reference strain (M21970 - accession code). Simplot software (Lole et al., 1999) was used applying the bootscanning method. Neighbor-Joining (NJ) trees were constructed under Kimura two-parameter model with sliding windows of 100, 160 and 200 base pairs (bp) with incremental steps of 20 bases. Query sequences were compared against reference sequences for each lineage defined by Valastro et al. (2016). S1 gene regions showing patterns of recombination were used as a query in a BLAST search to identify the source of potential recombination fragments. The top 10 hits for each query were downloaded and added to the sequence dataset in case they had not been yet analyzed.

Phylodynamic analyses

The temporal signal of the sequences to be submitted to phylodynamic analysis was investigated with TempEst software (Rambaut et al., 2016). Sequences outliers in the regression of root-to-tip divergence versus sampling time were excluded. Time-scaled phylogenetic tree estimation was performed using BEAST/BEAGLE software (Ayres et al., 2012; Drummond et al., 2012) through the Cipres Science Gateway (https://www.phylo.org). BEAST software allows for the combination of different clock, substitution, and demographic models, demanding an appropriate model test approach. In the current study, marginal likelihood estimation (MLE) (Baele et al., 2012, Baele et al., 2013) was applied to compare alternative models in a Bayesian framework. Trees were reconstructed using SRD06 substitution model (Shapiro et al., 2005) and the uncorrelated gamma distributed (ucgd) relaxed molecular clock (Drummond et al., 2006), which outperformed alternative models. IBV demographic history in Brazil was investigated by applying the non-parametric Bayesian Skygrid coalescent model, which estimates the product of viral effective population size (N) and generation time throughout evolutionary history (Gill et al., 2013). To avoid making assumptions regarding IBV generation time, here we refer to estimates of effective population size as relative genetic diversity. In addition, we tested with MLE the best demographic parameter that described the IBV population history.

Monte Carlo Markov Chains (MCMC) were run sufficiently long to ensure stationary and adequate effective sample size (ESS) for the main parameters. Tracer software (available at: http://beast.bio.ed.ac.uk/Tracer) was used to diagnose MCMC, adjust initial burn-in and to perform the Skygrid demographic reconstruction. TreeAnnotator was used to summarize the maximum clade credibility (MCC) tree from the posterior distribution of trees and the MCC tree was visualized and edited in FigTree (available at: http://tree.bio.ed.ac.uk/software/figtree/).

Results

The molecular diversity of IBV in Brazil

In order to classify the IBV lineages circulating in Brazil, a preliminary analysis was performed including all Brazilian IBV sequences available in GenBank and the reference dataset provided by Valastro et al. (2016). A total of 192 IBV Brazilian sequences for S1 gene was obtained. These sequences varied in size and covered different regions of S1 gene. A total of 140 sequences (73%) covered the first portion of S1 region, approximately between positions 190 and 700 (including hypervariable regions 1 and 2 located between nucleotides 112 and 423); 27 sequences (14%) covered the middle of the gene between positions 740 and 1090 (including hypervariable region 3 located between nucleotides 820 and 1161); and 25 sequences (13%) spanned the whole S1 gene. ML tree reconstruction using Valastro et al. (2016) reference sequences classified 143 (74.5%) Brazilian sequences as GI-11 lineage (former SA-I group), 46 (24%) as GI-1 (former Mass serotype), 2 (1%) as GI-13 (former 4/91 serotype) and 1 (0.5%) as GI-9 (former Ark serotype).

Phylogenetic analyses with sequences from worldwide were performed for the first portion of S1 gene, which was observed to be the most common sequenced region. After excluding sequences too short or identical, a ML phylogenetic tree with 1634 taxa (199 reference sequences and 1435 worldwide) was constructed (Fig. 1 ). GI-11 clade was composed by Brazilian, Argentinian and Uruguayan sequences and basal to that, a group of sequences from Argentina was observed (Fig. 2 ). These sequences did not group with other reference lineages and were analyzed together with GI-11 lineage aiming to assess the role of recombination in this basal group.

Fig. 1

Circular ML phylogenetic tree of IBV partial S1 gene. Sequences from worldwide were selected based on the availability of date and country of sampling information. Brazilian sequences are the green dots in external nodes and international sequences are the purple dots in external nodes. IBV Genotype I (GI) lineages are indicated in different colors. The classification was performed as proposed by Valastro et al. (2016) by using reference sequences provided in the study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

ML phylogenetic tree of partial IBV S1 gene for GI-11 lineage and related sequences and bootscanning analysis. The left panel shows a phylogenetic tree of GI-11 lineage and basally grouped sequences. The red star marks the most recent common ancestor of GI-11 lineage. Branch support (SH-like) is shown for the relevant clades and sequences' labels are colored according to the country of origin of the sample. The right panel shows bootscanning plots of three clusters of sequences. Y-axis represents the percentage of permutated trees in which the query sequence grouped with one of the reference lineages. X-axis represents different positions of the analyzed genomic region according to the sliding-window procedure. Each line in the plot represents one IBV lineage and colors are according to the legend in the right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Bootscanning analyses were applied to all sequences that belonged to the GI-11 lineage as well as for the six Argentinian sequences that were grouped basally to it. A possible recombination was observed in the six sequences from Argentina. However, it was not possible to determine which lineage gave rise to it. Analyses showed that this basal group has some similarities to the GI-11 lineage, but the sequences mostly did not resemble with high support (above 70% bootstrap) any of the lineages already described (Fig. 2).

Two other divergent clusters of sequences within GI-11 lineage showed distinct patterns in the bootscanning analysis. The first group is composed of sequences isolated in Mato Grosso state (MT) in Brazil, previously described as a subcluster of the Brazilian lineage (Fraga et al., 2013). Our analyses showed that these sequences are mostly GI-11 but have a fragment of around 250 bp that does not match with any other described lineage. Aiming to identify the origin of this fragment, a BLAST search was performed using the unassigned portion of S1 gene as a query. This procedure gave rise to two sequences that were not present in our compiled dataset and although these sequences were isolated in a distinct state, they also clustered with the group isolated in MT state and showed the same recombination pattern (Fig. 2). The second group with high divergence within GI-11 lineage is composed only of sequences from Argentina and Uruguay. Bootscanning revealed that the analyzed S1 region was mostly related to GI-11 cluster, with a small region with very low similarity to any of the described lineages. Analysis performed with bigger window sizes (Supplementary Fig. 1), however, shows small support for recombination in this region. BLAST search using the fragment between 100 and 150 bp did not reveal any new sequences that could indicate the source of this potential recombinant fragment.

In order to understand the origin of the GI-11 cluster and its onset in Brazil and South America, Bayesian analyses were performed with GI-11 clade sequences. TempEst analysis showed that good temporal signal (R2 = 0.48) was only present in the GI-11 dataset when the subclusters of potential recombinant sequences were removed. Phylodynamic analysis of the “pure” GI-11 lineage sequences estimated the time of the most recent common ancestor (tMRCA) of this group to 1951 (1917 to 1975, 95% highest posterior density [HPD]) (Fig. 3 ). The mean rate of evolution of the analyzed fragment of S1 gene was 4.1 × 10−3 substitutions/site/year (2.3 × 10−3 to 6.1 × 10−3, 95% HPD). Skygrid plot reconstructed a relative genetic diversity trajectory with a small variation since the introduction of IBV GI-11 lineage in Brazil in the beginning of the 1950s (Fig. 3). Among the parametric demographic models, the constant population was estimated as the best-fitted parameter. Bayes Factor (BF) calculations supported this model over expansion or exponential growth with BF = 2.4 and 3.7 (respectively), revealing a weak support. The analysis testing the logistic growth model could not be completed due to very poor model fit.

Fig. 3

Time-scaled tree and virus relative genetic diversity history of IBV GI-11 lineage. The Bayesian phylogenetic tree was estimated for sequences with no signal of recombination in the bootscanning analysis. Skygrid plot is superimposed to the tree in the same time-scale. Relative genetic diversity is represented here as a surrogate for the product of effective population size and generation time (see Materials and methods for details). The solid red line represents the mean relative genetic diversity and the red area around the line represents the 95% HPD of the estimates. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

IB is an endemic disease in Brazil and it is responsible for significant and measurable economic losses in Brazilian poultry flocks (Colvero et al., 2015). Despite the massive use of a vaccine derived from the Mass IBV serotype, outbreaks have occurred frequently in Brazil (Carranza et al., 2017; Fernando et al., 2017; Mendonça et al., 2009), as well as in some neighboring countries (Alvarado et al., 2005; Marandino et al., 2015; Rimondi et al., 2009). Due to the wide range of antigenically and genetically different viral types of IBV, the control and prevention of these outbreaks are challenging, making necessary a constant surveillance of the lineages circulating in the country.

In the present study, we have re-assessed the IBV molecular diversity in Brazilian isolates based on a recent update of the genetic classification for this virus (Valastro et al., 2016). Our analyses showed that 46 sequences (24%) isolated in Brazil were classified as GI-1 lineage, which corresponds to the Mass serotype. Since the end of the 1970s, a vaccine derived from this serotype has been extensively used in the field in Brazil. The Mass vaccine is available as inactivated or attenuated-live, but the second one is commonly used in the field poultry flocks. Previous studies reported that this vaccine strain could spread in the flock for few weeks, due to the intermittent shedding of virions from trachea and cloaca (Matthijs et al., 2008; Naqi et al., 2003). This ability of the vaccine strain to disseminate could be thought as risk to turn it into a new field strain. However, different studies demonstrated a very high identity among all GI-1 lineage sequences obtained from field flocks and Mass commercial vaccines distributed in Brazil suggesting they are the vaccine strains (Balestrin et al., 2014; Chacón et al., 2011; Fraga et al., 2013; Villarreal et al., 2010). Therefore, the high frequency of Mass serotype observed might be due to a local viral spread after the intensive vaccination practices.

A small percentage of Brazilian sequences were also identified as being GI-9 (0.5%) and GI-13 (1%) lineages, corresponding to the Arkansas and 4/91 (also known as 793B and CR88) serotypes, respectively. 4/91 strain is spread worldwide and was first reported in the United Kingdom in 1991 (Gough et al., 1992; Parsons et al., 1992), while Arkansas strain has been reported in the United States since 1972 (Johnson et al., 1973). In Brazil, the lineages Ark and 4/91 were isolated from vaccinated flocks (Mass serotype vaccine) in 1999 and 2007, 2008, respectively (Montassier et al., 2008; Villarreal et al., 2010). The Mass serotype was the only licensed attenuated strain to be used in field vaccination programs in Brazil until 2016 when the Brazilian GI-11 strain started to be used as an attenuated-live vaccine. However, both lineages Ark and 4/91 are used in vaccination protocols in other countries and the frequencies observed in our study could be related to the use of imported vaccines or for the importation of birds from places where those strains are present (Felippe et al., 2010; Villarreal et al., 2010). On the other hand, the fact that these strains were not found in samples isolated recently suggests that these are self-limiting cases of less epidemiological importance.

Most of the Brazilian sequences analyzed here (74.5%) were classified as GI-11 lineage. These findings are similar to the previous studies performed in order to classify the field isolates (Carranza et al., 2017; Chacón et al., 2011; Fraga et al., 2013). Our study also observed that a subclade in the GI-11 lineage presented possible recombination breakpoints in the bootscanning plots. This subclade is composed only of sequences from Brazil and was first named as BR-II (Fraga et al., 2013). The current study suggests it also belongs to GI-11 lineage. The bootscanning procedure also revealed that a subclade of sequences from Argentina and Uruguay might present recombinant fragments, however, in this case, the support of 160 and 200 bp window sizes were weak and new analyses with longer sequences should be performed to confirm our findings. A third potential recombinant group of sequences grouped basally to the GI-11 lineage and our results show that these sequences have indeed a fragment or two of GI-11 lineage in their genome. The backbone genome of these sequences was mostly unclassified, or at least not represented among the list of reference sequences used in this study.

Mutation and recombination events are very common in coronaviruses, such as IBV. These events could be responsible for the emergence of this diverse number of variant strains in the field (Jackwood, 2012). Due to the advance of the bioinformatics programs for recombination screening, several studies have been performed in order to evaluate if the recombinant events could be the source of these variant isolates described worldwide (Kiss et al., 2016; Naguib et al., 2016; Quinteros et al., 2016). Moreover, the high prevalence of avian coronaviruses circulating in birds other than chicken suggests that wild and synanthropic birds are potential disseminators and could play a role in the emergence of recombinant strains and the dissemination of IBV among distant regions (Durães-Carvalho et al., 2015). Therefore, it might be possible that wild birds are the source of the unclassified sequence fragments observed in the recombinant strains identified in the current study. Notwithstanding, our BLAST search approach would be able to identify the recombinant parental strains if their sequences were present in the public databases. Lastly, is important to note that regions in the bootscanning plots that do not resemble any of the reference IBV lineages might be simply originated by single nucleotide mutations and are not necessarily resulted from events of recombination.

The phylodynamic analysis estimated the introduction of IBV GI-11 lineage in Brazil in the beginning of the 1950s. The first reported case of IB in poultry flocks in the country was performed by Hipólito in 1956 (Hipólito, 1957). According to this author, the introduction of the IBV in Brazil, and probably into South America, could be explained by the importation of chickens from North American and European countries, where the disease was endemic in that time. The poultry production underwent major changes during the 1950s in Brazil, with the beginning of the industrial and intensive broiler rearing systems, which allowed an important expansion of the poultry production (Belusso and Hespanhol, 2010; Sorj et al., 2008). On the other hand, the intensive rearing systems and the high population of the flocks may allow the quick dissemination of respiratory diseases, which was common at the time (Hipólito, 1957).

The MLE analysis revealed a constant model as best describing the demographic history of IBV GI-11 lineage in Brazil. However, BF support was only moderate (BF = 3.7) over the exponential model and weak (BF = 2.4) over the expansion model, that might suggest a more complex scenario where fluctuations in the population size do not follow any of the tested parameters. By using the non-parametric coalescent Skygrid model, a smooth increase in relative genetic diversity can be observed since its introduction in Brazil until the late 1970s, when the vaccination based on the Mass serotype started. From this point onward, the Skygrid plot fluctuates around the same value, but we cannot observe a trend of fall in relative genetic diversity, suggesting a partially effective vaccination program. It is not uncommon for some IBV vaccines to induce poor protection against different serotypes of the virus (Ladman et al., 2006), which allows some viruses to scape of the immune system response provided by vaccination protocols. Most of the sequences analyzed here were originated from samples collected from vaccinated birds, that might support the idea that the Mass serotype has none or poor efficiency to control the GI-11 replication in the field, as indicated by studies of IBV experimental infection (Fernando et al., 2013, Fernando et al., 2017). Moreover, the oldest available Brazilian GI-11 sequence was originated from a sample collected in 1975. This, together with the GI-11 lineage estimated tMRCA (1951 [1917 to 1975, 95% HPD]), indicates that this lineage was already circulating in the country by the time when vaccination started. A very limited or no protection of Mass vaccine against these field variants may have allowed the spread and maintained of GI-11 lineage in Brazilian flocks, reflecting in an approximately constant viral population from 1980 to 2014. The constant viral population scenario revealed by our analyses, however with not strong BF support, also provide evidence that the vaccination based on the Mass serotype did not impact the GI-11 clade demographic history in Brazil. Franzo et al. (2016) have reported that a vaccination program based in a combination of Mass and 4/91 strains was more effective in controlling the spread of the GI-19 lineage when compared to protocols using only Mass vaccines. In this sense, since 2016, a new attenuated vaccine based on GI-11 (BR-I) lineage is available to be used in Brazilian commercial flocks and future studies are needed to evaluate changes in the viral molecular diversity and demographic dynamics of IBV in Brazil. It is important to highlight that the changes in viral relative genetic diversity estimated by the Skygrid model are subtle and the confidence interval of the analysis is wide, so only studies with more sequences could precisely describe the Brazilian IBV demographics' history and whether it is in expansion or controlled.

Conclusions

In summary, our study brings new insights about IBV genetic diversity in Brazil and unveils the evolutionary history of the GI-11 lineage in the country. Our bootscanning analysis highlights the importance of the recombination process in creating new variants of IBV in the region and the phylodynamic analyses reveal, for the first time, the virus demographic history in Brazil, which broadly agree with the country's history of poultry production and vaccination programs. Finally, a continuous surveillance of IBV infection would be valuable to ascertain long-term effects of GI-11 lineage vaccine introduction.

Conflicts of interest

The authors declare that they have no conflicts of interest.

Funding

This work was supported by the National Council of Technological and Scientific Development – CNPq (313564/2014-0; 313304/2014-9), Coordination of Improvement of Higher Educational Personnel – CAPES (181 – 18/12/2012) and the Rio Grande do Sul Research Foundation – FAPERGS (001/2013 – PqG, 16/2013 – TECNOVA). Tiago Gräf in 2016/2017 received funding from the VIROGENESIS project as part of the European Union's Horizon 2020 Research and Innovation Programme under grant agreement 634650 and from a Flagship grant from the (MRC) of the Republic of South Africa (MRC-RFA-FSP-01-2013/UKZN HIVEPI).