Molecular epidemiology of infectious bronchitis virus in Poland from 1980 to 2017.

The presence of infectious bronchitis virus (IBV) was identified for the first time in the poultry population in Poland at the end of the 1960s. From this time a few waves of epidemics caused by different IBV variants spread across the country. In order to gain more insight into the molecular epidemiology of IBV in Poland, in the present study the S1 coding region of 34 IBV isolates and nearly whole genome of 10 strains collected over a period of 38 years was characterized. Phylogenetic analysis showed that these strains belonged to five recently established IBV lineages: GI-1, GI-12, GI-13, GI-19 and GI-23. Additionally, two strains from 1989 and 1997 formed a separate branch of the phylogenetic tree categorized as unique early Polish variants, and one strain was revealed to be the recombinant of these and GI-1 lineage viruses. Irrespective of year of isolation and S1-dependent genotype, the genome sequences of Polish IBV strains showed the presence of six genes and 13 ORFs: 5'UTR-1a-1b-S-3a-3b-E-M-4b-4c-5a-5b-N-6b-3'UTR, however their individual genes and putative proteins had different lengths. The phylogenetic analyses performed on the genome of ten Polish IBV strains revealed that they cluster into different groups. The Polish GI-1, GI-19 and GI-23 strains cluster with other similar viruses of these lineages, with the exception of the two strains from 1989 and 1997 which are different. It seems that in Poland in the 1980s and 1990s IBV strains with a unique genome backbone circulated in the field, which were then replaced by other strains belonging to other IBV lineages with a genome backbone specific to these lineages. The recombination analysis showed that some Polish strains resulted from a recombination event involving different IBV lineages, most frequently GI-13 and GI-19.

Introduction

Infectious bronchitis virus (IBV) is the etiological agent of a highly contagious disease of chickens known as infectious bronchitis, but the virus can replicate in epithelial cells of different organs, also affecting the urogenital or digestive tracts beside the respiratory tract (Cavanagh, 2005, Cavanagh, 2007). Together with genetically similar viruses isolated from other domesticated Galliformes, IBV belongs to the Igacovirus subgenus within the Gammacoronavirus genus (Nidovirales order, Cornidovirinae suborder, Coronaviridae family, Orthocoronavirinae subfamily). The non-avian SW1 gammacoronavirus isolated from beluga whales was recently assigned to the separate Cegacovirus subgenus (Dong et al., 2007; King et al., 2018). The virus genome is an approximately 27 kb long single-stranded, positive-sense RNA consisting of several open reading frames (ORFs). Two thirds of the genome in the 5′ end are occupied by two overlapping ORFs encoding viral RNA-dependent RNA polymerase. The 1a and 1b ORFs encode 15 non-structural polypeptides (nsp2–16) which are associated with RNA replication and transcription. In the 3′ end are genes that among other products encode the four major structural proteins: spike (S), envelope (E), matrix (M), and nucleocapsid (N). The S glycoprotein is post-translationally cleaved into S1 and S2 subunits of about 534 and 627 amino acids during viral maturation. The S2 subunit anchors the spike into the virus membrane whereas S1 forms the extracellular part of the spike and plays a major role in tissue tropism and induction of protective immunity (Cavanagh and Gelb, 2008).

IBV undergoes many genetic changes generated both by recombinations and mutations such as substitutions, deletions and insertions, which could lead to the emergence of new variants. Among factors that create favorable conditions for such events are characteristic features of coronaviruses in the genome structure (large single-stranded RNA) and virus biology (minimal proofreading activity of viral polymerase) and modern poultry-rearing habits and immunological pressure caused by the worldwide use of vaccines (Ovchinnikova et al., 2011; Woo et al., 2009). Mutations within the S1 gene particularly result in new geno- or serotypes, and currently there are many such types around the world (de Wit et al., 2011). Their number, diversity and naming and the plurality of methods used for their determination for years have caused much confusion. To avoid it, new classification rules based on the whole S1 gene phylogeny (about 1600 nt) and new nomenclature have been proposed. This system distinguished and named 32 lineages, aggregating into 6 genotypes (GI to GVI) (Valastro et al., 2016). However, in the last three years, two more lineages (GI-28 and 29) and even one more genotype (GVII) have been described in China (Chen et al., 2017; Jiang et al., 2017; Ma et al., 2019).

In Poland, the first suspicion of IB was based on clinical observations, as respiratory symptoms incurable with antibiotics in some flocks and/or misshapen eggs from commercial flocks came to notice. Laboratory confirmation of IBV infection was obtained at the end of the 1960s. Between 1964 and 1965, 4165 sera from two hundred ten chicken flocks at the age of 6–18 months were examined in an agar gel precipitation test and only 13% of them were positive although 79% of flocks contained birds with positive serum (Karczewski and Cakala, 1967). Outbreaks of IB with respiratory signs and a drop in egg production and egg quality in non-vaccinated breeding and particularly laying chicken flocks were recorded in the mid-1980s (Bugajak et al., 1997). Since the mid-1990s, outbreaks of IB-nephritis have been reported in broiler flocks (Minta et al., 2000). A multiplex-PCR testing strains isolated between 1997 and 1998 revealed that the most of them belonged to the 793B type (Capua et al., 1999). The emergence of QX IBV was detected in 2004 (Domanska-Blicharz et al., 2006; Domanska-Blicharz et al., 2007). More recently, the next variant of IBV called Var2 which had been circulating only in the Middle-East region for the previous 20 years was also detected in Poland (Lisowska et al., 2017).

