Microevolution of Virulence-Related Genes in Helicobacter pylori Familial Infection

Helicobacter pylori, a bacterial pathogen that can infect human stomach causing gastritis, ulcers and cancer, is known to have a high degree of genome/epigenome diversity as the result of mutation and recombination. The bacteria often infect in childhood and persist for the life of the host. One of the reasons of the rapid evolution of H. pylori is that it changes its genome drastically for adaptation to a new host. To investigate microevolution and adaptation of the H. pylori genome, we undertook whole genome sequencing of the same or very similar sequence type in multi-locus sequence typing (MLST) with seven genes in members of the same family consisting of parents and children in Japan. Detection of nucleotide substitutions revealed likely transmission pathways involving children. Nonsynonymous (amino acid changing) mutations were found in virulence-related genes (cag genes, vacA, hcpDX, tnfα, ggt, htrA and the collagenase gene), outer membrane protein (OMP) genes and other cell surface-related protein genes, signal transduction genes and restriction-modification genes. We reconstructed various pathways by which H. pylori can adapt to a new human host, and our results raised the possibility that the mutational changes in virulence-related genes have a role in adaptation to a child host. Changes in restriction-modification genes might remodel the methylome and transcriptome to help adaptation. This study has provided insights into H. pylori transmission and virulence and has implications for basic research as well as clinical practice.

The significance of such extreme genome/epigenome diversity is not understood in detail, although the ability to adapt to the human host has been often suggested. Evidence for effects of DNA modification on gene transcription in bacteria has been limited primarily because of lack of techniques to sensitively and accurately detect methylated DNA bases and to accurately measure transcripts. The Single-Molecule Real-Time (SMRT) sequencing technology now allows methylome decoding at the single base resolution [13,14]. On the other hand, RNA-seq method in the next-generation sequencers allows accurate and sensitive measurements of transcriptome changes [15]. As a result, there are now an increasing number of reports on relationship between DNA methylation by modification enzymes and transcriptome in bacteria. For example, Type II restriction-modification (RM) systems affect global gene expression in Escherichia coli [16] and in H. pylori [17]. Type III restriction-modification systems affect expression of various genes in Neisseria, Haemophilus and H. pylori [10,18,19]. Type I restriction-modification systems affect gene expression in H. pylori [7] and Streptococcus pneumoniae [20]. There is also evidence that Type II modification enzymes affect H. pylori physiology [21][22][23].
There have been many studies of human familial transmission in attempts to identify transmission pathways. Early studies used electrophoresis-based methods or multi-locus sequence typing (MLST) based on seven conserved genes (~3 kb in total) [24,25]. Recently, more sensitive whole genome (~1.6 Mb) sequencing was used to analyze the process of familial infection [26][27][28].
In order to gain insight into H. pylori adaptive evolution in a new human host, especially in a child who was supposed to be newly infected, we used whole genome sequencing of very closely related H. pylori strains from the same family, which are difficult to distinguish by standard MLST analysis [29] and estimated the order of intrafamilial infection and then investigated which genes have experienced amino acid changes.

Results
Families infected with H. pylori of the same or very similar MLST sequence type In all, 19 H. pylori strains were isolated from five families (one strain per person) visiting a pediatrician in a hospital in Hokkaido, Japan ( Table 1). Each family member was infected with H. pylori strains of the same or very similar sequence type as judged by standard MLST analysis (S1 Table): the mother and the offspring in families K-1 and K-2, the parents and offspring in families K-3 and K-4, and the children in family K-5 [29]. The likely direction of familial transmission in families K-1 and K-2 was from mother to offspring. The allele type numbers of some of the seven MLST genes were different in families K-3 and K-4; this is based on a single nucleotide substitution, however, so we regard them as very similar sequences. The MLST results did not make it clear whether the father or the mother was the H. pylori donor to the offspring. In family K-5, transmission was likely from outside the family, which might be independent among children from the same reservoir or might have been followed by interchild transmission.
All these strains were shotgun sequenced by an Illumina high-throughput sequencer and mapped onto the known genome sequence of H. pylori strain F30 (Table 1), which belongs to the hspEAsia population (as do most isolates in Japan) [4,30]. About 90% of nucleotides were mapped with coverage of more than five reads, with detection of~35,000 nucleotide substitutions per strain. A lower mapping rate in the strains found in family K-4 can be explained by the lack of a cag pathogenicity island in their genomes. Little difference was observed when other hspEAsia genomes (F16, F32 and F57) [4] were used as a reference genome for mapping (S2 Table), so we used the mapping on H. pylori strain F30 for further analysis.

