Full-length genome reveals genetic diversity and extensive recombination patterns of Saudi GI-1 and GI-23 genotypes of infectious bronchitis virus

Abstract

Background

Despite numerous genetic studies on Infectious Bronchitis Virus (IBV), many strains from the Middle East remain misclassified or unclassified. Genotype 1 (GI-1) is found globally, while genotype 23 (GI-23) has emerged as the predominant genotype in the Middle East region, evolving continuously through inter- and intra-genotypic recombination. The GI-23 genotype is now enzootic in Europe and Asia.

Methods

Over a 24-month period from May 2022 to June 2024, 360 samples were collected from 19 layer and 3 broiler poultry farms in central Saudi Arabia. The chickens exhibited reduced laying rates and symptoms such as weakness and respiratory distress, while broilers showed respiratory issues. Samples, including tracheal swabs and various tissue specimens, were pooled, homogenized, and stored at -20 °C prior to PCR analysis. The samples underwent virus isolation in embryonated chicken eggs, RNA extraction using automated systems, and detection of IBV through real-time RT-PCR targeting a conserved 5’-UTR fragment. Full-length genome sequencing was performed, and recombination analysis was conducted using RDP 4.6.

Results

Saudi IBV strains were found to cluster into genotypes GI-1 and GI-23.1. The study identified critical amino acid substitutions in the hypervariable regions of the spike protein and detected recombination events in the ORF1ab, N, M, 3ab, and 5ab genes, with nsp3 of the ORF1ab showing the greatest number of recombination events.

Conclusion

The multiple inter- and intra-genotypic recombination events that were detected in different genes indicate that the circulating IBV strains do not share a single ancestor but have emerged through successive recombination events.

Introduction

Avian infectious bronchitis (IB) is a viral disease caused by avian infectious bronchitis virus (IBV), resulting in substantial losses within the poultry industry worldwide. IBV has been prevalent in the Middle East and North Africa, with limited success in controlling its spread or evolution. Moreover, the continuous emergence of new IBV variants adds to this challenge [1].

IBV belongs to the subgenus Igacovirus within the genus Gammacoronavirus, Orthocoronavirinae subfamily, Coronaviridae family, and Nidovirales order [2]. It is a single-stranded, non-segmented, positive-sense RNA virus with a large genome size of 27.6 kb [3]. The IBV genome encodes structural and non-structural proteins essential for viral replication, with a genome organization of 5’UTR-ORF1a/b-S-3a-3b-E-M-4b-4c-5a-5b-N-6b-3’UTR. ORF1ab constitute two-thirds of the IBV genome and encode polyproteins, which proteolytically cleave papain-like and 3 C-like proteinases into 15 non-structural proteins (nsp2–nsp16) [4].

The S protein, M protein, E protein, and N protein are structural proteins crucial for virus assembly and entry [5]. The S protein is highly glycosylated, and its actual molecular weight reaches approximately 200 kDa [6]. Cleavage of the spike protein occurs at the highly basic furin consensus motif RRFRR, resulting in the formation of two subunits, S1 and S2. The S1 subunit carries the N-terminal group and is primarily responsible for attaching the virus to the cellular membrane through interactions with cellular receptors. Notably, it encompasses three hypervariable regions (HVRs), located at positions 38–67, 91–141, and 274–387, which contribute to its variability and aid in evading host immune responses [6].

Infectious bronchitis virus (IBV) exhibits remarkable genetic and serological variation. This genetic variability arises from mechanisms such as amino acid substitution, insertion and deletion, as well as genetic recombination. These genetic alterations occur during replication because of the error-prone nature and random template switching of the coronavirus RNA-dependent RNA polymerase [7]. As a result, new variants and genotypes of IBV are constantly emerging worldwide. Genotyping relies on genetic sequencing to classify IBV strains into genotypes on the basis of S1 gene diversity. Phylogenetic analysis revealed seven genotypes (GI-GVII) that contain 35 genetic lineages. The GI-23 lineage, which is prevalent in the Middle East, poses a challenge globally [8,9,10,11,12]. GI containing the greatest number of genetic lineages, including the Massachusetts (Mass.) type (GI-1) and Egyptian variants I and II, as well as many Middle East strains (GI-23). The GI-23 lineage has spread to several countries in Africa, Asia, and Europe [8,9,10,11,12]. However, many IBV strains remain unclassified due to incomplete S1 sequences [8,9,10,11,12].

