Comparative genomics of Helicobacter pullorum from different countries

Background Helicobacter pullorum commonly colonized in the gastrointestinal tract of poultry and caused gastroenteritis. This bacterium could be transmitted to humans through contaminated food and caused colitis and hepatitis. Currently, the genetic characteristics of the H. pullorum were not recognized enough. In this study, the genomes of 23 H. pullorum strains from different counties were comparatively analyzed. Among them, H. pullorum 2013BJHL was the first isolated and reported in China. Results The genomes of the studied strains were estimated to vary from 1.55 to 2.03 Mb, with a GC content of ~ 34%. 4064 pan genes and 1267 core genes were obtained from the core-pan genome analysis using the Roary pipeline. Core genome SNPs (cg-SNPs) were obtained using Snippy4 software. Two groups were identified with the phylogenetic analysis based on the cg-SNPs. Some adhesion-related, immune regulation, motility-related, antiphagocytosis-related, toxin-related and quorum sensing related genes were identified as virulence factors. APH(3′)-IIIa, APH(2′’)-If, and AAC(6′)-Ie-APH(2′’)-Ia were identified as antibiotic resistance genes among the H. pullorum genomes. cat, SAT-4 and tetO genes were only identified in 2013BJHL, and tet(C) was identified in MIT98-5489. MIC determination revealed that the 2013BJHL showed acquired resistance to ciprofloxacin, nalidixic acid, tetracycline, gentamicin, streptomycin and erythromycin, only sensitive to ampicillin. The antibiotic resistance genetic determinants on the 2013BJHL genome correlate well with observed antimicrobial susceptibility patterns. Two types of VI secretion system (T6SS) were identified in 52.2% (12/23) the studied strains. Conclusion In this study, we obtained the genetic characteristics of H. pullorum from different sources in the world. The comprehensive genetic characteristics of H. pullorum were first described. H. pullorum showed highly genetic diversity and two sub-types of T6SSs were first identified in H. pullorum. 2013BJHL was found to be multidrug resistant as it was resistant to at least three different antibiotic classes.


Background
Helicobacter pullorum (H. pullorum) was a kind of Gram-negative and urease-negative bacterium, which initially isolated from the liver and intestinal of asymptomatic poultry [21]. The H. pullorum could be transmitted to humans through contaminated meat and it had been recognized as human pathogen which was associated with human colitis and hepatitis [5,21,22]. In these days, human recurrent diarrheal, Crohn's disease and

Open Access
Gut Pathogens  [3,12,21,23,24]. With the increased consumption of the chicken, H. pullorum might be one important emerging food-borne pathogen in the world. As the same time, isolates of poultry origin showed resistance to ciprofloxacin, gentamycin, erythromycin and tetracycline and was susceptible to colistin sulfate and ampicillin, and phenotypic resistance was almost the same as genotype resistance [4,8,9], however, there was no antibiotic recommendation for this organism. There were reports about the prevalence of H. pullorum in poultry but the positive ratios were variant (from 23.5 to 100%), which might be due to differences in detection methods and the isolation samples [4,8,9,14,25]. The isolation in poultry caecum seemed higher than in liver [1,6]. Currently, only few of genetic characteristics of H. pullorum were described, and showed diversity or similarity in different sources [4,18], but it was unclear in China. In this study, we attempted to characterize the H. pullorum strain isolated from retail chicken liver in China in the aspect of core-pan genome, the antibiotic resistance potential, virulence factor genes and population structure based on phylogeny. The comparative genomics of 22 H. pullorum strains from different sources and different countries released from NCBI to date also were investigated.

Genomic and antimicrobial susceptibility characteristics of strain 2013BJHL
The genome size of strain 2013BJHL was 1,875,659 bp, 34 scaffolds and GC content was 34.30%. 1939 genes were predicted and the total length of genes was 1,705,674 bp, which made up 90.94% of genome. A 7.3 kb incomplete phage-associated region was identified using online software PHASTER, which consisted of genes encoding metabolism, biosynthesis and hypothetical proteins. MIC determination revealed that 2013BJHL was resistant to ciprofloxacin, nalidixic acid, tetracycline, gentamicin, streptomycin and erythromycin but only sensitive to ampicillin.

Comparative genomic and genetic population structure
The genomes of the studied strains were estimated to vary from 1.55 Mb to 2.03 Mb for each, with a GC content of ~ 34%. 4064 pan genes and 1267 core genes were obtained from these 23 H. pullorum genomes. From the pan-genome analysis, we found that the strain (ID:229334/12) from fresh chicken meat in Portugal contained the most specific genes (191 genes), followed by strain NCTC12824 isolated from chicken faeces in Switzerland (100 genes). However, human strain MIT98-5489 had the fewest specific genes (0 gene). 76 specific genes were detected in strain 2013BLHL and 71.1% (54/76) of them were predicted to be hypothetical proteins. The distribution of the specific genes in different strains was shown in Fig. 1. We also found that free-range chickens had more specific genes than broilers in Indian. The individual genes for each strain were presented in the Additional file 1.
According to the isolation sources, these 23 strains were classified into three groups. These three groups share 1821 genes, and the human source group had the least number of specific genes (169 genes, Fig. 2). The specific genes in each group were present in the Additional file 2. Most of the human sources known functions specific genes were transferases associated genes.
A phylogenetic tree was built based on the cg-SNPs (Fig. 3). Two groups were identified. Most of the strains in Group2 were isolated from European. However, in Group1, the strains were isolated from Indian, China, Portugal and Canada.

