The Application of Mixed Organic and Inorganic Fertilizers Drives Soil Nutrient and Bacterial Community Changes in Teak Plantations

Appropriate fertilization can enhance forest productivity by maintaining soil fertility and improving the structure of the bacterial community. However, there is still uncertainty surrounding the effects of combined application of organic and inorganic fertilizers on soil nutrient status and bacterial community structure. A fertilization experiment was set up in an eight-year-old teak plantation with five treatments involved: mixed organic and NPK compound fertilizers (OCF), mixed organic and phosphorus fertilizers (OPF), mixed organic, NPK and phosphorus fertilizers (OCPF), mixed NPK and phosphorus fertilizers (CPF) and no fertilization (CK). Soil chemical properties and bacterial communities were investigated, and the co-occurrence pattern of the bacterial community under different fertilization treatments was compared. The results showed that the contents of soil organic matter and nitrate nitrogen, and the soil pH values were the highest after OCPF treatment, which were 20.39%, 90.91% and 8.16% higher than CK, respectively. The richness and diversity of bacteria underwent no obvious changes, but the structure of the soil’s bacterial community was significantly altered by fertilization. Of the dominant bacteria taxa, the relative abundance increased for Gemmatimonadetes, Myxococcota, ADurb.Bin063-13 and Candidatus_Koribacter, and decreased for Chloroflexi, Proteobacteria, JG30-KF-AS9 and Acidothermus under OCPF treatment in comparison to CK. The number of nodes and edges, the average degree and the network density of bacterial community co-occurrence networks were the greatest in OCPF treatment, indicating that application of OCPF could make the network structure of soil bacteria more stable and complex. Moreover, soil pH and organic matter were significantly correlated with bacterial community structure and were considered the main influencing factors. These findings highlight that the combined application of organic, NPK and phosphorus fertilizers is highly beneficial for improving soil quality and optimizing bacterial community structure in teak plantations.


Introduction
Fertilization is one of the pivotal management measures in the timber forest cultivation process. It directly affects soil fertility and forest productivity, as well as the diversity and composition of microorganisms by changing soil aggregate structure and nutrient levels [1,2]. Applying inorganic fertilizers is regarded as the fastest and most direct way to effectively increase soil nutrient content and accelerate tree growth and timber output. However, it can also generate negative impacts on soil health, such as soil acidification and degradation [3,4]. Likewise, the soil microbial community is also highly sensitive to the addition of inorganic fertilizers. Observable consequences seem to be the decreasing population and biodiversity of microorganisms, as well as instability of the community structure [5][6][7]. In recent years, the application of organic fertilizers for soil property has been greatly emphasized [8,9]. Some studies have shown that the application of organic fertilizers could significantly enhance soil quality in terms of soil structure, physicochemical properties and biological characteristics [10,11]. The combined application of organic and inorganic fertilizers could mediate organic carbon (C) sequestration, as well as elevate the nitrogen (N) and phosphorus (P) balance in soil environments [12]. Moreover, the addition of organic fertilizers is the key to relieving soil acidification [13,14].
Soil microorganisms are an extremely abundant biotic group and are important for maintaining the diversity and stability of the terrestrial ecosystem [15,16]. Bacteria play a decisive role in forest soil environments by participating in biogeochemical cycles, such as litter and dead root decomposition and nitrogen and phosphorus mineralization [17]. As the major biological regulator of nutrient cycling [18], bacteria groups are also very sensitive to changes in the soil environment, especially with the increase of available nutrients due to external inputs. However, it is widely acknowledged that mineral fertilizers generally decrease the diversity of soil bacterial communities [7,19] and may potentially decimate the ecosystem services provided by microorganisms [20,21]. In contrast, organic fertilizers directly provide labile carbon and a large number of energy substances, stimulating the activity of soil microorganisms [22]. The soil after the combined application of organic and inorganic fertilizers exhibited more abundant bacteria and a higher level of biodiversity [14]. The survival strategies of microbes also shifted to adapt to the changes in the soil environment after organic fertilization; in particular, the abundance of beneficial bacteria associated with the cycling of carbon, nitrogen and phosphorus in soil increased substantially [8,23]. For example, species of Firmicutes and Proteobacteria bacterial phyla prefer nutrient-rich environments and are enriched in the soil after organic fertilization [24]. Ye et al. also found that the addition of biofertilizers enriched the Actinobacteria phylum that was closely related to the decomposition of complex organic compounds [25]. However, some scholars considered that the diversity or richness had no notable changes after organic fertilizers [8,26], and the relative abundance of Firmicutes and Proteobacteria declined [27]. The response of microbes to fertilizers varies according to a series of factors and can be highly related to soil heterogeneity. Hence, more research is needed to understand the underlying changes in soil microorganisms after fertilization.
Bacterial communities are characterized not only by the composition and number of species, but also by the interactions and associations among different taxa [28]. The bacterial co-occurrence network was performed to reveal the interrelationships between species and the complexity of community structure and function [29][30][31]. Previous studies confirmed that alterations of community structure and bacterial network were closely related to soil properties, such as pH [32], soil matter organic [33,34] and available phosphorus [35], which have been considered major drivers of variation in biotopes. Although symbiotic network patterns were universally applied in forest ecosystems, the response of complex bacterial communities to the application of mixed organic and inorganic fertilizers in teak plantations is still unknown.
Teak (Tectona grandis L.f.) is one of the most precious hardwood species in the tropics, known for its excellent timber quality, and it is used extensively in the manufacture of high-grade furniture and wooden floors [36]. Fertilization is widely conducted in the practice of teak plantation forestry. Studies have shown that the application of chemical fertilizer can promote nitrogen and phosphorus absorption in trees and accelerate stem growth [37][38][39]. The addition of organic fertilizer can also promote the diameter growth of teak and improve the soil quality of its plantations [40]. There was little focus on the response of soil nutrients and bacterial communities to the combined application of organic−inorganic fertilizers in teak plantations. Thus, the objectives in the present study are to analyze the differences in soil chemical characteristics and bacterial community composition between different fertilization treatments, and to identify the key factors influencing the structure of bacterial communities.

