Eighteen sun-dried samples of kinema were aseptically collected from different places of India, Nepal, and Bhutan (Supplementary Table 1) in a pre-sterile container, sealed, labeled, and brought to laboratory for further analysis.

Ten grams of each sample were homogenized with 90 ml sterile 0.1 M phosphate buffer (pH 6.4) and kept undisturbed for 5 min to settle down the sample debris. Homogenate was collected, and genomic DNA was extracted using Nucleospin Food kit (Macherey-Nagel GmbH & Co., KG, Düren, Germany) following the manufacturer’s instructions. Quantification of DNA concentration was checked using a spectrophotometer (Eppendorf, United States), and quality of DNA was observed in 0.8% agarose gel electrophoresis followed by visualization under Gel Doc EZ imager (Bio-Rad, United States). Following the method described by Ray et al. (2019), the DNA from 18 samples of three countries were pooled in equal quantity (pooled from six initial samples of each country i.e., India, Nepal, and Bhutan), and the mixture was considered as one sample from each country and were used for metagenomic analysis.

The MinION gDNA sequencing Ligation Kit SQK-LSK109 (Oxford Nanopore Technologies, Oxford, United Kingdom) was used to prepare the metagenome libraries (Sevim et al., 2019) according to the manufacturer’s instruction. Genomic DNA (10 μg) was sheared using the Covaris g-tubes to produce the fragments of the size > 10 kb (Covaris Inc., MA, United States). The sheared genomic DNA with appropriate size was then selected using the BluePippin instrument (Sage Science, Beverly, MA, United States). End-repairing was performed on the sheared DNA using the NEBNext FFPE DNA Repair kit following the manufacturer’s protocol (New England BioLabs, Ipswich, MA, United States). Qubit HS DNA kit was used to quantify the gDNA. Ligation Sequencing Kit SQK-LSK109 was used for the clean-up and adapter ligation. The ligated sample was purified and eluted using the AMPure XP beads (Beckman Coulter Inc., United States) along with wash buffer and kit. Libraries were then sequenced (1D sequencing) on MinION MK1 device using R9 flow cell chemistry. The device was controlled through the MinKNOW software version 1.0.5 (Oxford Nanopore Technologies, Oxford, United Kingdom); Metrichor platform of ONT (Oxford Nanopore Technologies) was used to perform the 1D base calling.

LC-MS coupled to a Q-Exactive Orbitrap (Dionex Ultimate 3000, Thermo Fisher Scientific, United States) was used for the metabolite profiling of kinema samples (three samples each from India, Nepal, and Bhutan) following the method of Ramakrishnan et al. (2016). A 200-mg amount of each powdered sample was mixed with 1 ml of 80% methanol and sonicated in ice, the mixture was centrifuged at 14,000 rpm for 5 min at 4°C, and the supernatant was filtered using a 0.2-μm PTFE filter. Then, 10 μl of each sample was injected into the LC-MS system coupled to a Q-Exactive Orbitrap, and chromatography was performed on a hydrophilic interaction liquid chromatography (HILIC) (5 μ, 150 mm × 4.6 mm, Phenomenex Luna) with a flow rate of 0.4 ml/min at 40°C. The mobile phase A contained 10 mM ammonium acetate with acetonitrile in the ratio of 1:1 (0.1% FA), and phase B contained acetonitrile (0.1% FA) for positive mode. Whereas, for the negative mode, the mobile phase A contained 5 mM ammonium acetate in water, and the phase B contained 5 mM ammonium acetate in water with acetonitrile in the ratio of 1:9. The following gradient gave optimal resolution and was used in both ionization modes (0–2 min: 100% B, 2–15 min: 100–90% B, 15–25 min: 90–80% B, 25–30 min: 80–75% B, 30–35 min: 75–20% B, 35–40 min: 20–0% B, 40–45 min: 0% B, 45–45.1: 0–100% B, and 45.1–55 min: 100% B) at a 400-μl/min flow rate. The Q-Exactive Orbitrap (Thermo Fisher Scientific, United States) was set up for data acquisition in the full scan/data-dependent scan (FS/DDS) mode in a mass range of 70–1,050 m/z, alternating between MS and MS/MS scans. The MS operating conditions for all samples were as follows: spray voltage, 4,000 V (2,500 V for negative); vaporizer temperature, 320 °C; sheath gas flow rate of 30 arbitrary units (40 for negative); and auxiliary gas flow rate of 10 arbitrary units. Injector settings were as follows: 0–2 min: waste, 2–45 min: load, and 45–55 min: waste. The spiked reserpine and taurocholate were used as internal standards and later used for intensity normalization of data. The raw data files were imported into SIEVE 2.2 for the generation of peak list and component extraction (Hao et al., 2018), and the Human Metabolome DataBase (HMDB), KEGG (Zhang B. et al., 2018), and PlantCyc (Zhang et al., 2010) were used for possible identification of the compounds. The coefficient of variation (CV) was calculated from the pooled quality control sample data, and the data with CV_QC > 20% were removed.