In this study we attempted to molecularly characterize the field IBV strains detected in Poland during the period between 1980 and 2017. Strain determination was accomplished by phylogenetic analysis of the full S1 coding region sequences against reference strains representing all genotypes and lineages recently described (Valastro et al., 2016). Additionally, we also analyzed the complete genome sequences of ten Polish IB viruses. Recombination analysis was also performed using the obtained sequences of these strains.

Materials and methods

Polish IBV strains

Thirty four field IBV strains isolated between 1980 and 2017 in Poland were included in the study. These strains originated from poultry experiencing clinical forms of the disease as respiratory or enteric symptoms, nephritis, or problems with egg production.

Epidemiological information of the studied isolates is summarised in Table 1 . The samples were named to fulfill the previously described criteria, but to make it easier to follow the results of the analysis, in subsequent parts of the text they were shortened to the individual symbol given in the laboratory and the year of identification (Ducatez, 2016).

Table 1

Details of 34 Polish field IBV strains used in this study.

No	Isolate name	Chicken type	Age (days)	Symptoms	Lineage	Genbank Noa
1	gammaCoV/AvCoV/ck/Poland/01/1980	Layer	n/a	ovary dysfunction	GI-1	KT886443
2	gammaCoV/AvCoV/ck/Poland/78/1989	Broiler	21	respiratory	GI-1	MK581200a
3	gammaCoV/AvCoV/ck/Poland/79/1989	Broiler	n/a	respiratory	GI-1	MK581201a
4	gammaCoV/AvCoV/ck/Poland/80/1989	Broiler	49	respiratory	Unique variant	MK581202a
5	gammaCoV/AvCoV/ck/Poland/81/1989	Broiler	63	respiratory	Recombinant	KT886445
6	gammaCoV/AvCoV/ck/Poland/162/1997	Broiler	21	respiratory, enteric, nephritis	Unique variant	MK581203a
7	gammaCoV/AvCoV/ck/Poland/255/1997	Broiler	28	enteric, nephritis	GI-13	MK581204a
8	gammaCoV/AvCoV/ck/Poland/58/1998	Broiler	14	respiratory, enteric, nephritis	GI-13	KT886449
9	gammaCoV/AvCoV/ck/Poland/338/2004	Broiler	36	nephritis	GI-13	KT886450
10	gammaCoV/AvCoV/ck/Poland/548/2004	Broiler	10	enteric/respiratory	GI-19	MK581205a
11	gammaCoV/AvCoV/ck/Poland/29/2005	Broiler	24	nephritis	GI-19	KT886437
12	gammaCoV/AvCoV/ck/Poland/217/2005	Broiler	35	respiratory	GI-13	KT886451
13	gammaCoV/AvCoV/ck/Poland/14/2006	Broiler	28	nephritis	GI-19	KT886438
14	gammaCoV/AvCoV/ck/Poland/1387/2006	Broiler	42	nephritis	GI-19	KT886439
15	gammaCoV/AvCoV/ck/Poland/1612/2006	Broiler	21	nephritis	GI-19	KT886440
16	gammaCoV/AvCoV/ck/Poland/G074/2009	Broiler	n/a	respiratory, nephritis	GI-19	KT886454a
17	gammaCoV/AvCoV/ck/Poland/G018/2010	Layer	53	nephritis	GI-13	KT886452
18	gammaCoV/AvCoV/ck/Poland/G024/2010	Broiler	32	nephritis	GI-19	KT886441
19	gammaCoV/AvCoV/ck/Poland/G033/2011	Broiler	14	nephritis	GI-19	KT886442
20	gammaCoV/AvCoV/ck/Poland/G034/2011	Broiler	n/a	nephritis	GI-1	KT886446
21	gammaCoV/AvCoV/ck/Poland/G049/2012	Broiler	41	n/a	GI-13	MK576135
22	gammaCoV/AvCoV/ck/Poland/G195/2012	Broiler	40	nephritis	GI-19	MK581206a
23	gammaCoV/AvCoV/ck/Poland/G021/2013	Broiler	24	respiratory, poor growth	GI-13	MK576136
24	gammaCoV/AvCoV/ck/Poland/G019/2014	N/a	21	n/a	GI-19	MK576137
25	gammaCoV/AvCoV/ck/Poland/G193/2015	N/a	n/a	poor growth, nephrtitis	GI-13	MK576138
26	gammaCoV/AvCoV/ck/Poland/G029/2016	Broiler	42	respiratory, nephritis	GI-23	MK576140
27	gammaCoV/AvCoV/ck/Poland/G101/2016	Broiler	35	n/a	GI-23	MK576139
28	gammaCoV/AvCoV/ck/Poland/G103/2016	Broiler	39	nephritis	GI-23	MK581207a
29	gammaCoV/AvCoV/ck/Poland/G264/2016	Broiler	10	nephritis	GI-23	MK576141
30	gammaCoV/AvCoV/ck/Poland/G004/2017	Broiler	38	nephritis	GI-23	MK576142
31	gammaCoV/AvCoV/ck/Poland/G225/2017	Broiler	n/a	n/a	GI-23	MK581208a
32	gammaCoV/AvCoV/ck/Poland/G242/2017	N/a	70	respiratory, nephritis	GI-23	MK576143
33	gammaCoV/AvCoV/ck/Poland/G269/2017	Broiler	217	n/a	GI-23	MK576144
34	gammaCoV/AvCoV/ck/Poland/G326/2017	Broiler	42	n/a	GI-12	MK576145