Study Site
The study site is located in Luodian county, Qiannan Prefecture, Guizhou Province, China (106 • 46 E, 25 • 25 N, elevation 413 m to 449 m). This region has a subtropical monsoon humid climate with a mean annual temperature of 19.6 • C. The average annual rainfall is about 1400 mm. The daylight hours and annual frost-free period are 1350-1520 h and 335 days, respectively. The soil was identified as yellow−red earth and derived from sand shale. Teak seedlings of the fine clone 7544 were used for afforestation, which was received from the Research Institute of Tropical Forestry, Chinese Academy of Forestry in Guangzhou, China. The experimental plantation was planted in 2010 at a spacing of 2.5 m × 3.0 m.

Experimental Design
The fertilization experiment was established in a randomized block design and included five treatments: the combined application of organic and NPK compound fertilizers (OCF); the organic and CaMgP fertilizers (OPF); the mixed application of organic, NPK and CaMgP fertilizers (OCPF); the combined application of NPK and CaMgP fertilizers (CPF); and no fertilization as a control (CK). Each treatment was repeated four times. A total of 20 treatment plots of 270 m 2 each were set up. Fertilization was conducted in May 2018 and 2019, and annual fertilizer inputs are listed in Table S1.
The commercial organic fertilizer includes 45% organic matter content and more than 5% NPK content. NPK compound fertilizer is composed of N, P 2 O 5 and K 2 SO 4 in the ratio 14:16:15. CaMgP fertilizer contains 18% of P 2 O 5 , 45% of CaO and 12% of MgO. All inorganic and organic fertilizers were obtained from Hubei Xinyangfeng Fertilizer Co., Ltd., Jingmen, China.

Soil Sampling and Chemical Analysis
Soil samples of the surface layer (0-10 cm) were collected two years and four months after the final fertilization (in September 2021). Five sampling points in each plot were randomly selected, and soil samples taken from the five points were mixed together evenly. First, the plant roots, residues and break stones in the soil samples were removed. The samples were then passed through a 2 mm mesh sieve and divided into two sections: one was stored at −80 • C in a refrigerator for DNA extraction, and the other was used for the determination of soil chemical properties.
The pH value of each soil sample was determined with an acidimeter (PHS-3C, Shanghai, China) in a solution that was mixed with double-distilled water and soil (2.5:1). The content of soil organic matter (SOM) was measured using the potassium dichromate oxidation method [41]. Ammonium nitrogen (NH 4 + -N) was assayed using a spectrophotometric method, with Nassi reagent acting as the extraction agent. Nitrate (NO 3 -N) in the soil was extracted with 2M KCl and measured by ultraviolet spectrophotometry. Available phosphorus (AP) content was extracted from soil samples with a mixed solution of ammonium fluoride and hydrochloric acid and determined using the visible spectrophotometer [42,43].

