In Vivo Analysis of the Viable Microbiota and Helicobacter pylori Transcriptome in Gastric Infection and Early Stages of Carcinogenesis

16S rRNA reads were extracted from the RNA sequencing data using Metaxa2 (20) and normalized against the total read count in each library. Applying a cutoff of at least 1 rRNA read per million reads, 9 bacterial phyla and the eukaryotic kingdom Fungi were identified in at least one of the samples (Fig. 1A). Biopsy specimens from individuals initially scored as H. pylori positive by urea breath test (UBT), serology and culture contained on average 64 times (P < 0.01) more bacterial rRNA reads than those from individuals scored as H. pylori negative (see Fig. S2 in the supplemental material), and the majority of those were assigned to the phylum Proteobacteria (Fig. 1A). Most of these reads (4.8% to 99%; average, 82.1%) were assigned to the Helicobacter genus. Although the community compositions in clinically H. pylori-positive and -negative individuals seemed to be different (Fig. 1; see Fig. S3 in the supplemental material), those effects were not significant on the phylum level. Thus, if H. pylori infection does have an effect on the composition of other bacteria on the phylum level, it seems to be minor, as differences between individuals overshadow it.

Previous studies have shown that the gastric microbiota is complex, and conflicting evidence exists on whether there is a difference between the stomach microbial communities in H. pylori-positive and -negative individuals or at different disease stages (3, 5, 7, 16). Bik and colleagues found no significant differences in taxonomic composition between H. pylori-negative and -positive subjects, confirming our observations (5), while other studies have found that Helicobacter pylori-negative subjects had higher relative abundances of Firmicutes, Bacteroidetes, and Actinobacteria and that H. pylori-positive individuals had greater abundances of non-H. pyloriProteobacteria, Spirochetes, and Acidobacteria (21). Our data confirm the finding of Firmicutes, Bacteroidetes, and Actinobacteria in several of the subjects with low levels of H. pylori, with Actinobacteria being the dominating phylum in one case. These data support the findings of Maldonado-Contreras et al., who also reported higher levels of Actinobacteria in H. pylori-negative individuals (7).

Metaxa2 analysis revealed the presence of 350 bacterial genera across all samples, of which 33 had an abundance of higher than one per million reads in at least one sample (Fig. 2). The most abundant genera found in all biopsy specimens were Helicobacter, Pseudomonas, Acinetobacter, Escherichia, and unclassified Enterobacteriaceae (Fig. 2). Interestingly, Helicobacter was present in every sample, being the dominant genus also in 5 out of the 7 samples taken from individuals classified as uninfected using conventional methods. However, the seven samples from subjects classified as H. pylori negative contained approximately 100- to 10,000-fold fewer Helicobacter-associated reads than biopsy specimens scored as H. pylori positive. These findings were further supported by one-step PCR 16S rRNA amplicon sequencing of the same samples, showing highly similar results (see Fig. S4 in the supplemental material). In general, Helicobacter rRNA reads constituted 91.0% to 99.0% of all bacterial rRNA reads in the H. pylori-positive samples, except in sample Atr8 (56.7%), while in the individuals scored as H. pylori negative, the proportion of Helicobacter reads ranged from 4.7% to as many as 63.2% of all bacterial rRNA reads. Similar observations have been reported previously (5, 7, 8, 16), but in contrast to the case for previous studies, our findings represent actual expressed sequences and imply that these bacteria are viable and transcriptionally active, although they still could constitute a transient population. However, since these samples were derived from tissue biopsy specimens, to be detected by the method used in this paper, such transient bacteria would still need to be attached to the surface, speaking in favor of an established population rather than a transient one. These findings suggest that the frequency of low-level H. pylori infection might be higher than anticipated; it is even possible that H. pylori may be present to some degree in a vast majority of individuals. However, the numbers and/or activity of these bacteria must be very low in order to go undetected by conventional methods, and such low numbers of H. pylori are probably not contributing significantly to disease risk, since only very weak signs of inflammation were seen in 2 out of 5 of the H. pylori-negative controls (see Table S1 in the supplemental material).

