Genetic characteristics and potential pathogenic agents in Campylobacter upsaliensis based on genomic analysis

ABSTRACT

Campylobacter upsaliensis was the most common Campylobacter species in pets’ gastrointestinal tracts and has been isolated from patients with bacteremia, hemolytic-uremic syndrome, spontaneous abortion, and Guillain-Barré syndrome. However, the genetic characteristics and the full extent of its significance as a human pathogen remain to be fully understood. This study involved an investigation for genomic analysis of 154 strains from different sources and additional antimicrobial resistance profiles of 26 strains for this species. The genomes contained 1,558–1,971 CDS and the genome sizes were estimated to vary from 1.53 Mb to 1.86 Mb, with an average GC content of  34.71%. The entire analyzed genomes could be divided into three clades (A, B, and C) based on ANI and phylogenomic analysis. Significantly, nearly all strains in Clade B were isolated from patient samples, and the virulence-related sequences FlgD, GmhA, and CdtC might serve as determining factors for the classification of Clade B. Half of the tested isolates had MIC values over 64 μg mL⁻¹ for nalidixic acid, gentamicin, and streptomycin. Isolates from pets in China carried more resistant elements in the genomes. This study both provided a comprehensive profile of C. upsaliensis for its genomic features and suggested some pathogenic agents for human infection with this species.

KEYWORDS:

• Campylobacter upsaliensis
• pathogenic agent
• comparative genomics
• virulence factor
• antimicrobial resistance

Introduction

Campylobacter spp. are gram-negative, mobile curved rods, and emerging zoonotic pathogens responsible for various health issues, including foodborne and waterborne gastrointestinal (GI) illnesses, severe septicemia bloodstream infections, inflammatory bowel disease (IBD), reactive arthritis, and Guillain-Barré syndrome (GBS) in humans [1–4]. The most commonly isolated species were C. jejuni and C. coli [5,6]. Nevertheless, the growing awareness of other emerging Campylobacter species underscored their significance as pathogens in both human and animal contexts [3,7]. Among the emerging Campylobacter pathogens, C. upsaliensis stood out, being frequently reported in the gastrointestinal tracts of animals, particularly dogs and cats, with variable prevalence rates in different sources [7–10]. Despite its prevalence, C. upsaliensis has also been sporadically isolated from patients with conditions such as bacteremia, hemolytic-uremic syndrome, spontaneous abortion, and Guillain-Barré syndrome [11–16]. Consequently, research and monitoring efforts should persist to unravel the potential of C. upsaliensis as a human pathogen.

Antimicrobial resistance within the Campylobacter genus remains a pressing concern. Notably, for C. upsaliensis strains isolated from animals, resistance rates exceeding 40% have been observed for ciprofloxacin, nalidixic acid, gentamicin, streptomycin, and clindamycin. The resistance to macrolide antibiotics, which are typically the first-line treatment option, was equally concerning. Furthermore, a notable trend is the occurrence of multidrug resistance, which has surged to a rate of 50% [9,17].

Whole genome sequencing (WGS) exhibited a high discriminatory capability for analyzing species diversity [18,19]. However, reports on whole-genome sequencing and population genetics of C. upsaliensis have thus far been lacking. This study undertook a comprehensive comparative genomic analysis on the analyzed genomes which involved all genomes currently available in NCBI. Moreover, an assessment of C. upsaliensis prevalence across distinct animal species has been conducted. Additionally, the genetic and antimicrobial resistance analyses were performed on the isolates in this study. Simultaneously, genetic classification was conducted, and certain virulence markers potentially associated with human infection were identified.

Materials and methods

Sampling and bacteria isolation

