In order to conduct this research, samples of Sichuan dark tea were provided by the Sichuan Tea Industry Group Co., Ltd., and these were prepared from fresh leaves of tea plants collected in summer and autumn. During the fermentation process, the leaves of tea plants were mixed thoroughly to ensure homogeneity, and tap water was sprinkled in an adequate quantity to maintain 65–75% (w/v) solid contents and the temperature at 45–71°C. Samples were collected from tea piles at 15-day intervals and subjected to sensory evaluation as described by GB/T 23776-2009 (Gong et al., 2009). The fermentation process is over as the tea mass becomes reddish-brown and free of its stringent taste. Samples were collected during the fermentation process in triplicate and at the following time intervals; days 0 (YC), 8 (W1), 16 (W2), 24 (W3), and 32 (W4). The temperature of each fermentation tea pile was recorded from the center at the depth of 40 cm almost every day.

For DNA extraction, 5 g of each sample was suspended in 50 mL sterile Tween-NaCl buffer [0.9% NaCl (w/v), 0.05% Tween 20 (v/v) and 2% polyvinylpolypyrrolidone (w/v)] and homogenized by mixing thoroughly for 30 min via sonication at 4°C. The solution was further passed through a sterile gauze and subsequently centrifuged at 2,000 rpm for 2 min at 4°C to remove course material. To collect the microbial communities, the supernatant was shifted into a new tube and centrifuged at 10,000 rpm for 10 min at 4°C. Microbial genomic DNA was extracted from the precipitate using the E.Z.N.A.^TM HP Plant DNA Kit (Omega Bio-Tek Inc., GA, United States) according to the manufacturer’s protocol. Finally, DNA was purified using the E.Z.N.A.^TM Soil DNA Kit to remove impurities, such as polyphenol, that may have a negative effect on the PCR reaction (Yan et al., 2021).

The primer pair 338F: 5′-ACTCCTACGGGAGGCAGCAG-3′ and 806R: 5′-GGACTACHVGGGTWTCTAAT-3′ was employed to amplify the V3+V4 region of the bacterial 16S rRNA gene (Zeng et al., 2011; Apprill et al., 2015). The PCR reaction conditions were as follows: initial denaturation at 95°C for 3 min, denaturation at 95°C for 30 s, annealing at 55°C for 30 s, extension at 72°C for 45 s, 30 cycles in total; final extension was at 72°C for 10 min and held at 10°C for 10 min. The total PCR reaction system was 20 μL with the following ingredients; 4 μL of 5 × PCR buffer (with Mg²⁺), 2 μL of 2.5 mmol/L dNTP, 0.8 μL of 5 μmol/L P1 (338F), 0.8 μL of 5 μmol/L P2 (806R), 0.4 μL of 5 U/μL Taq Enzyme, 2 μL of DNA template, and 10 μL of ddH₂O. PCR products of all the same samples were mixed together and electrophoresed by running on a 2% agarose gel. Gel cutting was performed with a sterile sharp edge blade, and subsequently, the AxyPrepDNA Gel Recovery Kit (AXYGEN^®) was used for gel elution to recover the correct DNA band. Furthermore, PCR products with good quality were subjected to library construction and sequenced by using the robust Illumina MiSeq^TM PE300. Pair reads were merged into a single sequence by following the overlap relationship between PE reads and filtered for quality control to a avoid splicing effect. Finally, reads obtained by sequencing were checked for their quality, and sequence direction was corrected (Mukherjee et al., 2014; Yan et al., 2021). According to primer use, proper barcoding/labeling was performed on each step to avoid any error.

The USEARCH manual version 7.1¹ was deployed to extract non-repetitive sequences from optimized sequences, which is very useful to remove redundant calculations during data analysis², and to remove non-repetitive single sequences³. Similarly, UPARSE software (see text footnote 1) was deployed for clustering analysis at a 97% similarity index to get representative sequences of operational taxonomic units (OTUs) (Ye et al., 2017). The RDP classifier Bayesian algorithm was deployed to perform taxonomic analysis of OTUs at the 97% similarity index. Bacterial community composition present in every dark tea sample was classified and counted at each level of phylogeny, such as kingdom, phylum, class, order, family, genus, and species (Chen et al., 2016).