To verify these findings with deeper sequencing and complementary methods, we reverse transcribed the corpus RNA samples and amplified and sequenced the V3-V4 region of the 16S rRNA gene. This was also done for antrum samples from the same individuals. Although level were also very low using this method, we could verify the presence of Helicobacter pylori-specific amplicons in this analysis, after subtracting sequences found in no-template PCR control, in both antrum and corpus biopsy specimens.

It is possible that some of the reads classified only to the level of the Helicobacter genus derive from species other than H. pylori, some of which are not detected by traditional tests, as has been suggested previously (22). In the transcriptome sequencing (RNA-seq), the samples from some individuals contained reads classified as Helicobacter canis, although no such sequences could be detected through amplicon sequencing. However, the numbers of H. canis reads were very low, at most 0.44% of all Helicobacter reads and 0% in the clinically H. pylori-positive samples.

Finally, there is also the possibility that this is a technical artifact. Indeed, barcode leakage in Illumina sequencing is known to happen, but precautions to limit this were taken, as described previously (23). The samples were randomized in all steps of the workflow so that there would not be any batch effects from either the extraction, library preparation, or sequencing. However, this means that we cannot exclude contamination between sample groups, but we find this explanation unlikely in this case, since the transcripts detected in the clinically H. pylori-negative individuals were relatively diverse and did not represent only the most abundant transcripts overall and the 16S rRNA amplicon sequencing confirmed the finding of H. pylori rRNA in these samples. It should also be noted that Nicaragua is a country where the prevalence of H. pylori is high, around 70% (24), and it remains to be seen if these findings can be replicated in low-prevalence populations.

Testing for correlations between Helicobacter and the abundance of other genera, we found Campylobacter, Sulfurospirillum, and Deinococcus rRNA levels to be significantly associated with high abundance of Helicobacter (see Table S6 in the supplemental material). Deinococcus and Sulfurospirillum were also found exclusively in patients who were H. pylori positive by conventional methods. Campylobacter and Sulfurospirillum both belong to the family Campylobacteraceae and the families Helicobacteraceae and Campylobacteraceae together form the Epsilonproteobacteria subclass of Proteobacteria. It remains to be investigated whether these closely related species benefit from coexistence in the human stomach or whether the correlation is due to similar environmental preferences among Epsilonproteobacteria. Deinococcus, however, belongs to the phylum Deinococcus-Thermus and is regarded to be an environmental extremophile not expected to inhabit the human stomach, but the occurrence of Deinococcus in this study confirms one previous report (5).

Notably, we found that there seems to be a set of bacterial genera that are present in high, but similar, abundances across most of the study subjects, regardless of the levels of H. pylori, including Escherichia, Acinetobacter, Pseudomonas, and Streptococcus (Fig. 2). Contrary to previous reports (9), we did not see any significant effect of clinically determined Helicobacter status on the Shannon or Simpson diversity of other bacterial genera. However, the richness of genera was larger in Helicobacter-positive samples than in those scored negative by clinical tests (P = 0.01193). When Helicobacter was included in the analysis, genus richness was not significantly different (P = 0.4983), while the diversity and evenness estimates of the communities drastically changed (P < 10⁻¹⁰; specifically, for Shannon diversity P = 1.75e⁻¹⁵, for Simpson diversity P = 5.586e⁻¹⁴, and for evenness P = 1.436e⁻¹¹), resulting from lower diversity indices for the H. pylori-positive individuals than for the H. pylori-negative individuals due to the very high abundance of Helicobacter rRNA. Notably, earlier studies of the stomach microbiota based on sequencing of DNA as opposed to RNA have indicated substantially higher taxonomic diversity (5, 8, 9). However, many of these taxa exist in very limited quantities and are likely to be transcriptionally inactive, suggesting that their influence on the gastric ecosystem may be insignificant. Furthermore, this study also emphasizes the very small quantities of non-Helicobacter bacterial genera in stomach biopsy specimens (Fig. S4). Such small quantities of viable bacteria increase the impact of contaminations during nucleotide extraction and PCR amplification, which must be accounted for, for example, by comparison to negative controls.