Nucleotide diversity (Pi) per site of pooled samples was 0.72668 (India), 0.7243 (Nepal), and 0.72687 (Bhutan), respectively. The neutrality indices were performed using two tests: Tajima’ D with values of 1.363 (India), 1.347 (Nepal), and 1.364 (Bhutan), respectively, and Fu’s Fs statistic representing the values of -46.669, -31.869, and -43.373 of kinema from India, Nepal, and Bhutan, respectively. Pooled DNA samples were not significantly (p > 0.1) different in terms of sequence variation within the sample.

A total of 1,394,094,912 bp with an average of 464,698,304 ± 120,720,392 bp was generated from kinema metagenome. In comparison, among the samples, kinema from Nepal has the longest contig of 533,203 bp, whereas kinema from Bhutan has the shortest contig of 1,001 bp (Table 1). Raw reads from the taxonomic classification analysis were then normalized as relative percentages for visualization, and taxonomic abundances with > 1% and those with < 1% abundances were grouped as “others.” The overall taxonomy classification resulted in the identification of 47 phyla, 331 families, 709 genera, and 1,560 species. In the organismal diversity, the most abundant domain was bacteria (97.78%) in kinema metagenome, and the remaining domains were viruses, eukaryote, and archaea (Figure 2A). Firmicutes (93.36%) was the most abundant bacterial phylum (Figure 2B and Supplementary Table 2). Further classification led to the identification of abundant families represented by Bacillaceae (81.03%), followed by Paenibacillaceae, Planococcaceae, and Enterococcaceae (Figure 2C and Supplementary Table 3). Bacillus was the most abundant genus followed by Brevibacillus, Kurthia, Nit1virus (virus), Proteus, Acinetobacter, and other genera (Figure 2D and Supplementary Table 4). At the species level, B. subtilis was detected as the most abundant species in kinema, followed by B. glycinifermentans, B. cereus, Brevibacillus borstelensis, B. licheniformis, B. thermoamylovorans, Kurthia sp. 11kri321, B. coagulans, B. circulans, B. paralicheniformis, and others (Figure 2E and Supplementary Table 5). No eukaryote and archaea species were found with > 1% relative abundance in kinema metagenomes.

The result showed that B. subtilis was predominant in kinema samples of India and Nepal (Figure 3A). Contrastingly, B. cereus and B. borstelensis were significantly (p < 0.05) abundant in kinema samples of Bhutan (Figure 3A). In the overall distribution of species of all domains, bacterial species were the dominant domain in which 280 species of Bacillus (Supplementary Table 6) and 121 species of lactic acid bacteria (LAB) (Supplementary Table 7) were detected with < 1% abundances in kinema metagenomes. Apart from bacteria, the second abundant domain was viruses with a total of 42 phage(s), out of which 26 were Bacillus phage and 16 were others (Supplementary Table 8). Myoviridae was the predominant viral family, among which Bacillus phage phiNIT1 was the abundant phage observed in kinema. Among eukaryotes, a total number of 56 species were identified; 17 species each were classified as molds and yeasts, respectively, and 22 species were grouped as other microbial eukaryotes (Supplementary Table 9). Additionally, 11 species were identified under the domain archaea (Supplementary Table 10).

