Whole-genome characterization of Uruguayan strains of avian infectious bronchitis virus reveals extensive recombination between the two major South American lineages.

Infectious bronchitis virus (Gammacoronavirus, Coronaviridae) is a genetically variable RNA virus that causes one of the most persistent respiratory diseases in poultry. The virus is classified in genotypes and lineages with different epidemiological relevance. Two lineages of the GI genotype (11 and 16) have been widely circulating for decades in South America. GI-11 is an exclusive South American lineage while the GI-16 lineage is distributed in Asia, Europe and South America. Here, we obtained the whole genome of two Uruguayan strains of the GI-11 and GI-16 lineages using Illumina high-throughput sequencing. The strains here sequenced are the first obtained in South America for the infectious bronchitis virus and provide new insights into the origin, spreading and evolution of viral variants. The complete genome of the GI-11 and GI-16 strains have 27,621 and 27,638 nucleotides, respectively, and possess the same genomic organization. Phylogenetic incongruence analysis reveals that both strains have a mosaic genome that arose by recombination between Euro Asiatic strains of the GI-16 lineage and ancestral South American GI-11 viruses. The recombination occurred in South America and produced two viral variants that have retained the full-length S1 sequences of the parental lineages but are extremely similar in the rest of their genomes. These recombinant virus have been extraordinary successful, persisting in the continent for several years with a notorious wide geographic distribution. Our findings reveal a singular viral dynamics and emphasize the importance of complete genomic characterization to understand the emergence and evolutionary history of viral variants.

Introduction

Infectious bronchitis (IB) is still one of the most persistent respiratory diseases of chickens despite intensive control and research efforts made since its first description in 1931 (Shalk and Hawn, 1931). In regions where there is no highly pathogenic avian influenza virus or velogenic Newcastle disease virus, infectious bronchitis virus (IBV) is the cause of the major economic losses to the poultry industry (Cook et al., 2012).

IBV belongs to the genus Gammacoronavirus within the Coronaviridae family (de Groot, 2012). The positive-sense single-stranded RNA viral genome (27.6 kb) contains six genes and at least ten open reading frames (ORFs) in the order 1a-1b-S-3a-3b-E-M-5a-5b-N. The 1a and 1b ORFs code for the 1ab polyprotein using a ribosomal frameshift event at the end of ORF la; the large 1ab polyprotein is then post-translationally cleaved into 15 non-structural polypeptides, which are required for RNA replication and transcription. ORFs S, E, M and N encode four major structural proteins: the spike (S) glycoprotein, the small envelope (E) protein, the membrane (M) glycoprotein and the nucleocapsid (N) protein, respectively. ORFs 3a, 3b, 5a and 5b encode small non-structural proteins, which are not essential to viral replication (Cavanagh, 2007).

The IBV genome has high rates of mutation and recombination, leading to the continuous emergence of novel genetic and antigenic variants worldwide (Gough et al., 1992, Liu and Kong, 2004). Strain classification and evolutionary analysis of IBV variants are usually accomplished by the phylogenetic analysis of the complete coding region of the S1 subunit, a highly variable region of the S protein. The S1 subunit carries most of the virus-neutralizing epitopes and has a direct implication in protective immunity (Cavanagh and Davis, 1986, Cavanagh et al., 1992, Mockett et al., 1984). The most recent comprehensive classification identified six main genotypes (GI–GVI), 32 viral lineages (1 − 32) and a number of inter-lineage recombinants in global strains (Valastro et al., 2016).

Two lineages of the GI genotype (11 and 16), also called South America I (SAI) and Asia/South America II (A/SAII) genotypes, have been widely circulating for decades in South America (Alvarado et al., 2005, Marandino et al., 2015, Rimondi et al., 2009). These lineages are responsible for most outbreaks in commercial flocks and cause major economic losses in countries with high production levels, such as Brazil and Argentina, which are among the world's top producers of broiler meat (USDA, 2014). GI-11 is an exclusive South American lineage that emerged in the 1960s and is now apparently restricted to Argentina, Brazil and Uruguay (Chacon et al., 2011, Felippe et al., 2010, Marandino et al., 2015, Villarreal et al., 2007). The GI-16 lineage contains viruses collected in South America (Argentina, Chile, Colombia, Peru and Uruguay), Asia (China and Taiwan) and Europe (Italy, Russia and Slovenia) (Alvarado et al., 2005, Han et al., 2011, Huang et al., 2004, Liu et al., 2009, Marandino et al., 2015, Sesti et al., 2014, Tataje-Lavanda et al., 2016, Toffan et al., 2013, Valastro et al., 2016). Strains of this lineage were first reported in China in 1996, but an Italian strain from 1986 was later described (Valastro et al., 2016).

