Difficulties of Enterobacteriaceae genome annotation in deciphering gastrointestinal microbiome datasets obtained by 16S rRNA gene amplicon sequencing

Introduction

The human microbiome is a complex ecosystem consisting of the genetic material of more than 1 trillion microorganisms living inside the host [1]. Advances in massively parallel DNA sequencing technologies led to a significant increase in knowledge about microorganisms found in the natural environment, food systems, and the human body. In particular, sequencing of the 16S rRNA gene amplicon is a key method for studying the diversity and phylogeny of bacteria. This approach allows the simultaneous identification of most bacteria in complex microbial communities. Although the analysis of 16S rRNA gene diversity presents significant prospects for the study of bacteria from various habitats, the problems of standardizing approaches to sample preparation, DNA sequencing, and data analysis for obtaining reliable information on the composition, structure, and diversity of bacterial populations remain [2].

In studies of the microbiota of most ecosystems or habitats, identification at the species or strain level increases the ecological and/or clinical significance of the results, compared with identification at the genus level. For example, identification at the species level is often critical for host-associated microbial communities, as these communities often include commensal and pathogenic species of the same genus. In addition, some bacterial taxa include species that are specific to several localities and inhabit strictly defined niches of a given environment [3]. High-throughput sequencing of near-full-length 16S rRNA gene fragments (e.g., PacBio single-molecule circular consensus sequences, real-time sequencing, and whole genome sequencing) was expected to improve detection accuracy to species and strain levels. However, due to the greater availability and lower cost of ribosomal amplicon metasequencing, molecular epidemiological studies of the bacterial microbiota of humans, other animals, plants, and the environment are currently conducted on a population scale (i.e., thousands of samples) [4].

Bacterial reference databases with wide phylogenetic diversity, such as SILVA, RDP, and Greengenes, play a key role in data analysis [5-7]. Nevertheless, the taxonomic annotation of 16S rRNA gene sequences is incorrect in 10-17% of cases [8]. SILVA and RDP are regularly updated and represent extensive and complete libraries of 16S rRNA gene sequences from all studied habitats. In contrast, the Greengenes database was last updated in 2013 [7]. Taxonomic assignment using a reference database for large arrays of sequences and pipelines of metagenomic data processing platforms is associated with a certain percentage of misidentifications. Thus, correction of the taxonomic position is required using phylogenetic analysis and sequences of type strains.

Our previous studies described the frequency and structure of obesity-associated cardiometabolic risk factors in a cohort of children and adolescents [9], changes in the biochemical status of obese youths [10-12], and analyzed the intestinal microbiota [13, 14]. Species of Bifidobacterium are present in the gastrointestinal tract of a healthy person, and a change in the number and composition of their species is a sign of intestinal dysbiosis. Features of the gut microbiome associated with obesity include a decrease in Bifidobacterium counts and reduced phylotype diversity. Bifidobacterium species were also examined, for which accurate species identification could not be performed using V3-V4 variable regions or standards for amplicon sequencing [15]. In addition, gut microbiome dysbiosis in obese adolescents was associated with altered species spectrum of enteric bacteria. Members of this family include both normal intestinal microbes (Escherichia coli) and opportunistic pathogens, such as Klebsiella spp. In this regard, accurate species identification is essential for obtaining complete information about the composition of representatives of this family in the intestinal microbiocenosis [16]. The goal of our study was to evaluate the correctness of the taxonomic identification of enteric bacteria by means of using the SILVA 132 reference database. The present study was conducted on samples of patients with normal body weight and obesity, for whom the results of both bacteriological analysis [16] and amplicon sequencing data were available.

Material and Methods

When studying the educational preferences of students studying at the universities of the Stavropol Krai and practicing Islam, 1,500 questionnaires were evaluated, which was a sufficient reference sample.

Brief description of the research material