Estimation of transmission pathways from the number of nucleotide substitutions
Nucleotide substitutions between each strain pair within a family were counted and used as the distance between compared strains (S3 Table) to construct a phylogenetic tree for each family (Fig 1). As expected, the number of nucleotide substitutions between strains of distantly related MLST sequence types was greater compared to the same or very similar MLST sequence type by one to three orders of magnitude. From the whole genome sequence information, various evolutionary relations were inferred between multiple strains of the same or very similar MLST sequence type within a family.
In family K-2, comparison of three strains with the same MLST sequence type showed the distance between the two strains from children K36 and K37 was, on average, 1.7-fold smaller compared to K35 (from mother) and K36 (from child) and to K35 (from mother) and K37 (from child) (Fig 1B (ii)). The difference was statistically significant (P < 10 -5 ), supporting the hypothesis that the latest infection (transfer) occurred between the hosts of K36 and K37; i.e. transmission between children. (This inference is discussed in Discussion, below.) A comparison of three strains with the same MLST sequence type in family K-3 showed the distance between K27 (from mother) and K28 (from child) is about half compared to K26 (from father) and K27 (from mother) and to K26 (from father) and K28 (from child) (Fig 1C). The difference is significantly greater compared to assuming a random distribution of  substitution mutations (P < 10 -5 ), suggesting the latest infection occurred between the hosts of K27 (mother) and K28 (child); i.e. transmission between mother and child. The likely direction is from mother to child, considering the general pathway of H. pylori infection. (See Discussion, below.) No significant difference between distances was found for the four strains of very similar MLST sequence type in family K-4 (P > 0.01, Fig 1D).
There were three MLST sequence types corresponding to father, mother and offspring in family K-5 (see the preceding section). No significant difference between distances was detected for the three strains isolated from the offspring (P > 0.01, Fig 1E (ii)).

Detection and classification of strain-specific nucleotide substitutions
All the nucleotide variations in strains with the same or very similar MLST sequence type were compared within each family for the detection of strain-specific nucleotide substitutions. Bases corresponding to the same site were compared, and assigned as a strain-specific substitution when only a single strain had a base different from the other strains in the comparison. In the case of family K-1, there were only two strains with the same sequence type, so we could not distinguish K15-specific from K17-specific substitutions and treated all the substitutions as strain-specific mutations.
We classified base substitutions to single nucleotide variations (SNVs) or clusters of nucleotide polymorphisms (CNPs) in order to estimate the frequency of recombination ( Table 2). CNP is defined as a cluster of two or more substitutions separated by <200 bp and flanked by >200 bp of identical sequence on both sides and is considered as a sign of substitution by a recombination event [31]. Substitutions not clustered as CNPs were assigned as SNVs. A large number of CNPs was observed in family K-1, which accounted for about 85% of substitutions. The number of substitutions included in a CNP was, on average, about 11 for those of family K-1 but less than three in all the other strains. Strain-specific nucleotide substitutions were classified into nonsynonymous (changing amino acid), synonymous (retaining amino acid) and substitutions in noncoding regions. The mutation rate was calculated for isolates from children by dividing the number of SNVs (substitutions not included in CNPs) by their age in years, assuming H. pylori infection occurred soon after birth. The rate ranged from 0.7 x 10 -6 -5.2 x 10 -6 per year per site (S4 Table). The value is comparable to those reported for earlier estimations for intra-family evolution [26,31,32].