Virulence, antibiotic resistance and secretion system analysis
Numerous of virulence associated genes were identified in this study: the adhesion-related genes(wbp), immune regulation genes(napA), motility-related genes(fla), toxin-related genes(cdt), antiphagocytosis-related genes(wec) and quorum sensing genes(luxS). Some of these genes were homologous with the genes in C. jejuni, like immune evasion, motility related genes and glycosylation system related genes. The details of the virulence genes were listed in the Additional file 3 and presented in Fig. 4.
T6SS was identified in 52.2% (12/23) the studied strains (see Table 2 for detail). According to the identified gene arrangement, these T6SS could be classified into two types. Type I including 15-18 genes and had 18,520-20,025 bp long, Type II contained 19-22 genes and 20,056-22,572 bp long. The gene contents of each T6SS type were listed in Table 1. The gene arrangements and gene components of these two T6SS types were presented in Fig. 5.

Discussion
In this study, based on the comparative genomic analysis, the comprehensive characteristics of these H. pullorum genomes were obtained. The results of core-pan genomes analysis showed highly intraspecies genomic diversity in H. pullorum. Human strain MIT98-5489 did not contain any specific genes. However, another human strain NCTC13156 from Canada had 99 specific genes, which also indicated that H. pullorum had a high genetic diversity. For three different groups, we found these three groups share 1821 genes, and the human source group contained less number of specific genes. Since the number of strains from each group were different, we cannot exclude the bias caused by the number of strains. The cas1_2 gene (CRISPR-associated endonuclease Cas1) was the unique gene in two strains from human origin. In addition, there were 6 CRISPR-related genes in the pan-genome and found in all strains of human origin. The study of Marraffini et al. [15] showed that CRISPR loci counteracted multiple routes of HGT and could limit the spread of antibiotic resistance in pathogenic bacteria. Therefore, this might have an impact on the horizontal transfer of resistance genes.
2013BJHL was the first reported strain isolated from chicken liver in China. 76 specific genes were detected in this strain and four of them were related to antibiotic resistance. More drug resistance genes were found in the strains from Asia. For antimicrobial susceptibility testing, 2013BJHL showed the characteristics of multi-drug resistance. And antibiotic resistance genetic determinants on the 2013BJHL genome correlate well with observed antimicrobial susceptibility patterns. Resistance tests on 11 strains of H. pullorum in India have shown that all strains were multi-drug resistance [18]. Interestingly, in this study, the strains without individual missense mutations in the gyrA gene were isolated from European. We speculated that this might be related to the abuse of fluoroquinolone antibiotics in Asian settings. Moreover, all the antibiotic resistance genes in broiler chickens were AAC(6′)-Ie-APH(2′′)-Ia, and all the antibiotic resistance genes in free range chickens were APH(3′)-IIIa, and the number of antibiotic resistance genes in broiler chickens were more than that in free range chickens. This might be because in intensive farming, antibiotics were routinely fed to livestock as growth promoters and to prevent potential bacterial infections, which have contributed to increase in drug resistance worldwide, enabling re-emergence of zoonotic infections. At the same time, strains from different sources might have different patterns of drug resistance. These results also reflected the resistance situation of H. pullorum in Asia was more serious. All four strains were resistant to ciprofloxacin and two strains were resistant to erythromycin and tetracycline    respectively in Borges' study [4] and the perfect correlation between antimicrobial phenotypes and genotypes was found, which was also confirmed in our study. But our results were similar to those in India, showed characteristics of multi-drug resistance. Core genome analysis played a key role in determining population structure which in turn shed light on the evolutionary trajectories of the strains [13]. The phylogenetic tree (Fig. 3) indicated the Asian strains were concentrated in Group1 and the Europe strains were concentrated in Group2. In addition, no drug resistance gene or gene individual missense mutation in the gyrA gene were detected in Group2. The genetic population structure of H. pullorum might be influenced by geographical distribution but due to the limitation of the analyzed strains the ecological niche characteristics could not be well inferred.
There was previous study reported that H. pullorum was capable of host cell adhesion, and the adhesion rate was comparable to C. jejuni [20]. Some homologous genes in C. jejuni were also found in H. pullorum. T6SS played an important role in the virulence of pathogens such as Vibrio cholera O1, Salmonella and C. jejuni [2,10,16]. In this study two subtypes of T6SS were identified. The existence of T6SS-related genes in all strains were consistent. The tssA gene was truncated in strain MIT98-5489, which might be due to the quality of the sequences, as there were lots of gaps in the regions of this strain. Although the interspecies region of these two subtypes of T6SS were not totally identical, the constituent genes showed highly similarities with the T6SS in C. jejuni [4]. The pathogenetic characteristics of these two subtypes of the T6SSs need further study.