Good’s coverage ranging from 0.98 to 0.99 was observed in all samples. The alpha diversity indices revealed the diversity in the following order: Bhutan (0.92) > India (0.71) > Nepal (0.57) (data not shown). Similarly, a non-parametric Shannon index also followed the same pattern, Bhutan (5.45) > India (3.92) > Nepal (2.89). Intra-diversity showed no significant differences in pooled samples of kinema from India (p = 0.54), Nepal (p = 0.55), and Bhutan (p = 0.48). Based on the species distribution and its relative abundance of all domains, the inter-diversity, calculated by PCoA, was distinctly different among samples (Figure 3B). The hierarchical clustering UPGMA revealed two clusters: cluster I (kinema from India and Nepal) and cluster II (kinema from Bhutan). Cluster I was observed with two subclusters, cluster IA and cluster IB of kinema from India and Nepal, respectively (Figure 3C). Statistical variation in the distribution of species among samples showed that kinema from India was significantly different from Nepal (p = 0.01945) and Bhutan (p = 0.000127), respectively. Similarly, kinema from Nepal was also significantly (p = 0.00007398) different from Bhutan.

It was observed that 485 overall species including bacteria, viruses, eukaryotes, and archaea were common in all samples. We observed that 479 bacterial species were shared in all samples, whereas 56 species were unique to India, 95 species to Nepal, and 229 species to Bhutan, respectively (Figure 4 and Supplementary Table 11). Among viruses, four species were found common such as Bacillus phage Shbh1, Bacillus phage SP-10, Bacillus phage SPG24, and Bacillus phage SPP1, whereas 13, 5, and 10 viral species were unique to kinema from India, Nepal, and Bhutan, respectively (Figure 4 and Supplementary Table 12). Among eukaryotes, two species (Pichia kudriavzevii and Mucor ambiguousi) were the common species, whereas 27, 8, and 10 species were unique species to India, Nepal, and Bhutan, respectively (Figure 4 and Supplementary Table 13). No common and shared archaea species were detected in kinema samples, whereas two, six, and three species of archaea were found unique to India, Nepal, and Bhutan, respectively (Figure 4 and Supplementary Table 14).

Different enhanced functional pathways were observed after mapping metagenomic ORFs against eggNOG and KEGG databases. Around 56.41–56.84% of the total mapped ORFs were assigned to COG genes, and 43.16–43.59% were assigned as KEGG pathway genes. Metagenome ORFs were aligned with genes associated to 216 KEGG metabolic pathways and categorized into three different levels. In level 1, metabolism was the most abundant category followed by environmental information processing, genetic information processing, cellular processes, human diseases, and organismal systems (Figure 5A). At level 2, there were 17 super-pathways with > 1% abundance and minor super-pathways with a relative abundance < 1% (Figure 5B), and at level 3, there were 36 super-pathways with > 1% abundance (Figure 5C) with 180 minor super-pathways with a relative abundance of < 1% (Supplementary Table 15). The distribution of the functional features in kinema metagenome was compared, and significant differences among the functional features of kinema of India, Nepal (p = 0.04182), and Bhutan (p = 0.02698) were found. Similarly, significant differences between kinema from Nepal and Bhutan (0.00007833) were observed. Exploring the enzyme classification by predictive analysis revealed the presence of protease, glucosidase, galactosidase, amylase, lipase, and γ-PGA proteins and enzymes involved in its biosynthesis (Supplementary Table 16). Additionally, further analysis was carried out on enzymes involved in aminoacyl-tRNA biosynthesis (Figure 6 and Supplementary Table 17) and phosphotransferase system (PTS) (Figure 7 and Supplementary Table 18), which revealed the enzymes responsible in the biosynthesis of amino acids and sugar components present in kinema.

The analysis of metabolites associated with predictive functional pathways using MetaboAnalyst v5.0 network revealed that many compounds were found associated with different metabolic pathways (Figure 8 and Supplementary Table 19). Amino acids such as glutamine, asparagine, valine, alanine/sarcosine, arginine, lysine, phenylalanine, tyrosine, glutamate, leucine/isoleucine, tryptophan, ornithine, and 4-aminobutanoate (γ-aminobutyric acid, GABA) were associated with predictive pathways related to amino acid metabolism. Fatty acids were associated with fatty acid biosynthesis, biosynthesis of unsaturated fatty acids, and alpha-linolenic acid metabolism (Figure 8). Sugars were mostly found to be associated with carbohydrate metabolism. Whereas, in citrate cycle (TCA cycle) pathways of the carbohydrate metabolism, organic acids such as citrate and isocitrate were involved, which were also found to be associated with glyoxylate and dicarboxylate metabolism. Among metabolites, daidzein, genistein, isovitexin, genistin, maackiain, biochanin-A, (+)-pisatin, apigenin, chrysin, and chrysophanol were found to be associated with flavonoid and isoflavonoid biosynthesis. Biotin, nicotinamide, pyridoxamine, pyridoxine, pantothenate, and riboflavin were found to be associated with metabolism of vitamins and cofactors. (R)-lipoate was associated with lipoic acid metabolism. Additionally, sarpagine and catharanthine were associated with indole alkaloid biosynthesis; swainsonine with tropane, piperidine, and pyridine alkaloid biosynthesis; and solavetivone with sesquiterpenoid and triterpenoid biosynthesis (Supplementary Table 19).