Genes with an amino acid substitution
We investigated which genes had experienced an amino acid substitution using function-based gene grouping with clusters of orthologous groups (COG) [33] ( Table 3, S5 Table). The number of genes with at least one nonsynonymous substitution was counted. We found enrichment of amino acid substitutions in COG categories of M (cell wall/membrane/envelope biogenesis) Table 3. COG enrichment of genes with strain specific non-synonymous substitutions.  and T (signal transduction) genes in multiple strains. The category M genes included those for outer membrane proteins (OMPs), lipopolysaccharide biosynthesis, lipoprotein-related proteins, penicillin-binding proteins in the cell wall and cell division proteins (S6 Table). Category T genes included spoT (included also in category K) and chemotaxis-related genes (included also in category N) (S6 Table). The genes with nonsynonymous substitutions included those decayed in hspEastAsia, including molybdenum-related genes (mogA and modB), and those decayed in hpEurope, including ackA for acetate formation [4]. The substitution mutations might represent a step in their decay.
Genes not assigned to a COG category were also analyzed. Genes annotated as encoding OMPs [34] were the most frequently enriched with the amino acid changing mutations among the strains after transmission, as reported [27]. These included alpAB, hofABDEF, horD, hopADFGINQZ, homBD, babAB and other hypothetical OMPs ( Table 4 and S6 Table).
Various known virulence factors were identified among them ( Table 4). The cagACTWY genes reside on the cag pathogenicity island [35]; vacA is known to cause vacuolation in host cells [36] and virB2 and comH are included in Type IV secretion systems, which are used for DNA import from the surrounding environment [37]. Cysteine-rich proteins [38] and HtrA protease [39] are involved in host interaction. The tnfα gene induces tumor necrosis factor alpha in the host [40] and gamma-glutamyltranspeptidase promotes pathogenesis [41][42][43]. The fucT gene is used for Lewis antigen mimicry and is important for immunity avoidance [44]. We also found amino acid substitutions in several restriction-modification genes ( Table 4).

Common and different patterns of evolution in familial transmission
We focused on three families involving children to gain insight into the steps of H. pylori adaptation to a human child host ( Table 5 and S6 Table). We emphasize that each lineage has a unique pattern of amino acid changing mutations but they all showed a change in cag and other virulence-related genes.
(i) Family K-2: K36 and K37. K36 had specific nonsynonymous nucleotide substitutions in five motility genes, five signal transduction genes, cagA, tnfα, hcpD, fucT, rfaJ-1/2, and one Type III R (restriction) gene. K37 had these types of mutations in cagA, hcpD, the RNA   polymerase subunit gene (rpoB) and two restriction-modification genes (Type IIG and Type III R) as well as four OMP genes.
(ii) Family K-3: K28. Among the 49 genes with at least one nonsynonymous nucleotide substitution specific to the child strain, there are four cag genes (cagV, cagW, cagY and cag5), vacA, htrA, five OMP genes, three restriction-modification genes (two Type I R and one Type III M) and clpX/clpP protease genes ( Table 5).
(iii) Family K-5: K23, K24 and K25. Among the three strains isolated from the three children in family K-5, K24 had the fewest (n = 12) genes with an amino acid change, among which are a virulence-related gene (comH), a restriction-modification gene (Type I R), purB for nucleotide metabolism, a peroxidase gene and a heme biosynthesis gene. K23 had nonsynonymous substitutions in cagAC genes, four genes related to motility and chemotaxis and two genes related to signal transduction. K25 had nonsynonymous mutations in cagC, three genes (ppa, rfaJ-1 and cfa) related to lipids, atpB for membrane ATP synthase and pcm for protein repair.

