Shifts in Bacterial Community Associated with Green Manure Soybean Intercropping and Edaphic Properties in a Tea Plantation

: The continuous cultivation with excessive application of agrochemicals has led to the deterioration of soils. Incorporating leguminous green manure (GM) was found to improve the physicochemical and biological properties of soils. However, the inﬂuence of GM soybean intercropping on the temporal distribution of bacterial communities in strongly acidic soils is less explored. In this study, a nine-month ﬁeld trial of soybean intercropping was conducted in a deteriorated tea plantation. This was used to test the hypothesis that GM treatment ameliorates tea-cultivated environments through changing bacterial communities as well as edaphic properties. GM treatment was demonstrated to increase both functional and population diversity during all the months that were considered. An alteration of life strategies exhibited by bacterial communities in GM treatment was observed, which shifted from oligotrophy ( Acidobacteria , Chloroﬂexi , and the candidate phylum WPS-2) to copiotrophy ( Bacteroidetes and Proteobacteria ). This also contributed to the remarkable increase in metabolic potential of bacterial communities toward all six carbon source categories. The potentially versatile organic matter decomposers and/or plant growth-promoting bacteria, such as Burkholderiaceae , Chitinophagaceae , Sphingobacteriaceae , and Sphingomonadaceae bacteria, were identiﬁed as the most effective biomarkers in GM treatment. These bacterial groups showed strong correlation with soil pH; organic matter; and available K, Ca, and Mg. The increased diversity, metabolic potential, and copiotrophic taxa provided insight into the beneﬁts brought by soybean intercropping, with enhanced community stability, facilitated nutrient cycling, and microbe–plant interactions in the strongly acidic tea plantation. formal analysis, S.-H.L.; investigation, F.-T.S. and S.-H.L.; curation, preparation, F.-T.S.; writing—review F.-T.S.; supervision, F.-T.S.; project F.-T.S.; funding acquisition, F.-T.S.


Introduction
The microbes associated with plant roots is enormously diverse, and studies on plant-microbe interactions revealed that plants are able to shape their rhizosphere microbiome [1][2][3]. Several strategies have been highlighted to enrich microorganisms beneficial to plant growth and health [4][5][6]. These included changing the quantity and/or quality of root exudates via plant breeding or genetic engineering, or introducing and stimulating beneficial microorganisms by soil or crop management. Among various agricultural manipulations, the application of leguminous green manure (GM) has been found to restore soils through improving not only physicochemical but biological properties of soil [7,8].
This agricultural practice was also demonstrated to maintain microbial diversity, as well as increase the abundance of biologically and agronomically relevant groups and soil enzyme activities involved in N and P cycling [9,10].
The tea plant (Camellia sinensis L.) is a popular beverage and cash crop widely cultivated in Asian areas. In tea plantation soils, the most commonly found bacterial taxa at the phylum level were Acidobacteria, Actinobacteria, Chloroflexi, and Proteobacteria [11][12][13][14][15]. The dynamics of bacterial communities in the chronosequence of tea plantations have been studied, and beneficial bacteria were found to decrease after long-term tea cultivation [11,12]. Wang et al. [14] also indicated that the age of tea plantations served as an important factor changing the composition of soil aggregates, which further influenced bacterial metabolic activity and diversity. The tea plantation soil-associated bacterial community showed higher metabolic potential and diversity compared with other rhizo-compartments [15]. Besides, organic amendment was observed to improve microbial diversity, enhance community multi-functionality, and positively correlate with soil quality index [16,17].
The main purpose of microbial ecology studies within the rhizosphere is to facilitate the development of agricultural systems that deliver high levels of food security while reducing the environmental impacts [18]. Investigation on bacterial communities shaped by roots as well as soil properties helps to unravel microbe-plant interactions and develop more effective soil management strategies to sustain plant health. In Taiwan, there are 12,266 ha of tea plantation manipulated with conventional farming [19]. The continuous cultivation with excessive application of agrochemicals has led to the deterioration of soils such as the increase in acidity and decrease in organic matter content. Adopting soybean intercropping in the tea plantation may help to improve soil quality, which also serves as an alternative practice to suppress weed growth instead of using herbicides. However, the influence of GM treatment on bacterial communities in strongly acidic soils is still unclear. In this study, a field trial was conducted to investigate the effects of soybean intercropping on bacterial communities in a deteriorated tea plantation. Community-level physiological profiling and 16S rDNA amplicon sequencing were performed to determine the temporal distribution of bacterial community during a nine-month study period. These were used to test the hypothesis that GM treatment ameliorates tea-cultivated environments through changing bacterial communities as well as edaphic properties.

