Population Structure, Genomic Features, and Antibiotic Resistance of Avian Pathogenic Escherichia coli in Shandong Province and Adjacent Regions, China (2008-2023)

Abstract

Avian pathogenic Escherichia coli (APEC) poses a global threat to poultry health and public safety due to its high lethality, limited treatment options, and potential for zoonotic transmission via the food chain. However, long-term genomic surveillance remains limited, especially in countries like China where poultry farming is highly intensive. This study aimed to characterize the population structure, virulence traits, and antimicrobial resistance of 81 APEC isolates from diseased chickens collected over 16 years from Shandong and neighboring provinces in eastern China. The isolates were grouped into seven Clermont phylogroups, with A and B1 being dominant. MLST revealed 27 STs, and serotyping identified 29 O and 16 H antigens, showing high genetic diversity. The minor phylogroups (B2, C, D, E, G) encoded more virulence genes and had higher virulence-plasmid ColV carriage, with enrichment for iron-uptake, protectins, and extraintestinal toxins. In contrast, the dominant phylogroups A and B1 primarily carried adhesin and enterotoxin genes. Antimicrobial resistance was widespread: 76.5% of isolates were multidrug-resistant. The minor phylogroups exhibited higher tetracycline resistance (mediated by tet(A)), whereas the major phylogroups showed increased resistance to third- and fourth-generation cephalosporins (due to bla_CTX-M-type ESBL genes). These findings offer crucial data for APEC prevention and control, safeguarding the poultry industry and public health.

Keywords: antimicrobial resistance; avian pathogenic Escherichia coli; population structure; virulence; whole-genome sequencing.

Figure 1

Phylogenetic diversity and population structure of APEC isolates from eastern China. A neighbor-joining (NJ) phylogenetic tree was constructed based on the core-genome SNPs from 81 APEC strains isolated between 2008 and 2023 in Shandong Province and adjacent regions, China. Metadata—including the sample origin (liver, lung, spleen), isolation year, O:H serotypes, multilocus sequence types (STs), and Clermont phylogroups (A, B1, B2, C, D, E, G)—were annotated using iTOL to illustrate the genetic diversity and population structure.

Figure 2

Virulence characteristics and phylogroup-specific differences in 81 APEC isolates. (A) Detection rates of five major virulence factor categories—adhesins, invasins, iron-uptake proteins, protectins, and toxins—in 81 APEC strains. These factors were identified through the WGS data and annotated with Abricate using the Ecoli_VF database. The ExPEC Virulence Factors denote those associated with systemic infections in humans and animals, while the IPEC Virulence Factors correspond to genes linked to human diarrheagenic strains. (B) Comparison of the number of virulence genes per strain between the major phylogroups (A and B1, n = 58) and minor phylogroups (B2, C, D, E, G, n = 23). (C) Prevalence of the ColV virulence plasmid in the major vs. minor phylogroups. (D) Clade-specific differences in the detection rates of key virulence factors, where darker colors indicate higher prevalence. The p-values in (B) were calculated using the Mann–Whitney U test, and those in (C) were determined by Fisher’s exact test (GraphPad Prism 8.0). *** p < 0.001.

Figure 3

Antibiotic resistance profiles of APEC isolates (n = 81) from chickens with colibacillosis in eastern China over the 2008–2023 study period. (A) UpSet plots reflecting the multidrug resistance patterns of the APEC strains against 8 antibiotic classes. In the central dot matrix, each column represents a unique resistance phenotype profile. The colored dots within the column indicate resistance to at least one antimicrobial agent in the corresponding drug class, with the colors matching specific drug categories. The upper histogram shows the total number of strains with a specific resistance phenotype profile, while the left histogram displays the number of strains resistant to each drug class. Sulfamethoxazole/trimethoprim, a combination antibiotic, was classified under both the sulfonamides and trimethoprim categories for statistical analysis in this study. (B) Comparison of the number of antibiotic classes to which strains from the major phylogroups (A and B1, n = 58) and minor phylogroups (B2, C, D, E, G, n = 23) were resistant. (C) The prevalence of antimicrobial resistance in four common STs of ExPEC. The p-values in (B) were calculated using the Mann–Whitney U test, and those in (C) were determined by Fisher’s exact test (GraphPad Prism 8.0). ns, no significance; * p < 0.05, ** p < 0.01, *** p < 0.0001.

Figure 4

Distribution of the antibiotic resistance determinants across 81 APEC genomes. The central heatmap displays the presence (filled squares) or absence of 54 antibiotic resistance genes and three chromosomal resistance mutations, organized into ten antibiotic-class groups (columns color-coded by class). Above the heatmap, a bar chart quantifies the number of isolates carrying each determinant. To the right, a core-genome SNP phylogeny annotated with the seven Clermont phylogroups links each strain’s lineage to its resistance profile. The asterisk (*) indicates that the gene has undergone at least one point mutation associated with antibiotic resistance.

Figure 5

Comparison of the resistance gene detection rates for β-lactams (A), tetracyclines (B), fluoroquinolones (C), and aminoglycosides (D) across the APEC phylogroups. Notes: The major phylogroups include A and B1 (n = 58), while the minor phylogroups include B2, C, D, E, and G (n = 23). Extended-spectrum β-lactamase (ESBL)-encoding resistance genes are shown in the red font. the p-values in (A–D) were determined by Fisher’s exact test (GraphPad Prism 8.0). ns, no significance; * p < 0.05, ** p < 0.01.