Genome assembly, annotation, and archiving

Illumina raw reads were trimmed and filtered using TrimGalore! software v.0.3.7 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) applying the quality cutoff Q30 and only keeping reads longer than 30 bp. Filtered reads were de novo assembled using SPAdes v.3.9.0 [18]. The resulting draft and complete genomes were annotated as described previously using the online automatic pipeline Rast annotation v.2.0 [19]. Summary statistics from the sequencing and assembly were collected using MultiQC v.1.0 [20]. All genomes with accompanying metadata were submitted to NCBI under the BioProject number PRJNA601302 [21]. The 254 genomes we sequenced were analyzed with 469 H. pylori genome sequences published by others. The final study dataset comprised 723 genome sequences that were used for subsequent analyses. Of these, 337 isolates came from populations in the Americas and 386 from populations in other continents (Supplementary Table 1). All genomic sequences were submitted to the online automatic pipeline Rast annotation.

Core and accessory genome characterization

A pan-genome list was constructed from all 723 isolates in this study by automated annotation of all genomes in the dataset using RAST [19]. Duplicate genes were removed if they exceeded a BLAST threshold of 70% nucleotide identity [22]. This reference pan-genome list comprised 2545 unique genes that appeared in at least one of our 723 genomes. Consistent with previous studies and the whole-genome MLST principle [23], the gene complement and allelic variation of each isolate were determined by comparison with the pan-genome list, with gene presence recorded as a BLAST match exceeding 70% sequence identity over 50% of sequence length. Core genes were present in >90% of the genomes and accessory genes were present in at least one isolate. Each gene was aligned individually using Mafft [24] and concatenated into a single multiFASTA alignment file for each isolate.

Analysis of genome ancestry

Genome-wide haplotype data were prepared by calling SNPs from a core-genome alignment using BEAGLE v.3.3.2 [25]. A total of 498,461 SNPs were identified in 1385 genes. The genome-wide haplotype data was used in FineSTRUCTURE v.0.02 to define isolate populations based on the similarity of the haplotype copying profiles obtained by CHROMOPAINTER v.0.02 [26]. FineSTRUCTURE was run for 100,000 iterations of both the burn-in and Markov chain Monte Carlo (MCMC) method to cluster individuals based on the co-ancestry matrix as described [27]. The results were visualized as a heat map with each cell indicating the proportion of DNA “chunks” a recipient receives from each donor. The population structure was analyzed in more detail using principal component analysis (PCA) with the co-ancestry matrix obtained previously. Multiple PCAs were calculated from 1 to 11 principal components (PCs) using R.

To identify the proportion of ancestry of H. pylori isolates from the Americas, we conducted chromosome painting using ChromoPainterV2 [26], designating both the American and non-American populations as donors (723 genomes), and American isolates as recipients (337 genomes). The algorithm identifies important haplotype information from dense data such as SNP data, and efficiently describes shared ancestry within a recombining population. Each individual is painted using all the other individuals as donors and the result is visualized in a bar plot built with R [28]. Finally, we calculated the average and standard deviation of the proportion of ancestry from each population present in each American country, which were visualized in box plots using R.

Fst and GenomegaMap analyses to identify gene variants fixed in American populations

Local adaptation in American isolates compared to non-American ancestors was evaluated using the Fst test within the R package PopGenome v.2.2.4 [29]. The Fst calculation was performed by comparing all genes of the pan-genome that were present in >50% of the isolates (1649 genes) in American versus non-American isolates. The criteria for data inclusion in the multiple alignments were as follows. First, we selected genes present in >50% of the studied strains, then the sequences that did not constitute >90% of the length of the gene were removed; finally, each position that did not constitute >50% of alignment depth was also removed. The final alignments were then use to identify polymorphisms using SNP-sites. [30]. A list was obtained by ordering the 1649 genes according to the Fst value, and the top 35 genes (Table 1) were analyzed further. The consensus nucleotide and amino acid sequence alignments of these genes from all the American and non-American isolates were calculated using the SeaView v.4 software [31] and compared against the reference strain 26695 [32] to determine synonymous or non-synonymous mutations. From the top 35 genes, we selected virulence genes with available information on the structure of the proteins. In these genes, the structure of the protein was obtained from Protein Data Bank and the Fst values of each position were visualized in the structure with gradient color, using PyMOL software (The PyMOL Molecular Graphics System, v.2.3.2 Schrödinger, LLC), highlighting positions with an Fst value above the upper 99th percentile.