Site Descriptions and Sampling
This study was conducted during June 2020 and March 2021 in a tea plantation located in Mingchien (23 • 50 30.0 N, 120 • 38 33.5 E, elevation 369 m), Nantou, i.e., a county which occupied more than half of the tea-cultivated areas in Taiwan. The soil at the study site was classified as loam (Typic Paleudults) under the USDA soil classification system [20]. The tea cultivar grown in this field is Sijichun, a commonly found cultivar grown in this area and with an age of almost three years. The plant used in GM treatment is soybean (Glycine max L.) cultivar Tainan no. 7, a fast-growing variety with high biomass and superior drought and cold resistance [21]. During the study period, in November 2020, a total of 20 kg of vegetative meals (castor meal, peanut meal, and rapeseed meal) (N-P 2 O 5 -K 2 O: 6-2.1-2.1, pH 6.1) (consisting of 90% organic matter) was applied in the field with 0.1 ha in size (ca. 600 plants). Chemical pesticides, such as diafenthiuron and λ-cyhalothrin, were used in pest control.
As for GM treatment, the inter-row was manually sown with soybean seeds on 9 June 2020 at a seed rate of 30 kg ha −1 . A bare ground (BG) control which was mulched with polyethylene plastic film and without soybean intercropping was used for comparison. The experiment was set up on 6 experimental plots (both treatments in triplicates), and the size of the individual plots was 9 m × 0.9 m (Figure 1). After GM treatment, the soybean rhizosphere samples were collected with a trowel from 0-20-cm-deep soils within 30 cm from the base of the soybean plant. The BG control samples were collected in parallel from the same depth of soils. Three random sampling points within each experimental plot were mixed into one sample, and six samples were obtained from six experimental plots for each sampling time point. The collected samples were placed in ice boxes, brought back to the laboratory, and stored at 4 • C. For soil chemical properties analyses, a total of 18 samples which represent three sampling time points (after three-, six-, and nine-month soybean intercropping) were collected. Soils were immediately air-dried and sieved (20 mesh for pH, EC, and nutrients determination; 35 mesh for organic matter determination) before analyses. For bacterial community analyses, a total of 30 samples which represent five sampling time points (after one-, two-, three-, six-, and nine-month soybean intercropping) were collected and used immediately. During the study period, the monthly average temperature ranged from 15.4-28.7 • C and the monthly precipitation ranged from 0-210.5 mm (Supplementary Materials Figure S1). The average daily precipitation of 3 days before the sampling dates (2 July 2020; 13 August 2020; 15 September 2020; 8 December 2020; 29 March 2021) was 0, 21, 1, 0, and 0 mm, respectively.

Analyses of Soil Chemical Properties
The pH and EC were measured as described in [22]. The Walkley-Black wet oxidation method was used to determine organic matter content [23]. Total N was determined by the Kjeldahl digestion method [24], as described in [25]. Available P, K, Ca, Mg, Fe, Mn, Cu, and Zn were analyzed by the Mehlich No. 3 procedure [26], as described in [22].