Animal fecal or rectal douche fluid samples were collected from both pet hospitals and farms. The isolation of Campylobacter spp. was carried out using the Campylobacter isolation kit that incorporated a membrane filter method (ZC-CAMPY-002, Qingdao Sinova Biotechnology, China). To outline the procedure, a 0.4 mL suspension of stool specimen diluted with saline was transferred into a 4 mL enrichment buffer provided within the kit. This enriched suspension was then subjected to incubation at 37°C for 24 h in a microaerophilic atmosphere comprising 5% O₂, 10% CO₂, and 85% N₂. Subsequently, approximately 300 μL of the cultured enrichment suspension and stool specimen suspension diluted with saline were separately applied onto the surface of filters affixed to double medium plates. These plates, containing Karmali and Columbia agar supplemented with 5% defibrinated sheep blood, were positioned in a microaerophilic environment at 37°C for 48 h [20]. Suspected monoclonal colonies were singled out, purified, and stored for subsequent investigations. These isolates underwent preliminary characterization through PCR amplification and analysis of the 16S rRNA [21].

Antimicrobial susceptibility testing

The minimum inhibitory concentrations (MICs) for a total of eleven antimicrobials, including erythromycin, azithromycin, nalidixic acid, ciprofloxacin, gentamicin, streptomycin, chloramphenicol, florfenicol, tetracycline, telithromycin, and clindamycin, were determined for all isolates in this study. This was achieved using the agar dilution method (ZC-AST-001, Qingdao Sinova Biotechnology, China), as previously reported [22]. The type strain C. jejuni ATCC 33560^T was employed as a control for comparison.

DNA extraction and genome sequencing

The genomic DNA for genome sequencing was extracted using the QIAamp DNA Mini Kit (Qiagen, German), following the manufacturer’s instructions. The concentration and purity of the extracted DNAs were measured using the NanoDrop spectrophotometer (Thermo Scientific, United States). The genome sequencing procedures were conducted at the Novogene Corporation (Beijing, China). For sequencing, a total of twenty-five draft genomes were sequenced using the Illumina HiSeq 2500Xten platform. Additionally, one complete genome was sequenced by Single Molecule Real-Time (SMRT) technology and assembled using the Canu v2.0 [23] to generate contig without gaps.

Genomic analysis

The genomes underwent gene prediction and functional annotation using the Prokka pipeline v1.14.6 [24]. To detect prophage sequences, the phiSpy v4.2.21 [25] was employed. Protein sequences annotated by Prokka were subsequently assigned to the Clusters of Orthologous Groups (COG) database using eggNOG-mapper v2 [26].

Antimicrobial resistance genes were predicted using the Comprehensive Antibiotic Resistance Database (CARD) v3.2.6 with Resistance Gene Identifier 6.0.1 [27] and the MEGARES Database v2.00 [28] using the Abricate tool v1.0.1, respectively. In addition, virulence genes across all genomes were detected using VFanalyzer [29]. To determine the average nucleotide identity (ANI) values, pyani v0.2.10 [30] was utilized. Core-pan genome analysis was conducted using the Roary pipeline 3.13.0 [31], which utilizes annotated assemblies in GFF3 format obtained from Prokka results.

Multilocus sequence typing (MLST), which relies on the variation among seven housekeeping alleles (adk, aspA, atpA, glnA, glyA, pgi, and tkt), was employed for both phylogenetic and epidemiological analyses. The assignment of sequence types (STs) and clonal complexes (CCs) was carried out utilizing the Campylobacter pubMLST database [32]. Additionally, any newly identified alleles and profiles were submitted to this database for inclusion.

Phylogenetic and phylogenomic analysis

Homologous protein clusters were extracted using the CD-HIT 4.8.1 [33], utilizing the.faa files from the Prokka results and applying an amino acid identity threshold of 40%. Among these clusters, representative sequences present in all analyzed genomes were identified as core protein sequences.

Core Single Nucleotide Polymorphism (SNP) calling was conducted using the Snippy 4.6.0, utilizing the C. upsaliensis type strain (CCUG 14913^T) as a reference. The Gubbins 2.4.1 [34] was employed as the recombination-removal tool, enabling the extraction of pure SNP data by eliminating recombination events.

The virulence-related genes were aligned to the analyzed genomes and the subject sequences were extracted from the genomes using a Perl script. These subject sequences were then translated into amino acid sequences.