Means the nearly whole genome sequence deposited.

Sample processing

The earliest virus materials from the 1980s were available in the form of a lyophylizate of allantoic fluids from commercial chicken eggs. After propagation in SPF embryos, materials from the 1990s were in the form of allantoic fluids stored deep frozen. Field materials delivered to the Department of Poultry Diseases for diagnostic purposes between 2004 and 2017 were isolated in specific pathogen-free (SPF) chicken eggs as described previously (Gelb and Jackwood, 1998). Virus genome presence confirmation and genotype determination preceded the SPF egg isolation. Materials from 1980 to 1998 as referred to above were also refreshed using the virus isolation method on SPF embryonating eggs. Harvested allantoic fluids were processed using an RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's recommended procedure for RNA extraction, and isolated RNA was stored at −70 °C until analysis.

The S1 coding region and whole genome sequencing

The RT-PCRs were conducted on the one-step model using the One Step RT-PCR kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Various combinations of primer pairs described recently as well as additional primers specifically constructed for some strains (Appendix A Supplementary material) were applied for amplification and sequencing of the whole S1 coding region (Binns et al., 1986; Boursnell et al., 1987; Dolz et al., 2006, Dolz et al., 2008; Lisowska et al., 2017; Worthington et al., 2008). The reactions were run according to the recommended protocol for the kit with different annealing temperatures depending on the melting temperature of the primer pair used. Amplified PCR products were visualized by electrophoresis on a 2% agarose gel stained with ethidium bromide and then purified using a QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). Typically, for the S1 coding region of Polish IBV strains, 3–7 PCR products were sequenced in both directions using Sanger sequencing technology by Genomed (Warsaw, Poland). The complete genomes of these IBV strains were generated using Illumina MiSeq technology (Illumina, San Diego, USA) in several laboratories. The five IBV strains 78/1989, 79/1989, 548/2004, G195/2012 and G103/2016 were processed in the Department of Microbiology of the Swedish National Veterinary Institute (SVA, Uppsala, Sweden), four subsequent virus strains 80/1989, 162/1997, 255/1997 and G225/2017 were analyzed in the Department of Omics Analysis in our institute, and one IBV strain G074/2009 was sequenced by Genomed (Warsaw, Poland). Analyses in these organizations were made according to the standard procedure. Briefly, RNA extracted directly from the allantoid fluid was retrotranscribed into DNA using a Superscript IV First-Strand cDNA Synthesis Kit (Invitrogen, Waltham, USA) and the second strand was synthesized with the addition of Klenow polymerase (New England Biolabs, Ipswich, USA). A 300 bp-long paired-end DNA library was prepared using a Nextera XT sample preparation kit (Illumina, San Diego, USA) and sequencing was performed using a MiSeq Reagent kit v3 (Illumina, San Diego, USA).

Sequence analysis

Sequences of S1 coding region fragments obtained by Sanger sequencing were trimmed based on quality and assembled into consensus sequences using Geneious v11.1.3 (Biomatters, Auckland, New Zealand). Sequences of Polish viruses were searched with BLAST (Basic Local Alignment Search Tool) to find these ones with the highest similarity and include them in the phylogenetic analyses. Then the full S1 sequences were aligned with 199 sequences representing 32 lineages in their 6 IBV genotype groupings and 26 unique variants as Valastro et al. recommended (Valastro et al., 2016) using Clustal W. Additionally, eight sequences representing the two newly identified GI-28 and GI-29 lineages and one GVII-1 genotype were included in the analysis. Sequencing data from MiSeq technology obtained from the SVA (Uppsala, Sweden) and Genomed were processed with the CLC Genomics Workbench (Qiagen, Hilden, Germany). The reads obtained in the Department of Omics Analysis of our Institute were assembled into contigs with the SPAdes assembler using the website at http://spades.bioinf.spbau.ru (Bankevich et al., 2012). For phylogeny of the complete genomes of the Polish IBV strains, a preliminary analysis was carried out using all 362 gammacoronaviruses available in the NIAID Virus Pathogen Database and Analysis Resource (ViPR) through the website at http://www.viprbrc.org/ (Pickett et al., 2012). Next, 56 strains were selected for further analysis, taking into account their clustering in the ViPR analysis. Alignments of nucleotide sequences were performed using the multiple alignment using fast Fourier transform (MAFFT) method in Geneious software, v11.1.3 (Biomatters, Auckland, New Zealand). The alignments were then exported to the MEGA program, v7.0.26 (Tamura et al., 2013). Maximum likelihood (ML) phylogenetic analyses of the S1 coding region and of the complete genome were then conducted using the best-fitting nucleotide substitution models (the lowest Bayesian information criterion (BIC) scores in each analysis were for the general time reversible (GTR) model and a discrete gamma distribution (+G) with five rate categories, assuming that a certain fraction of sites are evolutionarily invariable (G + I)). Bootstrap analyses of the resultant trees were performed using 1000 replicates.

To detect any recombination events in the analyzed sequences, RDP4 software v4.97 was used (Martin et al., 2017). The full S1 coding region sequences of 241 IBV strains were screened to check if unusual clusters formed by Polish IBV strains are viruses representing real new IBV lineage or recombinants. Ten full genomes of Polish IBV strains were also analyzed for recombination events using the complete genomes of 56 representative viruses selected for analysis as described above. The RDP4 analysis was accomplished using different available methods with their default parameters, however recombination events were only considered proven if detected by at least seven programs (RDP, Geneconv, BootScan, Maxchi, Chimaera, SiScan and 3Seq) and the p-value was calculated at below 1.0 × 10E−30.