Determination of Functional Diversity and Metabolic Potential of Bacterial Communities
Community-level physiological profiling was performed based on the carbon source utilization patterns determined by Biolog EcoPlate TM (Biolog Inc., Hayward, CA, USA). Soil suspensions were prepared, 1000-fold diluted, and inoculated into each microplate well, as described in [22]. The optical density at both 590 and 750 nm was read on a microplate reader at 0, 24, 48, 72, and 96 h. The final values used to denote activity in each well were calculated, as described in [27]. Substrate richness (S) was calculated by counting total number of carbon substrates oxidized by individual treatment on Biolog EcoPlate TM (wells with OD higher than 0.25). The average well color development; substrate average well color development; and diversity parameters, including Shannon's diversity index (H) and evenness (E), were calculated, as described in [28].

Determination of Population Diversity and Structure of Bacterial Communities
Total DNAs of soils were isolated using a DNeasy PowerSoil Pro Kit (Qiagen Inc., Germantown, MD, USA) following the manufacturer's recommendations. Primer set 341F-805R with a barcode was used to amplify the V3-V4 region of 16S rDNA, which was further purified and used in sequencing libraries generation, as described in [22]. The libraries were sequenced on an Illumina Miseq platform and 300-bp paired-end reads were generated at Genomics BioSci&Tech Ltd., Taiwan. Raw sequence reads were submitted to the NCBI SRA database under BioProject PRJNA716948. Paired-end reads were trimmed using the Cutadapt program, merged using FLASH v1.2.11, and analyzed by FastQC v0.11.5 and MultiQC v0.9 for quality control. Operational taxonomic unit (OTU) picking was performed using mothur v.1.39.5 with 97% identity. Chimera sequences were identified using UCHIME v4.2 with Gold database, as reference data. The sequences derived from different OTU were used to perform species annotation after compared with representative sequences from the SILVA database [29]. According to species annotation, the statistical amount of sequences of every sample in each classification level was calculated. The relative abundance (percentage of each OTU in total OTUs) of dominant bacterial taxa (>2%) at the phylum and the family level was expressed. In order to evaluate bacterial diversity within a community from each treatment (alpha diversity), the Chao1, Shannon index, and InvSimpson index were calculated with QIIME software [30].

Statistical Analysis
All the analyses were performed in triplicates for the BG control and GM treatment during the study period, and the results were presented as mean values. One-way ANOVA (analysis of variance) and the Tukey HSD test (p < 0.05) were used to evaluate the significant differences between treatments using XLSTAT statistical software (New York, NY, USA). Canonical correspondence analysis (CCA) and Pearson correlation analysis were used to elucidate the relationships between dominant bacterial taxa and soil properties or metabolic potential. A hierarchically clustered hea -map was generated based on Euclidean dissimilarity metrics using the TBtools software [31]. Principal component analysis (PCA) based on Pearson correlation matrix was performed to determine the relationships between treatments and other variables, such as dominant bacterial phyla, dominant bacterial families, and color response data, based on substrate utilization. Substrate utilization assay data were analyzed after substrates were divided into six groups and the average absorbance per category was calculated [32]. All meaningful loadings (>0.5) were included and considered significant in the interpretation of principal components [33]. Biplots were constructed to interpret the analysis, with the original variables drawn as vectors to summarize the correlation between the variable and both illustrated axes [34]. Weighted UniFrac distance metric [35] was used on beta diversity analysis, to compare bacterial communities presenting in each treatment through principal coordinates analysis (PCoA) and non-metric multidimensional scaling (NMDS) with QIIME software. Linear discriminant analysis effect size (LEfSe) was determined to find biomarkers between treatments using relative abundances, and the graph which has bars representing the effect size (LDA) for particular taxa was generated. For analysis of similarities (ANOSIM), an ANOVA-like nonparametric statistic test which operates on a ranked dissimilarity matrix was performed using R packages.