In Saudi Arabia, IBV infection was initially identified in 1984 through the detection of an IBV isolate via RT‒PCR for the N gene, although its serotype was not specified [14]. Serological detection of IBV strains related to serotype 793/B (GI-13) occurred in Saudi Arabia in 1997 and 1998, with further confirmation in 2002 [15]. Subsequent studies in 2009 and 2010 characterized two IBV strains (IBV/CHICKEN/KSA/101/2010 and IBV/CHICKEN/KSA/102/2010) related to CH/LDL/011 (GI-16) and IBV/INDIA/TN/92/03 (H120, GI-1), respectively [16]. Recent reports from Saudi Arabia indicated that the cocirculation of IBV strains is related to Mass. (GI-1), 4/91 (GI-13), CK/CH/LDL/97I (GI-16), and Middle East IBV (GI-23) [11, 17].

Limited research utilizing full-length genomes has hindered the understanding of recombination mechanisms, particularly in Saudi strains. This study aimed to investigate the presence of IBV as a potential causative agent in layer and broiler farms and to conduct comparative phylogenetic analyses of IBV strains identified in Saudi Arabia, elucidating genetic diversity and recombination events that may aid in understanding the genetic evolution of IBV in the region.

Materials and methods

Sample collection

In the present study, 360 samples from 19-layer poultry farms and 3 broiler poultry farms in central Saudi Arabia were collected over a 24-month period from May 2022 to June 2024. All farms had a history of vaccination with IBV vaccines against the MA5, H120, and IB 4/91 strains, following standard vaccine regimens. Affected chickens presented reduced laying rates and symptoms such as weakness, nasal discharge, dyspnoea, oviduct prolapse, and peritonitis, which was observed during post-mortem examination, resulting in significant fatalities. Four tracheal swabs and six tissue samples, including those from the lungs, trachea, and kidneys, were collected from each farm. After collection, the samples were pooled (10 birds per pool), homogenized with 1 mL of phosphate-buffered saline containing 1,000 U/mL penicillin and 100 µg/mL streptomycin, and stored at -20 °C until processing. All analyses were conducted at the Central Veterinary Laboratory, Ministry of Water and Agriculture, Riyadh, Saudi Arabia.

Virus isolation

Tissue suspensions (10% w/v) were made in sterile phosphate-buffered saline (PBS), clarified via low-speed centrifugation and filtered through 0.45-µm bacteriological filters. The suspected clarified fluid was inoculated into the allantoic cavity of a 10-day-old embryonated chicken egg (0.2 ml/ECE). The inoculated embryos were incubated at 37 °C and 70% humidity and candled daily for 7 days. Embryo death within 24 h post inoculation was considered nonspecific death. Allantoic fluids from multiple eggs are pooled, diluted, and used for subsequent passages, with observable embryonic changes indicating strain adaptation.

RNA extraction

RNA extraction was conducted via automated Magna Pure Compact (Roche) RNA extraction by loading 400 µl of allantoic fluid from the 2^nd ECE passage into a MagNA Pure tube (2.0 ml), including at least one negative extraction control. The extraction process adhered strictly to the manufacturer’s guidelines for each system. The purified RNA was stored at -80 °C until use.

Detection of IBV via RT‒PCR

For the detection of IBV, a real-time RT‒PCR assay was used to target a 143-bp fragment of the 5’-UTR, which is conserved in the Ark, Beaudette, Conn, DE072, Florida, and M41 IBV strains. The primers and probes were designed on the basis of the M41 strain (AY851295) conserved sequence of the 5’ untranslated region (5’UTR). IBV5’GU391 (391-5’-GCT TTT GAGCCT AGC GTT-3’-408), the reverse primer IBV5-GL533 (533- 5’-GCC ATG TTG TCA CTG TCT ATT G-3’-512) and the dual-labelled probe IBV5 G probe (494-5’-FAM CACCAC CAG AAC CTG TCA CCT C-BHQ1-3’-473) [18]. The primers were synthesized by Integrated DNA Technologies (Coralville, IA), and the probe was synthesized by BioSearch Technologies (Novato, CA). RT‒qPCR was performed in a total of 25 µL reactions via one-step Quantitect master mix Probe RT‒PCR (Qiagen, Valencia, CA). The reaction was conducted in a SmartCycler (Cepheid, Sunnyvale, CA) and started with reverse transcription at 50 °C for 30 min; then, the reverse transcriptase enzyme was inactivated by incubation at 95 °C for 15 min, followed by 40 cycles of 94 °C for 1 s and 60 °C for 60 s. For each reaction, the cycle threshold (Ct) number was determined according to the number of PCR cycles at which the fluorescence of the reaction exceeded 30 U of fluorescence, which is the default value for the SmartCycler.

Full-length genome

