Evolution in an oncogenic bacterial species with extreme genome plasticity: Helicobacter pyloriEast Asian genomes

Japanese/Korean core genomes diverged from the European and then the Amerind

A phylogenetic tree was constructed from concatenated seven genes atpA, efp, mutY, ppa, trpC, ureI and yphC, which were used for multi-locus sequence typing (MLST) [18] and phylogenetic analyses [19, 20]) (Additional file 1 (= Figure S1)). The tree showed that the 6 East Asian strains, the 4 Japanese strains (F57, F32, F30 and F16) and the 2 Korean strains (strain 51 and strain 52), are close to the known subpopulation hspEAsia of hpEastAsia, whereas 4 strains (Shi470 [21], v225d [22], Sat464 and Cuz20) are close to another subpopulation of hpEastAsia, hspAmerind. Strains 26695, HPAG1, G27, P12, B38, B8 and SJM180 were assigned to hpEurope. Strains J99 and 908 were assigned to hspWAfrica of hpAfrica1. PeCan4 was tentatively assigned to hspAmerind although it appears to be separate from the above 4 hspAmerind strains and somewhat closer to other subgroups (a subgroup of hpEurope, hspMaori and a group of "unclassified Asia" in the HpyMLST database [18]).

We deduced the common core genome structure of these 20 genomes based on the conservation of gene order using CoreAligner [23] (Table 1). CoreAligner determines the set of core genes among the related genomes not by universal conservation of genes but by conservation of neighborhood relationships between orthologous gene pairs allowing some exceptions. As a result, CoreAligner identified different numbers of core genes among strains (1364-1424), which reflect deletion, duplication and split of the core genes in the individual strains.

For phylogenetic analysis among the strains, we further extracted 1079 well-defined core orthologous groups (OGs) as those that were universally conserved, non-domain-separated, and with one-to-one correspondence (see Methods). The concatenated sequence of all well-defined core OGs resulted in a well-resolved phylogenetic tree (Figure 1). The tree was composed of two clusters, one containing the Japanese, Korean and Amerind strains and the other containing the European and West African strains. The tree strongly supported a model in which the Japanese/Korean strains (hspEAsia) and the Amerind strains (hspAmerind) diverged from their common ancestor, which in turn diverged from the ancestor shared by the European strains (hpEurope) long before. This conclusion is robust, as shown by the high bootstrap values of the internal nodes, primarily because the tree is composed of a large quantity of sequence information with approximately 1400 genes. The Japanese and Korean strains were not separated into two clusters. PeCan4 appeared diverged from the other four hspAmerind strains as expected from the result of the phylogenetic analysis based on the 7 genes described above. SJM180 appeared diverged from the other hpEurope strains in the well-defined core gene-based tree.

Phylogenetic profiling to identify gene contents of hspEAsia

To thoroughly characterize the gene contents specific to the Japanese/Korean (hspEAsia) strains, we conducted phylogenetic profile analysis using the DomClust program [24]. This analysis determines the presence or absence of a domain, rather than a gene, and allows detection of split genes, partially deleted genes and partially duplicated genes (detailed in Methods). Their features will be explained in the next five sections.

Differences in outer membrane proteins and related proteins in the number of loci of gene families and in alleles at each locus

The notable exception was oipA, for which a secondary locus was found in hspEAsia (6/6 strains) and hspAmerind (5/5), but not in hpEurope (0/7) or hspWAfrica (0/2). This increase of the secondary locus can be explained by a novel DNA duplication mechanism associated with inversion [25]. The two hopMN loci in hpEurope (7/7 strains) and hspWAfrica (1/2) were reduced to one locus in the hspEAsia (6/6) and hspAmerind (5/5). This loss was likely caused by the same duplication mechanism [25].

For the babABC family, the babC locus [26] was empty in all the hpEastAsia strains (6/6 hspEAsia and 5/5 hspAmerind) as well as from all the hspWAfrica strains (2/2) and two hpEurope strains (B38 and B8). This is in contrast to the presence of three loci in the other (5/7) European strains (Table 2).

The strain J99 carried a sabA gene (jhp0662) at the sabA locus and a sabB gene (jhp0659) at the sabB locus [27]. All the hpEurope strains but the strain B38 (6/7) and this hspWAfrica strain (J99) had these two loci, whereas all the hpEastAsia strains but the strains 52 and PeCan4 (5/6 hspEAsia and 4/5 hspAmerind) lacked sabB locus (Table 2). These hpEastAsia strains all carried a sabA gene at the sabA locus. Genes of hpEurope differed among strains. Three strains (HPAG1, G27 and SJM180) carried a sabA gene at the sabA locus and a sabB gene at the sabB locus, as J99. The strain 26695 carried a sabA gene at both the sabA and sabB loci, whereas the strain P12 carried a sabB gene at both the loci. The strain B8 carried a sabA gene at the sabA locus and a hopQ gene at the sabB locus, along with another hopQ gene at the hopQ locus.

Some of these genes (oipA, babA and babB) and homAB genes were previously reported to diverge between the East Asian and Western strains [13, 14, 17]. Difference in the number of copies of homAB genes between East Asian and Western strains was reported [17].