Soil Chemical Properties
The effects of GM treatment on soil chemical properties were evaluated during the nine-month study period. In comparison with the BG control, soil pH was significantly (p < 0.05) higher in GM treatment after six months of soybean intercropping ( Table 1). The contents of available K, Ca, and Mg were 95%, 45%, and 103% higher in GM treatment than in the BG control. GM treatment was found to decrease EC and available Fe to a significant level after nine months. Organic matter content was slightly higher in GM treatment (1.40-2.43%) relative to the BG control (0.87-1.92%) during all the months that were considered. A higher amount of available P (256.65-331.97 mg kg −1 ) and low amount of Mn (6.83-9.42 mg kg −1 ) and Zn (1.21-2.59 mg kg −1 ) were observed in both treatments during the nine-month study period.
indicate the absolute amount of each characteristic. Different letters in the same row within the same sampling period indicate significant differences (p < 0.05) in one-way ANOVA.

Functional Diversity and Metabolic Potential of Bacterial Communities
The carbon source utilization of bacterial communities was determined in several incubation time, and the largest difference appeared at 96 h between treatments was demonstrated. During all the months that were considered, the substrate richness and H-index were significantly (p < 0.05) higher in GM treatment than in the BG control (Table 2). GM treatment also increased the metabolic potential of bacterial communities, as revealed by the remarkably higher average well color development values (0.796-1.077) comparing with that (0.231-0.58) in BG treatment. Higher utilization of all six carbon source categories was observed in GM treatment than in the BG control, regardless of the month used for comparison (Figure 2a). PCA showed distinct differences in carbon source utilization patterns among bacterial communities in the GM treatment and BG control samples ( Figure  2b). The proportion of variation explained by PC1 was 87.27%, while, by PC2, it was 4.59% in samples collected during nine months. Samples collected from GM treatment and BG control were widely separated from each other on the PC1 axis. Utilization of all six carbon source categories showed a significantly (p < 0.05) positive correlation with GM treatment while negatively correlated with the BG control (Supplementary Materials Table S1). Table 2. Functional diversity of bacterial communities in samples collected from the bare ground control and green manure treatment during the study period.

Bare Ground Control
Green Manure Treatment