To determine whether the composition of the stomach microbiota was influenced by changes in the stomach environment as a result of disease progression, we first analyzed the community composition in relation to the histopathological changes as determined based on four biopsy specimens in the corpus (see Table S1 in the supplemental material). The analysis did not reveal any segregation between the groups (see Fig. S3 and S5 in the supplemental material), but all groups of H. pylori-positive individuals overlapped with each other, with only a few individuals as outliers. However, as atrophic and metaplastic changes can occur very variably across the tissue and the overall histological scores were based on other, albeit adjacent, biopsy specimens, there is a substantial likelihood that these scores may not correspond to the atrophic status of the biopsy specimen that was used for RNA extraction. Therefore, to determine the level of atrophy in the tissue independently of the pathology scores, we took advantage of the fact that we also had information about the host gene expression in the RNA-seq data.

We used the expression levels of genes associated with atrophic gastritis, i.e., ATP4A, ATP4B, GHRL, GIF, CCKRB, PGC, PGA3, and PGA4, to determine the level of “molecular atrophy” in the tissues (see Materials and Methods for details). These genes are linked mainly to the presence of parietal cells, a cell type lost in atrophic gastritis (25, 26). Also, by using the gene expression of the proton pump-encoding ATP4A and ATP4B of the human hosts, we obtained a nonsubjective measure of the acid-secreting capabilities of the microenvironment of the bacteria analyzed. Regression analysis was performed on bacterial rRNA abundance relative to degree of atrophy and mean ATP4 expression. We found that the abundance of the bacterial genera did not change significantly in relation to the level of molecular atrophy in the tissue. Furthermore, we detected no association between corpus atrophy level or disease status and differences between the communities in the antrum and corpus sample within each individual. In line with our results, Dicksved et al. reported that gastric cancer patients were dominated by the genera Streptococcus, Lactobacillus, Veillonella, and Prevotella but that there were no significant differences from healthy controls (16). Aviles-Jimenez and colleagues did find significant differences in abundances of taxa when comparing nonatrophic gastritis to gastric cancer; however, subjects having intestinal metaplasia overlapped both the gastritis and cancer groups (6). Our results were based on samples ranging from no to high-grade atrophy, and in this set we could not identify any distinct set of bacterial genera that changed with more severe atrophy.

To study H. pylori gene transcription at the site of infection, we extracted the H. pylori transcriptome from the RNA-seq data. Due to the low proportion of bacterial RNA in the samples, we used only the expression values for the genes passing certain threshold criteria (see Materials and Methods), since expression estimates of genes with lower expression were at risk of be highly influenced by noise. The thresholds applied allowed differential expression analysis of 285 genes (see Table S5 in the supplemental material). In line with the rRNA classification, we detected H. pylori transcripts in all of the individuals, but the total amount varied substantially between the samples (Fig. S2), and for the clinically H. pylori-negative individuals, fewer than 100 reads per million reads mapped to the H. pylori genes. The number of reads mapping to Helicobacter pylori mRNA per million total reads was strongly correlated with the number of bacterial 16S rRNA sequences per million reads (r² = 0.76, P < 0.0001).

Similar to the results of previous studies (27), we found high expression of several of the classical toxins and adhesion factors of H. pylori (Table 1). The napA and flaA genes were among the most highly expressed genes, while, for example, vacA, babA, and sabA had lower expression levels. Interestingly, HP1192, a secreted protein involved in flagellar motility, was the most highly expressed gene after the rRNA genes. Among the H. pylori genes with highest overall expression were genes involved in pH regulation and ammonia production, including those encoding urease subunits A and B (ureA and ureB), a pH-regulated permease (ureI), and an aliphatic amidase responsible for ammonia production, acylamide amidohydrolase (amiE). Also, two genes encoding the recently characterized histidine-rich, metal-binding polypeptides Hpn (HP1427) and Hpn2 (HP1432) (28), involved in sequestration and intracellular trafficking of nickel, were among the most highly expressed genes.