For hopMN, two gene types (hopM and hopN) have been recognized [26, 27]. Phylogenetic network analysis revealed two variable regions within the hopMN family (region II and IV; Figure 2). Combining the two types of two variable regions defined four main gene types, of which two corresponded to hopM and hopN. The two types in region II were designated m1 and m2 (m for mid). The types in region IV were designated c1 and c2 (c for C-terminus); c3 was another variant type in region IV, composed of parts of c1 and c2. In this designation, previous hopM and hopN genes correspond to hopMNm1-c1 and hopMNm2-c1, respectively. All hpEastAsia strains except the strains 52 and PeCan4 (9/11) carry sequence type c2 at region IV. The c3 variant is observed in J99, PeCan4 and SJM180 (Figure 2A and 2F).

Three vacA paralogs and vacA itself were found in 26695 [28]. Those paralogs share the auto-transporter domain at the C-terminus with vacA [28]. A large deletion in vacA-2 (HP0289) (approximately 2400 amino acids) was found in all the hspEAsia strains except the strain 51 (5/6) (Table 2 and Additional file 2 (= Table S1)).

It was described earlier that horA OMP locus in 26695 is composed of two open reading frames (ORFs) (HP0078/HP0079) whereas that in J99 is composed of one ORF (jhp0073) [27]. The horA locus in all the hspEAsia strains shows apparent gene decay by fragmentation through various mutations (Figure 3). Whether the genes in the other strains are functional is not known.

A putative periplasmic endonuclease gene (nucG, HP1382) was split in all the hspEAsia strains examined (Table 2 and Additional file 2 (= Table S1)). Detailed analysis revealed that the split was mediated by recombination between short similar sequences [25].

Massive decay of molybdenum-related genes for two-electron reduction-oxidation reactions

In the 20 H. pylori genomes, the only gene for molybdo-enzymes identified was bisC. At least one gene in each of the three Mo-related functions, Mo transport, Mo cofactor synthesis and a Mo-containing enzyme, decayed in all hspEAsia strains (Table 3 and Figure 4). Detailed analysis of nucleotide sequences revealed a mutation in 10 of 12 Mo-related genes in some of the hspEAsia strains (Table 3 and Additional file 3 (= Table S2)). The occurrence of apparently independent multiple mutations (Additional file 3 (= Table S2)) suggests some selection against use of Mo in the hspEAsia strains. All other strains but P12 possessed all intact genes. The strain P12 had a truncation of moaD (Additional file 3 (= Table S2)). Tungsten sometimes substitutes for Mo, but genes for known tungstate/molybdate binding proteins (TupA and WtpA) were not found in the H. pylori genomes.

The sequences in the four Japanese strains were confirmed by polymerase chain reaction (PCR) with the primers listed in the Additional file 4 (= Table S3).

The Mo-related genes were in a list of "chronic gastritis-associated" genes [30], primarily because they are absent from three Amerind strains from the Athabaskan people [31]. The 5 Amerind strains analyzed in the present study are different from the three Amerind strains in this respect. This difference could reflect the later migration of the Athabaskans to the Americas [32].

Two pathways between acetyl~CoA and acetate in some Japanese strains

Our profiling revealed an important change at the center of energy and carbon metabolism related to acetyl~CoA. Two pathways connect acetyl~CoA and acetate (Figure 5A). In anaerobic fermentation, acetyl~CoA is converted into acetate by phosphoacetyl transferase (pta product) and acetyl kinase (ackA product) with generation of ATP (anaerobic pta-ackA pathway) [33]. The intermediate acetyl~P, a high-energy form of phosphate, likely serves as a global signal. Although these reactions are reversible, assimilation of acetate may be irreversibly mediated by acetyl~CoA synthetase (acoE product) by the generation of acetyl~CoA, which enters the TCA cycle to generate energy under aerobic conditions (aerobic acoE pathway).

It has been suggested that strain 26695 (hpEurope) carries a mutation in pta for the former pathway whereas strain J99 (hspWAfrica) lacks acoE for the latter [28, 34]. All European strains in this study (7/7) had at least one inactivated pta and ackA gene through a variety of mutations (Figure 5C). Two of five Amerind strains, PeCan4 and Cuz20, also had a mutated pta and ackA, whereas the other 3/5 Amerind, 2/2 African, and 3/6 hspEAsia strains had a pta and ackA intact but had a deletion of acoE. Exceptions to such apparent incompatibility between the two pathways were found for 3/4 of the Japanese strains (F16, F30 and F57), which had intact genes for both pathways (Figure 5BCD). The sequences in the four Japanese strains were confirmed (see Methods and Additional file 4 (= Table S3)).

A gene for an amino acid utilization

An ortholog of jhp0585 in J99 is absent from 26695 [2]. An ortholog is present in the six other hpEurope strains and both hspWAfrica strains, but absent from all hpEastAsia strains (hspEAsia and hspAmerind) (Additional file 2 (= Table S1)). It encodes a homolog of 3-hydroxy-isobutyrate dehydrogenase and the related beta-hydroxyacid dehydrogenase (COG2084). The 3-hydroxy-isobutyrate dehydrogenase degrades the branched-chain amino acid valine. H. pylori requires branched amino acids for growth. It is not known what the substrates or products of reactions catalyzed by this gene product are, or the biological relevance of its distribution.

Gene contents unique to other groups

Prophage-related genomic islands and other mobile elements