In order to understand IBV biology, in addition to the analysis of the S1 region, it is desirable to obtain complete genome sequences of different viral population. Genomic studies provide information about the origin and spread of IBV genetic variants (Naguib et al., 2016) and help to identify the genetic changes associated with inter and intra-continental spreading. So far, there are not genomes available for any South American strains of IBV. In the present study, we obtained full-length genome of two South American strains of the GI-11 and GI-16 lineages, and performed phylogenetic analyses to determine the evolutionary relationships between these strains and worldwide field strains.

Materials and methods

Strains

The UY/11/CA/18 and UY/09/CA/01 strains were isolated from IBV outbreaks in commercial broilers during 2009 and 2011 from Uruguay. These strains were previously classified within the GI-11 and GI-16 lineage by full-length S1 analysis using Sanger sequencing (Marandino et al., 2015, Valastro et al., 2016).

Propagation and purification of IBV particles

The strains were propagated in 10-day-old embryonated chicken eggs by allantoic route inoculation to obtain large quantity of IBV virions. Forty mL of allantoic fluid was harvested 72 h post inoculation and then concentrated through a 20% (wt/vol) sucrose cushion. Viral pellets were re-suspended in 1200 μl of phosphate buffered saline (PBS) and the virus suspension was layered over a continuous gradient of 30 to 55% (wt/vol) sucrose in PBS. The gradient was centrifuged at 100,000g for 4 h at 4 °C.

RNA extraction and Illumina sequencing

RNA extraction was performed using the Quick-RNA™ MiniPrep kit (Zymo Research, Irvine, CA, USA) and 200 μl of purified viral particles. The extracted RNA was eluted in 35 μl of RNase-free water. Reverse transcription was carried out using the Maxima H Minus Double-Stranded cDNA Synthesis kit (Fermentas Life Sciences Inc., Hanover, MD, USA) and 13 μl of extracted RNA. Viral RNA enrichment was determined by quantitative PCR (Callison et al., 2006). Nextera XT DNA Sample Preparation kit (Illumina, USA) was used from 1 ng of double strand cDNA. After 12 PCR cycles, the final library was purified with AMPure XP (Benchman, USA) and quantified with the Qubit dsDNA HS assay kit (Invitrogen, USA). Quality and length of the library were assessed with the Agilent high-sensitivity DNA kit (Agilent, USA) using an Angilent 2100 Bioanalyzer (Agilent, USA). Sequencing was performed on an Illumina MiSeq (Illumina, USA) platform at the Institute Pasteur, Montevideo, Uruguay.

Genome assembly and annotation

The raw reads generated by the Illumina sequencer were analyzed using the Galaxy web-based platform (Goecks et al., 2010). Paired-end reads were converted to FASTQ Sanger format using FASTQ Groomer (Blankenberg et al., 2010). The read quality was assessed by FastQC (Andrews, 2010) and adapter sequences were trimmed and low-quality reads were removed from the raw reads using Trim Galore (Krueger, 2012). Reads were mapped to the closest genomic reference IBV strains (accession numbers: EU418976 and KP780179) using Burrows–Wheeler Alignment tools with Galaxy's default settings. Alignments were visually inspected and manually optimized using Integrative Genomics Viewer (IGV) (Robinson et al., 2011). The sequences were deposited in the GenBank database under the accession numbers MF421319 and MF421320. The ORF prediction was carried out in http://covdb.microbiology.hku.hk and http://www.jcvi.org/vigor/ (Huang et al., 2008, Wang et al., 2010).

Genome sequence analysis

Sequences were aligned using MUSCLE algorithm implemented in MEGA 5.0 (Tamura et al., 2011). The best-fit model of nucleotide substitution was selected under the Akaike information criterion and Bayesian information criterion as implemented in jModelTest.

Identification of potential recombinant and parental sequences and localization of possible recombinant breakpoints were performed using the RDP4 program, which implements seven distinct algorithms for characterization of recombinant sequences (Martin et al., 2015).