Results

Accession numbers

Full S1 sequences of the 34 analyzed Polish IBV isolates as well as complete genomes of ten of them were submitted to the GenBank database and accession numbers were assigned as given in Table 1.

Genome organization

The nearly full genome sequences of ten IBV strains were obtained with the 5′ and 3′UTR fragments incomplete. In all genomes the analysis predicted 6 genes consisting of 13 open reading frames (ORFs) with a typical order for IBV of 5′UTR-1a-1b-S-3a-3b-E-M-4b-4c-5a-5b-N-6b-3′UTR, but with their individual genes and putative proteins having different lengths (Table 2 ). The ORFs/proteins with a constant conservative amount of nt and amino acids (aa) were 3a (174 nt/57 aa), 5b (249 nt/82 aa), and N (1230 nt/409 aa). The next ORFs/proteins of conservative length were 4b, 4c and 5a counting 285 nt/94 aa, 171 nt/ 56 aa and 198 nt/65 aa respectively in 8 strains, whereas in two IBVs each of these structures was different (longer or shorter by one nt codon/aa). The difference in ORF/protein length of accessory 3b protein was also slight as it fell within 192–195 nt/63–64 aa. The most diverse in terms of length was the ORF coding E protein, ranging from 282 to 285 nt (93–94 aa) in Polish GI-13 and GI-23 IBV strains to 330–333 nt (109–110 aa) in Polish GI-19 IBVs. A similar relationship was also observed in the case of the ORF encoding M protein, which was the shortest (627–681 nt/223–226 aa) in IBV strains of GI-13 and G-23 lineages and the longest (777 nt/258 aa) in IBVs of GI-19 lineage. The ORF of the S protein was of varying lengths from 3462 nt/1153 aa to 3510 nt/1169 aa and did not show any dependence for length on identified IBV lineages. The molecular relatedness extents of the compete genome and the individual ORFs between Polish and selected IBV strains were 81–100% (Appendix B1-B16 Supplementary material).

Table 2

Regions/genes positions and lengths (nt and aa) of the Polish IBV strains.

IBV strain	UTR/ORF	5′UTRa	1a	1ab	S	3a	3b	E	M	4b	4c	5a	5b	N	6b	3′UTRa
78/1989	Position	1-335	336–12,173	336–20,206	20,157–23,618	23,645–23,818	23,818–24,009	23,990–24,289	24,291–24,968	24,969–25,253	25,174–25,344	25,328–25,525	25,522–25,770	25,713–26,942	26,951–27,175	27,176–27,324
	Length (nt)	335	11,838	19,797	3462	174	192	300	678	285	171	198	249	1230	225	149
	Length (aa)		3945	6597	1153	57	63	99	225	94	56	65	82	409	74
79/1989	Position	1–532	533–12,370	533–20,403	20,354–23,815	23,842–24,015	24,015–24,206	24,187–24,486	24,488–25,165	25,166–25,450	25,371–25,541	25,525–25,722	25,719–25,967	25,910–27,139	27,148–27,372	27,373–27,640
	Length (nt)	532	11,838	19,797	3462	174	192	300	678	285	171	198	249	1230	225	268
	Length (aa)		3945	6597	1153	57	63	99	225	94	56	65	82	409	74
80/1989	Position	1–679	680–12,532	680–20,565	20,516–24,022	24,022–24,195	24,195–24,389	24,373–24,698	24,674–25,345	25,346–25,630	25,551–25,721	25,705–25,902	25,899–26,147	26,090–27,319	27,328–27,555	27,556–27,875
	Length (nt)	679	11,853	19,812	3507	174	195	324	672	285	171	198	249	1230	228	320
	Length (aa)		3950	6602	1168	57	64	107	223	94	56	65	82	409	75
162/1997	Position	1–514	515–12,364	515–20,397	20,348–23,857	23,857–24,030	24,030–24,224	24,205–24,489	24,482–25,162	25,163–25,447	25,368–25,538	25,522–25,719	25,716–25,964	25,907–27,136	27,145–27,366	27,367–27,644
	Length (nt)	514	11,850	19,809	3510	174	195	285	681	285	171	198	249	1230	222	285
	Length (aa)		3949	6601	1169	57	64	94	226	94	56	65	82	409	73
255/1997	Position	1–514	515–12,358	515–20,391	20,342–23,833	23,833–24,006	24,006–24,200	24,181–24,462	24,459–25,130	25,131–25,415	25,336–25,506	25,487–25,687	25,684–25,932	25,875–27,104	27,113–27,334	27,335–27,651
	Length (nt)	514	11,844	19,803	3492	174	195	282	672	285	171	201	249	1230	222	317
	Length (aa)		3947	6599	1163	57	64	93	223	94	56	66	82	409	73
548/2004	Position	1–716	717–12,569	717–20,602	20,533–24,050	24,050–24,223	24,223–24,417	24,401–24,730	24,603–25,379	25,380–25,667	25,585–25,749	25,730–25,930	25,927–26,175	26,118–27,347	27,356–27,577	27,578–27,844
	Length (nt)	716	11,853	19,812	3498	174	195	330	777	288	165	201	249	1230	222	267
	Length (aa)		3950	6602	1165	57	64	109	258	95	54	66	82	409	73