DNA Extraction and PCR Amplification
DNA extractions from 0.5 g fresh soil samples were performed using the Magnetic soil DNA kit (Tiangen Biotech Co., Ltd., Beijing, China), referring to the manufacturer's instructions. The extracted DNA concentration and purity were monitored on 1% agarose gels, and the concentration was diluted to 1 ng·µL −1 with sterile water to ensure amplification efficiency and accuracy. The V3-V4 hypervariable regions of bacterial 16S rRNA genes were amplified by polymerase chain reaction (PCR) with the forward primer (341F 5 -CCTAYGGGRBGCASCAG-3 ) and the reverse primer (806R 5 -GGACTACNNGGGTATCTAAT-3 ). The PCR amplification mixture contained 15 µL Phusion ® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 µM forward primer, 0.2 µM reverse primer and 10 ng template DNA. The thermal cycling program of the PCR reaction process included an initial denaturation step at 98 • C for 1 min, followed by 30 cycles of denaturation at 98 • C for 10 s, annealing at 50 • C for 30 s, extension at 72 • C for 30 s, and then a final 5-min extension at 72 • C. PCR reactions were subjected to three replicates [28,44].

Illumina NovaSeq Sequencing and Bioinformatics Analysis
PCR products were mixed with an equal volume of 1 X loading buffer and detected by electrophoresis on 2% agarose gel, which was then purified using a Qiagen Gel Extraction Kit (Qiagen Inc., Frank furt, Germany). NEBNext ® Ultra™ IIDNA Library Prep Kit (Cat No. E7645) was used to sequence the construction of libraries following the manufacturer's protocols, and the library quality was evaluated on the Qubit ® 2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). Finally, the library was paired-end sequenced on an Illumina Novaseq platform (Illumina, San Diego, CA, USA) according to the standard protocols of Novogene Co., Ltd. (Beijing, China). The raw data have been stored in the Sequence Read Archive (SRA) data of the National Center for Biotechnology Information (NCBI) database with accession number PRJNA820355.

Statistical Analysis
One-way analysis of variance (ANOVA) and Tukey's honestly significant difference (HSD) test in SPSS25.0 software (SPSS Inc., Chicago, USA) was used to determine the differences in soil chemical properties and alpha-diversity indexes that were normalized by using a log 10 (X + 1) conversion. Significant differences in the relative abundance of dominant bacteria at phylum and genus levels were analyzed by the Kruskal−Wallis test. Principal Coordinate Analysis (PCoA) of weighted unifrac distances was used to explore the beta-diversity changes between bacterial communities for each sample. Analysis of Similarities (ANOSIM, based on Weighted Unifrac distances) by Anosim function within QIIME2 software was used to test community structure discrepancies among different treatments. Linear discriminant analysis (LDA) and effect size (LEfSe) analysis were conducted to identify potential bacterial biomarkers of each treatment. Mantel tests (vegan package in R software, http://www.r-project.org/ 27 March 2022) were used to test the relationship between the soil environment and the distribution of bacterial communities (at the ASV level). Redundancy analysis (Canoco 5.0 software, http://www.microcomputerpower.com/ 27 March 2022) and Spearman's correlation analysis were applied to evaluate the influence of edaphic factors on the main bacterial phyla and genera.
Co-occurrence networks were used to explore the interactions between soil bacterial taxa (based on ASV levels). The network construction was based on Pearson's correlation analysis by using the psych package in Rstudio software (http://www.rstudio.com/ 27 March 2022) and the correlation coefficient and significance level (R > 0.06, p < 0.05) was set up. Analytical data were the ASVs with the top 100 absolute abundance in each sample [48,49]. Visualization of network relationships and calculation of topological properties were done using Gephi0.9.2 software (http//:gephi.org/ 27 March 2022).

Changes in Soil Chemical Properties in Teak Plantations
The combined fertilization enhanced the soil pH value and SOM content (Table 1, p < 0.05). The soil pH value in OCPF treatment was significantly higher than CK by 8.16%, and the SOM content also increased by 20.39% and 17.53% compared to CK and CPF, respectively. The average content of NH 4 + -N and NO 3 − -N in OCF, OCPF and CPF treatments were noticeably greater than those in CK. The content of NH 4 + -N and NO 3 − -N in OCF treatment were 152.94% and 90.59% greater than CK, respectively. The soil AP content showed an increasing trend with increased fertilization in the order: OCF > CPF > OPF > OCPF > CK.