A total of 361 metabolites were identified using the PlantCyc database. PCoA scores were plotted with a total cumulative variance of 77.1% (PCoA-1, 54.7%; PCoA-2, 22.5%) with distinct formation of metabolites (Figure 9A). The hierarchical clustering analysis of significant metabolites revealed the formation of two clusters (Figure 9B), which followed the same pattern as shown in the microbial distribution. Cluster I constituted kinema from India and Nepal and was further sub-divided into two sub-clusters, clusters IA and IB, respectively. Cluster II constituted kinema from Bhutan and was observed to be distinctly different from the other samples. In metabolomic analysis, 36 primary metabolites (Figure 10A and Supplementary Table 20A) and 37 secondary metabolites were identified (Figure 10B and Supplementary Table 20B). The presence of isoflavones and various bioactive compounds associated with immunomodulatory effects, immuno-stimulatory effects, phytoalexins, and anti-viral properties was detected (Table 2). The abundance of secondary metabolites and other bioactive compounds was compared among kinema samples based on log-transformed fold change (Figure 11). Among isoflavones identified in kinema, daidzein, chrysin, and chrysophanol were significantly (p = 0.014651) higher in Nepal than in India and Bhutan. Other than isoflavones, compounds such as p-coumaroyltyramine, melatonin, 3-hydroxy-L-kynurenine, and catharanthine were higher in India compared to Bhutan and Nepal (p < 0.05). Serotonin, solavetivone, swainsonine, N-acetylaspartic acid, and sarpagine were higher (p < 0.05) in Bhutan than in India and Nepal, and soyasaponin III was found to be higher in Nepal than Bhutan and India (p = 0.0000396). The metabolites with at least fourfold differences among samples from different locations were detected (Figure 11 and Supplementary Table 21).

The correlation between the abundant species and primary metabolites (amino acids, fatty acids, sugars, organic acids, vitamins, and nucleobases) was constructed with a bivariate result obtained from Spearman’s correlation (Figure 12). B. subtilis, B. paralicheniformis, and B. licheniformis were positively correlated (p < 0.01) with arginine, phenylalanine, glutamate, leucine, isoleucine, tryptophan, and 4-aminobutanoate (GABA) among amino acids; linoleate and linolenate among fatty acids; galactitol, sorbitol, mannitol, ribose, xylulose, arabinose, xylose, fructose, mannose, galactose, glucose, L-gulose, scyllo-inositol, and myo-inositol among sugars; and citrate, isocitrate, and salicylate among organic acids. B. borstelensis, B. cereus, and B. coagulans were positively correlated (p < 0.01) with alanine, valine, lysine, palmitoleate, benzoate, and pyridoxine. Furthermore, the correlation analysis observed the formation of galactitol, sorbitol, mannitol, ribose, xylulose, arabinose, xylose, leucine/isoleucine, and salicylate with B. amyloliquefaciens, B. thermoamylovorans, and Kurthia sp. 11kri321 (p < 0.01). B. circulans was found to be significant and positively correlated (p < 0.01) with ornithine, glutamine, cinnamate, terephthalate, riboflavin, thymine, GABA, and glutamate. B. glycinifermentans was the only species found to be positively correlated (p < 0.01) with the compounds glutamate, phenylalanine, arginine, tryptophan, linoleate, linolinate, fructose, mannose, galactose, glucose, scyllo-inositol, myo-inositol, citrate/isocitrate, biotin, nicotinamide, adenosine, and adenine.