The MAFFT 7.471 [35] was then employed to perform multiple sequence alignment of these sequences. Subsequently, the phylogenomic tree was constructed using FastTree 2.1.10 [36] with the maximum-likelihood (ML) algorithm. The phylogenomic tree was visualized using Dendroscope 3.8.3 [37] and further modified using the iTOL 6.7.2 [38] for enhanced visualization.

Statistical analysis and data visualization

Statistical analyses were performed using SPSS Statistics v27.0.1.0, and Chi-square tests were employed. A significance level of P < 0.05 was deemed statistically significant.

A variety of figures were generated using the R 4.2.2 [39] statistical language and environment. The design of bar charts, polyfrequency plots, and hierarchically clustered heatmaps was accomplished using the R packages ggplot2, ggvenn, RColorBrewer, reshape, reshape2, ggrepel, and dplyr.

Data  availability statement

In addition to the strains sequenced in this study, which have been deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/), other previously sequenced C. upsaliensis genomes were also obtained from the NCBI database. Comprehensive background information and biosample IDs for these genomes can be found in Supplementary Table S1.

Results

Isolation ratio from different animals

During the years 2019-2020, a total of 123 rectal swab samples were collected from dogs and cats undergoing medical examinations, while an additional 129 fecal samples were collected from captive sheep in Beijing, China. A total of 26 C. upsaliensis strains were successfully isolated from various animal sources in these samples. Specifically, 7 of these isolates (12.96%, 7/54) were obtained from cats, 18 (26.09%, 18/69) from dogs, and an additional isolate (0.78%, 1/129) was derived from sheep samples. The calculated isolation ratio from domestic pets was determined to be 20.33% (25/123).

This study encompassed the analysis of a comprehensive dataset consisting of 154 C. upsaliensis genomes. Of these, 26 genomes were sequenced within the scope of this study, while an additional 128 genomes were sourced from the NCBI database. The geographical distributions of these strains spanned across Asia, Europe, Africa, Oceania, and North America, with isolation dates ranging from 1979 to 2020. Notably, within this diverse distribution, the 112 strains from Asia, accounting for 72.73% (112/154) of the total, were exclusively obtained from the cities of Beijing and Shenzhen in China.

In terms of the source of these genomes, the majority were derived from dogs (58.44%, 90/154), followed by cats (14.29%, 22/154), and humans (11.69%, 18/154). A smaller proportion was isolated from coyotes (2.60%, 4/154) and sheep (0.64%, 1/154).

Genomic characteristics

Complete genome sequencing of strain XJK42-1 revealed the presence of a single chromosome and a circular plasmid. The chromosome exhibited a size of 1,644,360 bp with a GC content of 34.93%, while the plasmid measured 67,456 bp in size with a GC content of 29.57%. Genomic prediction and annotation results indicated that this strain encompassed 1,776 Coding DNA Sequences (CDS), along with 44 tRNAs, including 6 rRNAs (16S rRNA and 23S rRNA, with 3 copies of each), as well as a single tmRNA. Moreover, the strain harboured six prophages and notably lacked Type IV secretion system (T4SS).

Across the analyzed strains in this study, the genome sizes exhibited variation ranging from 1.53 Mb to 1.86 Mb each, while maintaining an average GC content of 34.71%. In terms of genomic content, C. upsaliensis genomes were comprised of 1,558–1,971 CDS, along with approximately 40.8 tRNA, 2.5 rRNA, and a single tmRNA.

All the new genomes sequenced in this study had 50 or fewer contigs, while parts of currently published genomes even had >100 contigs. The genomic characteristics of the complete genome of the XJK42-1 strain were consistent with all C. upsaliensis genomic characteristics and had the characteristics of C. upsaliensis [9]. The XJK42-1 strain held a 67 kb long plasmid which did not contain antimicrobial resistance genes. This plasmid was 96.51% identity to the 68 kb length plasmid (GenBank, CP059683.1) of human-isolated strain FDAARGOS_736 from clade C with 80% coverage.