A PCR was performed according to previous methods [19], and some changes were modified to cover the spike region via oligo Sense primer 5′- AGTBTCYACACAGTGTTAYAAGCG-3 and Anti-sense primer CATAACTAACATAAGGGCAA-3 to amplify approximately 3019 bp beginning from the end of the 1b gene and spanning two hypervariable regions (HVR1–2 & 3), with the remaining segment generated by overlapping genome fragments via 18 IBV consensus primer sets and segments [19]. A one-step RT‒PCR high-fidelity kit containing SuperScript™ III (Invitrogen) was used. Briefly, a 50 µL PCR mixture was prepared containing 7 µL of RNA, 25 µL of 2X reaction mixture, 2 µL of sense primer (10 µM), 2 µL of antisense primer (10 µM), 1 µL of Taq high-fidelity enzyme mixture, and 13 µL of PCR-grade water. The PCR conditions were as follows: reverse transcription at 50 °C for 30 min; one cycle of pre-denaturation at 94 °C for 2 min; 40 cycles of denaturation at 94 °C for 15 s, annealing for 30 s, and extension at 68 °C for 2 min; and a final extension at 68 °C for 10 min. The PCR products were first analysed via 1% gel electrophoresis. Equal volumes of both RT‒PCRs were pooled and used to obtain the specific target size. The PCR product was purified using the AMPure XP cleanup kit (AMPure XP, Beckman Coulter, Inc., Kraemer Blvd. Brea, CA 92821, USA). It was then quantified using a Qubit system (DNA kit, USA) and a Bioanalyzer (DNA assay kit, USA). Finally, the purified product was loaded into the MiSeq reagent kit v2 (2 × 250 bp, 500 cycles). All steps were conducted according to Illumina manufacturing in the sequencing unit of the Central Animal Health Laboratory, Weqaa, Riyadh, Saudi Arabia.

Sequence analysis

The multisequence analysis was conducted with MEGA 5.2. Phylogeny reconstruction was performed via the minimum evolution method, whereas phylogenetic validity was assessed via the interior-branch test with 1000 bootstrap replications. Nucleotide substitution is modelled via the maximum composite likelihood method, which incorporates both transitions and transversions. Uniform rates among sites and homogeneous patterns among lineages were assumed. Tree inference uses the Nearest-Neighbor-Interchange (NNI) as the ML heuristic method, with the initial tree derived from Tamura-Nei nucleotide substitution model. Deduced amino acid analysis of the S protein and nsp3 was conducted via MEGA 5.2.

Recombination analysis

Recombination analysis for Saudi IBV strains was conducted via RDP 4.6, a robust software tool designed for detecting and analysing recombination in nucleotide sequences. Initially, the nucleotide sequences of various genes (ORF1ab, S, M, N, 3ab, and 5ab) from the Saudi IBV strains were aligned via ClustalW to ensure accuracy. These sequences were then imported into RDP 4.6 for comprehensive analysis. The software employs several methods, including RDP, GENECONV, Bootscan, Maxchi, Chimaera, SiSscan, and 3Seq, each offering distinct approaches to identify recombination events. The RDP method detects recombination breakpoints on the basis of statistical models, whereas GENECONV assesses gene conversion probabilities. Bootscan and Maxchi use sliding window techniques to compare the likelihoods of parental sequences, and Chimaera evaluates sequence likelihoods to identify recombination. SiSscan analyses similarity scores across windows, and 3Seq examines three-segment models for parental sequence combinations. The analysis was conducted via default settings to ensure robustness and reliability. The detected recombination events were cross-verified with other methods and manually reviewed to confirm the accuracy of the breakpoints. The findings, including breakpoint positions, recombinant sequences, and parental strains, were systematically recorded and analysed [20].

3D structural modelling using swiss-model

The S1 protein sequences of various Saudi strains were modelled using Swiss-Model, an automated homology-based protein structure modelling platform. Swiss-Model identified homologous templates with the highest sequence identity to the query sequences, and 3D models of the S1 protein for each strain, as well as for the H120 and Egyptian G23 commercial vaccines, were constructed. Molsoft ICM 3.5 was used for further refinement of the generated models. After the initial modelling, energy minimization algorithms were applied to optimize the protein geometry, minimize steric clashes, and improve side-chain conformations. Structural differences between the Saudi strains and the commercial vaccines were assessed by comparing the 3D models of the S1 proteins, with particular focus on the highly variable regions (HVR-1 to HVR-3) and the amino acid substitutions in the spike protein.

Results

IBV detection