G074/2009	Position	1–516	517–12,342	517–20,375	20,326–23,823	23,823–23,996	23,996–24,187	24,168–24,500	24,373–25,149	25,150–25,434	25,355–25,525	25,509–25,706	25,703–25,951	25,894–27,123	27,132–27,373	27,374–27,418
	Length (nt)	516	11,826	19,785	3498	174	192	333	777	285	171	198	249	1230	222	45
	Length (aa)		3941	6593	1165	57	63	110	258	94	56	65	82	409	73
G195/2012	Position	1–516	517–12,369	517–20,402	20,353–23,835	23,835–24,008	24,008–24,202	24,183–24,515	24,388–25,164	25,165–25,449	25,370–25,540	25,524–25,721	25,718–25,966	25,909–27,138	27,147–27,368	27,369–27,433
	Length (nt)	516	11,853	19,812	3483	174	195	333	777	285	171	198	249	1230	222	65
	Length (aa)		3950	6602	1160	57	64	110	258	94	56	65	82	409	73
G101/2016	Position	1–606	607–12,459	607–20,507	20,443–23,934	23,934–24,107	24,107–24,298	24,279–24,563	24,566–25,236	25,237–25,521	25,442–25,612	25,596–25,793	25,790–26,038	25,981–27,210	27,219–27,443	27,444–27,791
	Length (nt)	606	11,853	19,827	3492	174	192	285	681	285	171	198	249	1230	225	348
	Length (aa)		3950	6607	1163	57	63	94	226	94	56	65	82	409	74
G225/2017	Position	1–497	498–12,347	498–20,395	20,331–23,822	23,822–23,995	23,995–24,186	24,167–24,451	24,444–25,124	25,125–25,406	25,330–25,497	25,481–25,678	25,675–25,923	25,866–27,095	27,104–27,325	27,326–27,613
	Length (nt)	497	11,850	19,824	3492	174	192	285	681	282	168	198	249	1230	222	288
	Length (aa)		3949	6606	1163	57	63	94	226	93	55	65	82	409	73

As one-sided PCR (5′ and 3′ RACE-PCRs) techniques were not applied, UTR fragments should be considered as incomplete.

Phylogenetic analysis

Genotyping based on phylogenetic analysis of full S1 coding region sequences of 34 Polish IBV strains from the years 1980–2017 grouped them into six groups: five distinct, previously known lineages and an additional new one (Fig. 1 ). Five isolates comprised of four early ones from the 1980s and one identified in 2011 were assigned to the GI-1 lineage. One strain identified in 2017 affiliated to the GI-12 lineage. The group of GI-13 lineage contained eight IBV strains: two from the late 1990s isolated between 1997 and 1998, two identified between 2004 and 2005, and four strains isolated after 2010. The GI-19 lineage comprised ten strains detected between 2004 and 2014 and the group of GI-23 lineage held eight isolates detected between 2015 and 2017. Two isolates, 162/1997 and 80/1989, were in the separate cluster designated early Polish on the phylogenetic tree.

Fig. 1

Phylogenetic tree of the S1 gene of 207 reference and 34 Polish IBV strains (bold underlined letters). The tree was constructed using MEGA 7 using the maximum likelihood method based on the GTR + G + I model and 1000 bootstrap replicates (bootstrap values shown on the tree). To make the tree clearer visually, branches with IBV lineages only distantly correlated with studied Polish strains are collapsed.

Sequence analysis revealed that five Polish GI-1 strains shared nucleotide identities of 88.4–100% and formed two clades. In one clade were the Dutch H120, North American, South African, Indian and four Polish isolates but the strain 81/1989 formed a distinct branch in this GI-1 subtree with low nucleotide identity of 86.2–88.1% to the rest of this group. The analysis of the S1 coding region sequences of eight Polish isolates of GI-13 lineage showed that four strains formed a common branch, sharing 98.8–99.4% nt identity with the 7/91 strain from the United Kingdom, whereas two earlier strains from 2004 and 2005 which were similar to each other with 99.3% identity form a sister group with the more recent viruses and had 97.9–98.2% nt similarity with the 7/91 strain. The earliest Polish GI-13 strains from 1997 and 1998 were visibly different and had nt identity of 93.9–96.1% to the rest of the G-13 IBV strains. One of them, strain 58/1998, occupied positions close to the Israeli Variant 1 strain from 1996 with identity of 96.1%, and its similarity to the Moroccan GI-13 lineage prototype G strain from 1983 was 93.9%. The similarity of the other, strain 255/1997, was 92.3% and 94% to strains from Israel and Morocco, respectively. The GI-12 lineage contained only one Polish strain, G326/2017, which shared 98.7% nt identity with the Dutch D274 virus. Eight Polish strains were in the GI-23 lineage and the similarity of their S1 coding region sequences was between 97.2 and 99.5%. Their identity with the pathogenic Israeli IS/1494/06 strain from 2006 showed as from 98.5 to 99.5%. The subtree of GI-19 lineage contained Polish QX strains in two branches, of which one contained nine strains with nt identity of 98.2–99.6% to the European QX prototype Dutch L-1449 K/04 IBV. One strain, G019/2014, with similarity to the previous ones of 96.3–97.5%, constituted the offshoot branch. The two Polish strains 80/1989 and 162/1997, isolated at an interval of 8 years from each other, formed a separate branch in the phylogenetic tree and they shared 85.1% nt identity.