Sequence Characteristics and Soil Bacterial Community Diversity
A total of 2,215,833 reads were obtained after sequencing all soil samples, and after denoising and filtering, 1,266,203 effective sequences were used for subsequent bacterial community analysis. The coverage index (99.95-100%) showed that the sequencing depth effectively reflected the actual richness of soil bacterial communities (Table S1). The average quantity of effective sequences decreased after fertilization when compared with no fertilization. The mean length of effective sequences under CPF was the longest, and was significantly higher than that of CK, OCF and OPF treatments (Table S2). A total of 10,048 ASVs were retained in 20 samples. The Venn diagrams showed that 837 bacterial ASVs were shared among different treatments and 993 ASVs were shared between four fertilization treatments. The number of specific ASVs in CPF and OCPF treatments was greater than in CK ( Figure S1). Based on the Shannon, Simpson and Chao1 indexes, no significant differences in soil bacterial diversity and richness were found between different fertilization treatments ( Figure 1A-C). However, the Pielou's evenness index in CPF and OCPF treatments was significantly higher than in OCF and OPF treatments ( Figure 1D).

Impact of Fertilization on Bacterial Community Structure
The community structure of soil bacteria, based on the ANOSIM analysis, varied significantly among different fertilization treatments (R = 0.372, p = 0.005). The PCoA analysis explained 60.94% of the total variation in bacterial communities, of which PC1 accounted for 42.32%. However, soil bacterial communities were not clearly clustered into different groups according to the different fertilization treatments in the ordination diagram ( Figure S2).

Impact of Fertilization on Bacterial Community Structure
The community structure of soil bacteria, based on the ANOSIM analysis, varied significantly among different fertilization treatments (R = 0.372, p = 0.005). The PCoA analysis explained 60.94% of the total variation in bacterial communities, of which PC1 accounted for 42.32%. However, soil bacterial communities were not clearly clustered into different groups according to the different fertilization treatments in the ordination diagram ( Figure S2).
A total of 36 bacteria phyla and 619 genera with definite groups were detected by species annotation of all ASVs. The mean relative abundance of 9 phyla and 11 genera exceeded 1% ( Figure 2B). At the level of phylum, the relative abundance of Acidobacteria (32.79-45.64%) was the highest, followed by Proteobacteria (18.24-23.27), Firmicutes (6.28-10.29%), Actinobacteria (6.97-8.68%) and Chloroflexi (4.33-10.10%). Other bacterial phyla occupied only a minor proportion ( Figure 2, Table S4). The relative abundance significantly declined for Chloroflexi, while it increased for Gemmatimonadetes and Myxococcota under OCPF and CPF treatments, compared to CK (Figure 3, p < 0.05). Verrucomicrobia and Bacteroidetes (except for OCF) showed higher relative abundances in fertilization treatments than those in the CK group. The relative abundance of Firmicutes and Actinobacteria was increased in OCP, OCPF and CPF treatments, whereas Acidobacteria and Proteobacteria in OCPF and CPF treatments were less abundant than that of CK (Table S4). At the genus level, Subgroup_2, Candidatus_Solibacter and ADurb.Bin063-1 were the three most dominant bacterial genera ( Figure 2B). The relative abundance of ADurb.Bin063-1, Acidothermus, Candidatus_Koribacter and The relative abundance significantly declined for Chloroflexi, while it increased for Gemmatimonadetes and Myxococcota under OCPF and CPF treatments, compared to CK (Figure 3, p < 0.05). Verrucomicrobia and Bacteroidetes (except for OCF) showed higher relative abundances in fertilization treatments than those in the CK group. The relative abundance of Firmicutes and Actinobacteria was increased in OCP, OCPF and CPF treat-ments, whereas Acidobacteria and Proteobacteria in OCPF and CPF treatments were less abundant than that of CK (Table S4). At the genus level, Subgroup_2, Candidatus_Solibacter and ADurb.Bin063-1 were the three most dominant bacterial genera ( Figure 2B). The relative abundance of ADurb.Bin063-1, Acidothermus, Candidatus_Koribacter and JG30-KF-AS9 presented significant differences among the fertilization treatments. Candidatus_Koribacter showed a prominent increase in OCF and OPF but a decrease for Acidothermus in comparison to CK. Application of OCPF markedly reduced the relative abundance of Acidothermus and JG30-KF-AS9 compared to CK (Figure 3, Table S4).

