Computational Identification of Specific Genes and Those Functions on the Gut Microbiota

1. Introduction

Humans coexist with a diverse array of bacteria, forming a complex microbial community that has co-evolved with its host over millennia [1]. These symbiotic bacteria are instrumental in maintaining human health, particularly within the digestive system, where they contribute to metabolic processes and the maturation of the immune system. The relationship between humans and their microbiota is characterized by intricate interactions that sustain homeostasis and confer physiological advantages. However, disturbances to this delicate equilibrium may result in a range of immune, inflammatory, and metabolic disorders[2]. The extensive genetic diversity observed among bacteria, which significantly surpasses the genetic variation found within the human population, carries important implications for human health and disease. An understanding of bacterial polymorphisms and their effects on the human host provides valuable insights into pathogens, cellular processes, and potential therapeutic interventions [3].

So far, the two principal methodologies, amplification sequencing, and shotgun sequencing, that have been employed to analyze the intestinal microbiota and all DNA present in the sample. On the other hand, there are some limitations to shotgun sequencing. High costs and data scalability are major issues[4]. When assigning genes to specific bacterial taxa, a significant barrier arises when the query sequence lacks a corresponding reference in existing databases. This absence of taxonomic context complicates the classification of metagenomic data [5]. Limitations of short read lengths, sequencing errors, and lack of taxonomically related reference sequences make metagenomic shotgun sequence classification unsuccessful. [6]

As indicated by prior research, PICRUSt2 serves as a predictive metagenomic analysis tool that estimates the constituent genes based on amplicon analysis data. This software predicts the functional composition of metagenomes utilizing 16S rRNA gene sequences and reference genome databases, allowing for accurate estimations of gene family abundances across various microbial communities[7]. Tax4Fun2 is an R package that predicts functional profiles and functional gene redundancy in prokaryotic communities from 16S rRNA gene sequences[8]. However, this method is unable to yield results for genes with unknown functions and lacks sensitivity regarding the reliable detection of functional changes associated with variations in bacterial populations. It is estimated that a significant proportion of proteins within prokaryotic genomes—ranging from 20% to 50%—remain uncharacterized.

Various computational methods have been developed to predict the functions of hypothetical proteins in microbial genomes. These approaches integrate multiple bioinformatics tools and databases to improve prediction accuracy. Common techniques include sequence similarity searches using BLAST, domain identification with InterProScan and Pfam, and ortholog detection using COGs[9,10,11]. Functional interaction networks have also been employed to predict protein functions and analyze their properties [12]. These methods have been applied to predict functions for hypothetical proteins in important pathogens like Mycobacterium tuberculosis and Vibrio parahaemolyticus, revealing potential roles in virulence, pathogenicity, and host adaptation[11,12].

In response to these limitations, we have developed a novel methodology for predicting the functions of bacterial proteins by integrating database searches for functional annotation with amplicon analysis. This approach enables the prediction of functions for proteins of unknown function by considering sequence homology and protein interactions. Furthermore, it can anticipate the implications of variations in bacterial populations by assessing the predicted functions of identified bacteria. Previous studies have suggested that exposure to ultrasound can increase and decrease bacteria[13]. It is also known that mice use ultrasonic vocalizations to communicate their emotional states and to engage in social interactions.[14] In this study, we also applied our novel methodology to analyze the changes in the intestinal flora of mice that exhibited heightened pain sensitivity due to exposure to ultrasonic frequencies.

2. Materials and Methods

2.1. Sample Preparation

2.1.1. Animals

Adult C57BL/6N mice (male, 7 weeks old) were obtained from Japan SLC (Shizuoka, Japan) and used for behavioral and biochemical studies. The mice were housed in the Laboratory Animal Research Center at Tokyo University of Science under standard conditions (23°C ± 1°C, 12-h light/dark cycle, lights on at 8:00 a.m.). The mice were housed in standard mouse cages with sawdust bedding and had free access to food and water. This study was conducted by the “Basic Guidelines for the Implementation of Animal Experiments at Academic Research Institutions, etc.” under the jurisdiction of the Ministry of Education, Culture, Sports, Science and Technology. All protocols were approved by the Tokyo University of Science Animal Experiment Committee (Approval Nos. Y21036 and Y22028).

