Four different small-scale batches were prepared to evaluate the effect of the initial treatment of sukumo (produced by A.S. Tokushima, Shikoku, Japan) on indigo reduction initiation. Sukumo (76 g/batch) was treated with 60°C tap water (heat treatment alone, batch 1), 25°C tap water [neither heat nor high alkaline treatment, batch 2 (control)], 25°C wood ash extract (pH ≥ 11.0; high alkaline treatment alone, batch 3), and 60°C wood ash extract (both heat and high alkaline treatments, batch 4) at half volume (500 ml) for preparation. The next day, another half volume of 25°C wood ash extract (500 ml) was added to each batch (day 1). The wood ash extract was prepared by immersing 380 g wood ash (Quercus phillyraeoides A. Gray, Nagomi Co., Gobo, Wakayama, Japan) into 5 l tap water. The solution was boiled for 10 min and then the supernatant was extracted and used for the aiming temperature after the wood ash settled in the water. Triplicate batches were prepared to confirm the difference in the initiation of indigo reduction and changes in the dyeing intensity based on the initial preparation.

The fermentation batches were kept at 26°C in a temperature-regulated room and gently stirred once a day using a laboratory spoon with a top width of 30 mm, to prevent the localization of acid produced by the bacteria in the fermentation fluid. Wheat bran (1 g) was added on days 5, 19, 51, 85, and 194. The pH and redox potential (ORP) of the fermentation fluids were measured with a D-71 pH meter (Horiba, Kyoto, Japan) and a D-75 pH/ORP/CO meter (Horiba), respectively. The pH of the fermentation fluid was maintained between 9.67 and 11.2 using Ca(OH)₂ to increase pH. The reducing state was tested by the dying ability of the fermentation fluid after dipping a small piece of cotton cloth into the fluid for 1 min and exposing it to air. After 2–5 min, the cloth was rinsed and dried. Dyeing intensity was measured by scanning the dyed cloth and analyzing the image using Mathematica (version 12.2). The intensity was expressed as L*a*b* value, which is the square root of L*² + a*² + b*² in International Commission on Illumination and its System Colorimetry (CIE) L*a*b* color space where L*a*b* color space, L* represents lightness, and a* and b* represent chromaticity, which indicates hue and saturation, respectively. The extension of a* and b* indicates color direction as follows: a*, red; −a*, green; b*, yellow; and −b*, the blue. Before sampling, the fluid was stirred using a laboratory spoon to homogenize the fluid, and 700 μL aliquots were collected and stored in 35% (w/v) glycerol at −80°C until use.

For DNA extraction, the FastDNA™ Spin kit for soil (MP Biomedicals, Santa Ana, CA, USA) was used according to the manufacturer’s instructions. The bacterial 16S rRNA gene sequence of the V3–V4 region (341F–805R) was amplified using the primer pair (5′-3′): V3V4f_MIX (ACACTCTTTCCCTACACGACGACGCTCTTCC-GATCT-NNNNN-CCTACGGG-NGGCWGCAG) and V3V4r_MIX (GTGACTGGAGTTCAGACGTG-TGCTCTTCCGATCT-NNNNN-GACTACHVGGGTATCTAATCC) purchased from Bioengineering Lab. Co. Ltd. (Sagamihara, Kanagawa, Japan). The primers consisted of an overhang necessary for the 2nd PCR, 0–5 bases of random sequences (described as N; adaptor), and a gene-specific sequence necessary for the amplification of 16S rRNA. The random sequences were used for quality control. The PCR solution (40 μL) consisted of 4 μL of 10× Ex buffer (TaKaRa Bio, Otsu, Shiga, Japan), 3.2 μL dNTPs (TaKaRa Bio), 2 μL each of the forward and reverse primers, 2 ng extracted DNA sample, and 0.4 μL of 5 U/ml Ex Taq polymerase (TaKaRa Bio). The cycling conditions were as follows: 94°C for 2 min; 25–35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s; extension of 72°C for 5 min. The amplification products were purified using the QIAquick PCR Purification Kit (Qiagen, Mississauga, ON, Canada) following the manufacturer’s instructions. Product identity and quality were confirmed by agarose gel electrophoresis.

For next generation sequencing (NGS) analysis, the first PCR products were submitted to Bioengineering Lab. Co. Ltd. (Sagamihara, Kanagawa, Japan). The second PCR was performed with index-adapted primers to generate paired-end libraries (2 × 300 bp) for NGS using the MiSeq platform (Illumina, San Diego, CA, United States) and the MiSeq reagent kit v3 (Illumina).