Except for the cag pathogenicity island (cagPAI), five genomic islands (GIs) were identified in the genomes of the four Japanese strains (Table 4, Figure 6 and Figure 7). In F32, the cagPAI was flanked by a 44-bp direct repeat, which extended the 22-bp sequence found in the other strains (Table 4). This length of sequence identity would allow homologous recombination [47] leading to the excision of cagPAI flanked by the repeat.

Strain	GI number	Type	Length	Start-end	Flanking repeat (bp)	Secretion system	Left gene (annotation)	Right gene (annotation)
F16	GI_HP_F16_1	prophage-like	12245	471964 - 484208 (HPF16_0465 - HPF16_0478)	N/D(a)	N/D(a)	HPF16_0464 (IS605 transposase (tnpB))	HPF16_0479 (typeIIR)
	GI_HP_F16_2	cagPAI	36761	871413 - 834651 (HPF16_0834 - HPF16_0810)	22(b)	Type IV	HPF16_0835 (hypothetical protein)	HPF16_0809 (glutamate racemase)
F30	GI_HP_F30_1 (left)	type 1b TnPZ partial	7246	1280406 - 1287651 (HPF30_1205 - HPF30_1211)	N/D(a)	tfs3b partial	HPF30_1204 (outer membrane protein horB)	HPF30_1212 (rpoD)
	GI_HP_F30_1 (right)	type 1b TnPZ partial	1655	1237267 - 1238921 (HPF30_1166 - HPF30_1167)	N/D(a)	N/D(a)	HPF30_1165 (hypothetical protein)	HPF30_1168 (5'-methylthioadenosine/S -adenosylhomocysteine nucleosidase)
	GI_HP_F30_2	cagPAI	37153	867993 - 830839 (HPF30_0803 - HPF30_0778)	22(c)	Type IV	HPF30_0804 (hypothetical protein)	HPF30_0777 (glutamate racemase)
F32	GI_HP_F32_1	type 2 TnPZ partial	24283	1058236 - 1082518 (HPF32_0988 - HPF32_1014)	N/D(a)	tfs3 partial	HPF32_0987 (hypothetical protein)	HPF32_1015 (hypothetical protein)
	GI_HP_F32_2	cagPAI	36609	534488 - 571096 (HPF32_0500 - HPF32_0524)	44(d)	Type IV	HPF32_0499 (hypothetical protein)	HPF32_0525 (glutamate racemase)
F57	GI_HP_F57_1 (left)	type 1b TnPZ partial	7246	103791 - 111036 (HPF57_0102 - HPF57_0109)	N/D(a)	tfs3b partial	HPF57_0101 (RNA polymerase sigma factor RpoD)	HPF57_0110 (hypothetical protein)
	GI_HP_F57_1 (right)	type 1b TnPZ partial	1625	152699 - 154323 (HPF57_0147 - HPF57_0148)	N/D(a)	N/D(a)	HPF57_0146 (5'-methylthioadenosine/S-adenosylhomocysteine nucleosidase)	HPF57_0149 (hypothetical protein)
	GI_HP_F57_2	type 1 TnPZ	38991	284353 - 323343 (HPF57_0279 - HPF57_0311)	8(e)	tfs3b partial	HPF57_0278 (typeIIM)	HPF57_0312 (type II DNA modification enzyme )
	GI_HP_F57_3	cagPAI	36797	562215 - 599011 (HPF57_0550 - HPF57_0575)	22(c)	Type IV	HPF57_0549 (hypothetical protein)	HPF57_0576 (glutamate racemase)

a) N/D, not detected

b) TTATAATTTGAGCCATTATTTA

c) TTTCAATTTGAGCCATTCTTTA

d) TTATAATTTGAGCCATTCTTTAGCTTGTTTTTCTAGCCAAACCA

e) ACATTCTT

A GI found in strain F16 lacked similarity to known GIs of H. pylori whereas the other four GIs were homologous to transposable elements TnPZs, as recently reported [48, 49]. The GI in F16 appears to be a remnant of a prophage inserted into a restriction-modification system (Figure 6A). It is homologous to the 5'-half of the Hac II prophage found in H. acinonychis Sheeba. The F16 GI appeared to have lost its 3'-half, presumably through deletion mediated by the inserted IS605 copy. The GI included putative phage integrase genes (HPF16_0475 and HPF16_0476) that suggest the mobility of this region, and a DNA primase gene (HPF16_0468). The gene (HPF16_0469) next to the DNA primase gene had weak sequence similarity to a putative phage helicase gene (ORF35 of bacteriophage phi3626, e-value 5e-5 by TBLASTN against phage nucleotide database), which can be assumed to be the primase-helicase system found in several bacteriophages such as T3, T4, T7 and P4 [50]. Recently, a partial Hac II prophage region was reported for another H. pylori strain [51].