2.1.2. Recordings of Cries

The mouse was fixed with a fixation device, and the mouse’s tail was pinched with a 5 cm long artery clip to record its cries (or vocalizations) during pain responses. Its cries were recorded with an UltraSoundGate (Avisoft Bioacoustics, Nordbahn, Germany). The ultrasonic regions (over 20 kHz) of the obtained sound sources were used in this experiment as psychological stress sound sources (sound stress).

2.1.3. Von Frey Test

Naive mice, different from those used for the recordings of cries, were exposed to sound stress at 80 dB for 4h (day 0) in a soundproof box, and tactile thresholds were measured after 1, 3, and 7 days. The von Frey filaments evaluated the sound stress-induced hyperalgesia.

2.1.4. 16S rRNA Amplicon Analysis

The contents of the large intestine of mice exposed to acoustic stress for 4 hours a day were collected and subjected to 16S rRNA amplicon analysis (Illumina 16S Metagenomic Sequencing Library Preparation Protocol[15]). 16S rRNA amplicon was conducted by Hokkaido System Science Co. (Hokkaido, Japan).

2.2. New Application Methods to Identify the Specific Genes from Amplicon Data

2.2.1. Getting Data from NCBI RefSeq and Creating a Database

Using the NCBI Datasets command line tool (CLI) v12[16], we obtained the FAA file of type strain data for all bacteria from NCBI RefSeq[17]. After removing duplicate data, we created a database via BLAST+[18]. We have created a database (BSDB0 to BSDB9) that divides these data into 10 parts. We used data obtained on December 6th, 2024. We accessed the NCBI FTP site to retrieve the FASTA file of proteins produced by bacteria in the RefSeq database and created a database (DB2) using BLAST+. We used data obtained on April 20, 2023.

2.2.2. Classification of Bacterial-Specific Amino Acid Sequences Using BLASTP Search

We focused on bacteria species whose presence rate changed more than 1% between healthy mice as control and mice in psychological stress conditions. Then we analyzed the proteins produced specifically by bacteria that showed significant changes in relative abundance. The genome data for the standard strains of the relevant bacteria were obtained from NCBI Genome in FASTA and TSV formats.

Standard Protein BLAST[19] (BLASTP) searches were performed on the acquired amino acid sequences via BLAST+, and sequences that matched only a single bacterial species were classified as bacterial-specific sequences. The BLASTP search parameters were as follows: an expected value threshold of 1e-10, BSDB0~9, and a maximum target sequence length of 20×10DB. All other parameters were left in their default settings.

We integrated the search results and extracted amino acid sequences only present in related bacteria. We also classified amino acid sequences that had not been characterized for function.

2.2.3. Functional Prediction by Sequence Analysis Using PSI-BLAST and InterPro Search

In Section 2.2.2, we set out to predict the functions of proteins with unknown roles by conducting sequence analysis using PSI-BLAST[20]. We performed PSI-BLAST searches focused on these uncharacterized proteins to forecast the functions that lacked annotations. The parameters for the search were as follows: Database: DB2, expected value threshold: 0.05, PSI-BLAST threshold: 0.0001, Number of iterations: 3, Maximum target sequences: 100, with all other settings maintained at their default values.

Furthermore, we employed InterProScan[21] to investigate domains by utilizing InterPro[22] to detect the presence of specific sequence domains. Subsequently, the results from both InterPro and PSI-BLAST were combined and visualized graphically. We filtered the outcomes of our sequence function predictions based on established criteria for acceptance.

((Positive > 40% AND Identity hits / Subject Total Length > 20%) OR (InterPro hits = “yes”)) (1)

2.2.4. Functional Prediction of Protein-Protein Interactions Using Clusters Based on Gene Neighborhood

We sought to predict the functions of proteins whose roles remained undetermined in Section 2.2.3, utilizing genetic information as a foundation. Prokaryotic genomes, particularly those of bacteria, typically exhibit a circular structure wherein adjacent genes can influence each other’s expression simultaneously. This configuration is referred to as the operon structure, and it is well established that proteins encoded by genes within an operon are likely to share similar functions.