Urease hydrolyzes urea to ammonia and carbon dioxide, which neutralizes the pH in the acidic stomach. Although our samples showed a progression toward corpus-predominant gastritis and intestinal metaplasia, conditions that usually are associated with an increase of pH, we could not observe a significant change in ureA and ureB expression at later stages of disease. This implies that H. pylori need to maintain high expression levels of urease even during progression toward gastric cancer. The high expression of hpn and hpn2 is also very interesting in this context, since nickel is essential for urease activity and these two proteins regulate the availability of nickel in the cell.

Ammonia is of great importance as a nitrogen source and as a cytotoxic molecule. In addition to urease, H. pylori possess two aliphatic amidases responsible for ammonia production: AmiE, a classical amidase, and AmiF, a new type of formamidase. The amiE and amiF genes were highly expressed in the biopsy specimens, confirming the importance of ammonia and pH regulation for H. pylori infection and suggesting that H. pylori, when associated with the corpus epithelium, faces extensive pH stress (29). In addition, genes involved in protection against oxidative damage, encoding catalase (katA), thioredoxin (trxA), alkyl hydroperoxide reductase (tsaA), flavodoxin (fldA), and superoxide dismutase (sodB), were highly expressed.

While the absolute levels of expression are important, they do not explain how H. pylori regulates its gene expression in response to the host environment and progressing disease stages, which change the microenvironment at the site of infection. To study the bacterial gene expression at sites of premalignant tissue changes, we compared the expression levels of the 285 identified H. pylori genes to the tissue atrophy levels determined as described above. Albeit not significant after correction for multiple testing, several H. pylori genes showed clear trends of expression related to the degree of atrophy in the tissue (Fig. 3; see Table S7 in the supplemental material). Several of these genes also showed differential expression between the disease stages gastritis and extensive atrophy, but apart from this, few changes were seen associated with the pathology groups. Furthermore, the genes encoding outer membrane protein AlpB, flavodoxin FldA, and beta-lactamase HcpD were more highly expressed at high atrophy than at low atrophy levels (Fig. 3). Several genes were also negatively correlated with the level of atrophy, such as the genes encoding transfer-mRNA SsrA and hypothetical proteins HP1289 and HP1223.

Few studies on Helicobacter pylori gene expression in vivo exist to date and only a few genes have been studied in relationship to pathological changes in the mucosa (30, 31). When the gastric corpus mucosa becomes atrophic, several factors change in the local environment for stomach bacteria due to changes in gland cell type composition, e.g., decreased acid secretion due to loss of parietal cells. Several studies have investigated the transcriptional response to acid and other environmental stressors in vitro, showing that katA as well as cagA, vacA, and alpB expression levels are affected by acidic stress. HP1289, which encode an H. pylori-specific protein of unknown function, also has been shown to be regulated by pH (32).

AlpB is a colonization factor and binds to laminin (33). Upregulation of alpB, but not other adhesion factors, at high levels of atrophy might indicate a specific role for this adhesion factor. As laminins are located in the extracellular matrix of the gastric basement membrane (34), the Alp adhesins might play more prominent roles in areas of tissue damage or in bacterial invasion of the tissue. Laminins have also been shown to be differentially distributed in different parts of the healthy stomach and differentially expressed in premalignant and malignant lesions (35), suggesting that it might be beneficial for the bacteria to adapt to these changes.

In the differential expression analysis, several genes showed a tendency of being influenced by tissue atrophy levels, albeit not significantly after correction for multiple testing. Our sample groups were small, particularly the group with low levels of atrophy, and thus the statistical power may not be sufficient for detection of more subtle changes. Additional studies using larger sample groups and sequencing depths will probably uncover additional biomarkers for severe atrophy and H. pylori-associated disease progression. It would also be valuable to study this in subjects from several geographical locations, since the microbiota has been shown to vary in different areas of the world (36). Furthermore, the low levels of bacterial RNA in the samples in relation to human RNA resulted in exclusion of many genes from the analysis due to the fact that their read counts were below the detection limit in many of the samples, which could be alleviated by increasing sequencing depth or by depletion of human RNA in further studies.