The phylogenetic analysis of the 66 analyzed full IBV genomes showed that ten Polish IBV strains grouped into four phylogenetic groups (Fig. 2 ). Two early strains, 78/1989 and 79/1989, clustered together with Massachusetts-like strains (Mass 41, Peafowl/GD/KQ6/2003 and SES 15AB-01) showing the highest nt identity of 90.5% with the sequence of the prototype GI-1 lineage Beaudette strain. The two most recent strains, G103/2016 and G225/2017, clustered with IBVs of GI-23 lineage. The sequence identities of G103/2016 and G225/2017 to the previously described Polish gammaCoV/Ck/Poland/G052/2016 strain were 99.8 and 95.5%, respectively. The three Polish strains 548/2004, G074/2009 and G195/2012 were in the same cluster as other QX strains from Europe and Africa and were distantly related to Chinese QX IBVs analyzed in this study (SDZB0808 and P100). They had nucleotide similarity to each other of 92.5–93.3% and were located in two subclades. A Polish strain from 2009 clustered together with SWE/0658946/10, the first described full-genome IBV strain of QX type in Europe, with similarity of 93.7%. The other two Polish QX strains were 93.3% similar to each other and formed a common branch on the phylogenetic tree. Three early Polish strains 80/1989, 162/1997 and 255/1997 were in a separate branch on the phylogenetic tree and showed nucleotide similarity with each other in the range of 90.5–92.5%.

Fig. 2

Phylogenetic tree of the complete genomes of 56 representative gammacoronaviruses and 10 Polish IBV strains marked with bold underlined letters. The tree was constructed using MEGA 7 using the maximum likelihood method based on the GTR + G + I model and 1000 bootstrap replicates (bootstrap values shown on the tree).

Recombination analysis

Recombination analysis of all 241 aligned full S1 sequences was performed to assess the existence of possible recombinants among the analyzed Polish IBV strains, especially those with less obvious membership to the lineage, i.e. 80/1989, 81/1998 and 162/1997. Our analysis identified only one S1 coding region which resulted from recombination events and it belongs to the 81/1989 strain; this event having taken place was supported by seven different methods with a very good global KA p-value of 1.277E−34. We confirmed this recombination breakpoint with phylogenetic trees. The region from 1 to 651 nt of the 81/1989 IBV strain clustered together with 80/1989 and 162/1997 isolates, the viruses which formed the separate early Polish cluster on the full S1 coding region phylogenetic tree (Fig. 3a). In turn, the region from 652 to 1585 nt clustered with viruses belonging to GI-1 lineages (strains 78/1989, 79/1989, IBV438 India 2012 and H120) (Fig. 3b). The relevant S1 coding region fragments of the other IBV strains analyzed in this study grouped in the same way as they did in the phylogenetic analysis of the entire, intact S1 coding region.

Fig. 3

Phylogenetic tree of the S1 gene fragment between potential recombination breakpoints 1 and 651 (a) and 652 and 1585 (b) among 241 IBVs included in the analysis. Sequences of Polish IBV strains are marked with bold underlined letters and recombinants with black dots. The tree was constructed using MEGA 7 using the maximum likelihood method based on the GTR + G + I model and 1000 bootstrap replicates (bootstrap values shown on the tree). To make the tree clearer visually, branches with IBV lineages only distantly correlated with studied Polish strains are collapsed.

To check if any of the analyzed genomes of Polish IBV strains result from recombination events, their sequences were thoroughly examined using the RDP4 program. Our analysis revealed many such events. However, we selected five of them identified in six strains and they were supported with seven different methods (RDP, Geneconv, BootScan, Maxchi, Chimaera, SiScan and 3Seq) and a very good global KA p-value below 1.0 x 10E−30 (Table 3 ). They were 255/1997, 548/2004, G074/2009, G195/2012 and two GI-23 IBV strains G103/2016 and G/225/2017. Most of the putative recombinant regions were relatively long (from 2193 to 3320 nt), however one was exceptionally long at 8207 nt and one atypically short at 561 nt, and they were contained within Genes 1 and 2 (pol 1a, 1b and S). Three of the Polish strains (255/1997, G074/2009 and G195/2012) have recombined with a 4/91 vaccine strain (Genes 1 and 2) and three further strains (548/2004, G103/2016 and G225/2017) with European and Chinese QX strains (Gene 1).

Table 3

Recombination breakpoints, fragments lengths, genes, potential parents in the genome of Polish IBV strains.

Strain	Start	End	Fragment length	Gene	Potential parents	Global KA p-value (1.0 × 10E-Na)
548/2004	11,393	19,599	8207	1b	SDZB0808	5.628E-221
					Ck/CH/LSD08-8
G074/2009	4768	7888	3121	1a	4/91	8.694E-129
					Ck/CH/LSD/110851
255/1997	20,411	23,730	3320	S	4/91	1.769E-85
					Ck/CH/LSD/110851
G195/2012	14,591	16,783	2193	1b	4/91	1.166E-44
					Ck/CH/LSD/110851
G103/2016	19,400	19,960	561	1b	ITA/90254/2005	1.241E-36
G225/2017					AR251-15

Exponent in scientific notation.

Discussion