Population Diversity and Structure of Bacterial Communities
Considering all the sampling periods, the Chao1 value and Shannon diversity index revealed in GM treatment ranged from 201 to 574 and 5.516 to 7.158, respectively (Table 3). Both these values were significantly (p < 0.05) higher than that (165-277 for Chao1 and 5.125-5.422 for Shannon diversity index) observed in the BG control. After three-month soybean intercropping, the InvSimpson diversity index was also significantly higher in GM treatment (45.450) relative to the BG control (19.357). GM treatment tended to increase the Chao1 values, Shannon, and InvSimpson diversity indices as the cultivation periods increased. Shifts in the bacterial community structure in relation to GM treatment was visualized through PCoA (Figure 3a). The proportion of variation explained by PC1 was 63.22%, while, by PC2, it was 9.41% in samples collected during nine months. The distribution of samples from the BG control was closer to each other, which was far from samples collected from GM treatment. Ordination by NMDS showed that overall variation in the bacterial community structure among samples was associated with GM treatment and the sampling periods (Figure 3b). GM treatment samples collected from five sampling periods plotted further from each other, whereas samples from the BG control were closer to each other. ANOSIM also revealed a significant (p < 0.005) difference between bacterial communities in the BG control and GM treatment samples.  Dominant bacterial phyla which had higher relative abundance (>2%) in either BG control or GM treatment samples were selected and used for comparison. A total of nine phyla dominating in the BG control or GM treatment were displayed, which accounted for 93.32-95.35% or 94.92-96.18% of the relative abundance in the whole bacterial community, respectively (Figure 4a). During all the months that were considered, the relative abundances of Acidobacteria, Chloroflexi, and the candidate phylum WPS-2 revealed in the BG control were 34.83-38.26%, 10.83-17.62%, and 1.72-4.39%, respectively, which were higher than in GM treatment (18.44-33.54%, 7.23-13.82%, and 1.36-2.86%). GM treatment was demonstrated to enrich Bacteroidetes, Gemmatimonadetes, Planctomycetes, Proteobacteria, and Verrucomicrobia regardless of the month used for comparison.
PCA showed distinct differences in bacterial phyla dominating in the BG control and GM treatment (Figure 4b). The proportion of variation explained by PC1 was 48.86%, while, by PC2, it was 18.46% in samples collected during nine months. In general, the distribution of samples from the BG control was closer to that from one-month and two-month GM treatments. Furthermore, they were widely separated from three-, six-, and nine-month GM treatment samples on the PC1 axis. In most of the cases, samples collected from the BG control showed a significantly (p < 0.05) positive correlation with Acidobacteria, Chloroflexi, and the candidate phylum WPS-2, while they significantly negatively correlated with Bacteroidetes, Gemmatimonadetes, Proteobacteria, and Verrucomicrobia (Supplementary  Materials Table S2). Samples collected from three-and six-month GM treatments were widely separated from nine-month GM treatment on the PC2 axis, and they showed significantly positive correlation with Actinobacteria and Planctomycetes. Besides, it was observed that nine-month GM treatment samples positively correlated with Bacteroidetes, Proteobacteria, and Verrucomicrobia. Dominant bacterial families which had higher relative abundance (>2%) in either BG control or GM treatment samples were compared. Considering all the sampling periods, the relative abundances of Burkholderiaceae, Chitinophagaceae, Microscillaceae, Sphingobacteriaceae, and Sphingomonadaceae were higher in GM treatment (2.12-8.45%, 0.41-4.1%, 0-2.85%, 0.20-2.12%, 1.52-5%) relative to the BG control (0.09-0.39%, 0.06-0.41%, 0-0.09%, 0-0.05%, 0.31-1.27%) (Figure 5a). The higher relative abundance of Acetobacteraceae, Acidothermaceae and Ktedonobacteraceae was recorded in the BG control (1.28-3.16%, 2.01-4.64%, 4.97-8.02%) than in GM treatment (0.8-2.99%, 1.26-4.03%, 2.25-4.63%) regardless of the month used for comparison. In the PCA biplot, distribution of samples collected from one-month GM treatment was closer to the BG control samples, while they were widely separated from other GM treatment samples on the PC1 axis (Figure 5b). The BG control samples showed a significantly (p < 0.05) positive correlation with Acetobacteraceae, Acidobacteriaceae, Acidothermaceae, Ktedonobacteraceae, and Solibacteraceae (Supplementary Materials Table S2). In addition to one-month GM treatment samples, most of the samples collected from GM treatment showed a significantly positive correlation with Burkholderiaceae, Chitinophagaceae, Microscillaceae, Sphingobacteriaceae, and Sphingomonadaceae. LEfSe was used to find biomarkers in the BG control and GM treatment, and the LDA scores for particular taxa were given. The most effective biomarkers (LDA score > 3.6) identified in GM treatment included the phylum Bacteroidetes, Gemmatimonadetes, Planctomycetes, and Proteobacteria; and the family Burkholderiaceae, Caulobacteraceae, Chitinophagaceae, Gemmatimonadaceae, Microscillaceae, Rhizobiaceae, Sphingobacteriaceae, Sphingomonadaceae, and Xanthomonadaceae (Figure 6a). The most effective biomarkers (LDA score < −3.6) found in the BG control included the phylum Acidobacteria and Chloroflexi, the candidate phylum WPS-2, and the family Ktedonobacteraceae (Figure 6b).