This study shows that metatranscriptomic analysis of the gastric microbiota is feasible and can provide new insights into how bacteria behave in vivo in different stages of disease progression and under different environmental conditions. Coupling metatranscriptomic data with host gene expression levels also gives the possibility to study host-bacterial interaction in an unprecedented manner. Here, we describe the viable stomach mucosa-associated microbiota in an unbiased fashion, and at the same time we in detail study a specific species, Helicobacter pylori, in its niche in individuals with different forms of disease-associated histological changes. We found that although H. pylori infection does not change the bacterial diversity, there are correlations between the presence of Helicobacter and the Campylobacter, Sulfurospirillum, and Deinococcus genera. Analysis of viable H. pylori at the site of infection showed that expression of several genes showed a tendency to vary with atrophy levels in the infected stomach mucosa. However, the low proportion of bacterial RNA in the biopsy specimens calls for further development of enrichment techniques for microbial RNA in complex samples. Better enrichment strategies would further contribute to the use of in vivo large-scale transcriptomics analyses and would shed essential light on the host-microbial interaction at different stages of disease.

Patients were selected from a cohort of 149 patients undergoing endoscopy during 2010 and 2011 at the Hospital Escuela Lenin Fonseca in Managua, Nicaragua, an area of high H. pylori and gastric cancer prevalence. The sample collection was approved by the Human Research Ethics Committee at the University of Gothenburg (approval number 176-10) as well as by Universidad National Autonóma de Nicaragua, Managua (UNAN). All individuals included in the cohort agreed to participate in the study, and informed oral and written consent was obtained from each patient before participation.

All patients were tested for Helicobacter pylori infection by the urea breath test (UBT) using an IR300 spectrophotometer (Otsuka, MD, USA), by H. pylori serology according to procedures described previously (37), and by H. pylori culturing. From all patients, several endoscopic biopsy specimens were collected from different anatomic locations of the stomach, together with serum samples and whole blood to enable a variety of analyses on each patient. From this cohort, 25 patients were selected based on pathology grading of 4 antrum and 4 corpus biopsy specimens from each individual. The biological groups selected were (i) patients with gastritis, i.e., no signs of atrophy in any of the corpus biopsy specimens (Gast), (ii) patients with low-grade atrophic gastritis in the corpus (Atr), (iii) patients with extensive corpus atrophy but no signs of metaplasia (EA), and (iv) patients with intestinal metaplasia in the corpus (Met). In addition to this, we also selected five subjects who had no Helicobacter pylori infection as determined by UBT, serology, and culturing and no or very low signs of inflammation in the corpus mucosa as a control group (Min). A brief description of the subjects is given in Table 2 and more detailed clinical data in Table S1 in the supplemental material.

Of the 30 subjects, the majority (22/30) were women, and the median age was 41.5 years, ranging from 23 to 66 years (Table 2). All the subjects with gastritis and corpus atrophy were H. pylori positive based on UBT and serology, but from one of the corpus atrophy patients (Atr2) it was not possible to render any viable H. pylori colonies. Two of the individuals with intestinal metaplasia (Met2 and Met4) were H. pylori negative in both serology and culture even though one of them displayed a positive UBT test (Met4).

Biopsy specimens for RNA extraction were obtained during endoscopy and were snap-frozen in RNAlater (Ambion) and shipped frozen on dry ice to Sweden. For this study, we used one corpus biopsy specimen per individual described in Table 2. The content of the vials was spun down to collect all cells and mucus, and the RNAlater was removed carefully. RNA was extracted using the RNeasy minikit (Qiagen) with the addition of lysozyme and proteinase K to the lysis buffer prior to disruption of the tissues using the TissueLyser (Qiagen). The RNA was subjected to double DNase treatment to ensure that the analysis would not be confounded by DNA contamination, as described previously (38). The integrity of the resulting RNA was assessed using a BioAnalyser (Agilent) with an RNA integrity number (RIN) of >6 as a cutoff, and the RNA concentration was determined by Qubit (Invitrogen) measurement just prior to the library preparation procedure.