The retrospective phylogenetic analysis of IBV included field strains collected over a period of 38 years, between 1980 and 2017. We investigated the S1 coding region of 34 IBV strains and the whole genome of ten strains. Polish IBV strains showed different molecular features of the S1 coding region allowing their genotype or lineage to be determined, and their appearance in time reflects the history of IBV epidemics in Europe (de Wit et al., 2011). The plot showing the timeline of various IBV lineages's detection and introduction of different vaccines to poultry population in Poland is given in Fig. 4 . The first identified IBV isolates in Europe belonged to the Mass type. In the Netherlands they were diagnosed in the middle of the 1950s and one of them was even attenuated for vaccine development purposes (Bijlenga et al., 2004). The first IBV material available in our laboratory originates from 1980 and its S1 sequence displayed the features of GI-1 lineage, although no information was provided about the disease symptoms observed in the chicken flock where it was identified. Later on, especially in the middle of the 1980s, numerous cases of a drop in egg production were recorded. The problem was so serious that Polish veterinary authorities decided to allow the first IB vaccine introduction, but only for immunization of commercial layer flocks. The health problems in layers were significantly mitigated, but at the end of the 1980s respiratory problems and mortality manifested in broiler chicken flocks, and so the vaccination of chickens of this production type was also started (Minta et al., 1990). We thoroughly examined four virus strains from that time, 78/1989, 79/1989, 80/1989 and 81/1989. Two of them, strains 78/1989 and 79/1989, have the S1 structure typical of the GI-1 lineage. Similarly, their entire genomes revealed the highest identity to Mass-like strains such as H120 or Beaudette. These two viruses came from broiler chickens on farms in the Silesia region separated by only a few kilometers. The two other viruses, 80/1989 and 81/1989, were identified in broilers delivered to the laboratory near the same period (in June 1989) but from farms about 30 km from the previous ones. However, our investigation revealed a distinction between them. Phylogenetic analysis of the full S1 coding region showed that strain 80/1989 forms a separate branch on the tree designated as early Polish and was classified as a unique variant of IBV within the GI genotype. Deep analysis of the 81/1989 strain strongly suggests that its S1 coding region was created as the result of a recombination event between the Mass-like strains and the unique early Polish variants circulating in the field at that time. It should be emphasized that the identified recombination breakpoint (651 nt) was in the intermediate region between highly variable regions (HVRs) 1 and 2 and HVR3 previously described as the most frequent locations of variations between IB viruses, and moreover, it exactly matches the breakpoint (550 and 652 nt) of recombinants between viruses of GI-19 and GI-22 lineages (Valastro et al., 2016).

Fig. 4

The timeline of various IBV lineages's detection and introduction of different vaccines to poultry population in Poland.

The introduction of vaccines based on the Mass-like strains for chicken immunization significantly reduced the economic losses caused by IB. This state of IB control lasted for about 10 years until 1997, when IB disease inducing kidney damage appeared, caused by 793B-like IBV strains. The first case of nephritis was in 3-week-old broilers in the south of Poland. The birds showed signs of severe enteritis and the observed gross lesions were congested tracheas and lungs and swollen and pale kidneys with the presence of urine. In subsequent months, further broilers with nephritis were provided for diagnostic purposes and the diseased flocks from which they came were located in all regions of Poland; most of them had not been vaccinated against IB, but some had been immunized with Mass-like vaccine in the first days of life. The strains identified at that time, 255/1997 and 58/1998, have an S1 sequence similar to IBV strains of GI-13 lineage. Surprisingly, one of the first isolates known to cause nephritis, strain 162/1997 together with strain 80/1989 inflicting respiratory disorders, were located in the early Polish branch of the phylogenetic tree with well supported uniformity to others (bootstrap value of 74) (Hillis and Bull, 1993). The next Polish GI-13 strains identified between 2004 and 2015 revealed the highest nt similarity to the S1 sequence of the 4/91 strain contained in the most commonly used vaccine in Poland at that time (Adzhar et al., 1997).

The next epidemic wave of IB in Poland was caused by QX strains. The first report of disease induced by this virus type was published in 2006, but our studies showed the presence of the virus in Poland in 2004 (Domanska-Blicharz et al., 2006). It was from this year that QX strains were first identified in Holland, Germany, Belgium and France, and next year they became dominant in some of these countries (Worthington et al., 2008). The analysis of the migration history of GI-19 strains suggests that most European ones came from a single introduction from China, which then spread in European countries, evolving in them separately, since they tend to cluster by country. However, the genetic variability of GI-19 IBVs sometimes identified within countries suggests subsequent introduction of the virus in epidemic waves (Franzo et al., 2017). The division of Polish GI-19 strains into four clusters could reflect a separate introduction or epidemic wave of this virus variant into the country.

The last large epidemic wave of IB was caused by GI-23 (Var2) strains. The first strain of this lineage was identified in December 2015 and in the following months it was the most common virus type detected in field samples delivered to our laboratory for diagnostic purposes apart from strains of 793B (Lisowska et al., 2017). The S1 coding region of most Polish GI-23 strains is in the same phylogenetic cluster, except for strain G269/2017, which constitutes a separate one and could result from a separate virus introduction or from its intensive evolution.

The single Polish G326/2017 virus strain of GI-12 lineage analyzed in our study was identified in a 6-week-old broiler flock vaccinated with Poulvac IB Primer so it is highly probable that the identified strain originated from vaccine virus. Although Ball et al. (Ball et al., 2017) showed that after vaccination of 1-day-old broilers with this vaccine only RNA of the Mass strain was detected in tissues and swabs and explained it through the higher replication potential of the Mass virus. It cannot be ruled out that, as some other IBV strains are, the D274 virus is deposited in the body of chickens (possibly in the cecal tonsils) and after some time it is shed with cloaca (Alexander and Gough, 1978; Naqi et al., 2003).