The RT‒PCR results for the IBV virus across various farms revealed that six farms tested positive, whereas the remaining farms tested negative. High positive samples were found on four farms (two broiler farms and two-layer farms), with Ct values ranging from 20 to 24 (Table 1). Additionally, two-layer farms presented low-positive samples, with Ct values of 33 and 34. The farms with positive results used a variety of vaccines, including 4/91, Ma5, Ma5 clone30, and H120, with poultry ages ranging from 4 to 50 weeks, and included both layer (seventeen farms) and broiler (four farms) types. In contrast, the farms with negative results, regardless of their vaccine combinations and poultry ages, did not show detectable levels of IBV.

Table 1 Real time RT-PCR detection of IBV 5’-UTR gene in pooled samples from chicken farms in different regions in Saudi Arabia

Phylogenetic analysis

The Saudi IBV in GenBank was found to be clustered in GI-1, GI-16, GI-13 and GI-23. The four IBV strains genotyped in the present study were related to GI-1 (Ck/SA/Kharj-1/2023 and Ck/SA/Shaqra-4/2023) and GI-23 (Ck/SA/Kharj-2/2023 and Ck/SA/Kharj-3/2023). GI-23 is subclustered into 23.1 and 23.2 (23.2.1, 23.2.2, 23.2.3). The Saudi strains in the current study, GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023, were found to be subclustered with GI-23.1, which contains the original Egypt/Beni-Suef/01/Var I strain that was first isolated in 1999. This subcluster also contains MRB02/2016- Iran, Var2-06, which also contains two strains from Saudi Arabia, SA/IH1/12 and IBV/CH/SA/5/2019. Another Saudi strain (IBV/CH/SA/6/2019) was found in the NCBI database and found to be related to GI-23.2.2 together with Egy/Var II. strains and a commercial vaccine strain against Egy Var II (Fig. 1).

Fig. 1

Phylogenetic tree of the full-length spike gene from infectious bronchitis virus isolates in Saudi Arabia. This figure illustrates the phylogenetic tree of the full-length spike gene of infectious bronchitis virus (IBV) isolates from Saudi Arabia. The tree was generated through maximum likelihood analysis with 1000 bootstrap replications to ensure robust statistical support. The isolates from the current study are represented in red font, whereas the vaccine strains are shown in blue font

ORF1ab showed that Ck/SA/Shaqra-4/2023, which was classified as GI-1, was clustered with ORF1ab from GI-13, including the 4/91 vaccine strain, whereas Ck/SA/Shaqra-1/2023, which was classified as GI-1, was clustered with other GI-1 strains. Saudi GI-23 strains did not cluster with other GI-23 strains but clustered with strains belonging to GI-26, GI-14, and GI-19 (Fig. 2).

Both the 3ab and E genes showed that GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023 gained their 3ab genes from GI-19 FN30414|ITA/90,254/2005 (Fig. 2).

The M gene of Ck/SA/Shaqra-4/2023, which was classified as GI-1, was clustered as a subtree from Saudi strains that belong to GI-23.1, whereas Ck/SA/Karj-1/2023 was clustered with other GI-1 strains. The Egyptian GI-23 strains were clustered in another subtree distant from that of the Saudi GI-23 strains. Moreover, the GI-23|Ck/SA/Kharj-2/2023 5ab gene recombined, and the 5ab gene was obtained from GI13 strains, including the 4/91 vaccine (Fig. 2).

N gene analysis revealed that the GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023 genotypes, which are GI-23.1 genotypes, are clustered with GI-13, including the 4/91 vaccine. In contrast, Ck/SA/Kharj-1/2023, the GI-1 genotype, is clustered with other GI strains.

Fig. 2

Phylogenetic tree of the full-length ORF 1ab, 3ab, E, M, 5ab, and N genes of infectious bronchitis virus isolates from Saudi Arabia constructed via maximum likelihood with 1000 bootstrap replications. Figure 3 shows the phylogenetic tree based on the full-length sequences of the ORF1ab, 3ab, E, M, 5ab, and N genes of infectious bronchitis virus (IBV) isolates from Saudi Arabia. The trees were constructed via the maximum likelihood method with 1000 bootstrap replications to ensure robust statistical support. Isolates from the current study are shown in red

Deduced amino acid sequence of the spike protein

An amino acid substitution at 9 V to an A was found in the signal peptide in the IBV-EG/1212B-SP1-2012/vaccine in comparison to the Saudi GI-23 strains (CK/SA/Kharj-2/2023 and CK/SA/Kharj-3/2023). The Saudi G1 genotype (CK/SA/Kharj-1/2023 and Ck/SA/Shaqra-4/2023) showed no amino acid substitutions from the H120 vaccine. The cleavage site between S1 and S2 was found to be RRFRR/S for the genotype 1 Saudi strains but RRTRR/S for Saudi genotype 23 (CK/SA/Kharj-2/2023 and GI-23|CK/SA/Kharj-3/2023) (Suppl. 1).