To construct a rarefaction curve, a specific number of individuals from each sample were randomly selected to evaluate their OTUs, and a number of individuals with their relevant OTUs were used to construct the curve. The rarefaction curve is used to compare the abundance of each specie in given samples with variable sequencing reads and to validate the amount of sequencing data. Rarefaction analysis of OTUs at the 97% similarity index obtained from 16S rRNA gene sequence reads via illumina MiSeq platform was performed at the mothur website⁴, and the curve was made using the R tool (Amato et al., 2013). To evaluate the sequencing depth index, α-diversity, such as Chao 1, ACE, Shannon, and Simpson of each sample, were individually calculated, and the abundance as well as the diversity of microbial communities in dark tea were analyzed at different fermentation stages (Rogers et al., 2016). A Venn diagram was constructed to count the number of common and unique OTUs in all samples (Fouts et al., 2012). According to the taxonomic analysis results, single or multiple samples were compared at each classification level, and the R tool was employed to construct community structure component diagrams, histogram combined analysis diagrams, and RDA principal coordinate analysis diagram (Lu et al., 2016; Zhou et al., 2017).

All data is explained in mean values of standard deviation (SD) and analyzed by one-way analysis of variance (ANOVA). A Duncan multiple-comparison test was applied to detect variations among means of all samples at a p-value < 0.05 level of significance. All correlation and path coefficient analyses were performed with SPSS Statistics 20.0 (SPSS Inc., Chicago, IL, United States) and Excel 2019.

We analyzed the diverse bacterial communities found in Sichuan dark tea at five different stages during the fermentation process by large-scale sequencing-based analysis of 16S rRNA gene sequences. The sequencing data of the bacterial colonies found in all samples collected at five different stages of fermentation are given in Table 1. The total number of sequencing reads of sample YC are 62,456, the total number of nitrogenous bases are 26,118,006 bp, and the average sequence length is 418.59 bp. Similarly, the total number of effective sequencing reads of sample W1 are 46,909, the total numbers of nitrogenous bases are 19,638,346 bp, and the average sequence length is 420.24 bp. Sequencing analysis revealed 40,399 full-length sequence reads in sample W2, the total number of nitrogenous bases is 17,285,142 bp, and the average read length is 427.75 bp. In sample W3; 43,362 full length sequencing reads were obtained with a total number of 18,361,869 bp nitrogenous bases, and the average read length is 423.51 bp. Finally, after elimination of raw and incomplete reads, 60,400 total reads were obtained in sample W4, containing 25,535,185 bp total number of nitrogenous bases, and the average length of sequence reads was 422.57 bp (Table 1).

The rarefaction curve was drawn to show the sampling depth of each sample, which can be used to evaluate either sequencing date and is sufficient to represent all bacterial strains present in each sample. The flatness of each rarefaction curve representing all five samples proved that the sampling process had adequately collected the experimental material (Figure 1). The bacterial diversity in the amount of each randomly collected sample was enough to construct an individual DNA library for each sample. The confidence level about the amount of each bacterial community structure in Sichuan dark tea samples was very high to accurately reflect the bacterial community.

The bioinformatics and statistical analyses were performed at the OTU level with a 97% similarity index. Furthermore, a Venn diagram was constructed to display compositional similarity among the different bacterial communities present in all samples, and the overlapping of different samples was intuitively displayed at different classification levels. A total number of OTUs obtained by clustering the valid data were 2,948, and these were further divided into 42 phyla, 98 classes, 247 orders, 461 families, 1,052 genera, and 1,888 species (Supplementary Table 1). Similarly, OTUs in each stage were 2,082, 2,116, 381, 399, and 580, respectively. Among all samples, YC and W1 samples shared the highest number of 1,211 OTUs, YC and W2 samples shared 15 OTUs, YC and W3 samples shared 25 OTUs, and YC and W4 samples shared 50 OTUs. The unique OTUs that were not common or overlapping in each sample were 360 in YC, 399 in W1, 70 in W2, 94 in W3, and 151 in W4, which represented 17.3, 18.9, 18.4, 23.6, and 26.0% of total OTUs, respectively (Figure 2). The total number of common OTUs among all five samples collected at five different stages of pile-fermentation were 103, which were 3.6% of the total number of OTUs in all samples. These results indicate that the bacterial community structure of Sichuan dark tea has undergone significant changes during pile-fermentation.