The 1st sequence primers and N-sequences were removed from the resulting gene sequences using Cutadapt version 1.18, the fastq files obtained from NGS were processed using QIIME2 ver. 2020.2 (Bolyen et al., 2019). Paired end read merge, quality control, and creation amplicon sequence variants (ASVs) were performed using DADA2 (Callahan et al., 2016). Taxonomic classification was conducted using the feature-classifier for the 16S rRNA gene sequence based on the primer pair of 341F–805R created based on the Silva database (Quast et al., 2013; Yilmaz et al., 2014). For updated taxonomic identification, representative sequences were used for a BLAST search in National Center for Biotechnology Information (NCBI¹) (Miura et al., 2019). Rarefaction curves of observed Operational Taxonomic Units (OTUs) and Shannon index of diversity (Shannon, 1948) were plotted according to the QIIME2 alpha diversity script. Linear discriminant analysis (LDA) effect size (LEfSe) was conducted to identify the biomarkers for each batch by the Galaxy of Huttenhower Lab.² Redundancy analysis (RDA) between the microbial community and environmental parameters was performed with R package (vegan ver. 2.6–4, ggrepel ver. 0.9.2, ggplot2 ver. 3.4.0, ggpubr ver. 0.5.0). Networks of the bacterial community based on the relative content change trend was analyzed as previously described (Tu et al., 2021). Predictive functions of the metagenomes were estimated using PICRUSt2 (Douglas et al., 2020). Based on the PICRUSt2 and QIIME2 results, taxa-functional relationships were analyzed using BURRITO (McNally et al., 2018). Correlation coefficient between changes in dyeing intensity of samples and corresponding changes in the intensity of the KEGG ortholog in “pred metagenome unstrat” file produced by the PICRUSt2 results was estimated.

From the 28 samples, 1,442,013 raw sequences were obtained, with an average of 51,500 sequences per sample (days 1–7) following fermentation preparation. After adapter trimming, and filtering for quality checking and chimeric exclusion, 687,015 sequences (average: 24,536 per sample) were obtained.

Triplicate experiments confirmed the reproducibility of initiation of indigo reduction based on the treatment of sukumo. Microbiota analysis was performed with one sample set. Although wood ash treatment produced only approximately one unit difference in pH (batches 3 and 4; pH 10.7) compared with tap water treatment (batches 1 and 2; pH 9.7–9.8), contents of “others” were much higher in batches 1 and 2 than those in batches 3 and 4 after day 2. “Others” represents the sum of taxa present in less than 4.1% of all the samples (miscellaneous taxa), indicating that initial exposure of the microbial community in sukumo at approximately one unit higher than pH 9.7–9.8 strongly converged the microflora. However, the differences in the content of “others” were not observed between batches 1 and 2 on day 2, suggesting that heat treatment of 60°C of tap water was not very effective in converging the microflora. The lowest percentage of “others” (5.8%) and the highest percentage of the relative abundance of obligate anaerobic Alkaliphilus oremalandii (similarity: 98.8%–99.3%; abundance: 82.8%) in batch 4 at day 4 were observed (Figure 1). Moreover, the highest percentage of “others” (41.4%) and the lowest percentage of the relative abundance of A. oremalandii in batch 2 on day 2 were observed. The higher abundance (>50%) of A. oremalandii in batch 4 continued longer than in the other batches. These results suggest a strong relationship between the selection intensity of existing microorganisms in sukumo and relative abundance of obligate anaerobes in the initial treatment and the successive environmental changes.

The effect of the initial high pH treatment that was maintained was more pronounced in batches 2 and 4 than 1 and 2. However, the abundance of the alkaliphilic obligate anaerobic bacterium A. oremalandii was increased until days 3 and 4 in batches 3 and 4, respectively. In addition, the ratio of “others” decreased until day 4 in batch 4. Considering the ratio of “others,” which constitute miscellaneous low ratio constituent microorganisms, the converging effect for the microbial community was further enhanced in batches 1 and 2 by the addition of wheat bran on day 5.

Considering the rapid decrease in ORP, the remaining microorganisms in batches 1 and 2 after the preparation consumed oxygen more rapidly than those in batches 3 and 4 (Figure 1, Supplementary Figure S1A). Meanwhile, a pH of 11.2 in batch 3 and 11.1 in batch 4 deterred the metabolism of oxygen metabolizable microorganisms (Supplementary Figure S1B). However, the selected microorganisms of the initial treatments in batches 3 and 4 produced acid, which reduced the pH and allowed them to consume oxygen.