The numbering of the amino acid sequences was based on the S1 sequence of IBV strain H120 and included signal sequences (18 amino acids, MLVTPLLLVT LLCALCSA). Amino acids 38, 43, 63, and 68 are critical for receptor binding. D 38 to S (GI-23|CK/SA/Kharj-2/2023, GI-23|CK/SA/Kharj-3/2023) was detected, whereas D 38 to T was detected in IBV-EG/1212B-SP1-2012|Vaccine (Suppl. 1). H43 was conserved in all the Saudi strains screened in the present study. G 63 to D was detected in GI-23|CK/SA/Kharj-2/2023, GI-23|CK/SA/Kharj-3/2023), whereas deletion of this residue was detected in the IBV-EG/1212B-SP1-2012|Vaccine (Supp. 1). Moreover, an I to A 68 amino acid substitution was detected in GI-23|CK/SA/Kharj-2/2023, GI-23|CK/SA/Kharj-3/2023) and the IBV-EG/1212B-SP1-2012|Vaccine (Fig. 3, Suppl. 1). GI-23|Ck/SA/Kharj-2/2023 and GI-1|Ck/KSA/Shaqra-1/2023 did not show any amino acid substitutions in HVR-I from the H120 or Ma5 vaccine. Among the 3 HVRs, only the H117Q aa substitution was detected in the GI-1|Ck/KSA/Shaqra-1/2023 Saudi GI-1 strain in comparison with the H120 vaccine. In contrast, 54 V to L, 57 S to A, Q63 to D, and S77 to T were detected in the Ck/KSA/Kharj_3/2023 and Ck/KSA/Kharj_2/2023 from the IS/1949/06 and Egy Var 2 vaccines. S63 deletion was also detected in both GI-23 vaccines (Fig. 3).

Fig. 3

Amino acid sequences of hypervariable regions (HVRs) HVR-I, HVR-II, and HVR-III of the spike protein from Saudi strains of infectious bronchitis virus compared with reference and vaccine strains. This figure illustrates the deduced amino acid sequences of the hypervariable regions (HVRs) HVR-I, HVR-II, and HVR-III of the spike protein from infectious bronchitis virus (IBV) isolates collected in Saudi Arabia. The sequences are compared to those of the reference and vaccine strains

In HVR-II, S105 to T, H/S 119 to N, V121S, S130 to H, and G131N/Y and deletion of the 142 aa residue were detected in IBV-EG/1212B-SP1-2012|vaccine|Var2 in comparison to Saudi GI-23 strains (Ck/KSA/Kharj-2/2023 and Ck/KSA/Kharj-3/2023). In HVR-III, T279Y, S/L285H, T288N, N292H, N295S and I296L were detected in the currently available vaccine of GI-23 (IBV-EG/1212B-SP1-2012|Vaccine||Var2). I609V amino acid substitution was detected in 3/4 of the Saudi strains, including GI-1: Ck/SA/Shaqra-4/2023 and genotypes GI-23|CK/SA/Kharj-2/2023 and GI-23|CK/SA/Kharj-3/2023 (Fig. 3). The sequence variation between Saudi IBVs is summarized in Suppl. 1. The 3D modelling of the S1 protein from different Saudi strains in the current study, along with closely related commercial vaccines, showed marked difference between the G1 and G23 strains. The HVR-1 to HVR-3 regions showed complete identity between PP840344|CK/SA/Kharj-1/2023|G1 and H120, while only H to Q117 was detected in PP840347|CK/SA/Shaqra-1/2023|G1 (Fig. 4a). Figure 4 also illustrates the amino acid differences between H120 and the G23 genotypes of Saudi strains, as well as the currently available commercial vaccine (KU979007|IBV-EG/1212B-SP1-2012). It highlights the amino acid substitutions in HVR-I to HVR-3 of G23 vaccine variant II compared to the Saudi strains of G23 (Figs. 3 and 4). Notably, significant 3D structural changes were observed in the S1 protein of the commercial G23 vaccine compared to the Saudi G23 strains.

Fig. 4