The other four GIs in the other three strains had sequence similarity to TnPZs [48]. One GI in F57 was entirely homologous to the type 1 TnPZ inserted into the coding region for a DNA methyltransferase with 8-bp target duplication (5' ACATTCTT) (Figure 6B). The GI in F32 appeared to have been deleted by a type 2 TnPZ (Figure 7B). Among the Korean strains, a Type 2 TnPZ was observed only in strain 51.

The plasmid in F32 (pHPF32) was similar in sequence to known theta replication plasmids with a RepB family (Rep_3 superfamily) replication protein and R3 iterons [52–54].

The plasmid in F30 (pHPF30) was similar to a group of previously characterized H. pylori plasmids such as pHel4 in H. pylori [52, 55]. This carries genes for microcin (7-aa peptide; MKLSYRN), MccB (microcin C7 biosynthesis protein), MccC (microcin C7 secretion protein), MobBCD (for plasmid mobilization), a replication initiator protein, and two relaxases. When compared to other related plasmids, a substitution in mobB and a deletion covering several small ORFs were seen. Homologous plasmids are found in G27 (pHPG27 [56]), P12 (HPP12 [49]), and v225d [22]. HPAG1 [30], B8 [57], PeCan4 and Sat464 carry a similar plasmid without the MccBC genes.

Insertion sequences (ISs) were searched for in the Japanese strains using GIB-IS [58]. An apparently intact known IS was detected in two strains: IS607 in F16; IS605 in F32.

Divergence of genes between the East Asian (hspEAsia) and the European (hpEurope) strains

We systematically examined the amino acid-based phylogenetic trees of the orthologous genes (gene families) common to the six hspEAsia genomes and the seven hpEurope genomes. Trees of 687 OGs were selected with genes of the hspEAsia strains forming a sub tree with no genes of the hpEurope strains and vice versa. Each of the orthologs was plotted according to two distance parameters: d _a for the hspEAsia-hpEurope divergence and d _b for intra-hspEAsia divergence (Figure 8A). An hspEAsia-hpEurope divergence greater than twice that of the well-defined core tree (d _a *) was seen in 47 gene families (Table 5 and 6; genes of those orthologs in each strain are listed in Additional file 5 (= Table S4)). These genes were further divided by the intra-hspEAsia divergence (d _b ) into zone 1 (lowest divergence), zone 2 (average divergence) and zone 3 (highest divergence) (Figure 8B). Six typical trees are depicted in Figure 8C. The cagA tree (e) (zone 3) has large d _a and d _b values and a low d _b /d _a value, primarily because of the divergence in a C-terminal region of the ORF. This region, including sequences known as EPIYA (Gln-Pro-Ile-Tyr-Ala) motif, is involved in host interaction [22, 59]. The tree here is consistent with previous results [22].

This tree-based analysis effectively extracted known pathogenesis-related genes (Table 5 and Table 6) as discussed below. The list also included several genes related to antibiotics. Amino acid alignments (Additional file 6) located the divergent sites. The distribution pattern of these sequences suggests a possible relationship between structure and function as detailed below for each protein. The divergence could be related to differential activity and adaptation.

The variable d _a for an orthologous group is expected to be sensitive to the presence of a member with an exceptional phylogeny. The strain B8, assigned to hpEurope in this work (Additional file 1 (= Figure S1)), has been adapted to a mongolian gerbil [57]. The strain SJM180, also assigned to hpEurope based on the tree of seven MLST genes (Additional file 1 (= Figure S1)), clustered with hspWAfrica strains rather than with hpEurope strains in the tree of the well-defined core genes (Figure 1). To examine robustness of the above classification into diverged genes, the same analysis was conducted using the 6 hspEAsia strains and 5 hpEurope strains excluding B8 and SJM180 (Additional file 7 (= Table S5)). These two analyses used all the 20 strains, because we expected inclusion of the hspAmerind and hspWAfrica strains may provide better classification of the sub trees. In addition to these two analyses, analysis with the 6 hspEAsia and 7 hpEurope strains or with the 6 hspEAsia and 5 hpEurope strains was carried out, which allowed assignment of a bootstrap value to the branch separating the hspEAsia and hpEurope strains. Comparison of these 4 analyses is summarized in Additional file 7 (= Table S5). The four sets of results agreed rather well, especially for those genes with larger d _a value: 34 among the 47 genes in Table 6 were extracted in all the 4 analyses. The bootstrap value supported the separation of hspEAsia and hpEurope well in most cases, with the bootstrap value ≥ 900 in 41 among the 47 genes.

Positively-selected amino-acid changes between the East Asian (hspEAsia) and European (hpEurope) strains

Divergence could be adaptive or neutral. We searched for sites where the hspEAsia-hpEurope changes in amino acids were positively selected [60] and found that 7 of 47 genes passed the likelihood test (Table 7; red dots in Figure 8B). These selected sites were mapped on the coding sequences (Figure 9A). For CagA, several sites were found outside the area of EPIYA segments.

Locus tag	Gene	Description	p-value(a)	Positively selected sites (b,c)
HP0547	cagA	Cag pathogenicity island protein	< 1E-21	V238R (0.994)
				A482Q (0.953)
HP0373	homC	Putative outer membrane protein	< 1E-14	E110N (0.978)
				K428H (0.986)
				T437D (0.979)
HP0492	hpaA-2	Hpa paralog	< 1E-5	S34V (0.970)
				A46Q (0.993)
				R122F (0.967)
				K127S (0.962)
HP1185	sotB	Sugar efflux transporter protein	0.00005	T50S (0.956)
				A57L (0.990)
				N134G (0.983)
				W186Y (0.980)
mHP0174		Hypothetical protein	0.0007	F144W (0.952)
mHP1415	miaA	General tRNA delta(2)-isopentenylpyrophosphate transferase	0.0002	H174A (0.992)
HP0887	vacA	Vacuolating cytotoxin A	0.002(d)	S793A (0.964) (d) N931A (0.960) (d)

a) Bonferonni adjusted.