It should be noted that the IBV strains discussed here in detail are those that caused the greatest losses in Polish poultry farming. During this period, strains of other genotypes and lineages also circulated in the field but their detection or type determination was not possible using available methods. The comprehensive studies of Polish IBV isolates from 1998 to 1999 using serological and molecular tests conducted in cooperation with Italian researchers showed that the 793B type was a major component of the IBV population in Poland during this period, but serologically the presence of 624/I isolates was also identified and one isolate even showed no serological cross-reaction in an HI test nor amplification in RT-PCR (Capua et al., 1999). Recently, the 624/I and Q1 types were determined to affiliate to the GI-16 lineage, which has been present in Europe (Italy) since 1963 and persists until now (Franzo et al., 2018). In the period 2011–2013 numerous cases of D1466 IBV (GII-1 lineage) were detected, however, in subsequent years (2014–2015), the number of D1466-positive samples dropped to 3.1%, and currently we do not detect these viruses at all (Domanska-Blicharz et al., 2012; Domanska-Blicharz et al., 2017).

The complete genome sequences of all Polish field IBV strains showed the presence of six genes and 13 ORFs in the order previously reported, irrespective of the year of first isolation (Abolnik, 2015; Gomaa et al., 2008; Hewson et al., 2011). Most accessory proteins are conservative in their lengths, in contrast to the structural ones which differ by even as much as 35 aa (M protein of GI-13/23 and GI-19). An interesting observation is clustering of Polish IBVs based on the complete genome sequences. The earliest strains 78/1989 and 79/1989 belonging to the GI-1 lineage cluster with other representatives of this lineage such as the Beaudette and Dutch H120 strains, which could indicate the common origin of Mass-like viruses. The viruses from the next epidemic wave are 162/1997 and 255/1997, and they were on the separate branch of the phylogenetic tree together with the strain 80/1989. The 80/1989 virus was located on the S1 coding region phylogenetic tree with the 162/1997 strain in the separate IBV branch of early Polish IBV. On the other hand, the third virus of this separate group, 255/97, found its place in the GI-13 lineage on the phylogenic tree based on the S1 coding region. Thorough analysis using RDP4 software revealed that the S gene of this virus was acquired from IBV strains of GI-13 lineage during a recombination event. The results suggest that in the late 1980s and 1990s two IBV variants circulated in the Polish poultry population: GI-1 and unique early Polish ones. These viruses differ not only in the S1 coding region, which is the basis for the differentiation of lineages, but also in the remaining part of the genome. The viruses with such a genome backbone recombined with other viruses that were donors of the S1 coding region. Grouping on a phylogenetic tree based on the complete genome of the other five Polish strains was as expected. Three G-19 strains, 548/2004, G074/2009 and G195/2012, took positions among other QX-like viruses from Europe (Sweden and Italy) and Africa (South Africa and Sudan) (Abolnik, 2015; Abro et al., 2012; Ducatez et al., 2009; Naguib et al., 2016). In turn, two strains of GI-23 lineage, the viruses G103/2016 and G225/2017, were in the branch with the previously characterized Polish G052/2016 and Iranian IS-1494 strains of GI-23 lineage isolated in 2015. It seems that in Poland in the 1980s and 1990s IBV strains with a unique genome backbone circulated in the field, which were then replaced by strains belonging to other IBV lineages with a genome backbone specific to these lineages. In addition to the aforementioned recombination, five such events were also identified in Polish IBV strains. Three strains of the GI-19 lineage had ORF1a and ORF1b which revealed a high frequency of recombination events with 4/91 and SDZB0808-like strains (QX type strains from China from 2008). Two strains of the GI-23 lineage exhibited recombination with the Italy/90254/2005-type IBV, and a similar recombination pattern was also previously indicated (Lisowska et al., 2017).

In conclusion, phylogenetic analysis performed on the S1 coding region of Polish IBV strains collected during a 38-year period (1980–2017) showed that these strains belonged to five recently established IBV lineages: GI-1, GI-12, GI-13, GI-19 and GI-23. Additionally, two strains formed a separate branch of the phylogenetic tree described as unique early Polish variants and one strain revealed itself to be the recombinant of GI-1 lineage viruses and these unique early Polish variants. The phylogenetic analyses performed on the complete genome of ten Polish IBV strains showed that they cluster into different groups. Polish GI-1, GI-19 and GI-23 strains cluster with other similar viruses of these lineages, with the exception of the strains from 1989 to 1997 which are different. The recombination analysis showed that Polish strains are a mosaic of different parental viruses most likely resulting from recombination events involving different IBV lineages, most frequently GI-13 and GI-19. It should be also stressed that the major epidemics of IB in Poland appeared every 8–11 years: GI-1 in 1988, GI-13 in 1997, GI-19 in 2004 and GI-23 in 2015. These subsequent IBV lineages could have reached chickens in Poland in various ways: carried by wild birds, or as a result of international trade, including uncontrolled movement of animals across borders (Domanska-Blicharz et al., 2014; Hussein et al., 2014; Kahya et al., 2013). Despite the apparent regularity in the appearance of subsequent IB epidemics, it is absolutely impossible to predict when the next one will appear. The most important impediment to prediction are the visible climate changes forcing changes in bird behavior, but another is the extraordinary intensification of the poultry industry in Poland.

Taken as a whole, the molecular characteristics of Polish IBVs presented here could help to understand the origin, spread and evolution of IB viruses in Europe and the rest of the world.

Declaration of Competing Interest

None.