Structural comparison of the S1 protein from different Saudi strains with closely related vaccine strains. (a) S1 protein of the H120 vaccine strain (G1), (b) PP840344|CK/SA/Kharj-1/2023 (G1), and (c) PP840347|CK/SA/Shaqra-1/2023 (G1), representing Saudi strains related to the G1 genotype. (d) S1 protein of the commercial Egyptian G23 strain (KU979007|IBV-EG/1212B-SP1-2012), (e) PP840345|CK/SA/Kharj-2/2023 (G23), and (f) PP840346|CK/SA/Kharj-3/2023 (G23), representing Saudi G23 strains. The amino acid differences in the highly variable regions (HVR-I to HVR-III) between the H120 vaccine strain and the G23 Saudi strains are highlighted. The variations in amino acids in the G23 vaccine strain (d) compared to the Saudi G23 strains (e, f) are also marked. Circles indicating structural differences in the loops of the 3D models

Recombination analysis

Our analysis identified multiple recombination breakpoints with varying parentages and detection methods in ORF1ab (Table 2). For GI-1|Ck/KSA/Shaqra-1/2023, recombination events were detected in the nsp3 (papain-like protease (PLpro)) protein at three breakthrough positions, as were nsp4, nsp5 (3 C-like protease (3CLpro)), nsp6, and nsp10. GI-1|Ck/KSA/Kharj-1/2023 showed recombination events in nsp3-nsp5 and nsp3. It also shared a similar recombination event at two positions in nsp3 with GI-1|Ck/KSA/Shaqra-1/2023. GI-23|Ck/SA/Kharj-2/2023 showed recombination events in nsp3 at four positions and showed recombination events at nsp2 (a replicase product important for proofreading viral replication) and nsp12 (an RNA-dependent RNA polymerase (Pol/RdRp)). GI-23|Ck/SA/Kharj-3/2023 showed similar recombination events as GI-23|Ck/SA/Kharj-2/2023 at 3 positions of nsp3 and a single position of nsp12.

For breakpoint positions 419 to 2254 (nsp2 and nsp3), the recombinant sequence GI-23|Ck/SA/Kharj-2/2023 was found to have recombined with Ck/KSA/Kharj-3/2023, with H120, 4/91, and Ma5 as major parents, detected via RDP, GENECONV, Bootscan, Maxchi, Chimaera, SiSscan, and 3Seq. Similarly, breakpoints from 4244 to 5200 and 6115–6476, both located in nsp3, revealed recombination involving GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023 with distinct minor and major parents, including NGA/A116E7/2006, 4/91, and IBV/ck/EGY-Monuf/NR725/16. Additional recombination events were detected at various positions, such as nsp12 (13756 to 14064), nsp3 (5593 to 5964), and nsp3-nsp5 (4286 to 8252), with diverse recombinant sequences and parental strains, including GI-1|Ck/KSA/Kharj-1/2023 and GI-1|Ck/KSA/Shaqra-1/2023, and parents, such as MF421319| UY/09/CA/01 and KR231009| B1648.

Table 2 Detection of recombination events of the ORF1ab of Saudi IBV strains

Table 3 Detection of recombination events of different ORFs of Saudi IBV strains

Notably, breakpoints of GI-1|Ck/KSA/Shaqra-1/2023 at nsp5-nsp6 (8794–9888) with H120 and Ma5 as major parents and nsp2-nsp16 (19876–791) highlighted recombination involving GI-1|Ck/KSA/Shaqra-1/2023 with a combination of minor and major parents detected by multiple methods, including RDP, GENECONV, and Chimaera (Table 2). In addition, the deduced amino acid sequence of nsp3 also revealed 133’ SAGAGECDTA’ 142 246‘EKSV’249, 295’AEE ’297 insertions in GI-23.1 Saudi strains (Kharj 2, and 3) from GI.19|FN430414|ITA/90254/2005 (Suppl. 2).

The detection of recombination events within various ORFs of Saudi IBV strains revealed distinct patterns of genetic exchange involving multiple strains and parent sequences. For the S gene, a recombination event spanning positions 828 to 4198 was identified involving GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023, with KU238176| D888/2/4/08_IR and JX173489| Eg/CLEVB-1/IBV/012 as major and minor parents, respectively. This event was detected via various methods, including GENECONV, Bootscan, Maxchi, Chimaera, SiSscan, and 3Seq (Table 3).

The M gene of GI-1|Ck/KSA/Shaqra-1/2023 had two recombination events, with GI-23| Ck/KSA/Kharj-2/2023 and H120 being the major parents for the two recombination events. The N gene GI-23|Ck/SA/Kharj-3/2023 underwent two recombination events, with MF421319| UY/09/CA/01 and B164 being the major parents for the two recombination events. GI-23|Ck/SA/Kharj-2/2023 exhibited recombination events, with KY805846|IBV/Ck/EG/CU/4/2014 as the major parent (Table 3).