Genomegamap was used for Bayesian estimation of d_N/d_S ratio (also denoted K_A/K_S or ω) of vacA, babA, and cagA genes. MCMC was run twice for each analysis with 1000,000 iterations and a burn-in of 22,000 iterations. Both runs were compared for convergence at several parameters and merged to obtain distributions; graphs of omega data were produced using R [33]. Differences in positions under positive selection between American and non-American strains were compared using the MCMC output from the previous analyses.

Novel American subpopulations

The addition of 254 isolates and genomes from 14 geographical locations (eight not previously represented in whole-genome analyses) greatly improved the resolution of H. pylori populations in the Americas (Fig. 1 and Supplementary Fig. 1). This permitted the identification of novel subpopulations across the continent and an improved definition of the ancestral contribution in each regional subpopulation. Three subpopulations with European ancestry were identified, all with a hspSWEurope predominance: hspSWEuropeColombia (previously reported as hspSEuropeColombia) and two novel subpopulations, the hspSWEuropeHonduras including isolates from Honduras and neighboring countries (Nicaragua, Guatemala, and El Salvador) and hspSWEuropeMexico with isolates from Mexico, but also from North, Central, and South America. The greater depth of sampling also allowed improved resolution of the hspAfrica1NAmerica subpopulation, previously described to be present in the US and Canada. This subpopulation was found to be of predominantly West African ancestry (Supplementary Fig. 1), and to include several newly sequenced isolates from Brazil and a few from Colombia. Two other Latin American subpopulations with African, specifically hspAfrica1SAfrica, ancestry were identified: the previously described hspAfrica1Nicaragua population [10, 11] (Supplementary Fig. 1) that includes isolates predominantly from Nicaragua and neighboring Honduras; and the previously reported hspAfrica1MiscAmerica [11], comprising isolates from Mexico and Colombia (Supplementary Fig. 1). These refined descriptions of subpopulations demonstrate ongoing admixture in this region.

Population structure was further studied using PCA analysis (Fig. 1), where the first PC reflects the divergence of Asian strains with hspEAsia at one extreme and African strains at the other. The second PC, on the other hand, documents the distance of the hpAfrica2 population from all others (Fig. 1a). We also observed a component that clearly shows the separation of African and European subpopulations (Fig. 1b). Whereas hspNEurope is closer to hpAsia2, hspSWEurope is closer to European subpopulations of the Americas. The components best illustrating the relationship and distribution of American and non-American subpopulations are depicted in Fig. 1c. This shows hspSWEuropeMexico to be closer to hspSWEurope, while hspSWEuropeHonduras is closer to hspAfrica1Nicaragua, and hspSWEuropeColombia is clearly separated from the others. Separation of African subpopulations in the Americas is also clear: hspAfricaNAmerica is closer to hspAfrica1Wafrica.

The only American isolates that grouped with the Asian populations HpAsia2 and hspEAsia were Indigenous American isolates (often referred to in the literature as Amerind, a term that may have ambiguous interpretations but is generally used to refer to Indigenous Peoples of the Americas) (Fig. 1b). This supports previous suggestions that the Indigenous American H. pylori ancestors were replaced by the conquistador H. pylori European strains [37] in American populations following European colonization. Furthermore, two Indigenous subpopulations were distinguished: one including strains from Canada (hspIndigenousNAmerica) and the other with strains from South America (hspIndigenousSAmerica), formerly jointly designated hspAmerind  (Supplementary Fig. 1). It has been suggested that Indigenous strains were outcompeted by European strains [38]. To address this we analyzed strains with >10% Indigenous ancestry (Supplementary Fig. 2). The figure shows strains ranging from 100% Indigenous ancestry in Canada (hspIndigenousNAmerica) and over 80% Indigenous ancestry, particularly in the hspIndigenousSAmerica subpopulation, down to around 10% Indigenous ancestry. These strains display a gradient of admixture, with chunks of European, African, and Asian ancestry consistent with in vivo competition; this pattern suggests that Indigenous strains recombined with European and African strains in mestizo populations and were gradually excluded from these populations.