b) Posterior probabilities of dN/dS > 1.

c) Positions are for H. pylori 26695. Residues were aligned at the same site by both Mafft [128] and PRANK [136].

d) Two vacA genes (in B38 and B8) were eliminated because they belonged to different subtypes of the gene.

Three-dimensional structure was available for mapping some of the selected sites for three of these genes (Figure 9B). The three-dimensional structure of part of VacA, the p55 fragment, is determined [61]. S793A mapped on the surface of the p55 at its C-terminal region (Figure 9B). Deletion of the p55 region reduces VacA binding to cells [62], so S793A might affect cell binding of the hspEAsia and hpEurope strains. Two selected residues of HpaA-2 were mapped (Figure 9B). The residue (H211) corresponding to the selected residue H174 of H. pylori MiaA mapped to the alpha helix 10 of E. coli MiaA [63, 64] (Figure 9B).

Diverged genes and possible biological significance

We explored the possible biological significance of the observed divergence in genes in Table 6 using gene and protein properties, as summarized in Table 5.

Host-interaction proteins

Many of the virulence factors show wide divergence between hspEAsia and hpEurope, most likely because of co-evolution with the host. We anticipate that the list of well-diverged genes (Table 6) is enriched for host-interaction and potential virulence genes. We detected positively-selected amino-acid changes in two virulence factors: cagA and vacA (Table 7).

Many OMP families showed loss of one of their resident loci (hopMN, babABC, sabAB), whereas one family (oipA) showed duplication of its locus. Some OMP genes showed internal deletions (vacA-2) or interallelic homologous recombination (hopMN). A group-specific repertoire was seen for other OMP genes (homB, hopZ and hopQ), for other criteria. We also found substantial hspEAsia-hpEurope divergence in many OMPs (Table 5). The OMPs play important roles in host interaction such as adhesion to the host cells and induction of immune responses [26]. For example, OipA induces IL-8 from host cells [70]. Systematic decay of OMP genes occurred during adaptation of H. pylori to a new host of large felines, generating the new species of H. acinonychis [36]. Hence, the above OMP changes might reflect selection and/or fine regulation in host interaction, and more specifically, may help avoid the host immune system. At least two OMPs show evidence for positive selection (Table 7). We do not yet know whether these OMP changes are related to immune response or adhesin activity.

Lewis antigen mimicry is important for gastric colonization and adhesion. The mimicry affects innate immune recognition, inflammatory response, and T-cell polarization. Long-term infection by H. pylori might induce autoreactive anti-Lewis antigen antibodies [107]. Divergence in transferase genes for LPS biosynthesis may have resulted from co-evolution with the host immune system and could be related to changes in Lewis antigens in human populations. For example, the Le(a+b+) phenotype is almost absent in Caucasian persons whereas it occurs with a higher frequency in the Asian population [108]. This might be related to differences in pathogenicity and adaptation [109].

Changes in transporter genes, the loss of a putative amino acid utilization gene, divergence in a branched chain amino acid metabolism gene, differences in acetate metabolism genes, and divergence in motility and chemotaxis genes could also be related to host interaction, because these are related to the stomach environment. An interesting question is if these changes are related to variation in human diets.

Electron transfer

Several key electron transfer components were diverged between hspEAsia and hpEurope. The multiple and drastic changes in redox metabolism were unexpected. The systematic decay of all Mo-related genes through mutations in all (6/6) hspEAsia strains was the most striking. We do not know whether our findings reflect the biased environmental occurrence of Mo or the dietary habits of human populations. The richest sources of Mo include legumes, cereal grains (and baked products), leafy vegetables, milk, beans, liver, and kidney, whereas fruits, stem and root vegetables, and muscle meats are poor Mo sources [110].

The BisC homolog, the only molybdoenzyme found in the H. pylori genome, is similar to a number of periplasmic reductases for alternative oxidants such as dimethylsulfoxide or trimethylamine N-oxide [87]. Western strains of H. pylori might be able to use N- and/or S-oxide as an electron acceptor in energy metabolism in addition to oxygen and fumarate. One hypothesis about decay of the Mo-related genes is that this anaerobic electron transport system became maladaptive in the East Asian lineage. One possibility is the radical reaction mediated by MoaA in molybdopterin synthesis is dangerous in the presence of oxygen. This could explain the observed changes in oxidative phosphorylation and acetate metabolism.

A candidate for the BisC substrate is an oxidized form of methionine, free or within a protein. Methionine is sensitive to oxidation, which converts it to a racemic mixture of methionine-S-sulfoxide (Met-S-SO) and methionine-R-sulfoxide (Met-R-SO) [111]. The reductive repair of oxidized methionine residues performed by methionine sulfoxide reductase is important in many pathogenic bacteria in general, and specifically for H. pylori to maintain persistent stomach colonization [112, 113]. H. pylori methionine sulfoxide reductase (Msr, HP0224 product) is induced under oxidative stress control and can repair methionine-R-sulfoxide but not the S isomer, even though it is a fusion of an R-specific and an S-specific enzyme [114]. BisC from other bacteria can reduce and repair the S but not the R form [111].