In 3ab, FN430414| ITA/90254/2005 was the major parent of the recombination event involving GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023. Moreover, at 5ab, ON419883|AR/13/BA/A255 was the major parent of the recombination event detected at GI-23|Ck/SA/Kharj-3/2023 (Table 3).

Discussion

In this study, we investigated IBV strains isolated from different Saudi farms to screen the circulating genotypes in Saudi Arabia. We identified two strains related to GI-1 and two isolates related to GI-23. The former strains were closely related to the M41 and H120 strains, whereas the latter were closely related to Egy. Var I, which was initially isolated from Egypt and other countries in the Middle East [21,22,23]. GI-23 is subclustered into 23.1 and 23.2 (23.2.1, 23.2.2, 23.2.3). The Saudi strains in the current study, GI-23|Ck/SA/Kharj-2/2023 and GI-23|Ck/SA/Kharj-3/2023, were found to be subclustered with GI-23.1. This subclade was previously assigned to the EGY-var1/ISR-var2 group, which is now clustered within subclade GI-23.1, a subclade that contains not only viruses from Egypt or the Middle East but also from Turkey, Poland and Romania [1, 24].

The full-genome sequence of IBV obtained via next-generation sequencing (NGS) enables further analysis of individual genes of the IBV genome. GI-23.1 is clearly not the only ancestor of the Saudi strains but is recombinant: both inter- and intragenotypic gene exchange also occurs in ORF1ab, S, N, 3ab and 5ab. in the GI-23|Ck/SA/Kharj-2/2023 and/or GI-23|Ck/SA/Kharj-3/2023 Saudi strains. ORF1ab revealed that FN430415| NGA/A116E7/2006 and H120 were the major parents for recombination. In addition, 4/91 is considered one of the minor parents of such a possible recombination event in ORF1ab. The spike gene also recombined with KU238176| D888/2/4/08_IR and JX173489| Eg/CLEVB-1/IBV/012 as possible parents. KY805846|IBV/Ck/EG/CU/4/2014, MF421319| UY/09/CA/01, and B164 are major parents for three different recombination events in the N gene. In contrast, FN430414| ITA/90254/2005 and ON419883|AR/13/BA/A255 are major parents for recombination events at 3ab and 5ab, respectively. 4/91, B164 and IBV/Ck/EG/CU/4/2014 were also found to be minor parents for 3ab recombination. On the other hand, GI-1|Ck/KSA/Shaqra-1/2023 showed evidence of recombination at many sites in ORF1ab, as did the M and N genes, with 4/91 Ma5, H120, UY/09/CA/01 ON419883|AR/13/BA/A255, and NGA/A116E7/2006 as potential major parents.

Previous studies reported intra- and interrecombination of Egyptian viruses at the S1 level [11, 21, 25]. These recombination events involving different viruses revealed a puzzling relationship between Saudi Arabian strains and IBV strains in the Middle East, Africa and Europe and revealed the possible role of live attenuated vaccines in the evolution of new recombinant strains of IBV.

Amino acids at positions 38, 43, 63, and 68 play crucial roles in the binding of the S1 spike protein to host tissues [26]. D 38 to S, G 63 to D, and I 68 to A amino acid substitutions were detected in GI-23|CK/SA/Kharj-2/2023 and GI-23|CK/SA/Kharj-3/2023. The IBV-EG/1212B-SP1-2012|Vaccine presented D 38 to T, whereas deletion of this residue was detected in the IBV-EG/1212B-SP1-2012|Vaccine, which also presented I 68 to A amino acid substitutions. The amino acid substitution V617I was found to be associated with the chicken kidney cell tropism of changes in the YZ120 strain [27]. Interestingly, this amino acid substitution (V609I) based on the current aa numbering was detected in GI-1: Ck/SA/Shaqra-4/2023 and genotype GI-23|CK/SA/Kharj-2/2023 and GI-23|CK/SA/Kharj-3/2023, which might reflect potential nephrogenic tendencies.

The current commercial GI-23 vaccine is related to GI-23.2.2, which is assumed to be not protective against Saudi strains in the current study, which belong to GI-23.1. This speculation was augmented by the presence of multiple amino acid substitutions, including 6 amino acid substitutions in HVR-I, 11 amino acid substitutions in HVR-II and 7 amino acid substitutions in HVR-III. Some serotypes present as few as 10 amino acid differences in S1 [13], indicating that a small number of epitopes might be responsible for most of the VN antibody response [28].