Maximum-likelihood trees, with approximate likelihood ratio tests for internal nodes support, were inferred using PhyML. Phylogenetic trees were visualized and edited with Figtree (Rambaut, 2012).

Results

Purification of IBV particles and genome sequencing

Purified virions were obtained from the UY/11/CA/18 and UY/09/CA/01 strains by gradient centrifugation; viral particles were observed as a single clear band and obtained by fractionation of the gradient. The virion-extracted RNA had a concentration of 9.83 × 106  RNA copies/μl for the UY/11/CA/18 strain and 1.35 × 106  RNA copies/μl for the UY/09/CA/01. The Illumina sequencer for these strains generated 2.3 × 106 and 3.7 × 106 reads, respectively; a low percentage of these reads matched with chicken genome (2.15 and 3.83%), while 86.4 and 86.8% of them matched with a reference IBV strain (quality control statistics is available as Supplementary material). The reads were assembled and the genome sequence of both strains were obtained; the sequences included the coding genome, the UTRs, and the poly-A tail. The complete genome sequence of UY/11/CA/18 and UY/09/CA/01 strains have a size of 27,621 and 27,638 nucleotides respectively, excluding the poly(A)-tail. Both South America strains have thirteen ORFs (1a-1b-S-3a-3b-E-M-4b-4c-5a-5b-N-6b). The ORFs have the same length, except for the ORFs 1a and S that differ in 24 and 36 nucleotides, respectively.

Strain classification

A full-length genome dataset was built using complete IBV genomes available in the GenBank database (n  = 158).

Strain classification in genotypes and lineages was performed by S1 phylogenetic analysis with non-recombinant S1 sequences retrieved from the complete genomes and prototype IBV strains following Valastro et al. (2016) (Supplementary Fig. 1). The UY/11/CA/18 and UY/09/CA/01 strains clustered within the GI-11 and GI-16 lineages, respectively.

Phylogenetic incongruence and recombination

The South American strains were analyzed using incongruence between phylogenetic trees derived from the essential ORFs present in all IBV strains (1a, 1b, S, 3a, 3b, E, M, 5a, 5b and N). Contiguous ORFs with the same phylogenetic relationship were considered as single units for the analysis, and denoted as regions I to V starting from the 5′ end.

In the ORF1a phylogeny (region I), the South American strain UY/09/CA/01 (lineage GI-16) clusters with the South American strain UY/11/CA/18 (lineage GI-11) with an average nucleotide identity of 96% (Fig. 1a). Other strains of the GI-16 lineage, the Asiatic strain CK/CH/LDL/97I (accession number: JX195177) and the European strain CoV/Ck/Italy/I2022/13 (accession number: KP780179), fall in different positions in the tree and have an average nucleotide identity ranging from 90 to 92%.

Fig. 1

Phylogenetic trees obtained with the maximum-likelihood method and JC model (a), GTR model with gamma distribution and invariant sites (b and c), and HKY model with gamma distribution and invariant sites (d and e). Phylogenetic reconstruction was carried out using the ORF1a (a), ORF1b (b), S1 sequence (c), S2 sequence and ORFs 3a, 3b, E, M, 5a and 5b (d) and ORF N (e). Mapping uncertainties for internal nodes are shown as approximate likelihood ratio test values. Strains of GI-16 and GI-11 lineages are labeled with an asterisk. This analysis was performed using a sub-sampled non-recombinant data set with representatives of the IBV variability.

In the ORF1b phylogeny (region II), strains of the GI-16 and GI-11 lineages form a monophyletic group (Fig. 1b). The nucleotide identity of this group ranges from 97 to 98%.

In the phylogenetic tree based on the S1 coding region of ORF S (region III), the UY/09/CA/01 strain clusters with Asiatic and European strains of the GI-16 lineage, while the UY/11/CA/18 strain appears separated as the only complete-genome representative of the GI-11 lineage (Fig. 1c).

The phylogenetic trees based on the S2 coding region and ORFs 3a, 3b, E, M, 5a and 5b (region IV) have the same topology for the GI-11 and GI-16 lineages and were then concatenated to produce a single phylogeny. As observed for the ORF1b phylogeny, the GI-11 and GI-16 strains cluster together in a monophyletic group with 94˗98% nucleotide identity (Fig. 1d).