If the sole function of BisC is to repair methionine-S-sulfoxide, another means to repair methionine-S-sulfoxide may have appeared in the East Asian H. pylori, for example by higher expression of Msr. In this case, BisC may have been inactivated because Mo-related reactions were no longer necessary. The substitution by a DNA element downstream of the msr gene in the hspEAsia strains (5/6, all but strain 52) could be involved in the hypothesized methionine-S-sulfoxide repair activity of its product.

Another possibility is decrease of oxidative stress generating methionine-S-sulfoxide in the East Asian H. pylori. Oxidative stress is induced by acid exposure, and msr is among the oxidative stress genes induced by acid [115]. H. pylori infection has different effects on acid secretion in Europe and Asia [116]. In Europe, antral-predominant gastritis with increased acid secretion is frequent, whereas in Asia, pan-gastritis and subsequent atrophic gastritis with decreased acid secretion are common. The decrease in acid experienced by East Asian H. pylori lineages may have decreased their methionine-S-sulfoxide and made its repair by BisC unnecessary.

Downregulation of some of the Mo-related genes in a European strain under acidic conditions may be related to their decay [30]. Downregulation may occur to avoid the possible toxic effects of Mo metabolism under conditions of acid adaptation.

Taken together, our results led us to predict that the East Asian H. pylori strains are different from the European strains in electron transfer reactions and responses to oxygen and acid. Possibly related to this alteration in redox is the presence of the two acetate-related pathways in 3 out of 4 Japanese strains. These are expected to be able to switch from acetate fermentation to acetate utilization under aerobic conditions, as seen for E. coli [117]. The European strains, some of the hspAmerind strains, and the other hspEAsia strains may be regarded as mutants that lack the pta-ackA pathway and the supposedly important acetyl~P signal. Global effects of these defects on chemotaxis, nitrogen and phosphate assimilation, osmo-regulation, flagellar biogenesis, biofilm development, and pathogenicity are expected, based on the various phenotypes of E. coli strains defective in these genes [33].

Translation fidelity

Translational proteins also diverged between hpEurope and hspEAsia strains. MiaA (tRNA delta(2)-isopentenylpyrophosphate transferase) and TilS (tRNA lysidine synthetase) affect accuracy in elongation. The amino-acid change in MiaA turned out to be adaptive (Table 7). TilS affects translation efficiency at various stages. Ambiguity in translation is proposed to be important in the evolution of novel proteins by generating phenotypic and genetic diversity in the proteome for selection [118]. This role of ambiguity is similar to the evolutionary role of genome-wide modulation of mutation rates by genes such as mutS [119].

Implications for medicine

East Asian (Japanese/Korean) H. pylori appear to be quite different from European H. pylori. Our results provide a solid starting point for understanding the biology, host interaction, and pathogenesis of the East Asian H. pylori, which in most previous works were inferred from a European strain. Divergences included virulence, cell surface-related, and drug target genes. These results will affect our strategy in developing effective therapies and drugs. Questions raised by our findings include whether East Asian VacA (Figure 9B) interacts with host cells in the same way as European VacA.

The diverged gene frxA is associated with resistance to antibiotics metronidazole [120], which is frequently used in H. pylori eradication. The divergence in the frxA could affect resistance to this group of drugs in various ways. More generally, if redox metabolism differs between hspEAsia and hpEurope strains, the same drug might produce different effects, depending on intra-bacterial redox reactions.

The diverged genes included two potential drug targets (def and ftsA), so drugs that target these proteins may have different effects in East Asian and European strains. We do not know, for example, whether anti-H. pylori drugs designed from structure of European Def [98] will be as effective against East Asian H. pylori.

Remaining questions

Clearly, many studies are needed to answer these and other questions raised by the genomics results presented here. Phylogenetic analysis in the present study used OGs where genes of hspEAsia were clustered separately from those of hpEurope. Some genes do not share this topology, as suggested above for acoE deletion and hopMN recombination. We plan to study the distortion in the tree. We focused on differences between a limited numbers of strains from each group. However, there are variations within East Asian strains (Table 5). Further experimental examination of the divergence within hspEAsia, and between hspEAsia and the other strains are necessary to understand their divergence in detail. Such examination might reveal complexity in evolution and will be the subject of a separate study. The mechanisms underlying the variation, such as mutations and rearrangements, will be a subject of a separate study [25].

H. pyloriculture

Four strains were isolated from patients with diffuse type gastric cancer, intestinal type gastric cancer, duodenal ulcer, and gastritis (F57 [121], F32, F30 and F16 [122]). The ABO blood groups of the hosts were: F57, B; F32, A; F30, O; F16, B. Studies were performed according to the principles of the Declaration of Helsinki, and consent obtained from each individual after a full description of the nature and protocol of the study.

