Introduction {#Sec1}
============

Infectious bronchitis virus (IBV), which is a gamma corona virus classified in the family *Coronaviridae*, causes respiratory signs, egg-drops, nephritis, and proventriculitis in domestic fowl, resulting in serious economic losses \[[@CR7]\]. The single-stranded positive-sense RNA genome is approximately 27 kb, with eight subgenomic mRNAs (sg-mRNAs) that are transcribed under the regulation of transcription regulatory sequences (TRSs). The leader TRS is located at the 5'-end of the viral genomic RNA and 8 body TRSs are located upstream of the coding regions of structural and nonstructural genes. The homology between the leader and body TRSs may affect the copy number of subgenomic mRNA transcripts and the level of gene expression \[[@CR4], [@CR44]\]. The viral genome encodes RNA polymerase/transcriptase (1ab), spike (S), envelope (E), membrane (M), nonstructural proteins (nsps), 3a, 3b, 4b, 4c, 5a, 5b and 6b, and nucleocapsid (N) \[[@CR7], [@CR35]\]. Large 1a and 1ab proteins are cleaved by viral proteases into approximately 15 nsps, which are involved in virus replication and pathogenicity \[[@CR2], [@CR13]\]. Spike is cleaved into S1 and S2 by proteases, and most epitopes for virus neutralization are within the hypervariable regions of S1 and S2 \[[@CR19], [@CR20]\]. Six epitopes, D (24-60), E (132-149), C/A/B (291-398) and F (497-543), in S1 and 1 epitope, G (548-574), in S2 have already been defined \[[@CR19]\]. Different IBV serotypes do not confer protective effects against each other. E and M proteins are important for virion assembly and are membrane-anchored glycoproteins \[[@CR9]\]. The short ectodomain of the M protein protrudes out from the viral envelope and contains a protective epitope \[[@CR14]\]. Viral variation and diversity are attributed to the long RNA genome and frequent nucleotide substitutions, deletions, insertions, and recombination events occur due to the poor proofreading capability of RNA polymerases \[[@CR16], [@CR27], [@CR31], [@CR43]\].

IBVs are distributed worldwide in poultry, and most countries have their own indigenous variants and common genotypes \[[@CR38]\]. According to a phylogenetic analysis based on the short variable region of S1, several indigenous and common IBV genotypes (e.g., K-I, KM91-like, New cluster 1, QX-like) have co-circulated, and in the past recombination events are suspected to have occurred in Korea \[[@CR16], [@CR24], [@CR25]\]. Nephropathogenic KM91-like viruses were first isolated in 1991 and became prevalent in the field. To prevent the economic losses caused by KM91-like virus infection, oil emulsion-inactivated and attenuated live vaccines have been used in Korea \[[@CR25]\]. QX-like IBVs are rampant throughout Asia and Europe. They were initially isolated in 1997 in China and might have been introduced to Korea in 2002--2003 \[[@CR1], [@CR8], [@CR18], [@CR24], [@CR25]\]. Subsequently, novel recombinant viruses (New Cluster 1, NC1), which were suspected to have originated through recombination between KM91-like and QX-like viruses, emerged in Korea in 2005 and became prevalent between 2007 and 2010 \[[@CR27]\]. Recent recurrences of QX-like viruses were observed in 2008 and 2011, and they have become prevalent in Korea \[[@CR31]\]. Phylogenetic and computational recombination analyses using genome sequence data have shed light on the occurrence of complex recombination events in the genomes of IBVs and the precise donor and acceptor relationships between concurrent field viruses during recombination \[[@CR6], [@CR11], [@CR15], [@CR16], [@CR35], [@CR42]\]. However, the minimum essential core replication and pathogenicity-related 1ab gene is not typically included in phylogenetic analyses. In the present study, we analyzed coding and non-coding IBV sequences and conducted computational recombination analyses as well as phylogenetic analyses with S1, 1ab, and full genome sequences of two early, three recombinant and two recurrent QX-like viruses to understand the molecular epidemiology and molecular evolution of QX-like viruses at the genome level. Consequently, we clarified the molecular epidemiology of infectious bronchitis caused by QX-like viruses in Korea. Furthermore, we classified IBVs into three distinct lineages, including North America-Europe, China/Taiwan, and China, by using a 1ab gene-based phylogenetic analysis. We found that the early Korean and European QX-like viruses were generated by recombination between the North America-Europe and China lineages.

Materials and methods {#Sec2}
=====================

Virus isolation and propagation {#Sec3}
-------------------------------

QIA-03342, QIA-KR/D79/05, and QIA-Q43/06 were isolated in 2003, 2005, and 2006, respectively, and form a cluster with the QX strain according to a partial S1 gene-based phylogenetic tree \[[@CR24]\]. SNU9106 and SNU10043 as well as SNU8065 and SNU11045 have been classified as NC1 and recurrent QX-like viruses, respectively \[[@CR31]\]. All isolates were inoculated into 10-day-old SPF embryonated chicken eggs (ECE; Valo BioMedia, Adel, IA, USA) via the allantoic cavity route and incubated for 48 h. Next, the eggs were chilled at 4°C overnight. The allantoic fluid was harvested, and the supernatant was stored at -70°C after centrifugation at 3,000 rpm for 10 min.

RNA extraction and RT-PCR {#Sec4}
-------------------------

Viral genomic RNA was extracted from infectious allantoic fluid using the Viral Gene-Spin Kit (iNtRON Biotechnology, Seongnam, Korea), and RT-PCR was performed using a One-step RT-PCR Kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's instructions. The cDNA was synthesized at 50°C for 30 min and then heated to 95°C for 15 min to inactivate reverse transcriptase. The PCR conditions were as follows: 40 repetitions of denaturation at 94°C for 30 s, annealing at 50°C for 30 s, extension at 72°C for 2 min, and a final extension at 72°C for 5 min. The primer sets for the amplification of overlapping fragments of the IBV genomes and for sequencing are listed in Supplementary data 1.

Sequencing and genome sequence analysis {#Sec5}
---------------------------------------

The amplicons were purified using the MEGA-quick-spinTM Total Fragment DNA Purification Kit (iNtRON Biotechnology) and sequenced with PCR primers using an ABI3711 automatic sequencer (Macrogen Co., Seoul, Korea). The overlapping gene fragments were assembled to obtain a single complete genome sequence using ChromasPro version 1.5 (Technelysium Pty Ltd., Brisbane, Australia). Nucleotide and amino acid identity estimates and amino acid translations were obtained using BioEdit (ver. 5.0.9.1.). All sequences were submitted to the GenBank database (accession numbers, KU900738--KU900744).

Phylogenetic analyses based on S1, 1ab, and the complete genomes sequences were conducted using MEGA (ver. 6.0.6) and the neighbor-joining method with Tamura-Nei distances and 1,000 bootstrap replicates. The GenBank accession numbers of IBV isolates included in the analysis are presented in Fig. [2](#Fig2){ref-type="fig"}.

Computational recombination analysis {#Sec6}
------------------------------------

Putative recombination events between the QIA-03342, QIA-KR/D79/05, QIA-Q43, SNU-9106, SNU-10043, KM91, and YX10 strains \[[@CR40]\] were analyzed with RDP4 (ver. 4.14; neighbor-joining, Kimura distance correction; P \< 0.05) \[[@CR30]\]. Seven methods in RDP 4.14, including RDP, GENECONV, Bootscan, MaxChi, Chimaera, SiScan, and 3Seq were used to detect recombination events. We inputted genomic nucleotide sequences of each Korean QX-like virus with the nucleotides of the KM91 and YX10 strains. The recombination detection methods sequentially tested every combination of the three sequences for evidence that one of the three sequences was a recombinant and the other two were its parents.

Results {#Sec7}
=======

Characterization of the genome sequences {#Sec8}
----------------------------------------

The genome sizes of QIA03342, QIA-KR/D79/05, QIA-Q43, SNU9106, SNU10043, SNU8065, and SNU11045 were 27684, 27682, 27675, 27684, 27683, 27679, and 27679, respectively, and the GC contents ranged from 38.08--38.23% (Table [1](#Tab1){ref-type="table"}). Most of the open reading frames overlapped, and the ribosomal frameshifting/slippery sequences were conserved in all of the QX-like viruses. The YX10 strain (CK/CH/Zhejiang/06/10), a recombinant nephropathogenic QX-like virus isolated in the Zhejiang province of China in 2010, was selected due to its high nucleotide identity with the recurrent QX-like viruses \[[@CR40]\]. The whole genomes of QIA-03342, QIA-D79/05, QIA-Q43/06, SNU9106, and SNU10043 had higher nucleotide identities for KM91 (95--96%) than for YX10 (88%) but had identities to SNU8065 and SNU11045 of 87% and 98%, respectively (Table [2](#Tab2){ref-type="table"}). The 5′- and 3′-UTRs of QIA-03342, QIA-D79/05, QIA-Q43/06, SNU9106, and SNU10043 were more similar to KM91 (98--99% and 95--98%, respectively) than to YX10 (94--95% and 80--83%, respectively), and SNU8065 and SNU11045 were more similar to YX10 (99% and 95--96%, respectively) than to KM91 (94--95% and 84--85%, respectively) (Table [2](#Tab2){ref-type="table"}).Table 1Genome size, GC content, and IBV proteins with differing amino acid lengths for the examined QX-like virusesQIA-03342QIA-KR/D79/05QIA-Q43SNU9106SNU-10043SNU-8065SNU11045Full length27,68427,68227,67527,68427,68327,67927,679G + C content38.08%38.2%38.08%38.23%38.13%38.16%38.14%S1,1651,1651,1631,1651,1651,1651,1653b65656565656262E109109109109109108108M226226226226226225225 Table 2Nucleotide and amino acid identities (%) of QX-like viruses, when compared to KM91 and YX10Nucleotide and amino acid identity (nt/aa; %)QIA-\
03342QIA-D79/05QIA-Q43/06SNU-\
9106SNU-10043SNU-\
8065SNU-11045GenomeKM9195959696958787YX10888888888898985′-UTRKM9199989999999594YX10959495949599991abKM9197/9796/9796/9796/9796/9787/9187/91YX1087/9187/9187/9187/9187/9198/9998/99S1KM9185/8685/8590/9092/9091/9185/8585/85YX1096/9595/9491/8988/8888/8998/9898/98S2KM9192/9494/9498/9798/9793/9393/9593/95YX1080/7180/7180/7179/7178/7099/10099/1003aKM9163/5263/5263/5263/5262/5072/6972/69YX1074/6374/6374/6375/6576/66100/100100/1003bKM9180/7280/7280/7280/7281/7474/6774/67YX1086/8587/8586/8585/8386/8499/9999/99EKM9198/10097/9899/10097/9694/9586/8586/85YX1086/8587/8586/8585/8386/8499/9999/99MKM9194/9696/9894/9595/9698/10089/9389/93YX1090/9491/9490/9389/9289/9499/9999/994bKM9185/8293/9485/8294/9385/8189/8789/87YX1085/8490/9185/8489/8785/83100/100100/1004cKM9188/7787/74887787/7288/7771/537153YX1075/5874/5374/5874/5875/58100/100100/1005aKM9196/9797/9895/9595/9294/9483/8383/83YX1083/8583/8581/8182/8180/8098/9898/985bKM9197/9399/9898/9598/9499/9891/8891/88YX1093/9192/9092/9092/9093/9193/8892/88NKM9197/9797/9797/9697/9797/9789/9289/92YX1088/9288/9288/9288/9288/9288/9288/926bKM9197/9696/9394/9297/95100/10084/7885/77YX1084/8183/7881/7784/8084/7896/9595/933′-UTRKM9198979695988485YX1083828080829596

Comparison of open reading frames (orfs) {#Sec9}
----------------------------------------

Except for the S, 3b, E, and M proteins, the amino acid lengths of the proteins did not differ among strains (Table [1](#Tab1){ref-type="table"}). The nucleotide and amino acid sequences of each virus were compared to those of the KM91 and YX10 strains (Table [2](#Tab2){ref-type="table"}). The 1ab, S2, E, M, 5a, 5b, N, and 6b genes of QIA-03342, QIA-D79/05, QIA-Q43/06, SNU9106, and SNU10043 were more similar to KM91 (nucleotide (nt) 92--99%, amino acid (aa) 92--100%) than to YX10 (nt 78--93%, aa 70--94%). S1 of QIA-03342 and QIA-D79/05 had higher identities to YX10 (nt 96% and 95%, aa 95% and 94%, respectively) than to KM91 (nt 85% and 85%, aa 86% and 85%, respectively), but QIA-Q43/06, SNU9106, and SNU10043 with respect to KM91 and YX10 were more similar (nt 90--92%, aa 90--91% and nt 88--91, aa 88--89%, respectively). The 3a, 3b, 4b, and 4c genes of QIA-03342, QIA-D79/05, QIA-Q43/06, SNU9106, and SNU10043 showed relatively low nucleotide and amino acid identities when compared to KM91 (nt 62--94%, aa 50--94%) and YX10 (nt 74--90%, aa 53--91%), but they showed relatively high identities when compared to other viruses, e.g., 3a (Connecticut/KF696629, M41/FJ904713, TCoV/EU022526, CK/CH/LDL/110931/KJ425485, 2994/02/GU386375; 94--95%), 3b (Delaware/GU393332, TCoV, and Gray/GU393334; 93--97%), 4b (CK/CH/LTJ/95I/EF602448; 98--100%), and 4c (CK/CH/LTJ/95I, 93--97%). The 1ab, S1, S2, 3a, 3b, E, M, and 5a genes of SNU8065 and SNU11045 were more similar to YX10 (nt 98--100%, aa 98--100%) than KM91 (nt 72--93%, aa 67--95%), but 5b and N had similarly low identities to KM91 and YX10 (nt 91% and 89%, aa 88% and 92%; nt 92--93% and 88%, aa 88% and 92%, respectively). SNU8065 and SNU11045 showed relatively high nucleotide identities when compared to other viruses; specifically for 5b (CK/CH/LHLJ/141105/KP790145 and CK/CH/IBYZ/2011/KF663561, 99%) and N (SD10/KF625030, CK/CH/SD09/006/HM230751, and GX-NN120084/KF996278, 99%). The amino acid identities of nsps were examined with respect to the nsps of KM91 and YX10 (Table [3](#Tab3){ref-type="table"}). The identities were different among nsps and were higher for nsp 2, 3, 4, 5, and 6 than for the other nsps, which indicated that the former were more variable.Table 3Nonstructural protein (nsp) amino acid identities between the QX-like viruses and KM91 and YX10NSPPutative function^b^Protease^\*^Size (aa)Amino acid identity with respect to KM91/YX10 (%)Percent differences^b^QIA-03342QIA-D79/05QIA-Q43/06SNU-9106SNU-10043SNU-8065SNU-110452UnknownPLP67395/8596/8696/8695/8695/8686/9886/989--12%3Papain-like viral proteasePLP1,59495/8695/8694/8594/8594/8584/9884/989--14%4UnknownPLP/3CLpro51498/8898/8998/8998/8897/8887/10087/1009--13%5Coronavirus endopeptidase C303CLpro30798/9198/9098/9197/9098/9089/9990/997--10%6Hydrophobic domain3CLpro29399/8699/8797/8698/8798/8687/9987/9911--13%7nsp7 superfamily3CLpro83100/9899/9699/9699/9696/9496/9996/992--4%8nsp8 superfamily3CLpro21098/9499/9399/9398/9397/9294/9894/984--6%9nsp9 superfamily3CLpro111100/100100/100100/100100/100100/10098/98100/1000--2%10nsp10 superfamily, RNA synthesis3CLpro14597/9499/9697/9498/9499/9595/10094/993--5%12RNA-dependent RNA polymerase3CLpro94099/9699/9698/9699/9698/9696/9996/992--3%13Viral RNA helicase3CLpro60098/9798/9799/9799/9899/9897/10097/991--3%14nsp11 superfamily; exoribonuclease3CLpro52199/9799/9799/9799/9798/9797/9997/992%15Nidoviral uridylate-specific endoribonuclease3CLpro33897/9597/9598/9596/9496/9394/9694/962--3%1623S rRNA methylase3CLpro30299/9699/9698/9798/9698/9696/9996/991--3%^\*^PLP, papain-like protease; 3CLpro, 3C-like protease; ^b^ difference between identities with respect to KM91 and YX10

Comparison of regulatory sequences {#Sec10}
----------------------------------

We examined the TRSs of each IBV and compared nucleotide sequences between the leader and body TRSs (Table [4](#Tab4){ref-type="table"}). The leader TRS sequences of all viruses were identical (CTTAACAA). The leader TRS sequence of each virus had 0--3 nucleotide mismatches with the body TRS sequences, and the body TRS sequences of the 4 and 6 sg-mRNAs were identical to the leader TRS sequence.Table 4Comparison of the transcription regulatory sequences (TRSs) in the full IBV genomesSg mRNATRSKM91Sequence and location of TRSQIA-03342QIA-KR/D79/05QIAQ-43/06SNU9106SNU10043YX10SNU8065SNU11045LeaderCTTAACAA(57-64)CTTAACAA(57-64)CTTAACAA(57-64)CTTAACAA(57-64)CTTAACAA(57-64)CTTAACAA((57-64)CTTAACAA(57-64)CTTAACAA(57-64)CTTAACAA(57-64)1(S)BodyAAG\-\-\-\--(20,311-20,318)AGG\-\-\-\--(20,314-20,321)AGG\-\-\-\--(20,314-20,321)AAG\-\-\-\--(20315-20322)T--T\-\-\-\--(20,314-20,321)AGG\-\-\-\--(20316-20323)\--G\-\-\-\--(20,311-20,318)\--G\-\-\-\--(20315-20322)\--G\-\-\-\--(20,314-20,321)2(3ab)Body\--G\-\-\-\--(23,831-23,838)\--G\-\-\-\--(23840-23847)\--G\-\-\-\--(23840-23847)\--G\-\-\-\--(23835-23842)\--G\-\-\-\--(23840-23847)\--G\-\-\-\--(20342-23849)\-\-\-\-\-\-\--(23837-23844)\--G\-\-\-\--(23835-23842)\--G\-\-\-\--(23,838-23,845)3(E)Body\--C\-\-\-\--(24,035-24,042)T-C\-\-\-\--(24088-24094)T-C\-\-\-\--(24087-24094)T-C\-\-\-\--(24082-24089)T-C\-\-\-\--(24087-24094)T-C\-\-\-\--(24089-24096)TAC\-\-\-\--(24084-24091)TAC\-\-\-\--(24089-24096)TAC\-\-\-\--(24,085-24,092)4(M)Body\-\-\-\-\-\-\--(24,377-24,384)\-\-\-\-\-\-\--(24,435-24,442)\-\-\-\-\-\-\--(24,435-24,442)\-\-\-\-\-\-\--(24430-24437)\-\-\-\-\-\-\--(24,435-24,442)\-\-\-\-\-\-\--(24437-24444)\-\-\-\-\-\-\--(24429-24436)\-\-\-\-\-\-\--(24434-24441)\-\-\-\-\-\-\--(24,435-24,442)IR(4bc)BodyG-C\-\-\-\--(24850-24857)G-C\-\-\-\--(24905-24912)G-C\-\-\-\--(24905-24912G-C\-\-\-\--(24899-24906)G-C\-\-\-\--(24934-24941)G-C\-\-\-\--(24905-24912)G-C\-\-\-\--(24896-24903)G-C\-\-\-\--(24896-24903)G-C\-\-\-\--(24896-24903)5(5ab)Body\-\-\-\-\-\-\--(25,485-25,492)\-\-\-\-\-\-\--(25543-25550)\-\-\-\-\-\-\--(25540-25547)\-\-\-\-\-\-\--(25538-25545)\-\-\-\--TAA(25543-25550)\-\-\-\-\-\-\--(25545-25552)\-\-\-\--G-G(25534-25541)\-\-\-\--G-G(25539-25546)\-\-\-\--G-G(25,535-25,542)6(N)Body\-\-\-\-\-\-\--(25,786-25,793)\-\-\-\-\-\-\--(25844-25851)\-\-\-\-\-\-\--(25841-25848)\-\-\-\-\-\-\--(25839-25846)\-\-\-\-\-\-\--(25844-25851)\-\-\-\-\-\-\--(25846-25854)\-\-\-\-\-\-\--(25839-25846)\-\-\-\-\-\-\--(25844-25851)\-\-\-\-\-\-\--(25,840-25,847)6bBodyAGG\-\-\-\--(27,004-27,011)AAG\-\-\-\--(27,059-27,066)AGG\-\-\-\--(27,056-27063)AAG\-\-\-\--(27,053-27060)AGG\-\-\-\--(27,059-27,066)AGG\-\-\-\--(27059-27066)--GG\-\-\-\--(27,054-27,061)AGG\-\-\-\--(27,054-27061)AGG\-\-\-\--(27,054-27,061)

Computational recombination analysis {#Sec11}
------------------------------------

Recombination events were detected in the S1 gene of QIA-03342, QIA-KR/79/05, QIA-Q43/06, SNU9106, and SNU10043, with high significance (Fig. [1](#Fig1){ref-type="fig"}). The genomic backgrounds of QIA-03342, QIA-KR/79/05, QIA-Q43/06, SNU9106, and SNU10043 were KM91, and the partial or complete S1 gene might have been acquired from QX-like viruses. However, the patterns of S1 recombination differed among the viruses. Additional recombination events were suspected in S2 (QIA-03342 and SNU10043), 3a/3b (all), and 4b/4c (QIA-03342), but statistical significance was not obtained for these analyses.Fig. 1Using the full genome sequences computational recombination analyses were performed with RDP (ver. 4.14; neighbor-joining, Kimura-distance measure method; P \< 0.05). Seven algorithms in RDP 4.14, including RDP, GENECONV, Bootscan, MaxChi, Chimaera, SiScan and 3Seq were used to evaluate the recombination events. P values for the Bootscan are represented. The locations of the S1 and S2 genes are represented by the red boxes

Phylogenetic analysis {#Sec12}
---------------------

Based on a phylogenetic analysis using the nucleotide sequences of S1, QIA-03342, QIA-KR/D79/05, SNU8065, and SNU11045 these isolates were classified as QX-like while QIA-Q43/06, SNU9106, and SNU10043 were classified as NC1 sub-genotypes, as previously reported \[[@CR24], [@CR31]\]. Interestingly, the European QX-like viruses ITA/90254/2005 and CK/SWE/0658946/10 clustered with the early QX-like viruses QIA-03342 and QIA-KR/D79/05 (Fig. [2](#Fig2){ref-type="fig"}a) \[[@CR1]\]. QIA-03342, QIA-KR/D79/05, QIA-Q43/06, SNU-9106, and SNU-10043 clustered with KM91, and SNU8065 and SNU-11045 clustered with YX10 in the phylogenetic tree constructed with full genome sequences. ITA/90254/2005 and CK/SWE/0658946/10 clustered together, but they did not form a cluster with the Belgian nephropathogenic reference strain B1648 (Fig. [2](#Fig2){ref-type="fig"}b). The Europe and American strains formed a cluster, but the European QX-like viruses clustered with strains from China, including QX-like viruses. To remove the effects of highly variable and frequently recombining genes, i.e., S, 3a/3b, and 4b/4c, on the phylogeny, we constructed a phylogenetic tree using the minimum essential core replication gene 1ab (Fig. [2](#Fig2){ref-type="fig"}c). The topology and clustering pattern of the 1ab-based tree were similar but not identical to those of the whole genome tree. Based on the clustering patterns, we classified IBVs into North America-Europe, China/Taiwan, and China lineages. Unexpectedly, KM91, the early QX-like viruses, and the NC1 viruses formed a cluster with the North America-Europe lineage viruses. CK/SWE/0658946/10 clustered with the French strain 4/91, but ITA/90254 clustered with B1648. Thus, viruses with the 1ab genes of the North America-Europe lineage might have acquired Spike genes directly or indirectly from QX-like viruses in the China lineage. SNU-8065 and SNU-11045 clustered with viruses from the China lineage, including LX4, DY07, and YX10. Ck/CH/LSD/11235, GX-YL9, and Ck/CH/LGD/120723 clustered with viruses from the China lineage but did not form a cluster with QX-like viruses in the S1 tree. Thus, they possessed QX-like 1ab but non-QX-like spike genes.Fig. 2Phylogenetic trees based on: 1.5 kb nucleotide sequences of S1 (a), full genome (b), and 1ab genes (c), were constructed with the neighbor-joining method using MEGA 5.05. The bootstrap values were determined from 1000 replicates of the original data. The branch number represents the percentage of times that the branch appeared in the tree. Bootstrap values greater than 90% are shown. The p-distance is indicated by the bar at the bottom of the figure. The QX-like viruses characterized in the present study are marked with closed circles

Molecular characterization of structural proteins {#Sec13}
-------------------------------------------------

The amino acid sequences of cleavage sites between S1, S2, and the fusion peptide were compared (Table [5](#Tab5){ref-type="table"}) \[[@CR41]\]. The NC1 viruses QIA-043/06, SNU9106, and SNU10043 possessed the same amino acid sequence, RRFRR/S, as KM91 at the cleavage site between S1 and S2, but the other viruses had HRRRR/S. The early QX-like and NC1 viruses possessed PSGR/S at the cleavage site of the fusion peptide except Q43/06 (PRRR/S); however, the recurrent QX-like viruses had PRGR/S. The partial or whole amino acid sequences of B cell epitopes, D (24--61), E (132--149), C/A/B (291--398), and F (497--543) in S1 and G (548--574) in S2 were compared, as summarized in Fig. [3](#Fig3){ref-type="fig"} \[[@CR19], [@CR20]\]. The early QX-like viruses QIA-03342 and QIA-KR/D79/05 exhibited amino acids differences in all of the epitopes, but QIA-043/06, SNU9106, and SNU10043 acquired similar amino acids to KM91 for the C/A/B, F, and G epitopes. However, the C/A/B epitopes of SNU9106 and SNU10043 were more similar at the amino acid level to KM91 than QIA-043/06 was, and they acquired additional amino acid changes in the F and G epitopes. The recurrent QX-like viruses SNU8065 and SNU11045 had different amino acid sequences in all of the epitopes. Similar to KM91, QIA-Q43 had two amino acid deletions, 24G and 25N, in epitope D.Table 5Comparison of functional motifs in the spike (S), envelope (E) and membrane (M) proteinsSEMProteolytic cleavage siteER retention signalVariable amino acid residue in ectodomainS1/S2S2'(fusion peptide)561112131617KM91~534~RRFRR/S~539~^a^~688~ PSGR/S~692~^a^~103~RHGKLHS~109~^b^GESTQAEQIA-03342~536~HRRRR/S~541690~ PSGR/S~694~-^c^-N-EVQQIA-KR/D79/05D\-\-\-\-\--QIA-Q43/06~534~RRFRR/S~539~^a^~688~ PRRR/S~692~^a^\--N-EVQSNU9106~536~RRFRR/S~541690~ PSGR/S~694~\-\-\-\-\-\--SNU10043\-\-\-\-\-\--YX10~536~HRRRR/S~541690~ PRGR/S~694103~RDKLHP~108~E.^d^DSEILSNU8065E.DSEILSNU11045E.DSEIL^a^The different location numbers of KM91 and QIA-Q43/06 are due to the 24G and 25N deletion^b^The different location number is due to an insertion at 105G^c^Same amino acid sequence as KM91^d^Non-comparable amino acid residue due to 6E insertion Fig. 3Comparison of the amino acid sequences of each B cell epitope in the IBV spike proteins. Identical amino acids are represented with dots, while deleted amino acids are represented with hyphens. The amino acids unique to KM91 are represented with red letters, and the amino acids unique to each NC1 virus are represented with blue letters. The putative N-linked glycosylation sites are underlined

The endoplasmic reticulum (ER) retention signal E protein localizes to the ER \[[@CR26]\]. We compared the amino acid sequences of QX-like viruses with KM91 and YX10 (Table [5](#Tab5){ref-type="table"}). The early QX-like and NC1 viruses were one amino acid longer and had a different ER retention signal motif, RHGKLHS (105G insertion), when compared to the recurrent QX-like viruses, RDKLHP. The ectodomain amino acid sequences (1-23 amino acid residues) of M proteins were also compared (summarized in Table [5](#Tab5){ref-type="table"}). Compared to KM91 the early QX-like and NC1 viruses shared a 6E insertion and were one amino acid longer than the recurrent QX-like and YX10 viruses. The amino acid sequences of the recurrent QX-like viruses were more different from KM91 than the early QX-like and NC1 viruses.

Discussion {#Sec14}
==========

In contrast to our expectations, our results revealed that early QX-like viruses were recombinant KM91-like and QX-like viruses. The absence of pure QX-like viruses before the appearance of the early QX-like viruses suggests that recombinant viruses were introduced from other countries. A computational recombination analysis did not support recombination events in the S2, 3a, 3b, 4b, and 4c genes of the early QX-like viruses. However, their low nucleotide and amino acid identities, to both KM91 and YX10, and relatively high identities to other viruses in the North America group suggest potential recombinant events with other unknown viruses. The NC1 viruses may be products of recombination between the early QX-like and KM91-like viruses, as previously reported \[[@CR27]\]. Furthermore, the low nucleotide and amino acid identities of the recurrent QX-like virus 5b and N genes to both KM91 and YX10, but high identities to other viruses in China also suggest frequent recombinant events in the 5b and N genes. A high number of recombination breakpoints were identified in the 1a region, but the probabilities for these hotspots were not significant \[[@CR39]\]. However, recombinant events among IBVs may occur more frequently in the latter fourth of the viral genome than the remainder of the genome containing the 1ab gene \[[@CR22]\]. Thus, phylogenetic analyses based on 1ab may provide more information on the evolution of IBVs. According to the results of the present study, IBVs were classified into the North America-Europe, China/Taiwan, and China genotypes. The clustering of KM91, the early QX-like viruses, and the NC1 viruses with the viruses of the North America-Europe lineage was unexpected, but recent reports on the spread of H5N8 highly pathogenic avian influenza viruses from Korea to North America and Europe as well as their reintroduction to Korea by migratory birds may provide a hint \[[@CR21], [@CR23]\]. Siberia and Beringia are breeding sites for migratory birds from Eurasia and parts of North America, China, Korea and Japan are wintering sites \[[@CR23]\]. The high prevalence and diversity of gamma coronaviruses in Beringia wild birds may contribute to the intercontinental spread and evolution of these viruses \[[@CR32]\]. The presence of the European QX-like viruses ITA/90254/2005 and CK/SWE/0658946/10 in a cluster with the early QX-like viruses in the S1 gene tree may also reflect recurrent recombination events between viruses of the North America-Europe and China lineages. The recent detection of IBV-like gene fragments in the wild bird population may support a possible role for wild birds as reservoirs or long distance carriers \[[@CR12], [@CR17]\]. However, it is not clear whether they play a role in the recombination events of IBV. In addition, it is not clear why the QX-like viruses provide the spike gene rather than the 1ab gene.

The 1ab and N genes are the minimum essential core elements for genomic RNA replication and persistence of coronaviruses in infected cells. The 1ab gene encodes nsps 2--16, which are essential for viral RNA replication and pathogenicity \[[@CR2], [@CR13], [@CR29]\]. Thus, 1ab genotypes that increase replication may be selected for in the competitive conditions of superinfected cells. The early QX-like viruses possessed the 1ab genes of KM91-like rather than QX-like viruses. The predominance of KM91-like 1ab genes supports the hypothesis that they conferred more favorable biological properties to the recombinants than the QX-like 1ab genes. For example, they may have conferred efficient virus replication and extended persistence in the trachea \[[@CR34]\]. Based on the nucleotide and amino acid sequence identities, nsp 2, 3, 4, 5, and 6 were more variable than the other nsps among the viruses (Table [4](#Tab4){ref-type="table"}). Nsp 2 is a papain-like viral protease and nsp5 is a coronavirus endopeptidase. They are important for the post-translational cleavage of nsps \[[@CR13]\]. Nsp 3, 4, and 6 possess transmembrane domains and form an internal membrane scaffolding complex that facilitates the assembly of the membrane-associated replication complex \[[@CR33], [@CR37]\]. The interactions between these molecules and their role in complex formation may facilitate the coevolution of nsp 3, 4, and 6. The deubiquitinating activity of nsp 3 and its frequent mutations during the attenuation of virulent strains as well as the autophagosome-generating activity of nsp 6 may be directly, or indirectly, related to viral competence. This in turn may influence how recombinant viruses overcome host responses and promote their own growth \[[@CR3], [@CR10], [@CR34]\].

The compatibility of the leader TRS sequence with the body TRS sequences is a determinant of viral gene expression levels and viral replication efficiency \[[@CR4]\]. The 2--3 nucleotide mismatches of the body TRS sequences of the early and the NC1 sg-mRNA 1 encoding Spike, which is the major target of humoral immunity, were comparable to the single nucleotide mismatch in recurrent QX-like viruses. The 1--2 mismatches in the sg-mRNA 3 (encoding the E protein) and lack of mismatch in the sg-mRNA 5 (encoding the 5a and 5b proteins) of the early QX-like and the NC1 viruses, except SNU9106, are also consistent with 3 and 2 mismatches, respectively, in the corresponding sg-mRNAs of recurrent QX-like viruses. Differences in the expression of the E, 5a, and 5b proteins may affect virus replication efficiency \[[@CR36]\]. Thus, these mismatches between the leader and the body TRSs may be related to the predominance of certain genomic backgrounds. Furthermore, the higher copy number of the 1(S) sg-mRNA in recurrent QX-like viruses than KM91-like viruses may increase the probability of RNA-dependent RNA polymerase mediated spike gene recombination \[[@CR22]\]. Therefore, further studies to investigate these mechanisms may be valuable in the future.

The proteolytic cleavages of the spike protein are essential for virus infection of host cells. The cleavage site between S1 and S2 contained the furin-recognition motif R-X-R(K)-R, and was shared by all of the compared viruses \[[@CR5]\]. The second proteolytic cleavage site near the fusion peptide of the spike protein is important for virus entry, syncytium formation and virus infectivity \[[@CR41]\]. The second cleavage sites of all the compared viruses contained only 1 to 3 arginines, and they may be cleaved by other serine proteases \[[@CR41]\]. The effects of the single amino acid insertions in the ER motif of E and ectodomain of M proteins are unknown. However, the differing amino acid sequences of the M protein ectodomains, identified when examining recurrent QX-like viruses, encourages further study of a role for IBV antigenic variation \[[@CR14]\].

The prevalence of the early QX-like viruses may be explained by amino acid variation in the Spike protein's epitopes, brought about by immune-selection pressure induced by KM91 vaccination. However, host immune pressure induced by infections with QX-like viruses between 2002 and 2005, and vaccination, might have resulted in the selection of NC1 viruses. Compared to the early NC1 viruses, QIA-Q43/06, SNU9106, and SNU10043 possessed the C/A/B epitopes, which were more similar to KM91 than QX-like viruses; however, they also acquired additional mutations in the F and G epitopes, which were similar, but not identical, to KM91. Despite pressure for the NC1 viruses to evolve new epitopes, a live-attenuated vaccine strain, K2, was successful in the field, and the NC1 has rarely been isolated since 2011 \[[@CR28], [@CR31]\]. However, the prevalence of the recurrent QX-like viruses and their epitope structures may reflect incomplete protective effects or vaccine breaks when conventional vaccines are applied in the field.

Conclusion {#Sec15}
==========

In conclusion, the early QX-like viruses are recombinant viruses between KM91-like and QX-like viruses, and NC1 viruses subsequently evolved via recombination and the accumulation of missense mutations in spike genes. The predominance of KM91-like 1ab and QX-like spike genes in the recombinant viruses suggests that these loci are beneficial with respect to replication efficiency and the capacity to evade immunity, respectively.

Electronic supplementary material
=================================

 {#Sec16}

Below is the link to the electronic supplementary material. Supplementary material 1 (PDF 37 kb)

This study was supported by the BK 21 plus Program for Veterinary Science Research.

This work was partially supported by a QIA grant (No. B-1543084-2016-18-1)

The authors declare no conflicts of interest.

This article does not contain any studies with human participants or animals performed by any of the authors.