Comparative genomic characteristics

Using the C. upsaliensis type strain (CCUG 14913^T) as a reference, the ANI values within each pair of strains exceeded 95%, the established threshold for species classification [30]. This confirmed their identification as C. upsaliensis. Based on the ANI values, the investigated strains were stratified into three discrete groups labelled group A, group B, and group C. Notably, the ANI values within each pair of strains greater than 97% were observed within each group, as depicted in Figure 1 (Table S2).

Figure 1. Heatmap of pairwise ANI values. The dendrogram reflects the degree of identity between genomes.

The orthologous groups of 1,061 homologous clusters (Table S3) shared by all analyzed strains were extracted and used to build the phylogenomic tree. According to the topology structure of the phylogenomic tree, these strains could be divided into three clades, named clade A, clade B, and clade C. The strains in each clade completely correspond to the ANI grouping. Clade A, B, and C contained 25, 14, and 115 strains, respectively. Except for 4 strains with unknown origins, all the other strains in clade A were isolated from animals. Clade C also contained most strains from animals (83.48%, 96/115), except only a tiny percentage were from the human origin (5.22%, 6/115). Interestingly, the strains in clade B were fully human origin except for two strains with unknown origin. The strains of animal origin isolated from China were all distributed in clade A and C. The strains in clade B were isolated from Europe and Africa, and none of them were from Asia or North America. In addition to strains of unknown origin, most of the strains were isolated from fecal samples (94.81%, 146/154), and three strains isolated from blood and intestinal samples were in clade C (Figure 2).

Figure 2. Phylogenomic tree based on core orthologous groups. The first ring represents the strains and clades, the second ring represents the isolated host, the third ring represents the isolated source, the fourth ring represents the isolated location, the fifth ring represents the ST, and the outside ring represents the CC.

A total of 93,610 core SNP loci were identified. The phylogenomic tree constructed by core SNP loci was also divided into three clades as well as that was constructed by the homologous clusters, and the strains in each clade in these two different phylogenomic trees were the same (Figure S1).

In this study, these strains were widely matched to 98 STs, and only six strains could match to 3 CCs. Altogether, 60.20% (59/98) of the STs were novel, while 98.12% (110/112) of Chinese isolates belonged to all these novel STs. ST-244 (9.74%, 15/154), ST-213 (4.55%, 7/154), ST-227 (4.55%, 7/154), and ST-242 (3.90%, 6/154) were the dominant STs (Figure 2).

Core-pan genome analysis

With the growing number of analyzed C. upsaliensis genomes, there was a sharp increase in the count of pan genes, while the pool of conserved genes showed a relative reduction. A total of 5,626 genes were contained in the pan-genome. In combination with the phylogenomic tree, which was constructed by the core homologous clusters, these three clades had significantly different gene family structures from each other, especially in clade B (Figure 3).

Figure 3. Pangenome analysis of C. upsaliensis genomes. Presence, blue; absence, white. Genomes were ordered based on a phylogenomic tree constructed by the core orthologous groups.

There were 1,157, 1,235, and 1,165 homologous protein clusters contained by clade A, B, and C, respectively. A total of 972 proteins were identified as core homologous protein clusters of these three clades. Clade A, B, and C contained 79, 82, and 16 clade-specific homologous protein clusters, respectively (Figure 4A). Clade C, encompassing the largest number of analyzed strains, exhibited the fewest specific homologous clusters. CDSs of these specific clusters were analyzed and assigned to COG by EggNOG. A total of 62, 66, and 11 proteins from clade A, B, and C, respectively, were assigned to COG. The largest proportion of proteins belonged to the S category, which represents “unknown function”. Clade C had no specific core clusters for eleven functions, which both belonged to the “CELLULAR PROCESSES AND SIGNALING” and “METABOLISM” classes of functional proteins. Clade B had significantly fewer proteins than clade A in “Nucleotide transport and metabolism”, “Cell motility” and “Secondary metabolites biosynthesis, transport and catabolism”. Only clade A had proteins in “Nucleotide transport and metabolism”. In “Transcription”, clade A had no proteins, but clade C held significantly more proteins than clade B. Similarly, in “Translation, ribosomal structure and biogenesis” and “Replication, recombination and repair”, clade C had significantly more numerous proteins than clade A and B. Clade B had fewer proteins in “Energy production and conversion” and “Amino acid transport and metabolism” than in other clades, but more in “Cell wall/membrane/envelope biogenesis” (Figure 4B).