Countries characterized by distinctive admixture profiles

To infer recent recombination in American H. pylori subpopulations we used a reference collection of American and non-American donor genomes to paint each American recipient genome, and CHROMOPAINTER v.0.02 to infer the number of DNA chunks donated to each recipient genome from each donor source (Fig. 2a). The three American subpopulations with European ancestry and the three with African ancestry were all painted with a mixture of European and African donors from both American and non-American populations, indicative of a high level of admixture occurring in the Americas. A large proportion of these genomes were painted by donors from their own population, particularly the subpopulations hspSWEuropeColombia, hspSWEuropeHonduras, and hspAfrica1Nicaragua, where 50% or more of the chromosome was painted by their own population (Fig. 2a). This effect was most pronounced in the central parts of the continent, specifically Honduras, Colombia, and Nicaragua (Fig. 2b). The proportion of self-identity (fraction of the genome painted by strains from their own population) is variable in each subpopulation with a gradient from low to high self-identity (Fig. 2a). Admixture was lower in the two Indigenous subpopulations, particularly in the hspIndigenousNAmerica, suggesting these populations have remained genetically isolated, and both carried traces of hspEAsia (Fig. 2b).

Open in a separate window

Fig. 2

Representation of the ancestry admixture in H. pylori strains from American populations, as analyzed with chromosome painting.

Each column represents one strain and the color indicates the proportion of the corresponding ancestry in that genome. The color code for each subpopulation is shown on the right of Fig. 2. a Representation by subpopulation, from left to right, African, European, and Indigenous American subpopulations; b admixture profiles by country along the American continent, from left to right, from Canada to Argentina.

Even though some isolates in different countries showed similar ancestry, the proportion of inheritance from mixed American and non-American ancestors varied across countries (Fig. 2b and Supplementary Fig. 3), consistent with ongoing admixture between regional subpopulations. Each country can be differentiated from their European and African ancestors as well as American neighbors. Figure 3 and Supplementary Fig. 3 illustrate the composition of genomic ancestry along the continent, showing a particularly high degree of admixture in Mexico—the gateway to the Americas for the Spanish. Figure 3 also summarizes probable migrations to the Americas from other continents.

Open in a separate window

Fig. 3

World map illustrating H. pylori ancestral contributions to the Americas (colored arrows) from other continents (black arrows).

Filled circles represent the dominant H. pylori populations at the different geographical locations and the left panel shows chromosome painting results for the isolates from the respective American countries studied.

Some American isolates retain European or African ancestry

Many American isolates retained strong ancestral roots from other continents and could be classified within European or African populations (Fig. 2a). The first group included isolates from North, Central, and South America, with high levels of hspSWEurope ancestry (Fig. 2b). A second group was classified as hspSEurope, containing mostly isolates from the US, Canada, and Argentina. In addition, a small group of isolates from North America (US and Canada) but none from Latin America was classified as hspNEurope. Finally, a small number of American strains were classified within African subpopulations: some from Honduras as hspAfrica1SAfrica and seven US isolates as hspAfrica1WAfrica (Fig. 2a).

Virulence genes are fixed in American populations, driving evolution of H. pylori

Fst analysis was used to identify genetic variants more common in the Americas than in populations in other continents. Whole-genome analysis (American vs non-American populations) was used with a relaxed threshold for gene presence to include genes present in >50% of the genomes, enabling us to include commonly known virulence genes, such as the cagPAI—which are only present in a subset of strains. The distribution of Fst values over the pan-genome (Fig. 4) was used to select 35 genes with the highest Fst values (labeled in red in Fig. 4 and Table 1). In total, there were 142 sites with significant Fst values in these 35 genes, of which 22 encode recognized virulence factors and membrane proteins, including 13 Cag pathogenicity island genes, vacA, babA, hofC, and ppiC. Furthermore, the strongest fixation index was identified in the gene encoding for CagF, a chaperone for CagA. Other genes with high numbers of significant sites were vacA (with 13 sites), cagA (with 9 sites), and the cagPAI ATPase (with 8 sites). Fifteen of the 35 genes were present in all strains studied including vacA, hofC, ppiC, three other membrane proteins, fecA and a Zn-dependent protease (Table 1). These results suggest that virulence plays a strong role in regional adaption to specific human populations.