Recombination is used by viruses for continuous evolution and as a mechanism of evasion of immune pressure. This may lead to the exchange of some parts of genes or alleles with more advantageous combinations. Recombination among different IBV subgroups can lead to altered pathogenicity and transmissibility. Previous research has highlighted how various recombination events between field and vaccine strains contribute to the pathogenic characteristics of IBV, particularly when complete genomes that result in increased pathogenicity and replication capacity in SPF chickens are examined [29]. Additionally, ORF 1ab and accessory proteins, including 3a, 3b, 5a, and 5b, are believed to play significant roles in viral replication and immune evasion, contributing to virus pathogenicity [30, 31]. In this study, we identified multiple unique recombinants resulting from recombination between field strains and vaccine strains, as well as among field strains. ORF1ab across several strains. GI-1|Ck/KSA/Shaqra-1/2023 exhibited notable recombination in nsp3, nsp4, nsp5, nsp6, and nsp10. GI-1|Ck/KSA/Kharj-1/2023 shared some nsp3 recombination events with the Shaqra-1 strain. GI-23|Ck/SA/Kharj-2/2023 exhibited recombination in nsp3, nsp2, and nsp12, whereas GI-23|Ck/SA/Kharj-3/2023 presented nsp3 and nsp12 recombination patterns similar to those of Kharj-2. Additionally, recombination events were detected at three positions in the nsp3 of the GI-23.1 Saudi strains (Kharj 2 and 3), which originated from GI.19|FN430414|ITA/90254/2005. The recombination breakpoints were predominantly located in nsp3, which aligns with findings from previous studies [32,33,34]. These multiple recombination events in nsp3 indicate that it is a hotspot for IBV evolution. This speculation aligns with previous findings that revealed that nsp3 is among the hot spots of coronavirus evolution [35].

Similarly, a study of 212 Chinese IBV isolates revealed that most breakpoints were in nsp3 and located both midstream and downstream of the S gene [33]. Notably, nsp3 plays a crucial role in the replication and transcription processes of coronaviruses (CoVs). It is characterized by eight distinct domains. These domains include ubiquitin-like domain 1 (Ubl1), a Glu-rich acidic domain often referred to as the hypervariable region, a macrodomain (commonly called the X domain), another ubiquitin-like domain (Ubl2), and papain-like protease 2 (PL2pro). Additionally, nsp3 contains an ectodomain known as 3Ecto or the zinc-finger domain, as well as two unidentified domains, Y1 and CoV-Y. Nsp3 also features two transmembrane regions, TM1 and TM2 [36].

We also detected recombination events in the S, M, N, 3a and 5a genes. Although the S gene is a hotspot for recombination [37], we detected only a single intergenotypic recombination event in this gene in the present study involving GI-23|Ck/SA/Kharj-2/2023, with the major parent KU238176|D888/2/4/08_IR.

Both inter- and intragenotypic recombination were also detected in the M, N, 3ab and 5ab genes in the present study. The M gene of GI-1|Ck/KSA/Shaqra-1/2023 underwent both intra- and intergenotypic recombination with the GI-1 and GI-23 strains. The N gene in GI-23|Ck/SA/Kharj-3/2023 presented two intergenotypic recombination events, with MF421319|UY/09/CA/01 and B164. GI-23|Ck/SA/Kharj-2/2023 exhibited intragenotypic recombination in the N gene. In the 3ab gene, both GI-23 strains underwent intergenotypic recombination, with FN430414|ITA/90254/2005 as the major parent, whereas GI-23|Ck/SA/Kharj-3/2023 exhibited intergenotypic recombination in the 5ab gene, with ON419883|AR/13/BA/A255.

In the present study, the cocirculation of various live attenuated vaccines and wild-type strains led to both inter- and intragenotypic recombination events across multiple genes, including the nsp3, nsp4, nsp5, nsp6, nsp10, nsp12, and nsp16 genes, as well as the S, M, N, 3ab, and 5ab genes. Previous studies have reported both inter- and intragenotypic recombination of Egyptian viruses at the S1 level [11, 21, 25].

Conclusion

In the present study, both the GI-1 and GI-23.1 IBV strains were detected on different poultry farms. Many point mutations of the spike gene in the HVRs of the spike protein have been detected, which may theoretically lead to immune pressure and the evolution of variant strains that enable viruses to infect IBV-vaccinated populations. There are many recombination events, with nsp3 resulting in the greatest number of recombination events. The multiple inter- and intragenotypic recombinations of different genes indicate that the circulating IBV strains do not share a single ancestor but have emerged through successive recombination events. Further investigations are needed to understand the pathological mechanisms associated with recombinant IBV. The currently available vaccines used to protect poultry populations against IBV need to be re-evaluated for potential cross-protection, and indigenous vaccines are highly recommended for amino acid diversity from IBV GI-1 vaccine strains.