Bacterial community richness and diversity was explored in Sichuan dark tea samples collected at different stages of pile-fermentation. Noticeably, the sequencing coverage of all samples was higher than 0.9984, which represents the actual conditions of bacterial communities present in Sichuan dark tea samples collected at different stages of fermentation (Table 2). In general, the higher Shannon and lower Simpson indices represents the higher bacterial community diversity index in a sample (Huang et al., 2018). We observed the highest bacterial community diversity in Sichuan dark tea samples collected at YC and W4 stages, followed by the W1 period, and the lowest bacterial community diversity was observed in samples collected at the W2 stage.

To investigate the abundance of the bacterial community structure in Sichuan dark tea collected at different stages of fermentation, the Chao and ACE indices were measured. We observed the highest bacterial community abundance in dark tea samples collected at the YC and W4 stages, followed by samples collected at the W1 stage, and the lowest bacterial community abundance was observed in samples collected at the W2 and W3 stages. Our results show that the abundance and diversity of the bacterial community present in Sichuan dark tea samples was quite different at different stages of fermentation.

All Sichuan dark tea samples collected at different stages of pile-fermentation were analyzed to disclose the structure of the bacterial communities. Bacterial communities with higher abundance were annotated and clearly classified into 18 genera, and bacterial communities with relatively low abundance were merged together into a single category expressed by others (Figure 3). For example, bacterial communities with relatively high abundance at the genus level were Pseudomonas, unclassified_f_Enterobacteriaceae, Bacillus, Kocuria, Ralstonia, Lactobacillus, Staphylococcus, etc. (Figure 3). At the YC stage of fermentation, the significantly higher bacterial communities were unclassified_Enterobacteriaceae (26.38%), and the relatively high abundant bacterial communities at the genus level were Kocuria (15.72%), Streptomyces (5.92%), Pseudomonas (5.02%), Bacillus (4.73%), and Staphylococcus (4.19%).

At stage W1, the dominantly higher bacterial communities were unclassified_Enterobacteriaceae (26.38%), and the relatively high abundant bacterial communities at the genus level were Pseudomonas (4.91%), Lactobacillus (3.18%), Bifidobacterium (2.90%), and Bacteroides (2.10%). Similarly, the highest abundance of bacterial community at stage W2 was unclassified_ Enterobacteriaceae (51.55%), followed by Pseudomonas (41.50%), unclassified_Bacteria (2.15%), and the relative abundance of other bacterial genera were all less than 1%. Pseudomonas (61.61%) was an absolutely dominant bacterial strain at the W3 stage, and with relatively high abundance were unclassified_ Enterobacteriaceae (8.98%), Kocuria (6.86%), and unclassified_Bacteria (6.97%). At stage W4, the highest relative abundance was of Bacillus (33.20%), followed by Kocuria (18.63%), Lactobacillus (13.74%), Staphylococcus (7.92%), and Ralstonia (7.39%). We also observed that Pseudomonas and unclassified_Enterobacteriaceae displayed higher abundance at the YC, W1, W2, and W3 stages, and Bacillus, Kocuria, and Lactobacillus have higher abundance in the W4 period (Figure 4). Our results show that the bacterial diversity has obvious differences with the different piling-fermentation stages of Sichuan dark tea.

Redundancy analysis (RDA) is a type of PCA analysis that is constrained by environmental factors. For the pictorial representation of the relationship between bacterial flora at the genus level and environmental factors, samples and environmental factors were drawn on the same 2-D sequence diagram (Figure 5; Huhe et al., 2017). Among the four active compounds of Sichuan dark tea, polysaccharides have a positive correlation with total flavonoids, and caffeine has a positive correlation with amino acids. Furthermore, correlation was analyzed among the active compounds of dark tea, bacterial taxa, and different stages of pile-fermentation. We observed that effective formation of caffeine and amino acids predominantly exists at the following stages of fermentation: YC, W1, and W2.

Noticeably, the polysaccharides and total flavonoid contents were increased significantly during pile-fermentation stage W4. Similarly, the following bacterial genus Bacillus, Ralstonia, and Saccharopolyspora were significantly correlated with total flavonoid contents, and Kocuria was closely related with the formation of polysaccharides. The dominant bacterial genus at stage W3 was Pseudomonas, which displayed no relationship with any of the abovementioned four active ingredients of dark tea. These results prove that the biosynthesis of active ingredients of Sichuan dark tea were significantly affected by the composition of the bacterial community found in pile-fermentation.