Gastric biopsy specimens from each patient were inoculated onto a trypticase soy agar (TSA)-II/5% sheep blood plate and cultured under microaerobic conditions (O₂, 5%; CO₂, 15%; N₂, 80%) at 37°C for 5 days. A single colony was picked from each primary culture plate, inoculated onto a fresh TSA-II plate, and cultured under the conditions described above. A few colonies were picked from each plate and transferred into 20 ml of Brucella broth liquid culture medium containing 10% fetal calf serum, and cultured for 3 days under the conditions as described above. A part of the liquid culture sample was stored at -80°C in 0.01 M phosphate-buffered saline (PBS) containing 20% glycerol. DNA from each H. pylori isolate was extracted from the culture pellet by the protease/phenol-chloroform method, suspended in 300 μl of TE buffer (10 mM Tris HCl, 1 mM EDTA) and stored at 4°C for PCR analysis and nucleotide sequencing.

Genome sequencing

The genome sequences of H. pylori strains F16, F30, F32 and F57 were determined by a whole-genome shotgun strategy. We constructed small-insert (2 kb) and large-insert (10 kb) plasmid libraries from genomic DNA, and sequenced both ends of the clones to obtain 26,112 (F16 and F57), 30,720 (F30) and 33,792 (F32) sequences using ABI 3730xl sequencers (Applied Biosystems), with coverage of 10.0 (F16)-, 11.5 (F30)-, 12.7 (F32)- and 10.0 (F57)-fold. Sequence reads were assembled with the Phred-Phrap-Consed program, and gaps were closed by direct sequencing of clones that spanned the gaps or with PCR products amplified using oligonucleotide primers designed against the ends of neighboring contigs. The overall accuracy of the finished sequence was estimated to have an error rate of less than 1 per 10,000 bases (Phrap score of ≥40). Sequences of the molybdenum-related genes and the genes in the acetate pathway of the four Japanese strains were verified by resequencing PCR fragments directly amplified from genomic DNA (primers are in Additional file 4 (= Table S3)). The genome sequences of other strains were obtained from National Center for Biotechnology Information (NCBI) [123]. Accession numbers are in Table 1.

Gene finding and annotation

We used the same protocol to identify genes in the four new strains and 16 other complete genomes (Table 1; gene assignment differences are in Additional file 8 (= Table 6)).

Protein-coding genes were identified by integrating predictions from programs GeneMarkS [124] and GLIMMER3 [125]. All ORFs longer than 10 amino acids were searched using BLASTP [126] against two databases, one composed of genes of 6 H. pylori genomes in RefSeq database at NCBI ("close" database), and the other composed of genes of 300 complete prokaryote genomes (one genome per one genus) available at the end of 2008, except for those in the Helicobacter genus ("distant" database). When the predicted start position differed in GeneMarkS and GLIMMER3, assignments were made by consensus of hits, with consensus against the "distant" database taking priority over the "close" one. The consensus start position among bidirectional best hits with 50% or more amino acid sequence identity for each matched region for each genome pair was determined by majority rule. Overlap of genes was resolved by comparing the results from four prediction programs. Genes encoding fewer than 100 amino acids and predicted only by Glimmer3 were dropped except for the microcin gene.

tRNA genes were detected using tRNAscan-SE [127]. rRNA genes were identified based on sequence conservation. Putative replication origins were predicted by GC-skew (window size 500 bp, window shift 250 bp).

Core genome analysis

The common core structure conserved among 20 H. pylori genomes was identified based on conservation of gene order among orthologs using the CoreAligner program [23] implemented in the RECOG system. Briefly, CoreAligner identifies the genomic core of the input genomes by taking the longest path of the neighborhood graph that consists of conserved neighborhood gene pairs, which are defined as pairs of OGs that are within a neighborhood of 20 genes in at least half of the genomes. For this analysis, we used as input a set of OGs generated by the DomClust program [24] (see "Phylogenetic profile analysis" section below for details about identification of OGs by DomClust). Absence of a gene in some genomes (at least half of the genomes) in each OGs among the core is allowed. In addition, as identified OGs are at the domain level, if a counterpart of a gene in one genome is split in another genome, different number of genes can participate in the OGs in different genomes. Thus, the number of core genes in each genome can vary. Still, the numbers of core genes varied less (1364-1424; SD = 13.5) than the total number of genes among the strains (1465-1593; SD = 33.9) (Table 1). Among those core OGs, 1079 OGs were universally conserved (conserved in the all genomes), non-domain-separated, with one-to-one correspondence, and designated "well-defined core OGs". Those 1079 OGs were used for phylogenetic analysis (Figure 1). Nucleotide sequences of genes in well-defined core OGs were aligned by the Mafft program [128], from which conserved blocks were extracted by the Gblocks program [129].

Phylogenetic profile analysis

Phylogenetic profiling was carried out using the set of OGs generated by DomClust [24]. We identified OGs with East Asian-specific features as those whose phylogenetic profiles were highly correlated to the template pattern (taking 1 for hspEAsia and 0 for hpEurope). The DomClust clustering program can identify OGs at the domain level, and was used to identify genes truncated in particular strains. Clustering was performed based on PAM (point accepted mutation) distance rather than score to ensure proper evaluation of evolutionary distances, even if one gene was truncated; in the latter case, scores may underestimate evolutionary relatedness. To clarify differences in gene-splitting patterns among strains, we did not use DomClust options to suppress domain splitting.

