Decoding the microbiome and metabolome of the Panchagavya —An indigenous fermented bio‐formulation

For the first time, updated molecular techniques were used to validate and elucidate the effect of the Panchagavya. Metagenomics was used to decipher the bacterial microbiome structure, which showed promising results for their existence and abundance in the Panchagavya. The word “Panchagavya” is derived from Sanskrit, meaning “five ingredients from a cow,” which includes dung, milk, urine, curd, and clarified butter (ghee) [1]. The “Panchagavya” application has been reported to enhance plant growth and nutrient uptake, reducing the biotic and abiotic stresses in crop plants. Preparation of the Panchagavya formulation involves several fermentation processes with a consortium of beneficial bacteria [2]. In such bio-formulations, cow dung is an important constituent that acts as a source of bacterial inoculum [3]. The Panchagavya formulations can be applied to crop plants using many methods, among which foliar spray is widely used. The bio-formulation can be suspended or dissolved in water and, therefore, can be used along with irrigation water. In organic farming, it is used as a foliar spray at a concentration of 3% to supplement plant nutrient requirements and suppress pests and diseases [4]. The Panchagavya is also used as a biofertilizer at 1:100 (Panchagavya: soil; v/wt) to promote plant growth [5]. It was reported previously that the Panchagavya formulation contains numerous plant growth-regulatory molecules, such as gibberellins (GAs) and indole acetic acid (IAA) [3, 6]. In addition, the beneficial microorganism, namely, the Plant Growth-Promoting Rhizobacteria (PGPR), synthesizes many amino acids and other essential enzymes. It is evident from the previous reports that the foliar application of the Panchagavya enhances the yield and quality of cereals and vegetables [7]. Several studies on the biodynamics of bacterial species in the Panchagavya exists that are linked to the solubilization of various carbon compounds and minerals using a polyphasic approach resulting in the most reported bacterial species to produce IAA and ammonia [3]. Despite its wide use in crop production, especially in organic agriculture, the metabolic profiles (hormones, proteins, nutrients, etc.) and microbial profiles have not been explored completely. Moreover, the profile and dynamics of the entire microbial community associated with the Panchagavya are unavailable due to limitations in culturing several genera/species using conventional microbiological techniques. Also, the complete enzymatic profile of the Panchagavya is still unknown. Not many efforts have been made to explore the complete biology involved in the plant-growth-promoting and ISR activity of the Panchagavya. It is possible to understand the fundamental aspect of a complex system only if different omics approaches are combined [8]. Therefore, we have used different omics approaches to decode the microbial community (metagenomics) and metabolite profile (metabolomics) in the Panchagavya. These approaches can reveal the microbial composition and their abundance, the functional annotation of genes and important protein compounds, hormones, and so forth, in the Panchagavya formulation. With the recent advancement in metagenomic sequencing, construction of the whole microbiome is easier and more comprehensive than the conventional 16S rRNA sequencing of culturable bacteria [9, 10]. However, metagenomics reveals only the microbial populations; therefore, metabolomics is required to study the in situ microbial metabolic activity and its capabilities [11]. Therefore, with the help of metabolomics, we can identify the plant growth-promoting compounds, including hormones, which can be further screened with high-throughput metabolite identification technologies such as liquid chromatography mass spectrometry (LCMS) [12]. In this study, through an integrated approach, we report the microbial population assemblages of the Panchagavya by metagenomics and identify its potential plant-growth-promoting molecules through LCMS attached to the “metlin” library. Furthermore, by combining different omics technologies, we have revealed the overall function of these microbes in the Panchagavya. This is the first comprehensive report of the microbial and metabolome profiles of the Panchagavya. Metagenomic DNA from the Panchagavya formulation was used to decode its microbiome structure. The total sequence data size obtained was 135.6 Mb with 160,945 assembled contigs (Supporting Information: Table S1A). The operational taxonomic unit (OTU) picking through Eukaryotic and prokaryotic diversity was estimated using CCMetagen v1.2.5, which revealed that the domain bacteria were the most abundant OTUs (99.58%), with a high abundance of species belonging to the phylum Proteobacteria (40.40%), Firmicutes (26.90%), and Bacteroidetes (14.20%). Rhodospirillales (24%), Lactobacillales (22.40%), Bacteroidales (11.30%), Enterobacteriales (6.60%), Pseudomonadales (2.40%), Eubacterials, and Burkholderiales (1.00%) were the most prevalent orders. Among the class, Alphaproteobacteria was highly abundant (24.30%), followed by Bacilli (23.00%), Bacteroidia (11.30%), Gammaproteobacteria (11.20%), and Betaproteobacteria (3.40%). The proportions of bacterial families such as Acetobacteraceae, Streptococcaceae Prevotellaceae, Enterobacteriaceae, Bacteroidaceae, and Pseudomonadaceae were 23.70%, 21.30%, 6.8%, 6.60%, 3.7%, and 2.40%, respectively, in the Panchagavya bio-formulation. The most predominant genera in the decreasing order of their abundance were Streptococcus (20.80%), Acetobacter (14.60%), Prevotella (5.60%), Bacteroides (3.50%), Gluconobacter (3.00%), Klebsiella (2.60), Enterobacter (2.30%), Azospira (2.10%), Azotobacter (1.70%), Gluconoacetobacter (1.60%), Acinetobacter (1.40%), Selenomonas (0.90%), and Pseudomonas (0.40%). The Panchagavya bio-formulation was rich in S. lutetiensis (9.90%), A. peroxydans (8.20%), S. equinus (5.00%), S. infantarius (3.30%), G. oxydans (3.00%), B. fragilis (2.70%), uncultured B. bacterium (2.50%), and A. chroococcum (1.20%) species. The relative distributions of bacteria belonging to the different phyla, classes, order, families, and genera are summarized in Figure 1A and Supporting Information: Table S2. Based on previous reports, the different genera/species of bacteria identified in the Panchagavya formulation were grouped into different functional groups. Several bacteria that promote plant growth were identified in the Panchagavya bio-formulation. Among them, Accumulibacter phosphatis, Achromobacter spp. Root565, Acidobacteria bacterium, Acinetobacter lactucae, A. kooki, Aquitalea magnusonii, Bacillus cereus, Candidatus Accumulibacter spp. 66-26, Citrobacter amalonaticus, Clostridium pasteurianum BC1, Enterobacter cloacae, E. asburiae, E. cancerogenus, E. kobei, E. xiangfangensis, Gluconacetobacter diazotrophicus, Rhizobium leguminosarum, Pseudomanas putida, Acinetobacter kooki, Azospirillium brasilense, A. lipoferum, P. fluorescens, Azotobacter beijerinckii, B. subtilis Lactobacillus mucosae, and L. harbinensis were found to be predominant (Figure 1B,C and Supporting Information: Table S2). Several bacteria reported as an ISR activators in plants against biotic stress were found in the Panchagavya bio-formulation. The unclassified Burkholderiales showed the highest reads (79 reads), followed by Pseudomonas fluorescens (12 reads), Bacillus cereus (12 reads), and Pseudomonas chlororaphis (3 reads), which are involved in imparting ISR in the crop plants (Figure 1 and Supporting Information: Table S2). Various bacterial genera such as Butyrivibrio spp. Caloramator spp. Cellulosilyticum lentocellum. Chloroflexi bacterium, Clostridium clariflavum, C. cellulosi, C. disporicum, C. homopropionicum, C. innocuum, C. kluyveri, C. leptum, C. magnum, C. neopropionicum, Herbinix hemicellulosilytica, Leclercia sp. Marvinbryantia formatexigens, Dyella jiangningensis, Pseudobacteroides cellulosolvens, Ruminococcus spp., Selenomonas spp., and so forth, which decompose cellulose, hemicellulose, and other fibers, were found in the Panchagavya formulation (Figure 1 and Supporting Information: Table S2). Several bacteria (Acidovorax ebreus, Acinetobacter tandoii, A. venetianus, Alicycliphilus spp., Alkaliphilus metalliredigens, Chloroflexi sp., Dechloromonas spp., Desulfovibrio spp., Lysinibacillus spp., Paenirhodobacter enshiensis, Pandoraea thiooxydans, Propionispora spp., Pseudomonas aeruginosa, P. mendocina, P. alcaligenes, P. balearica, P. knackmussii, etc.) reported to be associated with soil bioremediation by degrading hydrocarbons, polyesters, phenols, and so forth, were predominantly present in the Panchagavya formulation (Figure 1 and Supporting Information: Table S2). Nutrient recyclers such as Accumulibacter phosphatis, Acidobacteria bacterium, Acidovorax caeni, Candidatus Accumulibacter sp. Citrobacter amalonaticus, Dechlorosoma suillum, Denitrobacterium detoxificans, Desulfosporosinus orientis, Desulfovibrio spp. Gallionella sp. Hylemonella gracilis, Mitsuokella jalaludinii, Novosphingobium nitrogenifigens, Pandoraea thiooxydans, Pseudoxanthomonas spp., Sebaldella termitidis, Thiobacillus sp., and so forth, were also part of the Panchagavya microbiome (Figure 1 and Supporting Information: Table S2). As auxin transporters, Azotobacter chroococcum, Arcobacter butzleri, uncultured Ruminococcus sp., Sutterella sp., Bifidobacterium bohemicum, Prevotella sp. CAG:891 and Pedobacter sp. V48 play a pivotal role. Paenirhodobacter enshiensis, Dialister sp. CAG:357, Acidovorax sp. 202149, Comamonas aquatica, Lactobacillus mucosae, Streptococcus equinus, Enterobacter asburiae, and Lactobacillus harbinensis aid cytokinin transportation. In addition, several bacterial genera were identified. However, their function in the soil and plant system is unknown (Figure 1 and Supporting Information: Table S2). The genes involved in promoting plant growth, metabolism, nutrition recycling, and so forth, were predicted. The number of genes predicted by Prodigal was 122,976. Of these, 87.9% (n = 108,096) of the genes were annotated (without filtering). When genes were filtered with 50% sequence identity and 30% query coverage, 76.6% (n = 94,359) of the genes were annotated. Similarly, genes with 70% identity and 50% query showed 49.7% (n = 61,175) annotated genes. The predicted genes matched Uniprot protein sequences of bacteria, fungi, archaea, and viruses. The sequence annotations, gene names, and organisms with other details of alignments are listed in Supporting Information: Table S3. The data also revealed the presence of important regulatory proteins involved in the biosynthesis and transport of auxin (in Arcobacter butzleri and six other bacteria) and cytokinin (in Paenirhodobacter enshiensis and also in other 35 genera) (Supporting Information: Table S4). The gene ontology classification revealed that all predicted gene sets were involved in molecular functions (79522 enriched genes onto logs), biological processes (70463 enriched genes onto logs), and cellular components (43,369 enriched genes onto logs). The GO analysis revealed that the percentages of genes involved in cellular components, viz., the cell part, the cell, and the membrane were 40.20%, 40.20%, and 17.90%, respectively. In the molecular function component, gene ontology responsible for catalytic activity (67.30%), binding (52.30%), and transporter activity (6.00%) was abundant. Higher gene percentages of GO terms in the biological process were the cellular process (62.30%), the metabolic process (60.20%), localization (10.70%), and response to stimuli (9.70%) (Figure 1D and Supporting Information: Table S5). The Clusters of Orthologous Groups of proteins (COGs) analysis distributed the annotated proteins into three major functional sections, that is, (a) metabolism, (b) cellular process and signaling, and (c) information storage and processing. Further subsystem analysis predicted the proteins related to carbohydrate metabolism (12.44%), amino acid metabolism (8.90%), protein metabolism (7.40%), RNA metabolism (6.09%), DNA metabolism (3.89%), membrane transport (2.52%), and stress response (2.44%) (Figure 1D). The associated annotated genes of different pathways were predicted using the DAVID program. Some of the significant pathways obtained were Amino Acid Biosynthesis (AAB), Amino Acid Degradation (AAD), Amine and Polyamine Biosynthesis (APAB), amino–sugar metabolism (ASM), Bacterial Outer Membrane Biogenesis (BOMB), Cell Wall Biogenesis (CWB), Co-Factor Biogenesis (CFB), Glycan Degradation (GD), Glycan Metabolism (GM) and Metabolic Intermediate Biosynthesis (MIB), Nitrogen metabolism (NM), one-carbon metabolism (OCM), Protein modification (PM), Pyrimidine metabolism (PYM), Sulfur metabolism (SM) pathways, and so forth (Figure 1 and Supporting Information: Table S6). The microbial community composition of the Panchagavya formulation assessed using CCMetagen revealed the reads of eukaryotic (17.50%) and prokaryotic (57.10%) origins (Supporting Information: Figure S1). The assembled metagenome by MEGAHIT was found to have 189,455 contigs, out of which 188,989 contigs were revealed to be of prokaryotic origin by EukRep, which was retained for downstream analysis. Read coverage and the relative abundance-based binning of the prokaryotic contigs using MetaBAT2 resulted in 54 genome bins. De-replication of the genome bins with dRep with minimum completeness of >75% and maximum contamination ≤ 25% resulted in 19 genome bins. CheckM analysis revealed that of the 19 bins, 3 bins (bin 4, bin 31, and bin 54) were resolved to the genus level, 3 bins (bin 2, bin 7, and bin 40) were resolved to the family level, and 4 bins (bin 1, bin 24, bin 32, and bin 43) were resolved up to order level (Supporting Information: Table S1). Taxonomy was assigned to these 19 bins using GTDB-Tk against the Genome Taxonomy Database, as mentioned in Supporting Information: Table S1C. The de-replicated bins were annotated using PROKKA, and the annotated genomes in gbk format were used as input for AntiSMASH to predict secondary metabolites. A total of 29 clusters were predicted from 11 bins by identifying genes encoding putative biosynthetic enzymes (e.g., polyketide synthases, nonribosomal peptide synthases, terpene, molybdenum cofactor), and associated genes were identified and are tabulated in Supporting Information: Table S1D. The XCMS was used to decipher the metabolites of the Panchagavya formulation. A total of 5737 features were identified in XCMS metabolite analysis (2035 in the aqueous extract method and 1783 in the multisolvent extract method), of which 2970 were identified in the Metlin database (1086 and 985 features of aqueous and multisolvent methods, respectively) (Supporting Information: Table S7A). A plant database was used to determine the metabolites extracted from both aqueous and multisolvent extractions. The essential plant hormones, like gibberellic acid, and its eight other derivatives, such as Gibberellin A70, Gibberellin A45, Gibberellin A61, Gibberellin A51, Gibberellin A4, Gibberellin A20, Gibberellin A29-catabolite, and Gibberellin A14 aldehyde, were recorded. In addition, several compounds such as Kaurenol, Strigol, Brassinin, Brassicasterol,cis-Zeatin 9-glucoside IAA, and so forth, belonging to different chemical groups that are associated with plant-growth-promoting functions and biotic stress tolerance, were also recorded in the Panchagavya formulation (Supporting Information: Table S7B). Organic farming has been gaining more significance worldwide due to the alarming impact of synthetic-chemical-based agriculture on human health and the environment. In organic farming, locally made organic manure prepared from the on-farm organic waste of the farm is the main input [13, 14]. One such important producer of raw materials is the cow. It has been reported that several traditional bio-formulations like the Panchagavya made from cow-milk-derived products and cow excreta (dung and urine) have shown plant-growth-promoting functions and confer disease tolerance ability to several crops [4, 5]. The Panchagavya is a mixture of five cow products obtained from an indigenous cow that include dung, urine, milk, curd, and ghee. However, the microbial and metabolic profiles of the Panchagavya, which are responsible for its crop-benefitting property, have not been described completely. Moreover, the potentialities of the Panchagavya microbiome and their products in bio-organic metabolism in soil have still not been mapped. We prepared Panchagavya's indigenous bio-formulation and analyzed its microbiome and plant-growth-promoting metabolite profile in this work. The metagenome sequencing approach was used in this study to understand the diversity and function of microbial niches, which could not be addressed using conventional methods [15, 16]. We have used a Metagenomic sequencing approach for profiling all the associated microbiota of the Panchagavya. The analysis revealed that the Panchagavya harbors a wide range of bacterial genera with the domain bacteria as the most abundant OTUs of 99.58% with a plethora of species belonging to the phylum Proteobacteria, Firmicutes, and Bacteroidetes. Further, functional classification of identified bacterial genera revealed the many bacteria involved in PGPR, ISR, biocontrol of plant pathogens, nutrient recycling in soil, and plant hormone transporters. Among the PGPRs, Enterobacter cloacae, E. asburiae, Rhizobium leguminosarum, Pseudomanas putida, Acinetobacter kooki, Azospirillium brasilense, A. lipoferum, P. fluorescens, Azotobacter beijerinkii, B. subtilis, and Lactobacillus mucosae were found to be predominant, whereas L. mucosae was found in abundance in the bio-formulation. The Pseudomonas fluorescens and Bacillus cereus were involved in imparting ISR in the crop plants, showing 12 reads each. At the same time, the unclassified Burkholderiales showed the highest reads of 79. The exploitation of microflora from cow excreta for plant growth promotion, phosphate solubilization, nitrogen fixing, cellulase-producing bacteria, enzymatic activities, methanogenic bacteria, and ammonia-producing bacteria, and plant hormone transporters is well reported [17]. Inhabiting microorganisms in the Panchagavya play an important role in agriculture by solubilizing the available minerals for plant growth. This study thus provides microbiological evidence for enhanced plant growth after applying the Panchagavya to the crop plants, as reported previously [5]. Various other bacterial genera were found in the Panchagavya formulation, which decompose cellulose and hemicellulose. Several bacteria associated with soil bioremediation by degrading hydrocarbons, polyesters, phenols, and so forth, were predominantly present in the Panchagavya formulation. Bacteria associated with nutrient recycling in soil were predominant, and bacteria responsible for transporting the plant hormones like IAA and cytokinin were also found. Our study also identified several bacterial genera associated with ISR induction in crop plants. In addition, several bacterial genera that can degrade complex carbohydrates such as cellulose, hemicellulose, and lignin were identified, indicating the usefulness of the Panchagavya in recycling/maintaining organic carbon in the soil. Recently, soil pollution with several toxic hydrocarbons has been a common problem, especially in agricultural fields, due to excessive use of pesticides, machinery, and sewage water for crop production. The waste generated by the industries has become an environmental issue due to the lack of technologies available for effective waste management practices involving physico-chemical and biological treatment. Unutilized metals that are not essential for plant growth also accumulate in various parts of the plant body. Higher contents of these metals may tend to inhibit growth and can be toxic to plants, animals, and humans. We have identified several bacteria (Acidovorax ebreus, A. venetianus, Alicycliphilus denitrificans, Alkaliphilus metalliredigens, Chloroflexi bacterium, D. aromatica, D. suillum, P. enshiensis, etc.) in the Panchagavya formulation that can degrade the hydrocarbon pollutants of soil and bioremediation of heavy-metal contamination. Therefore, use of the Panchagavya represents a potential strategy for rejuvenating contaminated/polluted soils to produce healthy crops. Plant hormones are among the most important growth regulators. They have a major impact on plant metabolism and play a vital role in stimulating plant defense response mechanisms against stresses [18]. In the current study, the Acinetobacter spp. were found to be more abundant in the bio-formulation. Several species of this genera are known to solubilize the phosphate in the soil more effectively and suppress the soil-borne fungal pathogens [6]. In addition, Acinetobacter spp. are known to produce several growth-producing hormones like IAA, gibberellic acid, and their derivatives [19]. Our metabolome analysis revealed that many plant growth-promoting hormones were the major components, and possibly, the IAA would have been produced by Acinetobacter sp. Another important plant hormone identified in our bio-formulation through metabolomics was gibberellic acid and its eight derivatives. Gibberellic acid and its isoforms are essential hormones for plant growth. The bacteria genera such as Azotobacter, Azospirillum, Pseudomonas, Acetobacter, Burkholderia, and Bacillus were abundant in the Panchagavya. They are the species of bacteria reported to produce gibberellic acid and its isoforms [20]. Swain and Ray [21] (2009) reported that bacteria such as Bacillus, Corynebacterium, Lactobacillus, Leuconostoc, Bifidobacterium, Enterococcus, and Streptococcus are present in cow dung, a core ingredient for preparing the Panchagavya. A cytokinin riboside 5'-monophosphate phosphoribohydrolase (EC 3.2.2.n1) was produced by 36 organisms in the Panchagavya; this enzyme converts inactive cytokinin nucleotides into the biologically active free-base forms by phosphoribosyl hydrolase activity. It is directly involved in activating cytokinin, an important growth-promoting hormone [22]. Sreenivasa et al. [23] have reported that gibberellic acid (as high as 3.01 µg/25 ml) was produced by the bacteria present in the cow urine, a constituent of the Panchagavya. In our study, from our metagenomic data, we identified several species of Rhizobium (Bradyrhizobium japonicum, R. miluonense, R. leguminosarum), which are known to produce Gibberellic acid and its derivatives [24]. In addition to the hormones, we also identified several essential proteins and enzymes in the bio-formulation involving many metabolic pathways in the soil, contributing to the overall soil health and crop growth. The predominant genes involved in the molecular functions, cellular process, and biological process were the most abundant in gene ontology analysis. We observed that the catalytic activity, cellular and metabolic process, binding, and cell part-related genes relapsed in the Panchagavya bioformulation. These genes play a dominant role in increasing the rate of chemical reactions, transcription, DNA replication, and expression of genes in the cell, which play a significant role in the overall growth of the plants. Metagenomic analysis exploring the functional diversity suggests that microbial communities can degrade carbohydrates, amino acids [25], cell division, and the cell cycle [26]. By decomposing organic materials and litter, microbial populations appear to play an essential role in the carbon cycle by reducing carbohydrates [27]. 0.2% of functional genes were related to nitrogen utilization. This is consistent with previous studies, which included habitats with nitrogen-fixing bacteria, such as soybean crops [28]. This is a preliminary study unraveling the microbes and metabolites responsible for the beneficial effects of the Panchagavya formulation when applied to crop plants and to increase the abundance of beneficial microbes and metabolites, which can be recommended alternatives to chemical fertilizers to enhance crop productivity. In addition, our study has identified several bacterial genera in the Panchagavya which have not been characterized so far for any role in the soil and plant system, which necessitates a systematic study to characterize such potential bacteria. For the first time, updated molecular techniques were used to validate and elucidate the effect of the Panchagavya. Metagenomics was used to decipher the bacterial microbiome structure, which showed promising results for their existence and abundance in the Panchagavya (number of reads of particular taxa; the gene annotations and functional annotations were predicted). On the other hand, the metabolome predicted the byproducts of numerous plant growth promoters and other compounds released and fermented by these bacteria. According to these findings, the Panchagavya is a promising manure that contains useful microorganisms and plant development promoters. As a result, the Panchagavya is a bio-formulation that is both ecologically benign and effective in improving soil and crop health, Figure 2. The Panchagavya was prepared following the previously described protocol [27]. The ingredients for the bio-formulation were obtained from the native cow maintained at the Organic Research Institute, University of Agricultural Sciences, Bangalore, India. A wide-mouthed container was used to prepare the Panchagavya formulation. First, a mixture of cow dung and cow ghee (clarified butter) was allowed to ferment for 4 days. On the 5th day, the remaining ingredients, cow urine, cow milk, and cow curd, were added and incubated for another 8 days. The solution was frequently stirred to promote aerobic microbial activity for another 10 days. Finally, the fermented bio-formulation was used for analysis. A 2 ml aliquot (mixed and filtered) of the Panchagavya was used for DNA extraction using the DNASure® Plant Mini Kit (Genetix Biotech) following the manufacturer's instructions. The quality and quantity of DNA were assessed using the Qubit™ dsDNA BR Assay Kit (Thermo Fisher Scientific) and Nanodrop® (DS-11 FX+, DeNovix® Spectrophotometer/Fluorometer; DeNovix Inc.) for metagenome sequencing. The sequencing was performed at a commercial facility (Sandor Life Sciences Pvt. Ltd., Hyderabad, India) using Illumina NextSeq. 500. One microgram of DNA concentration was used to prepare the whole-genome library (Library Prep Kit for Illumina®; New England BioLabs® Inc.). Using Covaris, 151 bp size was fragmented, end-repaired, and immediately ligated with the adaptor; then, the library was enriched with 10 cycles of PCR. Using Agencourt®AMPure® XP beads (Beckman Coulter), purification was performed following the manufacturer's instruction, and their quality was checked using the Agilent 4200 TapeStation System (Agilent Technologies Inc.) before sequencing in Illumina NextSeq. 500. Sequences of paired-end reads were obtained as FASTQ, and their quality was checked using the FastQC program (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The quality of the raw reads was checked for the ambiguous bases and Phred score of Q ≤ 30 (base call accuracy 99.9%). The trimmed reads were then assembled using MetaSPAdes Assembler, a variant of SPAdes assembler designed explicitly for whole metagenome studies [28]. The contigs of length <150 bp were discarded from the analysis to reduce false-positive gene prediction and annotation. The filtered sequences with an accuracy of 99.9% were used for choosing the operational taxonomic unit (OTU) in “MEtaGenome Analyzer” (MEGAN 4). The metagenome assembly was used for gene prediction using Meta GeneMark [29]. The read depth of each predicted open reading frame (ORFs) was quantified to show the abundance of a particular taxonomic group. Identification of eukaryotic and prokaryotic lineages in processed metagenomic reads was performed using CCMetagen [v.1.2.5] [30]. The MEGAHIT assembler [31] implemented in the SqueezeMeta pipeline [32] was used to assemble the processed metagenomic reads. Separation of prokaryotic from eukaryotic contigs was performed using EukRep [33]. Read coverage and relative abundance were estimated using CoverM (https://github.com/wwood/CoverM), which was provided as input to MetaBAT2 [34] for binning the prokaryotic contigs to genome bins. dRep was used to “de-replicate” the genome bins by identifying groups of highly similar genomes and choosing the best representative genome for each genome set. CheckM [35] was used to estimate the genome completeness and the rate of contamination of the de-replicated genome bins. Taxonomy was assigned to the genome bins using GTDBTk [36]. The taxonomy assigned bins were further annotated using Prokka [v.1.14.6] [37]. antiSMASH v6 [38] was used to predict the gene clusters that have the potential to produce secondary metabolites (detection strictness parameter kept “relaxed”). The functional annotation was performed by aligning ORF sequences against the NCBI-nr protein database using the DIAMOND tool version 0.7.9.58 [39]. The coding genes for predicted proteins were identified using Prodigital (version 2.6.1) [40]. The annotated gene result was then used for taxonomy assignments. Further, as described previously, bacteria involved in PGPR, phytohormone production, and activating the ISR in plants were identified [5]. A circular dendrogram was plotted using the iTOL online tool (https://itol.embl.de/), depicting the various species. Finally, the InterPro and SEED database (http://pubseed.theseed.org/) was used to annotate proteins into a set of functional clusters. The protein sequences of the metagenome-assembled genome were subjected to diamond blastp against the UniProt database, and the output was annotated using BLAST2GO [41]. The BLAST2GO annotation was exported in WEGO annotation format, and the proteins without GO annotation were removed. The filtered annotated file was subjected to WEGO for GO Analysis [42]. The gene IDs were obtained from the BLAST results for PGPRs, and pathway analysis was done for the extracted genes. Several pathways were identified using DAVID (https://david.ncifcrf.gov/), where the count threshold value is set to default two and EASE as 0.1. The extraction of metabolites from the Panchagavya was carried out using multisolvent and aqueous methods. In the multisolvent method, methanol, acetonitrile, and water were used in the ratios of 2:2:1 as the extraction solvent. About 10 ml of the Panchagavya sample was incubated with 50 ml of extraction solvent for 1 min and sonicated (Misonix Sonicators) for 10 min. After sonication, the sample was centrifuged at 12,000g for 15 min at 4°C. The resultant supernatant was collected and stored at −20°C until LC-MS/MS analysis. However, in the aqueous method, ddH2O water was used as the solvent. About 10 ml of the Panchagavya sample was incubated for 1 min with 50 ml of extraction solvent and then sonicated (Misonix Sonicators) for 10 min. Later, the sample was centrifuged at 12,000g for 15 min at 4°C. The resultant supernatant was filtered using a Whatman filter (125 mm); the filtrate was then collected and stored at −20°C until LC-MS/MS analysis. Metabolite analysis was carried out using Liquid Chromatography (LC), followed by MS/MS analysis using a QTRAP 6500 mass spectrometer (ABSciex) coupled with an Agilent 1290 Infinity II LC system with a C18 RRHD Zorbax column (20 × 150 mm, 1.8 μm particle size; Agilent). Analyst software version 1.6.3 was used for data acquisition, and the Analyst Device Driver for setting the parameters for the analysis. The solvents used were 0.1% formic acid in MilliQ water as solvent A and 0.1% formic acid in 90% acetonitrile as solvent B; the flow rate was 0.3 ml/min. The LC method was performed for 30 min with the following gradient: 2% B for 0–1 min, 30% B for 10 min, 60% B for 11 min, 95% B for 13–17 min, and 2% B for 21–30 min. The Information Dependent Acquisition method (IDA) was carried out for the mass spectrometry data acquisition, which was built using Enhanced Mass Spectra (EMS) to Enhanced Product Ion (EPI) modes. The top five spectra from the EMS mode were used for analysis in the EPI (MS/MS) mode using high-energy Collision-Induced Dissociation (CID). The metabolite data acquisition was carried out in positive and negative modes at 4500 and −4500V, respectively, at a probe temperature of 450°C. The compound parameters were −100V of de-clustering potential (DP) and 40 V of collision energy (CE). All the LC-MS/MS acquisitions were carried out in triplicate. The XCMS online tool was utilized against the ‘Metlin’ metabolite database to analyze metabolite data. The LC-MS/MS result files uploaded for the analysis were in their original format, that is, wiff and raw format. The retention time correction was carried out using the Peak groups method with a step size of 1 m/z [43]. For the alignment of peaks, the allowable retention time deviation of 10 s, 0.5 min/frac, and 1 min samp for the minimum fraction and samples in a peak group were set. The matched filter algorithm was performed for the feature detection, with an m/z deviation of 5 ppm, 5–20 s of peak width, 0.1 Da of m/z difference, and a signal-to-noise ratio (S/N) of 3. For identifying and annotating metabolites, 5ppm tolerance and 0.15 absolute m/z error were set with isotopes and adducts against the Plant Database. The Adducts used were [M + H] and [M + Na] for the positive mode and [M − H] for the negative mode. Muthukapalli Krishnareddy Prasannakumar: Designed the study, analyzed the results, drafted and revised the manuscript. Hirehally Basavarajegowda Mahesh: Designed the study. Perumal Buela Parivallal: Performed the experiment, drafted and revised the manuscript Pramesh Devanna: Analyzed the results, drafted and revised the manuscript. Puneeth Makali Eregowda: Performed the experiment. Aditya Narayan Sarangi: Analyzed the results and drafted the manuscript. Manasa Kodihalli Govindaraju: Revised the manuscript. Sushil Kumar Middha: Revised the manuscript. Sahana Nagaraj Banakar: Performed the experiment and revised the manuscript. All authors have approved the final version of the manuscript. We are extremely grateful to Molsys Pvt. Ltd., Bangalore, India, for their guidance during metagenomic and metabolomic analyses. Financial support was provided by the Karnataka Watershed Development Project- II [Grant No. WBB/ED/KWDP-II (SUJALA-III)/REC-V-III/2-16-17]. The authors declare no conflict of interest. The Whole Metagenome Sequencing of the Panchagavya is deposited in the NCBI SRA database with the SRA Accession Number: SRR13235684 under the BioProject PRJNA683981 and BioSample SAMN17049059: PGMET (TaxID: 256318). The supplementary table is deposited in the FigShare database with an accessible link https://figshare.com/s/433ec5dfc5276c862f43. Supporting Information. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.