Relationships between Dominant Bacterial Taxa and Soil Properties or Metabolic Potential of Bacterial Communities
Among all the measured soil properties, pH has a profound influence on the relative abundance of dominant bacterial taxa, as revealed by CCA and Pearson correlation analysis (Figure 7 and Supplementary Materials Table S3). Acidobacteria, Chloroflexi, Gemmataceae, Ktedonobacteraceae, and Solibacteraceae showed a significantly (p < 0.05) negative correlation with soil pH, while Bacteroidetes, Burkholderiaceae, Gemmatimonadetes, Microscillaceae, Proteobacteria, Sphingobacteriaceae, Sphingomonadaceae, Verrucomicrobia, and Xanthobacteraceae positively correlated. Available K, Ca, and Mg showed the similar trends as pH did on these bacterial taxa. Organic matter contents significantly positively correlated with Burkholderiaceae, Gemmatimonadetes, Pedosphaeraceae, Sphingobacteriaceae, Sphingomonadaceae, and Verrucomicrobia. Both Planctomycetes and its affiliated family Gemmataceae showed a significantly negative correlation with total N, as well as available Cu and Zn, and the latter also negatively correlated with organic matter content, as well as available Ca and Mg. Available P and Fe significantly positively correlated with Acetobacteraceae, Chloroflexi, and Ktedonobacteraceae, while they showed a negative correlation with Burkholderiaceae, Chitinophagaceae, Proteobacteria, and Sphingomonadaceae. The utilization of all six carbon source categories showed a significantly (p < 0.05) positive correlation with the Bacteroidetes, Gemmatimonadetes, Planctomycetes, Proteobacteria, and Verrucomicrobia phyla, as well as the Burkholderiaceae, Chitinophagaceae, Microscillaceae, Sphingobacteriaceae, Sphingomonadaceae, and Xanthobacteraceae families (Supplementary Materials Table S4). The relative abundances of Acidobacteria, Chloroflexi, Ktedonobacteraceae, and the candidate phylum WPS-2 significantly negatively correlated with the metabolic potential of bacterial communities toward most of the carbon source categories.