Recent studies have introduced innovative strategies for predicting protein functions and interactions through protein-protein interaction (PPI) networks[23]. Prokaryotic genomes, particularly those of bacteria, typically exhibit a circular structure wherein adjacent genes can influence each other’s expression simultaneously[24]. The STRING database is a comprehensive resource for protein-protein interactions, integrating experimental data, computational predictions, and literature-based information[25]. The STRING database uses gene neighborhoods of up to 300 bases as a factor in predicting protein interactions[25].

We identified the genes associated with each protein to utilize this information for inferring protein function. We then analyzed the locations of these genes, creating gene clusters with known functions that extended 100 bases before and after each gene. Because only the gene neighborhood is used as a predictor, the value is set more strictly than that used in STRING. We assessed whether the protein-coding genes that could not be inferred in Section 3.2.4 were assigned to the established gene clusters. The protein with the closest inter-gene distance was chosen as the representative sequence for those transferred.

2.2.5. Gene Ontology Enrichment Analysis by DAVID

We created a list of the extracted proteins identified in 2.2.2, plus the adopted sequences identified in 2.2.3~2.2.4. We used DAVID (The Database for Annotation, Visualization, and Integrated Discovery) [26] to perform gene ontology enrichment analysis on these lists. The list of extracted proteins identified in 2.2.2 and the adopted sequences identified in 2.2.3-2.2.4 have been added. The list was submitted using the RefSeq protein ID.

3. Result

3.1. 16S rRNA AmpliconAnalysis Result

We obtained data on the percentage of bacterial composition within each sample. Bacterial data from healthy mice were treated as Ctrl, and bacterial data from ultrasonically irradiated mice and evaluated for hyperalgesia with the von Frey test were treated as USV (Figure 1Figure S1).

We identified seven species of bacteria at the species level. We then tested for differences in the sample proportion for these species and found that Mucispirillum schaedleri (M. schaedleri) had increased significantly. (Table 1Table S1)

We focused on M. schaedleri, which was present in more than 1% of the samples in the control and ultrasound-exposed groups and significantly increased by ultrasound exposure.

3.2. Analysis of Bacterial Genome Data Result

For the M. schaedleri genome data, we applied the developed analysis method and obtained the results as follows.

3.2.1. Classification of Bacterial-Specific Amino Acid Sequences Using BLASTP Search

We obtained the type strain genome annotation of M. schaedleri in TSV format. We received the amino acid sequence in fasta format. 2188 genes were annotated, of which 2118 were protein-coding genes. After eliminating duplicates, 2065 types of proteins were reported. As a result of performing a BLASTP query on the 2065 proteins produced by M. schaedleri, it was found that 303 of the proteins had amino acid sequences characteristic of M.schaedleri. As a result of performing a gene ontology enrichment analysis on these using DAVID, it was found that there are many proteins with transmembrane functions present in the outer membrane. Of the 303 proteins, 253 were “hypothetical proteins” with unknown functions. RefSeq Protein IDs were obtained for the 50 proteins with known functions. Functional prediction was performed for the remaining 253 proteins using method 2.2.2 (Table S2).

3.2.2. Functional Prediction by Sequence Analysis Using PSI-BLAST and InterPro Search

After analysis of Method 2.2.2, 46 of the 253 sequences met the filtering criteria. These were added to the analysis, and as in Result 3.2.2, many terms related to membrane proteins were found. (Table S3)

3.2.3. Functional Prediction of Protein-Protein Interactions Using Clusters Based on Gene Neighborhood

After the analysis was run, 439 clusters were generated. As a result of mapping the target protein-coding genes to the clusters in 2.2.3, it was found that 114/144 genes belonged to known functional clusters. Of these, 109 clusters contained the protein-coding gene under analysis. In total, 136 protein-coding genes could belong to clusters of known functions. The most proximal genes in the clusters were extracted from these 136 genes, and 116 genes were found to be in the clusters, excluding duplicates (Table S4).

3.2.4. Gene Ontology Enrichment Analysis by DAVID

In addition to the 303 proteins classified in 3.2.1, we created a list that included the predicted results of proteins that were considered to have similar functions in 3.2.2 and 3.2.3 (Table S5). We performed a Gene Ontology enrichment analysis using DAVID for these. As a result, we found that in addition to membrane proteins, the functions of Peptidase M23, Ankyrin repeat, and Tetratricopeptide domain were enriched (Figure 2).