The study was approved by the Ethics Committee of the Research Center for Family Health and Problems of Human Reproduction (RC FHPHR), Protocol No. 6 of 21 December 2015. The RC FHPHR previously studied the intestinal microbiome of adolescents with normal body weight and obesity [13, 15]. The general characteristics of the patients are presented in Table 1. Laboratory studies were carried out using standard operating procedures (SOP), IHMS_SOP 03 V2 and IHMS_SOP 06 V2, developed in the course of implementing the project of the international consortium, International Human Microbiome Standards. Amplicon analysis of the V3-V4 variable regions in the 16S rRNA gene was performed at Novogene (China). Primary data were deposited in the NCBI Sequence Read Archive (SRA) under accession numbers SRR11006336-SRR11006339, SRR11006343, and SRR11006351-SRR11006388 (PRJNA604466) [14]. The amplicon libraries were processed using the algorithms of the QIIME2 2019.4 platform [17].

Table 1. General characteristics of study participants

Bioinformatics data processing, phylogenetic and statistical analysis

We used SILVA 132 reference database for taxonomic assignment. To elucidate the phylogeny of amplicon sequence variants (ASVs), identified as Enterobacteriaceae sequences, sequences of the complete 16S rRNA gene of type strains of all Enterobacteriaceae family species were used.

A total of 63 nucleotide sequences were included in the tree, identified by comparison with the SILVA 132 reference database as belonging to the family Enterobacteriaceae, along with typical bacterial strains of this family. Multiple alignment and phylogenetic tree construction were performed using MEGA X software [15]. DNA sequence alignment was originally performed using the MUSCLE algorithm with default settings. The alignment was then visually checked to correct obvious alignment errors and remove areas of questionable alignment. The maximum likelihood method was employed to construct the phylogenetic tree. Statistical support for phylogeny was implemented using bootstrap (1,000 iterations). Bootstrap values ​​≥85% were considered highly supported, values ​​of 75-84% were classified as moderately supported, and values of ​​50-74% were categorized as poorly supported. Values ​​<50% were not specified [19].

Results

Basic statistics for library analysis

Molecular genetic analysis performed 2,590,453 reads. The number of reads per sample ranged from 52,945 to 77,290. A total of 2,890 phylotypes (ASVs) were identified, and the range per sample was 342-564. Depth of sequencing evaluation via the Michaelis-Menten approximation showed that the composition of the microbiome at the ASV level was underestimated by an average of 2.04%.

General characteristics of the representation of Enterobacteriaceae

The Enterobacteriaceae content in the total microbiome ranged from 0.76 to 23.45%, and there were no significant differences between the control and obesity groups.

A total of 63 ASVs were assigned to the Enterobacteriaceae family. Sensu the taxonomy of the SILVA reference database, this family is represented by the genera Citrobacter, Enterobacter, Klebsiella, Proteus, Raoultella, Serratia and two undifferentiated groups (Escherichia-Shigella and Hafnia-Obesumbacterium) in the adolescent gut microbiome. In addition to these genera, ASVs have also been examined that could not be assigned to any genus (unidentified Enterobacteriaceae). Among all ASVs, only 65be was present in all samples, and this ASV was identified by SILVA as an Escherichia-Shigella phylotype (Figure 1). None of the ASVs exhibited significant differences in size between the obesity and control groups (Supplementary materials).

Figure 1. Frequency heatmap of isolated ASVs in gastrointestinal microbiomes in obesity and control groups of youths.

Taxonomy of Enterobacteriaceae