Discussion
In this study, the potential of soybean intercropping to modify soil chemical properties and bacterial community, which further restored the deteriorated soil, was demonstrated. Numerous studies have clearly shown that both the plant genotype and the soil type drive the rhizosphere microbial community [36][37][38]. The exudates, such as sugars and organic acids derived from soybean root [39,40], may serve as alternative carbon sources, which changed bacterial community within the inter-row of the tea plantation. The soil pH is considered as the best predictor of community composition [41][42][43] may also play a crucial role affecting bacterial community in this study. In addition to soil pH, organic carbon and available K were found as the key characteristics significantly correlated with the variation of the soil fungal community in a tea plantation [17]. We also demonstrated that available K, Ca, and Mg significantly positively correlated with the most effective biomarkers identified in GM treatment. In this tea plantation, the ecological benefit brought by GM treatment was addressed, since the increase in biodiversity facilitated the stability of microbial communities against biotic and abiotic environmental perturbations [44]. The soybean intercropping was assumed to enrich nitrogen-fixing taxa in the tea plantation. As expected, the Rhizobiaceae family was identified as one of the most effective biomarkers in GM treatment, which accounted for higher relative abundance (0.08-1.54%) than in the BG control (0-0.02%). Within Rhizobiaceae, the relative abundances of two rhizobial taxa identified at the genus level (Allorhizobium_Neorhizobium_Pararhizobium_Rhizobium and Mesorhizobium) also increased in response to GM treatment (Supplementary Materials Figure S2). Besides, soybean intercropping was demonstrated to increase the relative abundance of Bradyrhizobium, a genus which is commonly found in soybean nodules [45]. However, no symbiotic nitrogen fixation occurred in this strongly acidic soil during the study period. This may be overcome after conducting agricultural manipulations, such as liming and inoculating effective rhizobia, which have been used to increase nodulation and nitrogen fixation of soybeans grown in strongly acidic soils [46].
Smit et al. [47] hypothesized that the ratio between Proteobacteria and Acidobacteria (P/A) may provide insight into the general nutrient status of soils. Low P/A ratios would be indicative of oligotrophic soils, while high ratios were displayed under copiotrophic conditions. In this study, the average P/A ratio in the BG control and GM treatment was 0.74 and 1.36, respectively, which suggested an increase in the nutritional level of soils after soybean intercropping. Based on the experimental and meta-analytical results, Fierer et al. [48] further proposed that bacteria belonging to the phylum Acidobacteria were categorized as oligotrophs, while Bacteroidetes and β-Proteobacteria exhibited copiotrophic attributes. The life strategy concepts regarding copiotrophy-oligotrophy continuum [49] and the competitor-stress tolerator-ruderal (C-S-R) framework [50] were considered in this study to explain the prevalence and distribution of bacterial communities affected by GM treatment. Samples collected from the BG control showed a significantly (P < 0.05) positive correlation with Acidobacteria, a phylum associated with oligotrophy or stress tolerator, while they were negatively correlated with Bacteroidetes and Proteobacteria, which were associated with copiotrophy or competitor [51]. Based on the concept of three potential life history strategies further interpreted by Ramin and Allison [52], we demonstrated that bacterial taxa dominating in the tea plantation shifted from maintenance strategists (Acidobacteria) to growth strategists (Proteobacteria) in response to soybean intercropping. In addition, bacterial phyla Actinobacteria and Bacteroidetes, which were referred to as maintenance and resource acquisition strategists, also increased after GM treatment.
During all the months that were considered, higher relative abundances of Acidobacteria, Chloroflexi, and the candidate phylum WPS-2 were observed in the BG control than in GM treatment. Acidobacteria and Chloroflexi were abundant among very slow-growing and mini-colony-forming soil bacteria, and both of them have been found to increase continuously with increasing age of tea plantations [14,53]. The relative abundance of Acidobacteria was shown to increase toward lower soil pH [54][55][56]. Members belonging to this phylum can metabolize nutrient-poor and recalcitrant carbon substrates, and the abundance was generally higher in soils with very low resource availability [48,57]. Chloroflexi was demonstrated to account for 33-54% of the bacterial sequences in low-pH soils in a study of chemoautotrophic carbon dioxide fixation in drained paddy soils [58]. Members of this phylum are likely oligotrophic taxa adapted to survive under the resource-limited conditions found in deeper horizons, and they have been shown to decrease upon soil fertilization [57,59]. The candidate phylum WPS-2 (Eremiobacteraeota) was observed as a dominant bacterial group in oligotrophic Antarctic soils, which might proliferate through chemosynthetic carbon fixation [60]. A metagenomic study also revealed that members belonging to this phylum are likely anoxygenic phototrophs capable of carbon fixation [61]. The strongly acidic and low nutritional status of this tea plantation soils led to the higher relative abundances of these oligotrophic taxa as well as the lower metabolic potential of bacterial communities toward six carbon source categories.
In comparison with the BG control, notable increases in relative abundances of phyla Bacteroidetes, Gemmatimonadetes, and Proteobacteria were observed in GM treatment. Both Bacteroidetes and Proteobacteria were considered as copiotrophic taxa that exhibited relatively rapid growth rates in elevated C, N, or S conditions [48,57,62]. The life strategy of Gemmatimonadetes has not been clarified in soils, with only a few studies dealing with their response to re-wetting and thawing events [63,64].
In this study, GM treatment was demonstrated to enrich Bacteroidetes members, such as Chitinophagaceae, Microscillaceae, and Sphingobacteriaceae families. Chitinophagaceae have been found to positively correlate with β-glucosidase activity and organic matter degradation indicators [65,66]. Members belonging to this family were inhabitants of rhizosphere or roots. They possess the abilities to solubilize tricalcium phosphate; produce indole-3-acetic acid, 1-aminocyclopropane, and 1-carboxylate deaminase; and stimulate plant growth [67,68]. The Microscillaceae family is less explored, and currently only two chemoorganotrophic genera showing gliding motility and gelatin-degrading ability were reported [69]. Members belonging to Sphingobacteriaceae possessed very divergent genome fractions responsible for secondary metabolism, dormancy and sporulation, and metabolism of aromatic compounds under different environmental pressures [70]. In this study, Mucilaginibacter was found as the most dominant genus affiliated with Sphingobacteriaceae, whose relative abundance accounted for 31-100% of Sphingobacteriaceae in GM treatment samples collected from five sampling periods. This genus was demonstrated to degrade a variety of plant residues, such as cellulose, hemicellulose, pectin, and xylan through its highly diverse carbohydrate-active enzymes [71,72].
Within Proteobacteria, members belonging to families Burkholderiaceae and Sphingomonadaceae were identified as the most effective biomarkers in GM treatment. The genus Burkholderia affiliated with Burkholderiaceae was proposed as the most diverse and environmentally adaptable plant-associated bacteria [73]. They display genomic and metabolic plasticity, showing a broad range of host and possess a variety of plant growth promoting traits [74,75]. Our previous study also indicated that Burkholderiaceae dominated not only the root surface but root interior, which might form intimate associations with its host tea plants [15]. The family Sphingomonadaceae has been intensively studied, owing to its extraordinary catabolic flexibility to degrade natural and xenobiotic compounds [76]. The abilities of this family to induce innate plant growth promoting traits have been investigated, which include the formation of phytohormones, siderophores, and chelators [77]. In this study, the genus Sphingomonas affiliated with Sphingomonadaceae accounted for higher average relative abundance in GM treatment (1.83%) than in the BG control (0.5%). Members belonging to this genus have been noted to improve plant growth during stress conditions, such as drought, salinity, and heavy metals in agricultural soil [78].
Ho et al. [51] suggested that tracking microorganisms at finer phylogenetic and taxonomic resolution (e.g. family level or lower) may be more effective to capture changes in community response. At the family level, we demonstrated that GM treatment decreased the relative abundances of oligotrophic taxa, such as Acidobacteriaceae and Solibacteraceae (belonging to the phylum Acidobacteria), and Ktedonobacteraceae (belonging to the phylum Chloroflexi). In contrast, the relative abundances of families, such as Burkholderiaceae, Chitinophagaceae, Microscillaceae, Sphingobacteriaceae, and Sphingomonadaceae, increased in response to soybean intercropping. The potential copiotrophic attributes were displayed in these families, since they showed significantly positive correlation with the metabolic potential toward all six carbon source categories. Many bacterial families, such as Acidothermaceae, Ktedonobacteraceae, Microscillaceae, Pedosphaeraceae, and Solibacteraceae, have not been extensively studied owing to their rare cultured representatives. They might be obtained using diluted media or through a long period of incubation instead of conventional cultivation methods [48,79]. Efforts should be made to obtain more representatives to clarify their eco-physiological roles, especially in nutrient cycling and plant growth promotion in a near future.