4. Discussion

M. schaedleri, which increased in the ultrasound irradiation group, is a Gram-negative anaerobic bacterium with a spiral flagellum that is mainly observed in the intestines of mice. It has been reported that it is more abundant in the intestinal mucosa than in feces. [27,28] Although it is present in small quantities in human feces, the presence of 16S rRNA genes has been confirmed in colon mucus, suggesting its presence. M. schaedleri suppresses the growth of other bacteria by competing with pathogenic microorganisms such as Salmonella enterica serovar Typhimurium. [29] On the other hand, M. schaedleri can hurt the host by proliferating excessively. M. schaedleri has been reported to increase inflammatory conditions in the intestine. [27] The mechanism by which M. schaedleri increases in these inflammatory conditions is thought to be that M. schaedleri can inhabit environments rich in nitrate in inflammatory conditions, and that it has a rich supply of genes that adapt to oxidative stress. [27,30,31,32,33]It has also been reported that M. schaedleri increases due to social stress.[34] In this study, no outward signs of inflammatory conditions were observed. The results of the brain DNA microarray suggest that inflammation is occurring in the brain. It has been reported that neuroinflammation in the brain causes chronic pain and hyperalgesia. [35]In this study, outer membrane proteins were found to be abundant in proteins specific to M. schaedleri.[36] Gram-negative bacteria can release membrane proteins in the form of outer membrane vesicles as carriers of toxins, DNA, etc. OMVs can cause inflammatory diseases due to the presence of lipopolysaccharide (LPS) on their surface.[37] LPS causes inflammatory diseases via the innate immune system, such as TLR4[38]. Peptidase M23 can destroy the bacterial cell wall by hydrolyzing peptidoglycan and eliminating the bacteria.[39] This process has been demonstrated to cause the diffusion of LPS in the bacterial outer membrane into the surrounding environment, resulting in inflammation of the host’s immune cells and nerve cells. Furthermore, the other report showed that inflammation caused by LPS may increase pain sensitivity.[40] In addition, it has been confirmed that ultrasound can increase or decrease the number of bacteria. The results of this study revealed that ultrasound irradiation can directly or indirectly increase the level of M. schaedleri.[41] The results suggest that the action of M. schaedleri peptidase M23 increases the concentration of LPS in the host, causing inflammation in both the host’s immune cells and nerve cells. This LPS-induced inflammation likely contributes in part to the increased pain sensitivity observed in the host. Recent studies have shown that some bacteria contain Peptidase M23 in outer membrane vesicles. [42] This also suggests that peptidase M23 may contribute to LPS secretion via outer membrane vesicles. ANK repeat proteins are translocated into host cells via the type IV secretion system in prokaryotes. The translocated proteins are associated with evading the immune system and host cell-mediated gene expression. Class II chaperones of the type III secretion system, which are TPR-containing proteins, inject bacterial effectors into host cells by creating pores in the cell membrane of eukaryotic proteins.[43] These results also suggest that M. schaedleri affects host cells.

Our method succeeded in predicting and assigning annotations to 162 protein genes of the 253 genes whose functions were previously unknown. The most distinctive feature of our method focused on unique genes which only exist on that bacteria genome sequence. In our method, it was very helpful to predict the features for unknown or hypothetical genes by use of genes in neighborhood locus with annotations.

DNA microarray data of brain tissues for those mice also suggested that deferentially expression genes were related to inflammatory diseases, such as Ptgs2, Lcn2, and Cxcl1. We then performed a Gene Ontology enrichment analysis using DAVID for deferentially expression genes and extracted the term GO:0032496 response to lipopolysaccharide. This result suggests that there is a causal relationship between the data obtained from gut microbes and the data in the host brain via LPS.

We found that the proteins produced by M. schaedleri are registered in the AlphaFold Protein Structure Database.[44] Although our proposed methods only use primary structures of proteins, our method can be expanded to utilize 3D structure information.

We believe that further insights can be gained by predicting functions based on these predicted structures using structural comparisons.

We applied our method to analyze M. schaedleri. However, we also identified more than 20 types of bacteria classified only at the genus or family level. Therefore, as future trials, we will apply this method to unique genes on genus and family-level taxa all at once.