Conclusion
In a summary, this study we obtained the genetic characteristics of H. pullorum from different sources in the world. H. pullorum showed highly genetic diversity and two sub-types of T6SSs were first identified in H. pullorum.

Sampling, isolation and identification of strain
2013BJHL was isolated from broiler retail chicken liver from Beijing, China in 2013. The isolate was grown on Campylobacter agar (OXIOD, UK) with 5% sheep blood in a microaerobic atmosphere (85% N 2 , 5% O 2 and 10% CO 2 ). After culturing, the DNA for genome sequence was extracted using the QIAamp DNA Mini Kit (Qiagen, German) according to the manufacturer's instructions for sequencing. Then the NanoDrop spectrophotometer was used to measure concentration and purity of DNAs.  The quality requirements were concentration ≥ 50 ng/μL and total amount > 20 μg. The purity requirement was as follows: OD260/OD280 value should be between 1.6 and 1.8.

Antimicrobial susceptibility testing
Antimicrobial susceptibility to seven antibiotics (ciprofloxacin, nalidixic acid, tetracycline, gentamicin, streptomycin, erythromycin and ampicillin) was determined using the gradient strip diffusion method (E-test ™ , bio Mérieux, Nürtingen, Germany) following the manufacturer's instruction. The bacterial suspension for the E-test was adjusted to 2 McFarland in 0.85% normal saline. 75 μl was evenly spread on a Campylobacter agar supplemented with 5% defibrinated sheep blood and a single strip was put on each plate. After 72 h of incubation at 42 °C under a microaerobic conditions, the minimum inhibitory concentration (MIC) was determined. The type strain of ATCC33560 was used as control. Thereafter, both the results were interpreted according to the Clinical and Laboratory Standards Institute (CLSI) recommendations.

Whole genome sequencing (WGS) for strain 2013BJHL
The DNA sequencing was performed by an Illumina HiSeq2500Xten platform (Illumina Inc., San Diego, CA, USA) at the Beijing Genomics Institute (BGI) with a depth of 500× coverage. To sequence the genomes, a 270 bp paired-end library was constructed and then 150 bp reads were generated. FastQC (http://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) and fastp (https ://githu b.com/OpenG ene/fastp ) software tools were applied to evaluate and improve the quality of the raw sequence data, respectively. Low quality reads were removed if the quality scores of ≥ 3 consecutive bases were ≤ Q30. The clean reads were assembled by SOAPdenovo (http://soap.genom ics.org.cn/soapd enovo .html). The assembled sequences were predicted to genes and annotated the function using Prokka pipeline [19] and glimmer software (http://ccb.jhu.edu/softw are/glimm er/index .shtml ). Phage Search Tool (PHAST) web server(http://phast er.ca/) was used to search for phage sequences. The genome of strain 2013BJHL was submitted to NCBI, and the accession numbers was assigned as JXTX01000000.

Core-pan genome analysis
Core-pan genome analysis was deduced using the Roary pipeline [17] with the.gff files from Prokka results. The parameters were chosen "-e (create a multi-FASTA alignment of core genes using PRANK), -n (fast core gene alignment with MAFFT) and -v (verbose output to STDOUT)". The genes contained in all strains were called core genes, and the genes contained in only one strain were called specific genes. The 3D map of the specific genes ( Fig. 1) was drawn by using the scatter-plot3d package. Then, the 23 strains were classified into three groups (human:2 strains, chicken feces:17 strains and chicken_meat:4 strains) according to the isolation source, and core-pan analysis was performed. Genes contained in all three groups were called core genes, and those that existed only in one group were considered specific genes. Venn diagram (Fig. 2) was drawn by using VennDiagram package to show the overlapping area of different element sets.

Phylogenetic analysis
Core SNPs were called using Snippy4.3.6 software (https ://githu b.com/tseem ann/snipp y) from the same reference (the complete genome NCTC13154). Gubbins software [7] was used as the recombination-removal tool to gain the pure SNPs without recombination. Phylogeny reconstruction was performed with the Maximum likelihood (ML) method using MEGA 7 software [11] with 1000 bootstraps.

Virulence factors, antibiotic resistance genes and secretion system analysis
The virulence genes of all the genomes were detected on VFanalyzer (http://www.mgc.ac.cn/cgi-bin/VFs/ v5/main.cgi?func=VFana lyzer ). Resistance genes were predicted using the ResFinder (https ://cge.cbs.dtu.dk/ servi ces/ResFi nder/) and Comprehensive Antibiotic Resistance Database (CARD) (https ://card.mcmas ter. ca/?q=CARD/ontol ogy/35506 ), with an E-value of at least 1e−10 as the cutoff. The identity cut-off and query coverage values were kept > 80% and > 60%, respectively. Blastn was used to detect individual missense mutation in gyrA gene, which was responsible for conferring ciprofloxacin resistance. The protein sequences were functionally annotated and classified using databases such as KEGG, COG, SwissProt and PHI et al. The nucleotide sequences of all the annotated T6SS-related genes were extracted from genomes, and related gene clusters were visualized by genoPlotR package. Heatmap of virulence genes was generated using the pheatmap package.