The phylogenetic tree based on the ORF N (region V), encoding the N protein, showed that the South American strains form a monophyletic group, with a nucleotide identity of 97%. The Asiatic and European strains of GI-16 lineage form another monophyletic group clearly separated from the South American strains (Fig. 1e).

The full-length sequence of the UY/11/CA/18, UY/09/CA/01, CK/CH/LDL/97I and CoV/Ck/Italy/I2022/13 strains were further examined for the detection of recombinant events. Six RDP4 algorithms detected four recombinant breakpoints in South American strains, located at the end of ORF1a (nt 12,319), at both ends of the S1 coding region (nt 20,330 and 21,854) and the beginning of the ORF N (nt 25,900) (Fig. 2 ).

Fig. 2

Genomic organization of the GI-16 and GI-11 strains. Essential (1a, 1b, S, 3a. 3b, E, M, 5a, 5b and N) and non-essential ORFs (4b, 4c and 6b) are indicated. European strain does not have ORF6b. Recombinant points involved in the emergence of the South American strains are labeled with a head arrow. The five regions of the South American mosaic genomes are delimited by recombinant points (indicated with different colors). Divergence of the sequences are indicated by different shade of the same color.

Discussion

The traditional sequencing methodology for IBV genomes requires the producing of several overlapping PCR fragments followed by Sanger sequencing (Abro et al., 2012, Mondal and Cardona, 2007). This approach is particularly laborious because coronaviruses have the largest known RNA genomes (Lai and Cavanagh, 1997). High-throughput or deep sequencing approaches drastically reduced time and cost requirements for genomic sequencing. Additional advantages of this methodology are the significant increase in sequence coverage and depth that permit to study co-infection, recombination, and quasispecies diversity (Isakov et al., 2015, Pérez et al., 2014).

Relatively few studies have used high-throughput sequencing to obtain IBV genomes (Abolnik, 2015, Naguib et al., 2016, Quinteros et al., 2015, Reddy et al., 2015). All methods use viral enrichment in allantoic cavity of embryonated chicken eggs and different protocols to purify IBV particles, including RNase treatment of allantoic fluid, sucrose gradient or genome amplification with commercial kits. The purification of IBV particles with sucrose cushion and gradient along with Illumina sequencing technology were here used with good effectivity and extremely high coverage values. This standardized protocol allowed the straightforward obtainment of the full-length genome of the UY/11/CA/18 and UY/09/CA/01 strains, including the coding genome and UTR regions. These strains were classified as belonging to the GI-11 and GI-16 lineages, supporting previous analysis using S1 analysis with Sanger sequencing (Marandino et al., 2015).

The GI-11 and GI-16 lineages are good examples of successful genetic variants that have been widely circulating for decades in South America (Marandino et al., 2015). GI-16 is also one of the four most distributed lineages of IBV, together with GI-1 (Massachusetts), GI-13 (793B) and GI-19 (QX) (Valastro et al., 2016). Despite their wide distribution, extended field persistence, and economical relevance, the genomic evolution of these lineages is poorly understood. The strains here sequenced are the first obtained in South America and provide meaningful contributions toward our understanding of IBV variability in the continent.

South American strains have a genome with thirteen ORFs (Fig. 2), including the ORFs 4b, 4c and 6b which are quite variable in both presence and absence, and in their nucleotide sequence. The same organization was observed in the Asiatic strain of the GI-16 lineage but the European strain of this lineage lacks the ORF for the 6b putative protein. The role of ORF6b in IBV remain uncertain but the ORF6b homologue of SARS coronavirus was identified as an endoplasmic reticulum/Golgi membrane-localized protein that induces apoptosis (Ye et al., 2010). This protein accelerated the replication of murine coronavirus, increasing the virulence of the original attenuated virus (Tangudu et al., 2007).

The different clustering observed in the phylogenetic analysis shows that the UY/11/CA/18 and UY/09/CA/01 strains have a mosaic genome formed by regions with different evolutionary relationships. The region III (S1 sequence) is clearly divergent between them and lead to the classification of the strains in different lineages (GI-11 and GI-16) (Fig. 1c). Region II (ORF1b) and region IV (S2 sequence and ORFs 3a, 3b, E, M, 5a, 5b) of the South American genomes examined in this study are closely related with Euro-Asiatic strains of GI-16 lineage (Fig. 1b and d). Region I (ORF1a) and region V (ORF N) differ from other IBV strains and probably came from the ancestral South American GI-11 virus (Fig. 1a and e).

