A total of 94 IBV strains isolated from 10 provinces in China, including Shandong, Jiangsu, Liaoning, Heilongjiang, Jilin, Anhui, Hebei, Guangxi, Yunnan, and Shanxi, from 2020 to 2024 were preserved in our laboratory.

The S1 gene was amplified using IBV-degenerate primers (forward: 5′-ATGTTGGKRAMRYCWCTD-3′; reverse: 5′-AGTAGAACGTCTAVRACGAY-3′), with a product of approximately 1700 bp. The RT-PCR conditions were as follows: reverse transcription at 45°C for 30 min; initial denaturation at 94°C for 3 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 52°C for 30 s, and extension at 72°C for 90 s, with a final extension at 72°C for 10 min. The RT-PCR products were purified using the Biospin Gel Extraction Kit (BSC02M1A, BioFlux, China) and then were cloned into the pMD18-T vector (6011, Takara, China). The positive samples were sent to Shanghai Shenggong Biological Co., Ltd. for sequencing.

Nucleotide and amino acid sequence analyses and alignment were performed on the S1 gene of the isolated strain and reference strains of different IBV genotypes downloaded from GenBank using DNAStar software. Following the typing method established by Valastro et al., a phylogenetic tree was constructed using the maximum likelihood method in MEGA 7.0, and 1,000 bootstraps were calculated (Valastro et al., 2016). Genetic distances were calculated, and a neighbor-joining (N-J) evolutionary tree was generated. The potential N-glycosylation and O-glycosylation sites of the S1 gene were analyzed using the online bioinformatics tools NetNGlyc 1.0 Server and NetOGlyc 4.0 Server, respectively (Fan et al., 2019).

Potential recombination events were evaluated using a dataset of S1 gene sequences, including both isolates and reference strains. Analyses were performed with RDP 5.42 software using seven detection methods (RDP, CHIMAERA, BOOTSCAN, 3SEQ, GENECONV, MAXCHI, and SISCAN). The recombination events were only considered to be significant the following criteria were met: detection by at least five of the seven methods with P-values < 1×10⁻¹² (Mo et al., 2013). Subsequently, potential recombination events were further verified with Simplot 3.5.1 software (Feng et al., 2018; Jiang et al., 2018).

The JModelTest software was used to identify the best model of nucleotide substitution based on 88 candidate models (Darriba et al., 2012), and the optimal model result was GTR+γ+I. A Bayesian tree was constructed using BEAST v1.10.4 software. A Strict molecular clock and Bayesian skyline aggregation model were selected, and the S1 gene was analyzed using the Markov Chain Monte Carlo (MCMC) statistical algorithm. The final effective sample size (ESS) of the operation chain was >200, and the number of operations was 400 million (Zhao et al., 2016). After setting the parameters, BEAST v1.10.4 software was used to generate the xml file and Tracer v1.7.1 was employed to carry out various numerical analyses. The molecular evolution rate of S1 was also confirmed (Baele et al., 2013).

Representative strains of GI-19, GI-13, GI-22, and GI-1 genotypes were inactivated for vaccine preparation and five 4-week-old SPF chickens were immunized in each group. Booster immunization was conducted 14 days after primary vaccination. Blood samples were collected 14 days after the booster vaccination, and sera were separated and stored at -20°C. The IBV strains of the target genotypes were pre-titrated, and the serum was serially diluted with 100 EID₅₀ virus. The mixtures were inoculated into 9-11-day-old SPF embryo (n=5 per dilution). The SPF embryos were observed for 7 days to calculate neutralization titers (highest serum dilution inhibiting ≥50% egg death/dwarfing). Neutralization titers are expressed as log2 values.

The SPF chickens were purchased from Shandong Haotai Experimental Animal Breeding Co., Ltd. The experimental protocol was approved by the Ethical Committee for Animals in Shandong Academy of Agricultural Sciences, China.

To investigate the molecular characteristics of IBV circulating in China, 94 IBV strains were isolated from 10 provinces between 2020 and 2024. The nucleotide and deduced amino acid sequences were systematically aligned against representative reference strains spanning the major genotypes (per the Valastro classification). The homology among the different strains was presented in Table 1 . The phylogenetic tree was constructed using MEGA 5.0. The 94 IBV isolates were clustered into five GI genotypes: GI-19 (70.21%), GI-13(21.28%), GI-22 (4.26%), GI-28 (2.13%), and GI-1 (2.13%) ( Figure 1A ). Notably, GI-19 was the predominant lineage in all years ( Figure 1B ).