Figure 4. The homologous clusters of C. upsaliensis. (A) Different numbers of core homologous clusters from different clades. (B) COG functional classification of core specific clusters in each clade. The digit indicates the number of clusters.

Antibiotic resistance

The 26 C. upsaliensis isolates in this study showed high MIC values for most of these selected antibiotics. The percentage of isolates whose MIC values ≥64 μg mL⁻¹ for nalidixic acid, gentamicin, and streptomycin exceeded 50% and even reached 80% for nalidixic acid. Isolates with MIC values ≥32 μg mL⁻¹ for ciprofloxacin and tetracycline also exceeded 50%, while those for erythromycin, azithromycin, telithromycin, and clindamycin were also approached. Notably, only chloramphenicol and florfenicol had relatively low MIC values. For these 11 antibiotics, all isolates had one or more antibiotics with MIC values ≥64 μg mL⁻¹. There were 88.46% (23/26) of the isolates whose MIC values were ≥64 μg mL⁻¹ for two or more antibiotics, while 96.15% (25/26) of the isolates with MIC values ≥32 μg mL⁻¹ (Table 1).

Table 1. MICs of isolated C. upsaliensis in this study to antimicrobials.

	Erythromycin	Azithromycin	Nalidixic acid	Ciprofloxacin	Gentamicin	Streptomycin	Chloramphenicol	Florfenicol	Tetracycline	Telithromycin	Clindamycin
XJK1-1	1	<0.5	64	32	>64	>64	8	4	32	4	<0.5
XJK12-1	>64	>64	64	1	<0.5	32	8	4	16	>32	>32
XJK25-1	>64	>64	64	32	<0.5	1	2	2	32	>32	32
XJK42-1	2	<0.5	64	32	>64	32	8	4	32	8	0.5
XJK54-1	4	<0.5	>64	32	>64	1	8	4	8	16	0.5
XJK56-1	2	<0.5	>64	16	>64	1	8	2	16	8	<0.25
YQ10-1	2	<0.5	64	64	<0.5	1	<0.5	4	32	4	1
SYD1-1	>64	>64	64	2	<0.5	1	2	2	32	>32	>32
SYD2-2	4	<0.5	64	32	>64	>64	16	4	16	16	<0.25
SYD7-1	1	1	4	4	>64	>64	2	2	4	2	0.5
SYDD1-1	<0.5	1	64	32	>64	>64	8	2	32	8	8
SYDD18-1	1	1	64	32	>64	>64	8	2	<0.5	8	1
SYDD19-1	1	>64	16	16	>64	2	4	1	32	>32	16
SYDD4-1	32	4	>64	32	1	2	8	4	16	>32	0.5
SYS8-2	>64	>64	32	8	>64	>64	8	4	32	>32	16
XJK23-1	2	<0.5	64	32	<0.5	64	16	4	32	16	<0.25
XJK24-1	<0.5	>64	16	8	>64	>64	4	2	8	>32	8
XJK3-1	>64	>64	32	8	>64	>64	8	2	16	>32	>32
XJK47-1	1	1	64	>64	<0.5	1	8	8	32	8	<0.25
XJK6-1	<0.5	<0.5	64	16	<0.5	2	8	4	16	2	<0.25
XJK63-1	2	<0.5	>64	16	<0.5	32	8	2	8	4	0.5
XJK69-1	<0.5	<0.5	64	64	>64	>64	4	2	16	4	0.5
XJK79-1	>64	>64	64	32	>64	>64	8	8	64	>32	32
YQ17-1	>64	>64	64	64	>64	64	16	4	32	>32	>32
YQ8-1	>64	>64	>64	32	>64	64	8	2	32	>32	>32
YQ9-1	>64	>64	64	32	>64	>64	8	2	32	>32	>32