Similarly, the relationship of abundant species identified in kinema with secondary metabolites and other bioactive compounds was explored using a bivariate analysis results obtained from Spearman’s correlation (Figure 13). The correlation of B. glycinifermentans with the formation of daidzein, chrysin, chrysophanol, genistein, apigenin, isovitexin, and genistin was positively significant (p < 0.01). B. borstelensis, B. cereus, and B. coagulans were positively correlated (p < 0.01) with the formation of jasmonate, benzimidazole, demethylphylloquinone, 4-methylpyrazole, and (+)-pisatin. Phenylacetaldehyde, swainsonine, serotonin, solavetivone, sarpagine, maackiain, biochanin-A, 3-hydroxyanthranilate, and N-acetyl-L-aspartic acid showed a positive correlation (p < 0.01) with B. circulans, whereas 19 (R)-hydroxy-prostaglandin e2, 3-hydroxy-L-kynurenine, catharanthine, harmalol, p-coumaroyltyramine, melatonin, and salicylate showed a positive correlation (p < 0.01) with B. thermoamylovorans, Kurthia sp. 11kri321, and B. amyloliquefaciens. Among other metabolites, methyl jasmonate was the only metabolite showing a significant correlation (p < 0.01) with B. subtilis, B. licheniformis, and B. paralicheniformis.

Complex microbial community of major or minor domains exists or co-exists during uncontrolled natural fermentation of plant or animal substrates into fermented foods (Capozzi et al., 2017; Tamang et al., 2020). Kinema is also a fermented product of uncontrolled natural fermentation of soybean with existence of diverse microbial communities. The amplicon sequence of OTUs of kinema samples were earlier analyzed targeting the 16S rRNA gene for bacteria and ITS gene for fungi and yeasts (Kharnaior and Tamang, 2021). However, the profile of the microbial domains that appeared during kinema fermentation is important to know the entire microbial community structure of naturally fermented kinema. Hence, the shotgun metagenomic analysis was performed to explore the microbial community structure of kinema samples collected from India, Nepal, and Bhutan. DNA from the samples of the three countries were pooled in equal quantity and were used for metagenomic analysis. Pooled DNA samples were not significantly (p > 0.1) different in terms of sequence variation within the sample, since the pooled sequence could represent an actual sequence from each initial DNA sample (Anand et al., 2016).

In the organismal diversity, classified bacteria were the most abundant domain in kinema metagenomes, which was also observed in other naturally fermented soybean foods such as pe poke of Myanmar (Tamang et al., 2021b) and cheonggukjang of Korea (Tamang et al., 2022). In terms of relative abundance and occurrence, a colossal diversity of species of Bacillus was detected in kinema samples with 280 species including those detected at < 1% abundances. Previously, by amplicon sequence analysis, only 27 species of Bacillus were reported in kinema samples (Kharnaior and Tamang, 2021). Depth learning of genome mining of kinema has depicted huge microbial ecology, which is possible only by metagenomics using random shotgun sequence of DNA to highlight the higher resolution on taxonomic analyses at low abundance genera or species level (Tovo et al., 2020). Intra/inter-species diversity of Bacillus, both culturable and unculturable, in kinema exists and co-exists for the survival and competition among different species. The dominance of Bacillus in fermented soybean foods is probably due to the alkalinity of the product during fermentation (Gopikrishna et al., 2021) and also its proteolytic activity (Nguyen and Nguyen, 2020). Predominant Bacillus with different species has been reported in other Asian fermented soybean foods such as cheonggukjang (Nam et al., 2012; Tamang et al., 2022), natto (Kubo et al., 2011), pe poke (Tamang et al., 2021b), douchi (Hu et al., 2019), and thua nao (Pakwan et al., 2020). B. subtilis does not cause any health risk to consumers; hence, it is a safe and beneficial bacterium (Su et al., 2020). Besides B. subtilis, the next abundant species of Bacillus detected in kinema metagenome was B. glycinifermentans, which is known to produce many bioactive compounds (Stadermann et al., 2017). Some other major beneficial species of Bacillus, namely, B. thermoamylovorans, a thermostable enzyme producer (Yamada et al., 2017), B. coagulans, a probiotic bacterium (Cao et al., 2020), and B. licheniformis, functional bacterium (Muras et al., 2021), were also detected in kinema samples. However, the abundant species of Bacillus in samples of Bhutan was B. cereus, which was also detected by amplicon sequence analysis of kinema samples (Kharnaior and Tamang, 2021). B. cereus might have existed in fermenting soybeans during the uncontrolled spontaneous fermentation. Though B. cereus is a food-borne pathogen, its pathogenicity varies from harmless to lethal (Kavanaugh et al., 2022). B. borstelensis was also an abundant rod-shaped and Gram-positive bacterium in kinema samples, which is mostly thermoresistant, enzyme producing, and alkaline tolerant (Cotta et al., 2021). About 121 species of LAB (<1% relative abundance) were also detected from kinema samples, showing their co-existence with Bacillus spp. Although their abundance is low (<1%), their existence in naturally fermented soybean food may have fermentative or other functional roles, which may be investigated further.