Hypervariable regions (HVRs) may account for the antigenicity and serotype variation and associated with receptor binding. A comparative analysis of the amino acid sequences of the HVRs between the 94 IBV strains and their genotype-matched reference strains was performed using DNASTAR software. The results demonstrated different degrees of amino acid mutations within the HVR, along with sporadic mutations and deletions observed at other sites. Amino acid sequence alignment of the 66 GI-19 isolates with the reference strains LX4 and QX revealed that the amino acid mutations were predominantly located at positions 58, 65, 70, 121, 126, 129, 284, and 293 within S1 ( Figure 2A ). Amino acid sequence alignment between the 20 GI-13 isolates and 4/91 reference strain demonstrated that seven isolates exhibited similar variations, including deletions and point mutations. Notably, these mutations predominantly clustered in HVR I and HVR II ( Figure 2B ).

Analysis of the glycosylation sites of the 94 isolates indicated that only SDJN3/23 exhibited one potential O-glycosylation site at position 416 in the S1 gene and no O-glycosylation sites were predicted in the remaining isolates. All isolates harbored varying numbers of N-glycosylation sites, with significant differences in both quantity and positional distribution across the five genotypes. Compared to the reference strain LX4, all GI-19 isolates (66 isolates) contained 14–18 predicted N-glycosylation sites. Although 52 isolates showed nearly-identical patterns to LX4, 14 isolates displayed only minor positional variations ( Table 2 ). The GI-13 genotype (20 isolates) consistently carried 16–17 N-glycosylation sites. Notably, seven isolates shared a novel site (NGSL) at position 143 with reference strain 4/91 ( Table 3 ). The remaining 12 isolates largely matched the reference patterns. All 4 GI-22 genotype isolates possessed 16 N-glycosylation sites, compared to the CH/LSC/99I reference strain, uniformly lacking 3 sites: NKTN (position 53), NNST (167), and NVTG (518). Two GI-1 isolates retained 16 N-glycosylation sites, indicating complete conservation with the Mass reference strain. Both JSTZ-C/20 and SD/WF/202307 contained 16 N-glycosylation sites. A shared novel site (NVTN) was identified at position 279 compared with the LDT3 reference strain.

The cleavage sites of all isolates were identified, and the results showed that there were 5 cleavage site sequences: HRRRR (58 isolates, 61.7%), RRSRR (20 isolates, 21.28%), HRHRR (7 isolates, 7.45%), RRFRR (5 isolates, 5.3%), and HRRKR (4 isolates, 4.26%). The cleavage site of the S protein of GI-19 (56/66) isolates was the same as that of the reference LX4 (HRRRR), however, three variants were observed: HRRKR, HRHRR, and RRFRR. The cleavage sites of the 20 GI-13 isolates were the same as those of the reference strain 4/91 (RRSRR). The cleavage sites in the four GI-22 isolates were HRRRR and HRRKR. Both GI-1 and GI-28 exhibited uniform RRFRR. Details of the S1 sequence lengths and cleavage site distributions of the isolates are shown in Table 4 .

Recombination analysis revealed that 14 of the 94 IBV isolates (14.89%) exhibited statistically significant recombination events in S1 ( Table 5 ), as detected by RDP 5.42 (P < 0.05). The predominant recombinations involved the GI-13/GI-22 (7/14), GI-19/GI-7 (4/14), and GI-19/GI-22 (3/14) genotypes. To validate the recombination events identified by RDP 5.42, three representative strains (SDJN3/23, SDLY1/23, and JS3/20) exhibiting distinct genotype combinations were selected for in-depth analysis using Simplot v3.5.1 ( Table 6 ). These results indicated that SDJN3/23 was a recombinant strain derived from the parental strains GI-19 LX4 and GI-22 SAIBK, with a breakpoint at position 1102bp. SDLY1/23 was derived from parental strains GI-13 4/91 and GI-22 CK/CH/LSC/99I with a breakpoint at position 720 bp. JS3/20 was derived from the parental strains GI-19 LX4 and GI-7 TW2575/98, with breakpoints at position 604 bp and 1288 bp ( Figure 3 ). These findings are consistent with those of the RDP 5.42 analysis, confirming the occurrence of recombination events in these isolates.

The evolutionary rate of the remaining 80 isolates was analyzed after excluding 14 recombinant strains. The optimal model as GTR+γ+I was selected using JModelTest software, and estimated under a strict molecular clock with Bayesian skyline coalescent model. The mean nucleotide substitution rate of the isolated S1 gene was 1.98×10–⁴ substitutions/site/year, and the 95% HPD interval was 1.55-2.46×10^-4, as estimated using BEAST software ( Table 7 ).

Cross-neutralization tests were conducted between the main genotypic representative isolates SDJN3/23 (GI-19), JS3/20 (GI-22), and SDLY1/23 (GI-13), and the classical vaccine strain H120 (GI-1). The results showed that JS3/20 had the highest self-neutralization titer (1:256), but minimal cross-reactivity with GI-19/GI-13/GI-1 sera (<1:8). The SDJN3/23 strain, isolated in 2023, exhibited a strong self-neutralization titer (1:128) and demonstrated good cross-neutralization with GI-13, GI-1 and GI-22 sera (≥1:32) ( Table 8 ).