Our findings suggest that the mosaic genome of South American strains here analyzed arose by recombination between the GI-11 and GI-16 lineages. Recombination is frequently described in IBV field strains (Dolz et al., 2008, Kottier et al., 1995, Lee and Jackwood, 2000) as a result of the large genome size, a replication machinery that dissociates and reassociates from the template RNA (site-assisted copy choice recombination), and the availability of full-length and subgenomic-length strands for template switching (Fu and Baric, 1992). Most recombination studies in IBV have been conducted using the S1 coding region, but relatively few analysis have been performed in a genomic context (Ammayappan and Vakharia, 2009, Thor et al., 2011, Xu et al., 2016).

The recombination in these South American strains involved four recombinant points that delimit the five genomic regions and were then involved in the generation of the mosaic genomes (Fig. 2). Notable, these breakpoints do not disrupt protein coding regions and two of them are located within genes in inter-ORF regions. Intragenic breakpoints are not very frequent in coronaviruses as recombination generally occurs at intergenic regions by specific template-switching events during the discontinuous transcription of subgenomic mRNAs (Simon-Loriere and Holmes, 2011).

We propose the following evolutionary scenario to explain the pattern of variability in the strains of the GI-16 and GI-11 lineages. This proposal supports the emergence and spreading of the GI-16 lineage in Eurasia before being introduced in South America (Marandino et al., 2015, Yu et al., 2001). This Euro Asiatic spreading of the GI-16 lineage involved sequence divergence (Fig. 1a) (Franzo et al., 2015), and changes in the structure of the genome as evidenced by the differences in the presence of the 6b ORF between the Chinese and Italian strains (Fig. 2).

In South America, the strain of the GI-16 lineage would have experienced extensive recombination with GI-11 ancestral strains to produce two recombinant lineages that have retained the S1 sequences of the parental viruses. Notable, both strains remain extremely similar in the rest of their chimeric genome that is comprised by the ORF1b, S2, 3a, 3b, E, M, 4b, 4c, 5a and 5b of the GI-16 Euro-Asiatic strains, and the 1a, N and 6b of a South American ancestral strain of the GI-11 lineage.

Our findings underscore the fact that different S1 lineages can exist with virtually identical genomes originated by recombination between ancestral strains. It is notable that the S1 protein had persisted as a single unit without the detection of recombination within the S1 coding region in all the strain already analyzed (Marandino et al., 2015). The fact that recombinant breakpoints occur outside the S1 region suggests that the S1 coding region is positively selected and transmitted as an intact unit to avoid disrupting favorable intra-protein interactions (Lefeuvre et al., 2007). S1 is a predominant determinant of antigenicity of IBV and its variability could be related to the evasion of the immune response produced by the vaccine used for IBV control. Most South American countries apply vaccines with the Massachusetts serotype (GI-1 lineage) that is highly divergent from the field strains here analyzed (Fig. 1) (Marandino et al., 2015). In this context, the selection exerted by the continuous use of the Massachusetts vaccine, together with gene flow and local differentiation, have been driven forces for the successful persistence of the GI-11 and GI-16 lineages in South America.

Together, our results reveals that IBV dynamics in South America is singular as it involves two main lineages that have emerged by the transferring of genomic regions between viruses of different origin without disrupting the S1 coding region. These recombinant virus have been extraordinary successful, persisting in the continent for several years with a notorious wide geographic distribution. Our results reveal an interesting recombination pattern that emphasizes the importance of complete genomic characterization to understand viral emergence and evolutionary history.

The following are the supplementary data related to this article.

Phylogenetic tree obtained with the maximum-likelihood method and GTR model with gamma distribution and invariant sites. Phylogenetic reconstruction was carried out using all non-recombinant S1 sequences (n = 147) retrieved from the completely sequenced genome strains and prototype strains (PS) as described by Valastro et al. (2016). Mapping uncertainties for internal nodes are shown as approximate likelihood ratio test values. South American strains are labeled with an asterisk. Complete genome strains grouped in 18 of the 32 recognized lineages.

Supplementary material

Quality control of alignment sequencing data obtained with Qualimap 2 (Okonechnikov et al., 2016).