The earliest initiation of indigo reduction was observed in batch 3. The appearance of Alkalihalobacillus macyae (similarities: 98.1%–98.4%; content ratio: 5.3%) on day 3 in batch 3 was considered a contributing factor to the earlier indigo reduction. Initiation of indigo reduction was observed following the appearance of Alkalicella caleodiniensis (formerly Proteinivoraceae; 99.0%–100%; Tu et al., 2019b) in batches 1, 2, and 4. Indigo reduction was initiated on day 6 in batches 1 and 2. The appearance of Tissierella spp., Amphibacillus spp., and Atopostipes spp. may also contribute to the initiation of the indigo reduction. Although characterizations of these taxa were our unpublished results, these bacteria are indigo-reducing bacteria. This conversion effect by continuous high pH and low ORP was reflected in the ratio of “others” on day 7.

Differences in the microbiota due to treatment method of sukumo was most noticeable from day 2, when the effect of treatment began, prior to the effect of adding wheat bran, which appeared after day 6. Biomarker detection in each treatment batch was performed during days 2–5 (Figure 2). Heat treated (60°C) batch 1 compared with non-treated batch 2 showed that Actinobacteria (Brachybacterium), Bacillaceae (Mesobacillus, Pseudogracillibacillus, Paenibacillus, Salinicoccus, and Planifilum), and Lactobacillales (Atopostipes) were inhibited by treatment, while Amphibacillus (Bacillaceae) was enhanced (Figure 2A). High pH (pH ≥ 11.0) treated batch 3 compared with batch 2 showed that Actinobacteria (Yaniella), Bacillaceae (Bacillus, Mesobacillus, Oceanobacillus, Pseudogracillibacillus, Sinibacillus, Paenibacillus, Salinicoccus, Planifilum and Thermoactinomyce), and Caldicoprobacteraceae (Caldicoprobacter) were inhibited by the alkaline treatment, while Bacillaceae (Alkalihalobacillus and Sutcliffiella) and Eubacteriales were enhanced (Figure 2B). Heat and high pH treated batch 4 compared with batch 2 showed wider range of inhibited taxa compared with alkaline treatment (Figure 2C). Indeed, Actinobacteria (Brachybacterium and Nocardiopsis), Bacillaceae (Paenibacillus) and Lactobacillales (Atopostipes) were further inhibited compared with alkaline treatment alone; meanwhile, Mesobacillus, Sutcliffiella, and Eubacteriales (Alkaliphilus) were enhanced by the treatment. Therefore, high pH treatment had larger impact than heat treatment on microbiota, while heat combined with high pH treatment showed the biggest impact on microbiota among the three treatments.

Since batches 1 and 2 showed similar changes in their microbiota on 7, we focused on batches 2, 3 and 4 from days 10–209 for subsequent microbiota analysis. A total of 873,986 raw sequences were obtained from a total of 18 samples (average: 48,555/sample), one from each day (days 10–209) after fermentation preparation. After adapter trimming, filtering for quality checking, and chimeric exclusion, 353,688 sequences (average: 19,649/sample) were obtained.

The batches exhibited clear dyeing on day 7 and further strong dyeing on day 10. This may be attributed to the effect of the maintained high pH from day 2, low ORP, and the introduction of wheat bran on day 5 across all batches. Consequently, all three batches exhibited relatively strong dyeing intensity and the constituents in the microbiota were relatively similar after day 10, indicating that all the batches reached a stable state for indigo reduction in fermentation, regardless of the initial treatment.

Although no strong selection was observed in batch 2 on day 0 compared with the other batches, the microbiota converged appropriately on day 7 due to the continuously high pH starting on day 2 with the addition of Ca(OH)₂, low ORP and wheat bran. Furthermore, the establishment of desirable microbiota was observed on day 10, and relatively similar microbiota among the batches was observed on day 75 in this batch (Figure 3). The microbiota in batch 2 were mainly comprised of A. caledoniensis (similarity: 99.0%–100%; abundance ratio: 11.5%–17.8% on days 10–75), A. oremlandii (similarity: 98.8%–99.3%; abundance ratio: 8.2%–27.4% on days 10–75), Enterococcus gallinarum (similarity: 99.8%–100%; abundance ratio: 6.6%–16.9% on days 29–75), Atopostipes suicloacalis (similarity: 98.1%–99.1%; abundance ratio: 4.0%–12.1% on days 10–49) and Amphibacillus xylanus (similarity: 99.3%–100%; abundance ratio: 3.1–11.3 on days 10–49). The microbiota of batch 2 was different with regards to the lower abundance ratio of A. xylanus, the earlier appearance of Alkalibacterium iburiense and the higher ratio of “others.”