Biomarker Taxa of Soil Bacterial Communities
A total of 218 biomarkers of 17 phyla were identified in all treatments (Figures 4 and  S3). There were 28 taxa in OCPF, with the major enriched taxa to OCPF being Gemmatimonadetes and Nitrospirae (Figure 4). 18 taxa were significantly enriched in OPF treatment, and Acidobacteriales (order) had the highest relative abundance. The number of biomarkers enriched in OCF treatment was minimal, and the main significant groups were Crenarchaeota and Acidobacteria. There were more biomarkers in CK and CPF than in other treatments, and Chloroflexi (Ktedonobacterales) in CK and Verrucomicrobia and Bacteroidetes in CPF were more abundant ( Figure S3).

Biomarker Taxa of Soil Bacterial Communities
A total of 218 biomarkers of 17 phyla were identified in all treatments (Figures 4  and S3). There were 28 taxa in OCPF, with the major enriched taxa to OCPF being Gemmatimonadetes and Nitrospirae (Figure 4). 18 taxa were significantly enriched in OPF treatment, and Acidobacteriales (order) had the highest relative abundance. The number of biomarkers enriched in OCF treatment was minimal, and the main significant groups were Crenarchaeota and Acidobacteria. There were more biomarkers in CK and CPF than in other treatments, and Chloroflexi (Ktedonobacterales) in CK and Verrucomicrobia and Bacteroidetes in CPF were more abundant ( Figure S3). Species with no significant differences are colored uniformly yellow, and differential species biomarkers follow the group for coloring.

Co-Occurrence Network of Bacterial Communities
Co-occurrence networks were built based on Pearson's correlation analysis of soil bacterial communities. The topological attributes showed that fertilization changed the interactions among bacterial taxa ( Figure 5; Table 2). The number of edges in OCF, OPF, OCPF, CPF and CK was 285, 431, 500, 426 and 273, respectively. The proportion of positive edges tended to decrease and negative edges increased after fertilization ( Table 2). The network density and average degree also increased and reached the maximum in OCPF treatment, indicating that the application of OCPF increased the complexity of the bacterial community. In contrast, the number of network edges in OCF and CK was lower than in the others, showing an obvious modularization and a longer average path length. The hubs (higher connect nodes and edges) in OCPF were ASV204, ASV106 and ASV79, which belonged to Acidobacteria, Proteobacteria and Gemmatimonadetes, respectively. Acidobacteria and Proteobacteria in CPF and CK, Proteobacteria and Firmicutes in OCF, Acidobacteria, Proteobacteria and Gemmatimonadetes in OPF were the central taxa in the corresponding soil ( Figure 5). . LEfSe analysis of soil bacterial biomarkers after different fertilization treatments. The circles radiating from the innermost to the outermost represent the taxonomic levels from phylum to species in the cladogram. The size of the circles is proportional to the size of the relative abundance. Species with no significant differences are colored uniformly yellow, and differential species biomarkers follow the group for coloring.

Co-Occurrence Network of Bacterial Communities
Co-occurrence networks were built based on Pearson's correlation analysis of soil bacterial communities. The topological attributes showed that fertilization changed the interactions among bacterial taxa ( Figure 5; Table 2). The number of edges in OCF, OPF, OCPF, CPF and CK was 285, 431, 500, 426 and 273, respectively. The proportion of positive edges tended to decrease and negative edges increased after fertilization ( Table 2). The network density and average degree also increased and reached the maximum in OCPF treatment, indicating that the application of OCPF increased the complexity of the bacterial community. In contrast, the number of network edges in OCF and CK was lower than in the others, showing an obvious modularization and a longer average path length. The hubs (higher connect nodes and edges) in OCPF were ASV204, ASV106 and ASV79, which belonged to Acidobacteria, Proteobacteria and Gemmatimonadetes, respectively. Acidobacteria and Proteobacteria in CPF and CK, Proteobacteria and Firmicutes in OCF, Acidobacteria, Proteobacteria and Gemmatimonadetes in OPF were the central taxa in the corresponding soil ( Figure 5).