Discussion
H. pylori is known for sequence diversity between different strains, but strains from the same lineage can be difficult to distinguish by the standard MLST analysis using only seven genes. We undertook whole genome sequencing and distinguished the strain-specific nucleotide substitutions for isolates with the same or very similar MLST sequence type from the same family. On the basis of sequence difference, we revealed the likely pathway of evolution between these strains in some cases. Furthermore, after analyzing the nonsynonymous mutations, we suggested the strategy of H. pylori evolution during infection.
Following the construction of phylogenetic trees based on SNPs in the whole genome sequences. We inferred the direction of transfer between family members and other details (Fig  1). However, a recent study revealed a broad diversity in genome sequences in strains isolated from one specimen from one person [45]. The diversity represented in the above phylogenetic trees might well be accounted for by the diversity of lineages within an individual. In Fig 1B  (ii), for example, the mother might have transferred one lineage (the ancestor of K36) to one child and transferred another, well diversified strain (the ancestor of K37) to the other child. In order to clarify transmission pathways accurately, we need genome sequences of multiple strains from one host. The mutation rate was calculated on the basis of strain-specific SNVs. Most of the value was included with the mutation rate range in earlier calculations from whole genome information [26,31,32] and only the K24 strain was a little below the range (S4 Table). This difference might reflect the origin of H. pylori strains. Earlier work using whole genome sequences was with strains from South Africa [26], whereas strains isolated in Japan were used in the present study. Another possibility is recombination between lineages [45]. By contrast, work using 78 genes for analysis gave a much lower value for the mutation rate, using strains derived from USA, UK, Colombia, the Netherlands and South Korea [32]. This result might be related to the large difference in sequence diversity among genes: indeed, different genes can evolve at very different rates [4,5]. More analysis of whole genomes and individual genes in strains from various regions is required to fully understand the apparent variation of mutation rate.
Genes with an amino acid change might provide insight into the adaptation process. Many of these genes are related to the surface structure of the H. pylori cell, including OMP genes, lipoprotein-related genes and fucT, as found in other work on intrafamilial transmission and intrahost transmission [27,31]. H. pylori can attach to human gastric epithelial cells through various kinds of adhesion factors, including BabA (HopS), BabB, SabA (HopP), SabB, AlpA and AlpB. Protein BabA is one of the major adhesion molecules associated with severe pathogenesis in H. pylori infection, although babA expression was reported to disappear by six months after infection of Mongolian gerbils with nucleotide changes introducing a stop codon of the gene [46]. AlpA and AlpB were shown to contribute to laminin binding of H. pylori and to induction of inflammatory changes of gastric mucosa [47]. HopQ might be important in initial colonization and long-term persistence of H. pylori in the stomach by modulating the adherence to gastric epithelial cells [48]. Two different alleles of hopQ were shown to be associated significantly with the positivity of other virulence genes, including cagA and vacA [49,50].
An unexpected finding was the occurrence of amino acid substitutions in many virulence genes other than OMP genes and other surface-related genes. The genes of vacA and cagA are well known as important virulence-related genes of H. pylori. Mutations in these virulence genes were detected in H. pylori isolates of one or two family member(s) in each family but not, in general, in the isolates of all members of a family. One interesting possibility responsible for this observation is that these changes in virulence factors are related to adaptation to children in intrafamilial transmission. A related finding is that CagA is the most reactive antigen recognized by H. pylori-positive sera from children [51].
Amino acid changes in restriction-modification systems were detected in the three families involving children (see the last section in Results). A restriction-modification system consists of DNA methyltransferase, a modification enzyme, and a restriction endonuclease. DNA methyltransferase transfers a methyl group to a specific DNA sequence in the genome, which likely affects global gene expression among others. The restriction enzyme destroys DNA lacking such specific methylation resulting in genetic isolation. Recent work demonstrated the restriction-modification systems in H. pylori frequently change their presence/absence, sequence specificity and expression to remodel the methylome [7,10]. The mutations mentioned above could be related to adaptation to a new host through such epigenetics-driven adaptive evolution [52]. Substitution of target recognition domains of restriction-modification systems underlying drastic changes in recognition sequence [11,12], however, cannot be detected, in principle, by the present method based on mapping a genome sequence and SNP analysis.
Mutations were found in earlier work comparing whole genome sequences of closely related strains, especially in OMP genes, which is consistent with the results presented here [27,28,31,53]. Comparison of whole genome sequences of H. pylori isolated from grandfather, son and grandson of a family in England found amino acid changes in OMP genes [28]. Substitution mutations in OMP genes were found in inter-spouse transmission in Australia [27]. A mutation burst was found during the acute phase of H. pylori infection leading to mutational changes in OMPs and cag-related genes in humans and primates [53]. These results are consistent with the present study and we additionally found mutations commonly observed in isolates from children in two categories; i.e. virulence factors other than cag-related and restriction-modification enzymes. This difference might be caused because of isolation in young children, compared to isolation from adults in other studies, but we cannot exclude the possibility it is caused by differences in environment.
We assumed that the individual human hosts have driven the bacterial mutations in the above genes. Our procedures involved culturing the bacteria ex vivo as in almost all the works on bacterial variations. Have the ex vivo steps, as opposed to in vivo steps, induced or selected the mutations we observed here? We think that most of the strain-specific nonsynonymous substitutions were generated in vivo for the following reasons. First, many of the genes with those mutations also show rapid sequence changes in phylogeny giving long branches in their phylogenetic tree [4,5]. A simplest interpretation is the host adaptation through a nonsynonymous mutation is repeated for many generations to result in the rapid evolutionary rate. Second, comparable mutation rates per year were obtained in this and other studies based on strain culture ex vivo ( Table 2). This indicates that the number of mutations is approximately proportional to the years within a human body. This cannot be expected only through mutagenesis and selection in vivo. Third, the strain-specific nonsynonymous mutations are unique to each of the strains. For example, K23, but no other strains, carries such mutations in multiple cag genes and several motility/chemotaxis-related genes among 14 genes with annotation ( Table 5). Many strains carry those strain-specific nonsynonymous mutations in OMP genes, but the OMP repertoire is quite diverse among the strains. The difference is not likely a result of mutagenesis and/or selection during ex vivo culture, during which we used a medium of the same recipe, especially, the same batch of horse serum (see Materials and methods). From these considerations, it is, at present, natural to interpret that most of these strain-specific nonsynonymous mutations were introduced during the long-term (years of) growth in individual human stomachs although we cannot exclude some contribution of ex vivo growth.
In conclusion, our whole genome decoding of H. pylori strains from family members including children suggested adaptation of these bacteria to a new human host through mutations in virulence-related genes and restriction-modification genes in addition to OMP genes.