The microbiota of batch 3 was relatively similar to that of batch 4. However, the present ratio in Tissierella creatini (similarity: 94.5%) to A. oremalandii (similarity: 98.8%–99.3%) differed in the taxa comprising obligate anaerobes. Moreover, the ratio in the fluctuations of A. xylanus (similarity: 99.3%–100%) was also different.

On day 209, all three batches contained Proteobacteria in their microbiota, which comprised more than 20% in the microbiota ratio suggesting that although sukumo introduced Gram-positive bacteria into the fermentation system. However, Gram-negative bacteria exhibited superior adaptation ability compared with Gram-positive bacteria in the liquid environments. In addition, the ratio in Alkalibacterium spp. was increased, while other taxa existed categorized in Lactobacillales [e.g., Atopostipes suicloacalis and Enterococcus gallinarum (99.8%–100%)] were decreased. In the obligate anaerobes category, the ratio of preferred taxa on day 75 (e.g., Alkakucella caledoniensis and Alkaliphilus oremlandii) was decreased and facultative anaerobic Erysipelotrichaceae (93.2% similarity with Erysipelothrix inopinate) and Corynebacterium efficiens (similarity: 99.5%) appeared.

Changes in alpha diversity at a sampling depth of 8,825 reads from the rarefaction curves based on the observed OTUs and Shannon index of the four batches are shown in Supplementary Figure S2. High alkaline treated batches 3 and 4 showed significant decreases in both the observed OTUs and Shannon index from days 1 to 2. Moreover, the decrease in these values by heat (60°C) treated batch 1 was not as intense as a high alkaline treatment (batches 3 and 4). Although the high pH treatment (batch 3) and both high pH and heat treatment (batch 4) were performed at the preparation of fermentation (day 0), the decrease in the diversity of microbial constituents continued until days 2 and 4, respectively. After the decrease in diversity, recovery was observed on days 3 and 5 in batches 3 and 4, respectively.

Although the Shannon index of batches 1 and 2 did not changed the diversity from days 3–7, the observed OTUs gradually decreased. Furthermore, the introduction of wheat brane caused a further decrease in the observed OTUs. Surprisingly, all the alpha-diversity of the four batches, except for observed OTUs in batch 1, converged to almost the same values. These results are comparable with the similar dyeing intensity among the four batches and the ratio changes of the “others.”

The relationship between the changes in microbial communities among the three batches (batches 2–4) and the influence of environmental parameters were estimated using RDA (Figures 4A,B). The high impact of hot wood ash extract in batch 4 on microbiota was reflected in the distance between day 1 and after day 2 (Figure 4A). Meanwhile, the higher impact of the high pH in batch 3 compared with the heat treatment in batch 1 was reflected by the distances between days 1 and 2. Changes in microbiota of batch 1 was very similar that of the non-heat-treated batch 2. High pH treatment (batches 3 and 4) produced more rapid transitional changes from days 1 to 2. The microbiota treated with high pH was changed and influenced by low ORP from day 2. However, susceptibility to low ORP decreased in batches 1 and 2 compared with 3 and 4. The changing velocity of each microbiota was lower during days 2–5 compared with until day 2. However, the velocity became increased by the wheat bran introduced at day 5. Although the microbiota of each batch was differed, the wheat bran strongly shifted the microbiota towered indigo reducing state from day 5. The contribution to indigo reducing state remains unknown, however the appearance of Alkaliphilus oremlandii is strongly linked to the high pH and low ORP. Furthermore, Tissierellaceae, Atopostipes suicloacalis, Amphibacillus and Alkalicella caledoniensis are core members at days 6–7 (Figure 4A). Batch 3 started indigo reduction on day 3, which was attributed to the contribution of Bacillaceae-1. Although the transitional change from days 1 to 5 of each batch differed, all day 7 positions were similar. This is comparable with the changes in the Shannon index of alpha diversity in the 4 batches converged to nearly the same values on day 7.