Library preparation of cDNA was performed using the strand-specific ScriptSeq v2 RNA sample preparation kit (Epicentre Biotechnologies), including RiboZero eukaryotic rRNA depletion treatment (Epicentre). After library preparation, libraries were quantified, pooled for clustering, and sequenced using the Illumina HiScanSQ instrument, rendering for each sample a total of 20 to 30 million paired-end reads of 100 bp each. The RNA sequencing was performed at the Genomics Core Facility, Sahlgrenska Academy, University of Gothenburg.

To verify the presence of low-abundance Helicobacter rRNA in the clinically H. pylori-negative samples, RNA samples were reverse transcribed using the QuantiTect kit with accompanying genomic DNA (gDNA) depletion (Qiagen). In addition to the corpus samples, RNA was extracted from antral biopsy specimens from the same individuals using the same procedures as described above. For three individuals, Min1, Min3, and Min4, no antrum biopsy specimens were available for RNA extraction.

16S rRNA amplicon libraries were constructed from the cDNA samples using the universal primers 341F and 805R to amplify the V3-V4 region. The libraries were generated in a single PCR step with 25 cycles using a high-fidelity polymerase (Phusion Hot Start II high-fidelity PCR master mix; Thermo Scientific). Subsequent purification was carried out using a magnetic bead capture kit (AMPure; Agencourt) and quantified utilizing a fluorometric kit (QuantIT PicoGreen; Invitrogen). The purified amplicons were normalized to 4 nM and pooled equimolarly or, if the library was of lower concentration, using the same volume as if the library had been 4 nM. Sequencing was performed on the MiSeq platform using v3 chemistry, generating 2 × 300-bp pair-end reads. Negative PCR controls without template were also sequenced to monitor the background from the 16S rRNA PCR amplification.

Raw RNA-seq reads were trimmed using TrimGalore! version 0.3.5 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ ) with a quality cutoff of Q30, Illumina adapter trimming, and removal of reads shorter than 30 bp and reads left unpaired. In addition to this, low-complexity reads were filtered out using PrinSeq version 0.20.4 (39), applying the DUST algorithm (40) to avoid unspecific matches in the alignment analyses. After quality filtering, 17 to 36 Mreads (average 25) remained for each sample.

To analyze the microbial composition of the biopsy specimen tissues, we used the Metaxa2 software version 2.1.1 (20) to extract and classify small-subunit (SSU) rRNA transcripts from the RNA-seq data. Metaxa2 was applied to the quality trimmed and filtered reads using default settings, with the exception of using blast+ rather than the default blastall program and setting the minimum read length to be included in the analysis to 75 bp to ensure a high specificity. Genus annotation in Metaxa2 is done by an internal database, based on the high-quality SILVA database, release 111 (July 2012), which consists of SSU rRNA sequences from the Bacteria, Archaea, and Eukaryota domains (41). The resulting abundance files on the different taxonomic levels were filtered so that only taxa with more than one read per million total reads were used for downstream analysis. The raw genus counts of reads classified as 16S rRNA for each sample are given in Table S2 in the supplemental material. Principal-component analysis (PCA) and other statistical analyses were performed in the R software v3.3.0. Relationships between Helicobacter abundance and the abundance of other genera were assessed using Spearman's rank correlation, and P values were corrected for multiple testing using the Benjamini and Hochberg method (significance cutoff, 0.05). Significant effects of Helicobacter presence (as assessed by traditional tests) on overall community structure on the phylum and genus levels were assessed using metaxa2_uc, part of the Metaxa2 Diversity Tools (42) (default options). In this analysis, the within-group and between-group Bray-Curtis dissimilarities were tested using subsampling for 10,000 iterations (1,000 sampled reads per iteration). Simpson and Shannon diversities as well as richness of genera in the samples were assessed using the R package Vegan version 2.0-10 (43). Significant diversity differences between clinically classified H. pylori-positive and -negative individuals were detected using Student's t test. Regression to human ATP4 gene expression and atrophy scores (see below) was performed using fitting of linear models in R. P values were corrected for multiple testing using the Benjamini and Hochberg method, and corrected P values were considered significant when they were below 0.05.