Note: The unit of antibiotic concentration is μg/mL.

There totally seven antimicrobial resistance relative genes or mutations were identified in this study, namely 23S rRNA (A2075G), AAC(6”)-Ie-APH(2”)-Ia, APH(2”)-If, tet(O), gyrA (T86I), bla_OXA, and cmeABC. Among those, gyrA (T86I), bla_OXA, and cmeABC were found in all these twenty-six isolates. Among the isolates with MIC ≥32 μg mL⁻¹ for erythromycin, 50.00% (5/10) isolates had 23S rRNA (A2075G) mutation. There were 6 isolates containing 23S rRNA (A2075G) mutation but all with MIC ≤4 μg mL⁻¹ for erythromycin and azithromycin. All 26 isolates contained gyrA (T86I) mutations. Among them, 21 isolates (80.77%, 21/26) had MIC ≥64 μg mL⁻¹ for nalidixic acid and 23 isolates (88.46%, 23/26) had MIC ≥8 μg mL⁻¹ for ciprofloxacin. However, one isolate held MIC ≤4 μg mL⁻¹ for both nalidixic acid and ciprofloxacin. In isolates with MIC ≥64 μg mL⁻¹ for gentamicin, the APH(2”)-If were present in 41.18% (7/17) isolates, while AAC(6”)-Ie-APH(2”)-Ia genes were present in 29.41% (5/17) isolates. In isolates with MIC ≥32 μg mL⁻¹ for streptomycin, 35.29% (6/17) isolates contained AAC(6”)-Ie-APH(2”)-Ia gene. Although 80.77% (21/26) isolates with tetracycline MIC ≥16 μg mL⁻¹ were detected, only one isolate was detected with the tet(O) gene (Figure 5).

Figure 5. Heatmap of the distribution of antibiotic resistance relative genes or mutations. Orange indicates the presence and skyblue indicates the absence.

The XJK42-1 strain exhibited elevated MIC values for nalidixic acid, ciprofloxacin, gentamicin, streptomycin, and tetracycline. It encompassed four antimicrobial resistance-related genes or mutations, specifically 23S rRNA (A2075G), AAC(6”)-Ie-APH(2”)-Ia, gyrA (T86I), bla_OXA, and cmeABC. The resistance to quinolone and aminoglycoside antibiotics could be attributed to gyrA (T86I) and AAC(6”)-Ie-APH(2”)-Ia, respectively. However, the presence of 23S rRNA (A2075G) mutations in the genome did not lead to antimicrobial resistance for erythromycin (MIC, 2 μg mL⁻¹) and azithromycin (MIC, < 0.5 μg mL⁻¹), and the genome did not reveal resistance genes associated with tetracycline resistance (MIC, 32 μg mL⁻¹) (Figure 5).

More than 99% (99.35%, 153/154) of analyzed strains contained coding genes of bla_OXA beta-lactamases and CmeABC multidrug efflux pump. 70.78% (109/154) strains with gyrA (T86I) mutation that could be conferring resistance to fluoroquinolones. Notably, all gyrA (T86I) mutant strains were isolated from China, and the detection rate in China strains reached 97.32% (109/112). Similarly, other detected antibiotic-resistant genes or mutations were only found in isolates from China, and the detection rate of 23S rRNA (A2075G), AAC(6”)-Ie-APH(2”)-Ia, APH(2”)-If and tet(O) in isolates from China was 25.00% (28/112), 16.96% (19/112), 28.57% (32/112) and 15.18% (17/112), respectively (Figure 5).

Virulence genes