All batches exhibited similar changing tendencies from day 10 compared with the changes in days 1–7 (Figure 4B). The changing direction was the same with similar positioning of microbiota. Figure 4B shows that the change rate of microbiota was low based on the slow decomposing substrate, wheat bran, while the continuous introduction of the substrate is strongly associated with the maintenance of the dyeing intensity during the aging of each microbiota. However, the dyeing intensity gradually decreased with aging during fermentation. During the convergence of each microbiota toward an aging state, successive changes in several functional-redundant-indigo-reducing different taxa appeared and sustained indigo reduction for a long duration. Around day 29, Alkaliphilus oremlandii, Amphbacillus Atopstipes suicloacalis, and Alcalicela caledoniensis were the core bacteria that sustained this ecosystem. However, after day 29, in addition to Atopstipes suicloacalis, and Alcalicela caledoniensis, Tissierellaceae and Alkalibacterium were the core members in the microbiota. Meanwhile, after day 75, Eubacteriaceae, Eubacteriales, and Enterococcus gallinarum were added as the core members in the ecosystem, and at day 209, Gram-negative bacteria became the core members of the microbiota. The isolated positioning of the microbiota in this series of change to day 209 was in the opposite to the high pH, possibly due to local changes in pH of the fermentation fluid.

We analyzed the interaction of existing bacteria on days 1–7 (Figure 5A) and 10–209 (Figure 5B) and found insignificance differences in the detail between the batches, with similar tendencies. Therefore, the results from batch 4 are presented as the representative. Firstly, a relatively strong network was constructed based on the survived oxygen-metabolizing taxa such the 3 taxa of Actinobacteria and Sinibacillus, and 2 taxa of Oceanobacillus, Salinicoccus kunmingensis, Bacillus infernus, Planifilum composti and Caldicoprobacteraceae. This initial network was destroyed by the high pH and the low ORP. Subsequently, Msobacillus fermenti, Sutcliffiella cohnii and Alkaliphilus oremlandii appeared. Toward an enhanced indigo reducing state at day 7, network consisted of 2 taxa of Amphibacillus, Tissierellaceae, Atopstipes suicloacalis, and Alcalicela caledoniensis was constructed. After day 10, the similar network, constructed until day 7, was succeeded; however, the inter positive relationship had weakened. The gravity of the network gradually moved to the latter stage as shown in Figure 5B. Networks containing Gram-negative bacteria established after day 75 were relatively strong.

The predictive function of the metagenomes was estimated to clarify the prerequisite functional abundances for the initiation of indigo reduction. Since the initiation of indigo reduction was observed only in batch 3 on day 3, the functional abundance ratio between batch 3 and the other batches was estimated. The ratio in the functional subpathways containing a ratio > 1.05 in batch 3(B3)/batches 1(B1), B3/batch 2 (B2), or B3/batch 4 (B4) on day 3 are listed in Supplementary Figure S3. Among the presented subpathways, “phosphotransferase system (PTS)” and “starch and sucrose metabolisms” followed by “Prokaryotic defense system” and “nicotinate and nicotinamide metabolism” were identified. In addition, although their ratios are not very high, the subpathways related to the reconstitution of cellular functions and metabolism related to electron carriers and motility were also found.

Although the dyeing intensities during the initiation days were not always strong, batches 2–4 exhibited strong dyeing intensities compared with day 10. To understand the prerequisite functional subpathways that initiate dyeing, the functional abundance ratios of day 10 (beginning of intense dyeing)/day 2 (before initiation of dyeing) was estimated. In addition, ratio of day 16 (stable phase in intense dyeing)/day 2 in batches 2–4 were estimated (Supplementary Figure S4). To understand the important functions for enhancing dyeing intensity, the functional abundance ratios of day 10/day at the initiation of dyeing and day 16/ day at the initiation of dyeing in batches 2–4 were estimated (Supplementary Figure S4). The listed categories and super pathways were similar to the list in Supplementary Figure S3, except for the subpathways in super pathways of cell motility and metabolism (category of unclassified). Subpathways in the category of “unclassified” are important based on the subpathway in either batch 2 or 4. Increased items from Supplementary Figure S3 were reflected in the diversity change in the microorganisms related to indigo reduction. Although the listed subpathways were increased, it was considered that the “phosphotransferase system (PTS)” and “starch and sucrose metabolisms” remain important for the initiation of dyeing and enhancement of the dyeing intensity. In addition, although their ratios were lower in the two subpathways, “Prokaryotic defense system” and “nicotinate and nicotinamide metabolism” were relatively high.