Abundance of viruses, mostly bacteriophages, was only 1.85% in kinema metagenomes; their appearance might suppress the growth of some pathogenic bacterial species, which were detected in low abundance, probably combating their multiplication for food safety to extend the shelf life of the product (Endersen and Coffey, 2020). Few archaeal species were also detected in kinema; though their abundance was very low, their contributions to health and disease remain unknown (Adam et al., 2017).

Microbial communities in kinema samples of Bhutan showed the highest diversity index compared to those of India and Nepal, probably due to environmental factors (Lee S. et al., 2017), seasonal changes (Kumar et al., 2019), geographical locations (Lee et al., 2019), and inconsistency in traditional preparation practices adopted by the local people. We observed that a set of species was common and unique to each country, and the reason could be the unusual associations of microbial communities and environmental conditions including temperature and humidity. However, the predominance of beneficial microbes with various antimicrobial properties could be the reason for the reduction of the abundance of unwanted microorganisms.

The predictive functional analysis of kinema metagenomes by KEGG annotation revealed the abundance of several metabolic pathways, where carbohydrate metabolism was predominant. The abundance of genes related to galactose metabolism could involve in the degradation of complex sugars in soybeans during fermentation (Chun et al., 2021). In amino acid metabolism, the abundance of genes associated to sub-pathways alanine, aspartate, and glutamate metabolism in kinema metagenome could enhance the aroma and flavor of the product (Yi and Hong, 2021). Glutamate is known as the precursor for the biosynthesis of GABA (Dhakal et al., 2012), which has several bio-functional properties (Diez-Gutiérrez et al., 2020). The presence of genes related to branched-chain amino acids was identified, which could be involved in the metabolism of carbon or nitrogen by B. subtilis (Belitsky, 2015). β-glucosidase detected in kinema could be involved in the biosynthesis of isoflavones and hydrolysis of oligosaccharides (Khosravi and Razavi, 2021).

Metabolite profiled by LC-MS analysis revealed the presence of several bioactive compounds indicating kinema has several bio-functional properties, which may impart health-promoting benefits to the consumers. Although raw soybeans have various functional compounds, the presence of microorganisms can enhance the levels of production of amino acids and secondary metabolites (Hyeon et al., 2020). The abundance of B. subtilis during the natural fermentation processes of kinema could enhance the nutrient content because of the ability to release hydrolase enzymes that breaks down macromolecules into small molecules and increases the content of metabolites, as also observed in cheonggukjang (Lee et al., 2017; Ali et al., 2018). Some metabolites such as amino acid, namely, alanine, leucine, isoleucine, phenylalanine, and valine, and phenylacetaldehyde were identified in metabolomic analysis of kinema, which are known to contribute to taste and flavor development in the product (Lioe et al., 2018). Different types of isoflavones such as daidzein and genistein were detected in kinema, which are mostly present in soybeans (Choi et al., 2020a,b) and are associated with many health-promoting or immune-stimulating properties (Kim, 2021). The content of isoflavones in kinema could be indirectly influenced by the presence of microbes in the rhizosphere of soybean plants that secretes various secondary metabolites, since soybean is known to produce specialized metabolites such as isoflavones and soyasaponins (Sugiyama, 2019). Group B soyasaponin (soyasaponin III) was also detected in kinema samples, which has an antitumor effect, anti-inflammatory activity, and antimicrobial and cardiovascular-protective activities (Omar et al., 2020). Metabolomics result revealed the presence of vitamin B complex such as biotin, riboflavin, nicotinamide, pyridoxamine, pyridoxine, and (R)-pantothenate in kinema, which are the sources of micronutrients required for bio-functional activities including host immunity (Yoshii et al., 2019). Vitamin K, also known as demethylphylloquinone, was also detected in kinema, which has an anti-inflammatory effect (Harshman and Shea, 2016). The presence of metabolites such as lipoate, 3-hydroxyanthranilate, and melatonin in kinema may be related to antioxidant activities (Galano et al., 2016), and detection of chrysophanol in kinema may have immuno-stimulating effects, anti-inflammatory effects, and anticancer activities (Wen et al., 2018; Prateeksha et al., 2019). Also, p-coumaroyltyramine and 4-methylpyrazole were detected in kinema, and these bio-active compounds are used as anti-diabetic and anti-acetylcholinesterase (Masi et al., 2021). Benzimidazole, as an anti-inflammatory agent (Hashem and Bakri, 2021), is also detected in kinema samples. Among the samples, kinema samples of Bhutan depicted a higher diversity and secondary metabolites comparable to those of India and Bhutan. Therefore, we believe that microorganisms play a crucial role in secreting and changing the chemical diversity of various metabolites under certain environmental conditions including temperature and pH. We also observed that the metabolite abundance differs among kinema from different locations, which suggested that the differences in cultivars or soybean seeds may vary the composition of metabolites (Wang et al., 2019), and also, the formation of metabolites could enhance or affect the accumulation of secondary metabolites by various environmental factors such as regional variation, climates, and methods of preparation (Iqbal et al., 2018).