All the analyzed genomes were used to assess the presence of known virulence genes associated with Campylobacter adherence, immune modulation, invasion, motility, secretion system, and toxin. Adherence and toxin-related virulence genes were present in all strains. Among virulence genes related to immune modulation, glf was present in 55.84% (86/154) strains from all strains of clade B and other clades. Besides, Cj1138 associated with LOS was present in 35.06% (54/154) strains, and none of the strains in clade B contained this gene, while neuB1 and neuC1 were absent in clade B strains. The invasion-related virulence gene ciaA was present in all strains, while ciaB was in 55.19% (85/154) strains. In motility-related virulence genes, clade A and B strains contained pseD but not pseE. The T4SS existed in 21.43% (33/154) strains (Figure 6).

Figure 6. Heatmap of the distribution of virulence genes. Orange indicates the presence, and skyblue indicates the absence.

No unique virulence-related gene was identified in clade B. Using the virulence gene’s sequences of strain XJK420-1 as query sequences, a total of 87 complete genes were acquired from these strains (Table S4). Both these genes were found in over 95% of the strains and could be accurately translated. The results of these gene alignments revealed that flgD, gmhA, and cdtC exhibited the highest specific variability among strains in clade B, all exceeding 6%. The flgD gene product had a total length of 293 amino acids, and the entire strains in clade B displayed 11 specific amino acid variation sites and 14 amino acid deletions (AA31, AA275-287, Figure 7). The gmhA gene product, spanning 194 amino acids, exhibited 13 specific amino acid variation sites in strains from clade B. In addition, the cdtC gene product, comprising 189 amino acids, displayed 12 specific amino acid mutation sites in strains from clade B. Aligning the amino acid sequences, the mutation sites specific for the strains of clade B were most in the signal region in CdtC (8/12). Unlike CdtC, most of the specific amino acids were in the N-terminal, and the amino acid differences in both FlgD and GmhA between clade B and clade A/C strains were distributed throughout the protein (Figure 7). Meanwhile, the phylogenetic tree constructed using the cdtC, flgD, and gmhA sequences showed that strains in clade B still robustly clustered together, while some strains in clade C clustered together with strains in clade A (Figure S2, Figure S3, Figure S4).

Figure 7. Alignment of amino acids of FlgD, GmhA, and CdtC. (A) FlgD, (B) GmhA, (C) CdtC.

Discussion

In this study, we conducted a comprehensive comparative genomic analysis of 154 C. upsaliensis strains based on a wide geographic distribution and diverse sample sources from humans, animals, and others. Little studies have described C. upsaliensis isolates so comprehensively until now [9,14,15,40].

Among these analyzed strains, 26 C. upsaliensis were isolated in this research. A total of 112 strains (72.72%, 112/154) of C. upsaliensis were isolated using the same isolation method from Beijing and Shenzhen, China [9]. The main isolation sources of these strains were dogs and cats. An increasing number of C. upsaliensis strains have been isolated and cultured, indicating that Campylobacter's isolation and culture technology is becoming more sophisticated [41].

In this study, more than half of the ST types were novel. As the genomic gold standard for prokaryotic species definition [42], all the analyzed strains met the ANI standard of the same species and could be divided into three groups. Meanwhile, evolutionary analysis based on the core homologous clusters and the core SNPs found that strains could be included robustly in three clades and that strains in each clade were the same as those in the group formed using ANI. The ST types and genomic evolutionary analysis suggested a high genomic diversity within C. upsaliensis. The strains isolated from Chinese animals were distributed in clade A and C, while those in clade B were mainly isolated from humans (P < 0.01), and this might suggest that these strains isolated from animals in China might not possess the genetic characteristics of human colonization. Similarly, C. concisus consisted of two genetic species, GS1 and GS2, and GS2 strains were better adapted to the human gastrointestinal tract [43,44]. C. lari strains were distributed in several subspecies or groups of members, i.e. Cll, Clc, C. lari UPTC, and the other four main clades, and a significantly different distribution of strains to subspecies levels within the sources were highlighted [45]. C. ureolyticus was a complex of several Campylobacter species [19].