The functional abundance of the subpathways, which were picked up in Supplementary Figure S4 and the contribution of the subpathways by major constituted taxa were estimated (Supplementary Figure S5) on day 3. The initial fermentation stage in batch 3 could be explained by the functional abundances of the subpathways of “phosphotransferase system (PTS)” and “starch and sucrose metabolisms.” The specific abundancies in batch 3 were attributed to the genus Alkalihalobacterium, which specifically existed in this batch. This genus exhibited a higher contribution than other taxa considering the existing abundance ratio. The contribution for the listed subpathways by the genus Alkaliphilus was relatively low considering the existing ratio. However, this genus contributed to “amino acid metabolisms” rather than that of the carbohydrates.

The time from the preparation (day 0) to the initiation of indigo dyeing was different depending on the initial preparation. However, a moderate and similar intensity of dyeing was observed on day 7 in all batches (Figure 1). All four batches exhibited high functional abundance ratios in the three subpathways of “phosphotransferase system (PTS)” and “starch and sucrose metabolisms” (Figure 6). These findings consistent with the dyeing results. Taxa that occupied major constituents were increased as the abundance ratio of Alkaliphilus was decreased in all the batches. The functional abundance of subpathways by genera Amphibacullus and/or Atopstipes contributed to fulfilling the abundance of the important subpathways involved in indigo dyeing. The contributions of subpathways important for dyeing in obligate anaerobes (e.g., Alkaliphilus and Umitella) were relatively low considering the present ratio. However, the contribution by genus Tissierella in batch was high considering the existing ratio. These obligate anaerobes most likely did not significantly contribute to the metabolisms of carbohydrates in this ecosystem.

From days 7 to 10, the dyeing intensity increased; however, no significant change was found in the ratio of the functional abundance ratio, which was important for indigo reduction (Figure 6 and Supplementary Figure S6). The increase in staining intensity was likely due to the qualitatively improved metabolic function of the microbiota rather than a quantitative increase in abundance. The decrease in the ratio of Alkaliphilus and increase in the ratio of Amphibacillus and/or Alcalicella and Tissierella were improved during days 7–10. The former exhibited a relatively low contribution for the important subpathways for indigo reduction, while the latter exhibited relatively high contribution for the important subpathways for indigo reduction (Supplementary Figure S6). In this way, important functions for the reduction of indigo and material circulation in the system shifted to the main constituent taxa on day 10.

After day 10, the functional abundance ratio of the subpathways of “phosphotransferase system (PTS)” and “starch and sucrose metabolisms” remained relatively high. The established effective metabolic pathway was maintained after day 10. The contribution of the subpathways by major taxa was also estimated in the microbiota at day 49 (Figure 7). The highly contributed taxa for the important subpathways were different based on the batch. Alcalicella, Tissierella, Enterococcus, and Ampibacillus were the major contributors to the important subpathways in batch 2. Meanwhile, Alcalicella, Enterococcus, and Amphibacillus were the major contributors to the indigo reduction in batch 3.

Table 1 lists the KEGG orthology correlated with staining intensity, which shows the correlation coefficients of ≥ 0.69 selected up from KEGG orthology items of 5,458 described in “pred metagenome unstrat” file created with PICRUSt2. Although the functional category differed from the results of the BURRITO analysis, the functions described the “phosphotransferase system (PTS),” “starch and sucrose metabolisms,” “Prokaryotic defense system,” and “nicotinate and nicotinamide metabolism” were found in the list. Among the KEGG orthology items, seven NAD(P)-dependent oxidoreductases were selected and the related taxa in each orthology and ratio “CountContributedByOTU” at the initiation of indigo reduction in batches 3 [from days 2 to 3 (D3/D2)] and 4 was [from days 3 to 4 (D4/D3)] were estimated using the “pred metagenome contrib.legacy” produced by PICRUSt2 analysis (Table 2). Enhancement of H⁺-translocating NAD(P) transhydrogenase subunit beta [EC:1.6.1.2 7.1.1.1] (K00324, K00325) of Alkalihalobacillus macyae and NADP-dependent alcohol dehydrogenase [EC:1.1.1.2], uronate dehydrogenase [EC:1.1.1.203], and D-lactate dehydrogenase [EC:1.1.1.28] of Atopostipes suicloacalis caused the initiation of indigo reduction in batch 3 (Table 2). Meanwhile, enhancement of formylmethanofuran dehydrogenase subunit E [EC:1.2.7.12] and anaerobic carbon-monoxide dehydrogenase catalytic subunit [EC:1.2.7.4] of Alkalicella caledoniensis might have contributed to the initiation of indigo reduction in batch 4 (Table 2). In batches 1 and 2, changes toward the indigo reduction are attributed to the increase in the ratio of Atopostipes suicloacalis and Alkalicella caledoniensis (Figure 1).