Relationships between Soil Factors and Bacterial Community Structure
The Mantel tests suggested that soil pH and SOM have a strong effect on the structure of the bacterial community (Table S5, p < 0.05). Redundancy analysis between the relative abundance of bacteria and soil factors demonstrated that the soil chemical properties explained 21.36% and 40.23% of the total variation on phylum and genus levels, respectively ( Figure 6A,B). The contribution of soil pH to bacterial community composition was the greatest and the most significant (Table S5).

Relationships between Soil Factors and Bacterial Community Structure
The Mantel tests suggested that soil pH and SOM have a strong effect on the structure of the bacterial community (Table S5, p < 0.05). Redundancy analysis between the relative abundance of bacteria and soil factors demonstrated that the soil chemical properties explained 21.36% and 40.23% of the total variation on phylum and genus levels, respectively ( Figure 6A,B). The contribution of soil pH to bacterial community composition was the greatest and the most significant (Table S5). Spearman's correlation analysis of major bacterial phyla and soil properties showed that the relative abundance of Chloroflexi was significantly negatively related to soil pH Spearman's correlation analysis of major bacterial phyla and soil properties showed that the relative abundance of Chloroflexi was significantly negatively related to soil pH value, NH 4 + -N and NO 3 − -N contents. The abundance of Nitrospirae was positively correlated with soil pH value, SOM, NH 4 + -N and NO 3 − -N contents. Most phyla such as Verrucomicrobia, Gemmatimonadetes and Nitrospirae were positively linked with soil pH value. The abundance of Gemmatimonadetes and soil AP content, Patescibacteria and SOM content were significantly positive correlations ( Figure 6A). For the genera ( Figure 6B), Subgroup_2 and Acidothermus, which belonged to Acidobacteria and Actinobacteria, respectively, and their abundances were negatively related to soil pH value. HSB_OF53-F07 and Acidibacter showed strong negative correlations with most edaphic factors; they belong to Chloroflexi and Proteobacteria, respectively. The inverse relationship existed between the soil pH value and the abundances of Subgroup_2 (Acidobacteria) and Acidothermus (Actinobacteria).

Variations in Soil Chemical Properties under Different Fertilization Treatments
It is widely known teak grows better in rich soils with a neutral and slightly alkaline pH [50,51]. Teak is a species that shows a strong preference for calcium and is sensitive to soil acidity in its growth process [39]. In this study, we found that the soil pH value increased after the combined application of organic and inorganic fertilizers, and soil pH was significantly higher in OCPF treatment than no fertilization treatment (p < 0.05, Table 1). The main reason could be attributed to the high calcium oxide content in the soil due to the addition of CaMgP fertilizer [52].
The increase in the soil pH value might also be related to the decomposition of soil organic matter [53]. A significant improvement in soil nutrient status characterized by organic matter, nitrogen and phosphorus enrichment was observed as a result of the mixed fertilization ( Table 1). The application of organic−inorganic fertilizers could regulate soil aggregates and provide a source enriched in carbon and nitrogen for soil microorganisms, which promotes their activity and growth [22,54]. Correspondingly, bacterial taxa can mediate soil carbon cycling and effectively release more nutrients from organic matter [55,56]. The cumulation of soil nutrients was directly or indirectly influenced by bacterial communities, such as Gemmatimonadetes and Nitrospirae taxa, which are closely related to soil carbon degradation and nitrogen mineralization [57,58]. In this study, the increase of bacterial abundance in OCPF ( Figure S3) might promote the level of soil carbon and nitrogen. Gemmatimonadetes taxa were identified as one of the biomarkers and keystone taxa in OCPF treatment ( Figures 5 and 6), which is directly related to the highest level of organic matter in this soil.
Soil pH may play a major role in transforming soil nutrient contents by changing bacterial taxa. Microbial communities exhibit poor growth in lower pH environments, and a relatively acidic condition usually inhibits both enzyme activity and overall cellular metabolism [59,60]. A previous study also found that the application of organic fertilizer would strengthen the activity of phosphatase and invertase [61]. Moreover, bacterial networks and interactions among taxa exhibited higher complexity after OCPF treatment, which had a positive impact on nutrient cycling and accumulation [62].