Predictive functionality generated by KEGG annotation is a hypothetical claim, and prediction of genes is associated with functional profiles. Validation with real-time experiment may verify the predictive metabolic pathways to claim the presence or absence of some metabolites in the samples. Hence, we compared the KEGG database generated by metabolomics analysis with that of predictive metabolic pathways of the abundant metagenomes in kinema samples of India, Nepal, and Bhutan. Stickiness in kinema is due to γ-polyglutamic acid (γ-PGA) (Chettri et al., 2016), which determines the quality of the product preferred by the consumers. B. subtilis and B. licheniformis were both correlated with glutamate production in kinema, probably associated with γ-PGA production (Cai et al., 2017; Zhang C. et al., 2018). The positive correlation between B. amyloliquefaciens and many amino acids in kinema may enhance the antioxidant activity and formation of bioactive compounds and increase the levels of isoflavones (Shahzad et al., 2020). B. coagulans showed a positive correlation with tryptophan, tyrosine, and genistin in this study, which may improve the amino acid absorption (Stecker et al., 2020). The correlation between amino acids and B. borstelensis observed in kinema metagenomes might be the ability of thermostable D-amino acid amidase to produce amino acids at higher temperature (Baek, 2003). Unsaturated fatty acids (UFAs) such as oleate, linoleate, and linolenate were detected in kinema, the biosynthetic pathway of which were also predicted in kinema metagenomes. UFAs are important untargeted metabolites in fermented soybean foods (Park and Kim, 2021), which play an important role in biotic and abiotic stresses (He et al., 2020). The co-existence of B. subtilis with some fungi in kinema might be correlated with increase of unsaturated fatty acids in fermented soybeans (Kanghae et al., 2017). Pantothenate and pyridoxamine were positively correlated with B. subtilis strain, which may enhance the biosynthesis of vitamins and cofactors in kinema, imparting health-promoting benefits to consumers. Serotonin, a neurotransmitter (Sahu et al., 2018), was found to be correlated with the neuroactive ligand–receptor interaction via predictive function, and it could control the signaling pathways either intracellularly by serotonylation or extracellularly via membrane receptors (Bockaert et al., 2021). Swainsonine predicted in KEGG database was also detected in metabolomics analysis in kinema samples, which is involved in the biosynthesis of tropane, piperidine, and pyridine alkaloid (Li and Lu, 2019). A positive correlation between the predominant Bacillus species and the secondary metabolites such as aglycones and glycosides was observed, with the potential ability of Bacillus species in producing secondary metabolites that mediate antibiosis and exhibit as antimicrobial compounds (Caulier et al., 2019). A variation in microbial and metabolite formation in naturally fermented kinema samples of three different countries has been observed, and this may be due to different cultivars of soybean seeds used for fermentation, soil types, geographical and ago-climatic variation, and community- or region-specific conventional methods of preparation. However, the results obtained from the metagenome–metabolite interaction by correlation analysis may provide a better understanding in the selection of strains or targeting metabolites of interest for further studies.