Ethics statement
This study was undertaken with approval from the Ethics Committees of Kyorin University, Tokyo (No. 537) and Sapporo Kosei General Hospital (H24-104). Written informed consent was obtained from the patients (5 years old or older) and also from their parents when the patients are minor.

Strains
In all, 19 H. pylori strains were obtained from five families during April 2011-December 2012 in Sapporo Kosei General Hospital, Sapporo, Hokkaido, Japan. A single colony was isolated and subcultured on Brucella medium supplemented with 1.5% (w/v) agar and 7% (v/v) horse serum (BHS medium) under microaerobic conditions. The same batch of horse serum was used for the culture to minimize possible variation between cultures. Typing of strains was done initially by seven-gene MLST for all five families [29].

Genome sequencing and mapping
After incubation for 48 h under microaerobic conditions at 37°C on BHS medium, the culture of H. pylori (about 5×10 8 colony-forming units) was collected. A Wizard Genomic DNA purification kit (Promega, Madison, WI, USA) was used according to the manufacturer's instructions to isolate genomic DNA. A DNA library for genome sequencing was constructed by Nextera XT (Illumina, CA, USA) and sequenced by HiSeq2500 (Illumina, CA, USA). About 1.4×10 6 reads (~×200 coverage) with a length of 100 bp in the form of paired ends were selected from each read data (DRA accession no. 002504) and mapped against the genome sequence of H. pylori strain F30 (accession no. NC_017365) by BWA [54]. Nucleotide substitutions were detected by SAMtools software [55] without misalignment filtering to avoid pseudo-negative detection. Lists of nucleotide substitutions were compared by customized Perl scripts for calculation of the distance between strain pairs and for detection of strain-specific nucleotide substitutions. Nucleotide sites with coverage of more than five reads for all the members of a family with the same or very similar MLST sequence type were used for the detection of nucleotide substitution. Classification of nucleotide substitution to nonsynonymous, synonymous or substitutions at noncoding regions were done according to the gene annotation of H. pylori F30 [4]. For the calculation of strain-specific substitutions, the substitutions in strains K17 and K15 in family K-1 were counted together because the family has only two strains with the same MLST sequence type and it is not possible to assign substitutions to either of the strains.
The significance of differences between distances among strains with the same or very similar sequence type was analyzed by generating a matrix assuming the same probability of nucleotide substitution accumulation for all strain pairs. Matrices were constructed 1×10 6 times and the rank list of standard deviation was compared with the standard deviation of distances in the real data for calculation of the P value.

COG enrichment
COG of genes in the H. pylori F30 genome was annotated by rpsblast [56]. Genes with strainspecific nonsynonymous substitutions were counted and the significance of COG enrichment was tested by Fisher's exact test. A gene was counted only once even if it had more than two strain-specific nonsynonymous substitutions.
Supporting Information S1