Response of Bacterial Community to Soil Chemical Properties
Fertilization had a positive impact on soil pH and nutrient levels, directly or indirectly influencing changes in bacterial community composition and abundance. In the present study, soil pH and organic matter were the main impact factors after fertilization in the teak plantation (Table S5), which appears consistent with the findings of Wang et al. [63] and Zhang et al. [64], who researched soil bacterial communities in response to organic fertilizer application. However, Zhao et al. [65] reported that soil nutrients (e.g., organic matter, total N and total P), rather than pH, showed a significant correlation with the majority of abundant taxa. This may be caused by different fertilization strategies. Our study found that the abundance of Chloroflexi and WPS-2 decreased with increasing soil pH value ( Figure 6). This result suggested that acidic conditions are more suitable for the development of both types of bacteria. On the contrary, most bacterial phyla are obviously enriched as the soil pH value increases. For example, Verrucomicrobia, Latescibacteria, Nitrospirae and Gemmatimonadetes prefer neutral over acidic environments [66]. It may be possible that the degeneration or loss of metabolic function of these phyla in acidic soils occurs when the pH value is below a certain threshold [67]. Our study also found that the relative abundance of Gemmatimonadetes was positively associated with soil available phosphorus, which was also corroborated by the result that the eco-function of Gemmatimonadetes taxa was positively correlated with soil nutrients [57].
The application of organic fertilizer plays a substantial role in sustaining soil organic matter levels and simultaneously, the input of organic materials is also beneficial for the improvement of soil aggregate structure and available nutrients [68]. These changes further affect the distribution, abundance, activity and composition of the microbial community [2]. A clear negative relationship was observed between the relative abundance of Acidibacter (which belongs to Proteobacteria) and SOM content, which was inconsistent with the phenomenon that Proteobacteria taxa have a fast growth rate and are more likely to grow in nutrient-rich conditions [69]. This is potentially because of the site heterogeneity exhibited by bacterial taxa in soil, along with their substrate preferences, and the exact reasons for this inconsistency need further investigation.

Effects of Fertilization on Soil Bacterial Community Composition and Structure
Bacteria are the most abundant and diverse group of soil microorganisms, maintaining the stability and sustainability of the soil ecosystem in forests [70,71]. The decline in bacterial diversity may directly affect the stability of soil microbial communities [31]. Fertilization leads to changes in the structure of soil bacterial communities by changing the living environment of microorganisms. Previous studies found that urea or NPK fertilizer application inhibited the growth of microorganisms and showed a significant negative effect on bacterial diversity [72], while the addition of organic manure [73] or biological organic fertilizer [74] shaped a richer variety of soil bacteria. It was also reported that the diversity and richness of the bacterial community did not change significantly with either the application of organic fertilizer or mineral fertilizer [64,75]. In this study, the number of effective sequences showed a slight decline after fertilization (Table S3), and the alpha diversity of bacterial communities did not differ significantly between CK and fertilization treatments ( Figure 1). Possible reasons for this result are the short duration of the fertilizer's application or soil heterogeneity between treatments.
Determining the variations in soil bacterial composition between fertilization versus no fertilization treatment can effectively assess the effect of soil quality, health and fertilization on artificial forest ecosystems [76]. In the present study, although fertilization had no effect on the alpha diversity of bacteria (Figure 1), significant differences were revealed in bacterial community structure among the different treatments (p < 0.01, Figure S2). Acidobacteria and Proteobacteria were two of the most abundant phyla observed in the bacterial community, as shown in previous studies [77,78]. Both of them can properly enhance soil fertility and sustainability when they were enriched in soil [60]. It was reported that Acidobacteria were functionally diverse and contributed to the degradation of root exudates and litter to stabilize organic matter content in soil [79]. The activity of Acidobacteria taxa (such as Subgroup_2 and Candidatus_Solibacter) was more influenced by OCPF because the effect of soil pH may be stronger than that of organic matter, resulting in a decrease in the abundance of Acidobacteria (Table S4, Figure 6A). Firmicutes were the third most dominant bacterial phylum in the tested soils ( Figure 1, Table S4), and were considered to be the copiotrophic taxa in previous reports [63,80,81]. The addition of effective nutrients through OCPF and CPF treatments might create a more favorable environment for Firmicute taxa multiplication, thus expanding the abundance of Firmicutes. Additionally, Gemmatimonadetes, Myxococcota and Actinobacteria under OCPF treatment were significantly enriched, which may be attributed to their ecological function related to the degradation capacity of organic matter [82,83]. Contrarily, the application of OCPF, OCF and CPF had significant side effects on the relative abundance of Chloroflexi, and the growth and development of this phylum taxa are likely to be limited by the accumulation of soil nutrients. This is similar to the results of Liu et al. [76] and Durand et al. [84], who indicated that Chloroflexi is slow-growing taxa that is classified as an oligotroph.