After trimming, 16S rRNA amplicon reads were trimmed using Cutadapt (44), merged using vsearch –fastq_mergepairs (45), and thereafter processed using Metaxa2 as described above. Raw 16S rRNA read counts classified on the genus level are given in Table S2. PCA plots were constructed in the R Studio software using log-transformed raw counts from Metaxa2. To assess differences between corpus and antrum samples, metaxa2_uc was used to calculate Bray-Curtis dissimilarity between samples and within each corpus/antrum sample pair (additional options, “-g none -r 100 –table T –matrix T”). The within- and between-sample dissimilarities were related to disease stage as well as human ATP4 gene expression.

For identification of Helicobacter pylori transcripts, we created a custom mapping reference. For 20 of the individuals, we had previously whole-genome sequenced one Helicobacter pylori isolate from the antrum and one from the corpus (Table 2). In addition, from one individual we had sequenced two antrum isolates, from one subject two corpus isolates, and from one other subject an antrum isolate only (see Table S3 in the supplemental material). From one patient who was Helicobacter pylori positive based on UBT and serology, we could not obtain any viable colonies in culture (46). From these isolates, draft genomes were prepared by de novo assembly using SPAdes v3.7.0 (47) and annotated as described previously using the prokka annotation pipeline (48) with a recent annotation of strain 26695 (49) as the primary annotation source, further curated with respect to outer membrane proteins (50) and factors involved in DNA transfer and competence (51). A complete list of the annotation is given in Table S4 in the supplemental material. From the annotated draft genomes, FASTA files of nucleotide coding sequences (CDSs) with in total 82,661 CDSs were clustered on 100% identity using USearch v8.0.1623 (52) in cluster_fast mode, rendering a FASTA file of 46,515 centroid sequences. This file was concatenated with a file of all transcripts of the human genome annotation GRCh38, release 79. The combined FASTA file was used to generate a reference for kallisto (53) with a kmer length of 31 bp. Kallisto mapping was performed to the combined reference and generated transcript per million (tpm) counts for each transcript in the reference. To aggregate the CDSs of the different H. pylori strains, counts were first aggregated on gene name, i.e., that they had been annotated to the same source (e.g., a certain CDS in the 26695 file or a specific UniProt number). To further group the different H. pylori CDSs into orthologous groups of genes, we again performed USearch clustering, this time on the protein FASTA output files from prokka, with an identity cutoff of 80%. The H. pylori counts were then aggregated for each such orthologous group to collect all counts.

To determine the extent of atrophic changes in the tissue, we used eight tissue-specific human genes for areas of atrophic changes: ATP4A, ATP4B, GHRL, GIF, CCKRB, PGC, PGA4, and PGA3 (25, 26). Based on the RNA-seq counts per million reads for these genes, every sample received a score on a three-level scale for each gene, where score 1 corresponded to expression lower than the 95% confidence interval (CI) of the expression of that gene across all 30 individuals, 2 corresponded to expression within the 95% CI, and 3 corresponded to expression greater than the 95% CI. The overall molecular atrophy score was determined using the median value of the scores, where a median of 1 was considered high-grade atrophy, a median of between 1.5 and 2.5 intermediate atrophy, and a median of 3 low-grade atrophy. Atrophy scores, mean ATP4 (ATP4A and ATP4B) expression for each sample, and the correlation between the two are shown in Fig. S1 in the supplemental material.

On the kallisto counts, we applied a cutoff a median tpm higher than 1 for a gene to be included in further analysis, resulting in a list of 285 H. pylori genes (see Table S5 in the supplemental material). Since the reads mapping to H. pylori accounted for only a small fraction of the total reads, the kallisto tpm for each gene was normalized to the total tpm mapping to H. pylori genes in that sample. Samples from individuals clinically H. pylori negative were excluded from differential expression analysis due to radically lower total H. pylori counts which otherwise would confound the analysis. Differential expression analysis was performed using limma v3.28.20 (54, 55), comparing the normalized expression values (i) between pathology groups, (ii) between groups of molecular atrophy, and relative to (iii) mean ATP4 expression (ATP4A and ATP4B) and (iv) mean atrophy score.

The RNA-seq data set has been deposited to the ArrayExpress public repository (www.ebi.ac.uk/arrayexpress ) under accession number E-MTAB-3689 .