C. upsaliensis isolated and cultured in this study had high MIC values to antibiotics of quinolones, aminoglycosides, and tetracycline and had relatively serious antimicrobial resistance to macrolides, telithromycin, and clindamycin. These results were consistent with Shenzhen's, but the antimicrobial resistance to tetracycline was more serious. Moreover, the strains isolated from cattle in South Africa were more resistant to erythromycin, and multidrug resistance (MDR) was serious in different regions and sources [9,17]. No correlation found between the distribution of antibiotic-resistant genes and clades. The resistance genes mediating antibiotics of β-lactam and quinolone were detected in all used genomes. In addition, other resistance genes were only found in Chinese strains, indicating that the situation of antimicrobial resistance in China was serious and antimicrobial resistance gene transfer was common.

Except for two immune evasion-related genes specifically absent in clade B strains, there was a limited correlation between the presence of virulence genes and clades, similar to antimicrobial resistance genes. These results differed from C. lari, whose distribution was different in different subspecies [45]. However, sequence specifics in FlgD, GmhA, and CdtC were identified in clade B strains.

The assembly of the bacterial flagellar hook requires FlgD, a protein responsible for forming the hook cap. This protein might play key roles in reducing biofilm production, influencing motility and virulence, and retaining the ability to form a flagellar bundle behind a cell body [46,47]. The LPS inner-core biosynthesis protein GmhA supported their potential as novel vaccine candidate antigens, and its inactivation resulted in significant changes in the redox status of cells [48,49]. The specific amino acids in clade B strains might enable specific pathogenic functions of these two proteins.

In the cytolethal distending toxin, CdtA and CdtC constitute regulatory subunits, while CdtB acts as the catalytic subunit with phosphatase and DNase activities, leading to cell cycle arrest and cell death. The phylogenetic tree based on the cdtC sequence showed that clade B strains were consistently clustered together. By alignment of CdtC, we found that the mutation sites of CdtC in clade B were specific but conserved in clades A and C. CdtB is the active subunit, whereas CdtA and CdtC are required to bind and deliver CdtB within the target cell. CdtA/CdtC bound to the lipid rafts and cholesterol-rich microdomains within the host cell membranes [50], and there was also a report that CdtA was not essential in this process [51]. Inactivation of CdtC confirmed its role in establishing infection and ameliorates H. pylori-induced pathogenesis by hijacking cholesterol [52,53]. Interestingly, C. upsaliensis was reported to produce higher CDT titre than other Campylobacter [40]. Whether the type of CdtC in clade B determines CDT's higher affinity and pathogenicity to humans remained to be further studied and explored. In a way, however, the specific amino acids of CdtC in the clade B which was highly associated with human infection might be potential pathogenic factors in C. upsaliensis. In any case, further investigations were needed to explore these potential pathogenic factors.

Conclusions

In this study, a total of 154 available genomes were collected for genetic characterization and comparative genomic analysis. All analyzed strains can be divided into three clades, of which nearly all strains of Clade B were isolated from humans, and the CdtC, FlgD, and GmhA might serve as potential pathogenic factors for the Clade B strains associated with human infections. Meanwhile, strains isolated from pets in China carrying more drug-resistant elements needed to be paid more attention for antimicrobial resistance control.

Acknowledgments

We extend our gratitude to Academician George Fu Gao, Chinese Academy of Sciences, for his invaluable support for the project formulation and sample collections. We thank Dr. Lang Liu from the New Ruipeng Pet Healthcare Group for his kind help with the pet’s sample collections. We also should thank our colleagues from Nanshan Centre for Disease Control and Prevention, and Shunyi District Centre for Disease Control and Prevention.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Sponsored by the National Key Research and Development Program of China [grant no 2021YFC2301000], the Campylobacter standard strains applied to “National standards for food safety” [grant no TS202302], the Project for Research on detection and tracing of important pathogen combinations, identification and prevention of drug resistance hazards (29140) and Project for Gonorrhea vaccine development (33076).