To identify genes with characteristic patterns of hspEAsia strains, we constructed a phylogenetic profile for each OG as a vector of examined property values (e.g., number of domains or number of duplications). For surveying patterns of gene splitting and deletion, a phylogenetic profile was constructed for each OG using the number of domains for each gene that resulted from the clustering. For surveying patterns of gene duplication, a phylogenetic profile was constructed using the number of duplicated genes (in-paralogs). To find OGs with a characteristic hspEAsia pattern, equality of the medians among different populations was tested by Kruskal-Wallis test. Tests between East Asian and European strains used the six hspEAsia strains and the seven hpEurope strains. Tests among four subpopulations used six hspEAsia, five hspAmerind, seven hpEurope, and two hspWAfrica strains.

Phylogenetic network analysis of the hopM/N family was carried out using NeighborNet [130] implemented on SpritsTree [131].

Analyses of molybdenum-related genes

H. pylori protein sequences were searched against the CDD conserved protein domain database, by RPS-BLAST [132]. Protein families extracted from the search results for Mo-cofactor synthesis or binding domain were: PF03404 (Mo-co_dimer), PF03205 (MobB), PF02738 (Ald_Xan_dh_C2), PF01568 (Molydop_binding), PF02730 (AFOR_N), PF02597 (ThiS), PF03454 (MoeA_C), PF06463 (Mob_synth_C), PF03453 (MoeA_N), PF01315 (Ald_Xan_dh_C), PF01493 (GXGXG), PF02579 (Nitro_FeMo-Co, PF01967 (MoaC), PF03459 (TOBE), PF02391 (MoaE), PF00384 (Molybdopterin), PF04879 (Molybdop_Fe4S4), PF02665 (Nitrate_red_gam), PF00174 (Oxidored_molyb), PF00994 (MoCF_biosynth), PF03473 (MOSC), PF02625 (XdhC_CoxI), PF01314 (AFOR_C), PF01547 (SBP_bac_1) (pfam name in parentheses). Homologs of two molybdoproteins [133] that were not detected in the above protein families were absent in the H. pylori genomes.

bisC was the only molybdoenzyme gene in the 20 H. pylori genomes with detected domains PF01568 (Molydop_binding) and PF00384 (Molybdopterin). A multidomain TIGR00509 (bisC_fam) was also detected in bisC.

Analyses of horizontally transferred regions

GIs were detected by searching for regions that fulfilled the conditions of: (i) longer than 5 kb; (ii) continuous ORFs not perfectly conserved in all 20 H. pylori strains; and (iii) whole regions assumed as extrinsic by Alien Hunter [134]. Counterparts of detected GIs in Amerind strains were previously reported as TnPZ [48, 49].

Genes with a large distance between East Asian and European strains

OGs diverged between six hspEAsia and seven hpEurope strains were screened based on two values related to their phylogenetic tree. The d _a value was the distance between the last common ancestral (LCA) node of hspEAsia and the LCA node of hpEurope. The d _b value was the average distance of hspEAsia from its LCA node. OGs with hspEAsia-diverged genes were screened by introducing the following conditions (with hspAmerind omitted): (i) OGs in which all the hspEAsia genes of the OG formed a sub tree without any hpEurope genes in the phylogenetic tree; (ii) OGs universally conserved (not less than 12 of the 13 genomes; not less than 10 among 11 genomes for comparison of 6 hspEAsia and 5 hpEurope strains in Additional file 7 (= Table S5)); (iii) genes with no domain fusion/fission event among the 13 genomes (within ± 20% of the mean length of the OG, measured in amino acid residues); (iv) d _a value greater than twice the d _a value of the concatenated well-defined core tree (of amino-acid sequences) (denoted as d _a *; with the resulting cutoff of d _a > 0.02324; 1079 OGs; see "core genome analysis" section above). Among 1248 OGs that satisfied the criteria (ii) and (iii), 692 OGs satisfied the criteria (i), that is, complete separation of genes of hspEAsia from those of hpEurope. The d _b * ± sd values in logarithmic scale, corresponding to 0.00550 and 0.0231 (d _b * = 0.01128) in the original scale, were used as threshold values for the three zones (N = 687; five OGs with d _b = 0 were excluded from 692 OGs satisfying the above criteria (i)-(iii)).

Amino acid sequences of the genes were aligned by the einsi command of the MAFFT program [128], from which a neighbor-joining tree was constructed by the ClustalW program [135].

A branch-site likelihood ratio test of positive selection was carried out using PAML [60] based on the multiple alignment by the einsi command of MAFFT [128]. Only residues aligned at the same site by the einsi command and by PRANK (with codon option) [136] were considered. Positively-selected residues were mapped on the p55 structure of VacA using PyMol).

Statistics

The equality of means for phylogenetic profiling between East Asian and European strains was tested by Kruskal-Wallis one-way analysis of variance by ranks, a non-parametric method for testing equality of population medians among groups. The tests were conducted using the R statistics package [137].

Accession Numbers

The accession numbers of the H. pylori genome sequences reported in this paper are: F16 [DDBJ:AP011940.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011940.1], F30 [DDBJ:AP011941.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011941.1, DDBJ:AP011942.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011942.1], F32 [DDBJ:AP011943.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011943.1, DDBJ:AP011944.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011944.1] and F57 [DDBJ:AP011945.1 http://getentry.ddbj.nig.ac.jp/cgi-bin/get_entry2.pl?database=ver_ddbj&query=AP011945.1].