Network Patterns of Soil Bacterial Community
Co-occurrence network analysis provides evidence of the direct and indirect cooperation or composition relationships among microbial taxa [85]. The interaction among bacterial taxa is closely correlated with their ecological functions comprising nutrient metabolism associated with the C, N and P cycles [86,87]. The linkage density among species based on the network contributed to predicting ecosystem stability and complexity. This is because the network complexity is often dependent upon the number of edges and nodes within the network [28]. More network interactions and connectors (taxa) are fundamental to stabilize the microbial community structure [88].
Our study found that fertilization altered the relationship between bacterial species, and soil bacterial communities in different treatments exhibited different patterns of cooccurrence networks ( Figure 5). The number of linkage edges and network degrees was increased in all fertilization treatments. The number of nodes and edges and the network density after OCPF were greater than in CK and other fertilization treatments ( Table 2), indicating that interspecies relationships and community structures of bacteria were more complicated in the soil. The stronger interactions of bacterial taxa in OCPF treatment compared to CK possibly accelerated the metabolic activity of microorganisms and promoted the accumulation of soil nutrients. The higher content of soil organic matter and nitrate−nitrogen in OCPF also supported this result (Table 1). Yu et al. [49] also indicated that bacterial taxa within the complex network structure play a prominent role in promoting soil nutrient cycling in teak plantations. Moreover, the addition of fertilizers decreased the proportion of positive edges compared to no fertilization ( Figure 5, Table 2), indirectly suggesting that soil effective nutrients aggravated competition and niche separation of bacterial species [89,90]. The lower nutrient environment allowed a large number of bacterial taxa to coexist, while high nutrient environments caused more negative interactions among species [91]. These opposing relationships inhibited each other and excluded more species in the community composition, resulting in a diminution of species diversity [18,28]. Network diameter, modularity coefficient and average path length based on OCPF treatment were relatively smaller than others ( Figure 5). The result demonstrated that an intensive world network formed after OCPF treatment to rapidly respond to alterations in the soil microenvironment [92].

Conclusions
In this study, mixed fertilization with organic, NPK and phosphorus fertilizers significantly increased the content of soil organic matter and available nutrients and reduced soil acidity in young teak plantations. The changes in organic matter level and pH value in soil further altered the composition and network of the bacterial community. The relative abundance of Gemmatimonadetes and Myxococcota in OCPF treatment was clearly enriched, while it was decreased for Chloroflexi, Acidobacteria and Proteobacteria. The topological properties of the co-occurrence network showed that the nodes, edges, network density and average degree were the highest in OCPF, indicating that the bacterial community structure in this soil was more stable and complex. Moreover, soil pH and organic matter were closely associated with bacterial species and were regarded as the primary factors shaping bacterial community structure. Further research on the impact of keystone species within microorganisms is required to identify their contribution to nutrient cycling and to predict the coupling relationship between their function and forest productivity in plantation ecosystems.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10050958/s1, Table S1: Fertilizer input for each treatment; Table S2: Raw reads, effective tags, mean length and goods coverage of soil samples after sequencing; Table S3: The values (mean, n = 4) of raw reads, effective tags, mean length and goods coverage among the different groups; Table S4: The relative abundance (n = 4) of soil bacterial communities at dominant phylum and genus levels under different fertilization treatments; Table S5: Mantel tests of the correlation between bacterial abundance on ASV levels and soil properties, and the contributions of soil chemical properties to variation in soil bacterial taxa; Figure S1: The similarities and differences in the number of ASVs among the different treatments; Figure S2: Principal coordinate analysis (PCoA) and ANOSIM based on weighted unifrac of a bacterial community on ASV level for 20 soil samples from different fertilization treatments; Figure S3: The linear discriminant analysis (LDA) effect size (LEfSe) of soil bacterial biomarkers under different fertilization treatments. The values were significant (p < 0.05) when the LDA score was more than 2.

Data Availability Statement:
The data presented in this study are available in the NCBI Sequence Read Achieve database (accession number: PRJNA820355).