At the time of writing, this family was represented by 32 genera and 124 species. The following species had subspecies: Enterobacter cloacae, Enterobacter hormaechei, Klebsiella pneumoniae, Klebsiella quasipneumoniae, Klebsiella variicola, and Salmonella enterica. The December 2017 update of the SILVA reference database contains the genera Proteus, Serratia, Hafnia, and Obesumbacterium, which are not currently included in Enterobacteriaceae according to LPSN (https://lpsn.dsmz.de/). Proteus was moved to Morganellaceae [17], Serratia to Yersiniaceae, Hafnia and Obesumbacterium to Hafniaceae, and Pantoea to Erwiniaceae family.

Phylogenetic analysis revealed that the genera belonging to Enterobacteriaceae, according to the studied fragments, were polyphyletic, and they formed mixed clades (Figure 2). Only the genera Cedecea, Gibbsiella, Phytobacter, Pseudocitrobacter, Franconibacter, Mangrovibacter, Izhakiella, Rosenbergiella, Trabulsiella, and some groups of Citrobacter, Klebsiella, and Kosakonia were monophyletic and formed clades with good statistical support. Twenty-five ASVs were identified to the generic level, whereas 11 ASVs were assigned to the Escherichia-Shigella cluster. The identification matched for 22 ASVs, whereas the reclassification affected 14 ASVs. Twenty-seven ASVs remained identified only at the family level (Enterobacteriaceae).

Figure 2. Phylogenetic tree of the studied ASVs and sequences of type strains of Enterobacteriaceae. A – The outer group is marked with a diamond-shaped marker. Grey clusters denote monophyletic genera that do not include ASVs. Black markers denote ASVs that were assigned to a specific genus or taxonomic group. White markers denote unidentified ASVs. The scale is five substitutions per 100 base pairs (bp). B – Extended group 1. Gray clusters denote monophyletic genera not including ASV. Black markers identify ASVs that were assigned to a specific genus or taxonomic group. White markers denote unidentified ASVs. The scale is one substitution per 100 bp. C – Extended group 2. Grey clusters denote monophyletic genera that did not include ASVs. Black markers denote ASVs that were assigned to a specific genus or taxonomic group. White markers denote unidentified ASVs. The scale is two substitutions per 100 bp.

Sequences identified as taxa that do not currently belong to Enterobacteriaceae formed separate clades (Nos. 1-4, Figure 2). The identification of ASV 5535 as a member of Proteus was not confirmed. Clade #1 contains sequences identified as Hafnia-Obesumbacterium and unidentified (UI) Enterobacteriaceae. The remaining clades included both sequences characterized by SILVA as representatives of reclassified genera and re-identified by phylogeny. Clade #3 contained Escherichia-Shigella and Pantoea sequences, while clade #4 included Escherichia-Shigella and Serratia sequences. All of them featured medium to strong bootstrap support. Clade #2 contained sequences identified as Enterobacter, Pantoea, and UI Enterobacteriaceae, and branch nodes had weak bootstrap support. For these sequences, a search for the nearest homologs was performed using BLAST software (Table 2). SILVA identification was identical for 11 ASVs. Among the mismatched were representatives of the genera Erwinia, Pantoea, Serratia and Yersinia.

Table 2. Search for the closest homology using BLAST

The reclassified taxa accounted for 0.009-3.6% of the total microbiome. After the taxonomy correction, the content of several genera in the gut microbiome was changed, including Enterobacter, Klebsiella, and the Escherichia-Shigella group. The frequency distribution of other taxa was sporadic (Figure 3). Analysis of the overall frequency of ASVs with the same generic identification did not reveal significant differences between obese and control groups in counts for any genus.

Figure 3. Frequency distribution (% of the total count in the family) of Enterobacteriaceae family representatives in control and obesity groups.

Hence, the phylogenetic analysis of this family for the studied V3-V4 fragment was complicated by the polyphyly of some genera. The genus of half of the ASVs could not be specified.

Discussion

Dysbiosis is mainly associated with an increased number of pathobionts, such as Escherichia and Klebsiella spp. caused by a reduction in the number of taxa with useful metabolic activity, including lactobacilli and bifidobacteria. Dysbiosis is also associated with a decrease in biodiversity, i.e., a reduction in the number of microbial species present in the microbiome and lower complexity of the microbial community [21]. As for gastrointestinal microbiota in obese and normal-weight children, the former category has fewer counts of Bifidobacterium and higher counts of E. coli. Studies have shown that a high number of bifidobacteria in infancy and adulthood protects against obesity [22]. It was also revealed that with a decrease in body weight in children achieved by modifying their diet, the numbers of Bifidobacterium and Lactobacillus increases, while the number of enterobacteria decreases [23]. In inflammatory bowel disease, there is an increase in proteobacteria, viz. intestinal bacteria, including the opportunistic pathogens E. coli and K. pneumoniae, which increases mucosal inflammation and the risk of infections. Many studies have described a decrease in the number of bifidobacteria and lactobacilli and an increase in the numbers of Enterobacter in patients with irritable bowel syndrome (IBS) and diarrhea. Other researchers linked IBS with Campylobacter, Yersinia, Salmonella, Shigella and E. coli. The heterogeneity of the results is explained by the variety of methods used to determine the microbiota and the different criteria for enrolling patients [21]. While there is controversy as to which types of bacteria are associated with being overweight, some specific genera and species of bacteria appear important.

It is known that even the complete sequence of the 16S rRNA gene has a low low discriminatory power. Branching of genera and species within this family during phylogenesis, based on the 16S rRNA gene, has a significant stochasticity depending on the used algorithms and analyzed bacteria [20]. According to the results of some studies, it can be said that the entire Enterobacteriales order is generally characterized by polyphyletic branching and the absence of connected monophyletic groups [20, 24]. The used fragment is not optimal for phylogenetic analysis of this family. Different parts of the genome may have distinct phylogenetic similarities to other taxa. In other words, a group can be monophyletic for some parts of the genome and simultaneously paraphyletic for other parts. In analytical results, this may reflect either analytical ambiguity or actual phylogenetic inconsistencies.

The composition of the microbial community depends on such factors as the time of breastfeeding, family history (health status of the mother and other family members), and the dominant component of the microbial community that determines the human enterotype. The health status of the mother during pregnancy (past inflammatory and infectious diseases) can affect early (intrauterine) colonization of the child’s body with bacteria, such as Enterobacter, Enterococcus, Lactobacillus, Photorhabdus and Tannerella [25]. The main stage of colonization of the child’s body by symbiotic bacteria occurs at the time of birth. The mode of delivery largely determines the future composition of the microbiome. The microbiome of children born by caesarean section is very different from the microbiome of children born by vaginal delivery. Breastfeeding is the second step in the colonization of the baby’s intestines after birth. The method (artificial or natural) and the time of feeding strongly influence the composition of the intestinal microbiome, determining the dominant and minor components of the community. Feeding has the greatest impact on the diversity of representatives of the genus Bifidobacterium [25]. The dominant component of the community influences the minor components. There are several approaches to typing the intestinal microbial community according to the dominant component. The first approach involves the use of partitioning around medoids (PAM) and dividing them into three groups (Bacteroides, Prevotella and Ruminococcus). The second approach is based on the Dirichlet multinomial mixtures (DMM) and gives a division into 4 groups (Ruminococcaceae [R], Prevotella [P], Bacteroides 1 [B1] and Bacteroides 2 [B2]). In some cases, an additional enterotype is identified with a predominance of representatives of the Enterobacteriaceae (H) family – but, as a rule, it is associated with the presence of inflammatory diseases, alcohol dependence, or other ailments. Enterotypes do not depend on gender, age, ethnicity or geography. Rather, they depend on the characteristics of long-term nutrition [26].

Conclusion

Most taxa were characterized by the presence of a single sample; no dependence on division into groups was observed. However, it can be assumed that a more detailed study of taxonomic diversity, taking into account factors, such as enterotype, duration of breastfeeding and family history, may reveal differences in the frequency distribution. Future studies should also include the analysis of these samples using the whole genome sequencing technology due to such type of information retrieval analysis on a larger scale.

Conflict of interest

The authors declare that they have no conflicts of interest.

Funding

This study was conducted with the financial support by the Council for Presidential Grants of the Russian Federation (NSh- 3382.2022.1.4).

Ethical approval

All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with 1964 Declaration of Helsinki and its later amendments or comparable ethical standards.