Open in a separate window

Fig. 4

Whole-genome Fst analysis to identify genetic variants that are significantly more common in the Americas than in the rest of the world.

The X-axis indicates the nucleotide sites and the Y-axis shows the Fst value for each site. cagPAI genes include HP0522, HP0523, HP0524, HP0528, HP0525 HP0530, HP0534, HP0537, HP0538, HP0539, HP0541, HP0543, HP0544, and HP0547.

Fixed non-synonymous mutations occur in functionally important domains of virulence proteins

In order to relate mutations significantly fixed in American populations to specific protein domains we further characterized these mutations in virulence genes with corresponding protein structure data. In cagA six non-synonymous (NS) substitutions were identified, four in domain II and two in domain III (Fig. 5a and Supplementary Fig. 4a), both regions important for the diverse functions displayed by the protein, specifically the binding of CagA to the inner surface of the host cell and recruitment of PAR1, respectively. In vacA, nine of the significant Fst positions led to NS changes and all these amino acid changes were located in the mid region (referred to as the m region) of the protein (Fig. 5b and Supplementary Fig. 4b), the region that targets the protein to the host cell receptor. In BabA, three aa positions—aa 49, 53, and 153—showed a NS substitution (Fig. 5c and Supplementary Fig. 4c); this last residue has been implicated in sensitivity to acid. Finally, in ppiC, a peptidyl-prolyl isomerase, three NS substitutions occur in the chaperone domain (two substitutions in D101N and X97A) (Fig. 5d and Supplementary Fig. 4d). The position with the second-highest Fst value among all the genes was a synonymous substitution in aa position 100 of PpiC, also in the chaperone domain (Supplementary Fig. 4d).

Open in a separate window

Fig. 5

Localization of non-synonymous mutations fixed in the Americas in virulence proteins of H. pylori.

The available crystal structure of each protein was used to indicate the sites with amino acid changes. a CagA, b VacA, c BabA, and d PpiC. The color scale represents the Fst value for each position.

To further document specific differences in virulence genes between American and non-American strains, we studied positions under positive selection using GenomegaMap. As with Fst, results showed that major virulence genes (cagA, vacA, and babA) present several positions under strong positive selection (Supplementary Figs. 5–7). Although the regions of the genes with positions under selection were similar for American and non-American strains, the number of positions and strength of the ω values were higher for the American strains. This suggests stronger selective pressure in functionally active regions of virulence genes. In cagA most of the positions along the gene have as high as 35-fold ω values in American isolates, occurring mostly in regions representing extra EPIYA-C and CM motifs (Supplementary Fig. 5) [39]. Similar to Fst analyses, we found positions in the mid region of vacA that are under positive selection in the Americas, but not in strains from other continents (Supplementary Fig. 6). In babA, GenomegaMap showed positions in the head domain that were under positive selection (Supplementary Fig. 7), consistent with the Fst analyses.

This work was supported by “Consejo Nacional de Ciencia y Tecnología” (CONACYT; no. 576951) grant to ZYM-R, and by Magnus Bergvalls Foundation and Swedish Society for Medical Research (SSMF) to KT. SKS and BP are funded by the Medical Research Council (MR/L015080/1). Supported in part by the US National Cancer Institute National Cancer Institute to DRM (P01 CA028842, R01 CA190612, K07 CA125588, P30 CA068485). We are very grateful to Daniel Falush (Unit of Statistical Genetics of Bacteria, Institut Pasteur of Shanghai) for his critical input during the analysis of the presented data. We acknowledge the Aklavik H. pylori Project Planning Committee (Northwest Territories, Canada) for providing input on our presentation of this analysis.