Conclusions
In this study, the temporal distribution of bacterial communities in response to soybean intercropping in strongly acidic soils was revealed. GM treatment was observed to boost both functional and population diversity, which helps to enhance the stability of bacterial communities within the tea-cultivated ecosystem. The inter-row habitats stimulated by soybean intercropping further serve as hot spots for copiotrophs-mediated nutrient cycling and microbe-plant interactions. Besides, the dominant bacterial taxa provide opportunities to study their eco-physiological roles in the tea plantation.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/su132011478/s1, Figure S1: Monthly average temperature and monthly precipitation of the experimental site during June 2020 and March 2021. (CWB Observation Data Inquire System, https: //e-service.cwb.gov.tw/HistoryDataQuery/ (accessed on 20 September 2021)), Figure S2: Hierarchically clustered heat-map of relative abundances of rhizobial taxa in the bare ground (BG) control and green manure (GM) treatment samples collected after one month (_1M), two months (_2M), three months (_3M), six months (_6M), and nine months (_9M) of treatments. Heat-map was generated based on Euclidean dissimilarity metrics, normalized from −1 to 2.5 and color scale represented from blue to red through yellow, Table S1: Principal component loadings (correlation between original loading variable and principal component) as a measure of influence of a loading variable on treatment differences, Table S2: Principal component loadings (correlation between original loading variable and principal component) as a measure of influence of a loading variable on treatment differences, Table S3: Pearson correlation between dominant bacterial taxa and soil properties, Table S4: Pearson correlation between dominant bacterial taxa and metabolic potential of bacterial community toward six carbon source categories.

Data Availability Statement:
The raw sequence reads obtained in this study have been submitted to the NCBI SRA database under BioProject PRJNA716948.