Changes in the Microbial Community in Soybean Plots Treated with Biochar and Poultry Litter

: The application of organic materials that promote beneﬁcial microbial activity is vital to maintaining soil health and crop productivity. We investigated the effect on the soil microbiome of applying biochar (BC), poultry litter (PL), and a combination of biochar and poultry litter (BC/PL) in soybean cultivation at the Red River Research Station (Bossier City, LA, USA). We characterized the microbial proﬁles, community structure, and co-occurrence network from sequencing data to infer microbial interactions in the soil samples collected in the ﬁrst and second years of each soil treatment (2016 and 2017, respectively). Our results showed that soil treatments with BC, PL, and a combination of both moderately changed the microbial community composition and structure. In particular, genera signiﬁcantly affected by the different soil treatments were identiﬁed via differential abundance analysis. In addition, canonical correspondence analysis revealed that soil chemical properties, total N in the ﬁrst year, and total C and pH in the second year inﬂuenced the community variability. The differentially enriched bacterial ASVs and co-occurring taxa were linked to nutrient cycling. This study provides insights into the impact of soil carbon amendment on the soil microbiome, a process which favors beneﬁcial bacteria and promotes soybean growth.


Introduction
Intensive farming practices, including the usage of pesticides and synthetic fertilizers, hamper sustainable agriculture methods which provide a stable supply of food and feed without significant environmental damage [1]. The soybean (Glycine max. L.) is the most valuable legume crop, ranking as the ninth most-produced crop in the world [2]. As a leguminous plant, the soybean can acquire nitrogen through the N-fixation by symbiotic and asymbiotic bacteria in nodules and root surfaces [3][4][5].
Soil microbial communities are vital in the ecosystem, playing key roles in processes such as nutrient and soil carbon cycles [6][7][8]. Likewise, modifications in the soil microbial community through management practices such as the addition of soil amendments subsequently impact their assemblage and composition in the plant rhizospheres [6,[8][9][10]. In this regard, improving knowledge on the capacity of the soil microbiome is crucial for sustainable crop management, and will aid in understanding the influence of the soil

Soil Physicochemical Analyses
Soil pH was measured using a pH meter at a 1:1 soil:water ratio [36]. The total C and N contents in the soil were measured by dry combustion using a Thermo Finnigan FLASH EA 1112 CN analyzer (CE Elantech, Inc., Lakewo1od, NJ, USA) at the water quality laboratory at the Red River Research Station, LSU AgCenter [36]. Mehlich-3 extractable nutrients (P, K, Ca, Mg, and S) were determined with induced coupled plasma atomic emission spectroscopy (ICP-AES) (Spectro CIROSCCD, Mahwah, NJ, USA) from the Soil Testing Laboratory, LSU AgCenter [36]. Bulk density was measured by collecting multiple cores of known volume from different irrigation schemes [37].

Biochar and Poultry Litter Analysis
The chemical characteristics of the biochar and the poultry litter used in this study (Table 1) were analyzed with the following protocols. The pH of biochar samples was measured using a pH meter at a 1:5 solid:water ratio. The ash content of the biochar was determined by overnight combustion at 750 • C in a muffle furnace. The total C and total N were analyzed with the same method as soil sample analysis. The macro-and micro-element contents in biochar and poultry litter were analyzed from the Soil Testing Laboratory at the LSU AgCenter.

DNA Sample Processing
Total genomic DNA was extracted from 0.25 g of soil using the PowerSoil ® DNA Isolation Kit (Qiagen, Germantown, MD, USA) following the manufacturer's instructions. The extracted DNA samples were analyzed for quality and quantity using NanoDrop 1000 (Thermo Scientific, Wilmington, DE, USA) and were verified with agarose gel electrophoresis. The DNA samples were sent to Genomics Research Laboratory of the Biocomplexity Institute of Virginia Tech (Blacksburg, VA, USA) for 2 × 250 bp paired-end sequencing following the 16S Illumina Amplicon Protocol from the Earth Microbiome Project (https://earthmicrobiome.org/protocols-and-standards/16s/ accessed on 1 March 2017). The standard library preparation protocol as earlier described was followed [38]. Briefly, the V4 region of the 16S small subunit rRNA gene was amplified using the universal primer set, 515F (GTGYCAGCMGCCGCGGTAA) and 806R (GGACTACNVGGGTWTCTAAT). The PCR reaction mixture included 13.0 µL of PCR-grade water, 1.0 µL of template DNA, 10 µM of each primer, and 10.0 µL of 5PRIME HotMasterMix (2×) (Quantabio, Beverly, MA, USA). Samples were amplified in duplicate under the following thermocycler conditions: 94 • C for 3 min for initial denaturing, then 35 cycles of 94 • C for 45 s, 50 • C for 60 s, and 72 • C for 90 s. A final elongation step occurred at 72 • C for 10 min followed by a hold at 4 • C. After pooling duplicates, all amplicons were visualized on a 2% agarose gel and quantitated on a Qubit fluorometer (FisherScientific, Hampton, NH, USA). Normalization was performed based on Qubit results and amplicons were pooled. The pool was quantitated by qubit and tapestation assays to determine size (380 bp). The pool was then quantitated using qPCR. Amplicons were loaded at 9.5 pM along with a 25% Phix spike and sequenced using the MiSeq v2 500-cycle kit on the Miseq platform (Illumina, Inc., San Diego, CA, USA). Quality scores were within Illumina specifications. The Illumina MiSeq sequencing generated raw reads of 6,193,760 and 5,470,755 in 2016 and 2017, respectively (Table 2). Individual samples ranged from 250,643 to 478,815 sequences per sample (Table 2).

Microbial Community Analysis
The microbial communities were analyzed with the Quantitative Insights into Microbial Ecology (QIIME 2 2019.4) pipeline [39]. Demultiplexed paired-end fastq files and a metadata mapping file were used as input files. Individual samples ranged from 250,643 to 478,815 sequences per sample, with a sequence length of up to 250 bp (paired-end sequences, 2 × 250 bp), which were the V4 region of the 16S rRNA genes of the prokaryotic community in each sample ( Table 2). The Divisive Amplicon Denoising Algorithm (DADA2) in QIIME2 was used for sequence correction, removal of chimeras, and denoising. The sequences were trimmed using the -p-trim-left-f/r and -p-trunc-len-f/r function. The forward reads of all the sequences were truncated to 245 bases and the reverse reads at 128 bases and 122 bases for 2016 and 2017, respectively, based on the quality scores and the minimum length identified during subsampling from all the sequences. Paired sequences reads were joined and quality filtered using the paired-end DADA2 pipeline [40]. After quality filtering, we finally obtained 2,572,488 and 3,315,244 unique sequences for the samples from 2016 and 2017, respectively. Taxonomy was assigned to Amplicon Sequence Variants (ASVs) using the q2-feature-classifier [39] against the SILVA 128 database and BLAST+ consensus taxonomy classifier [41] directly after quality filtering and building of feature table and feature data. The unique sequences were aligned to a total of 51 phyla, 254 classes, 576 orders, 929 families, and 1763 genera (Table S1). All ASVs were aligned with MAFFT [42] and used to construct a phylogeny with fasttree2 [43].
In comparing microbial communities in 2016 and 2017, the resulting feature tables and representative sequences in each year were merged using feature-table merge and featuretable merge-sequence QIIME2 plugins for diversity analyses. The 2 years of sequencing data were merged for a direct comparison. Although differences in observed sample composition could arise due to different times of sequencing runs even with the same amplification protocol, the differences attributed to technical variation are expected to be relatively small compared to the biological variation, as previously reported [44,45].
The "q2-diversity" plugin was used to calculate the alpha diversity metrics and beta diversity metrics after samples were rarefied at an even sampling depth. The rarefaction curves showed a nice saturation pattern, indicating that the ASVs detected from this DNA sequencing represent most of the microbial community in each sample ( Figure S1). Diversity and community richness were analyzed with Kruskal-Wallis non-parametric tests to compare species richness among treatments in q2 diversity metrics using observed features "observed OTUs" [46], Faith's phylogenetic diversity [47], and Shannon's diversity indices [48]. The variation in community composition was calculated using weighted and unweighted UniFrac metric [49]. Differential abundance of features was analyzed using DESEq2 package in R [50]. Nonmetric multidimensional scaling (NMDS) was conducted with the "vegan" R package using Bray-Curtis dissimilarity, unweighted Unifrac, and weighted Unifrac values among samples, which were subjected to a Wisconsin double transformation before analysis [51]. Statistical significance corresponding to differences in dispersion among groups was calculated using "permdisp" homogeneity in group dispersions and differences in means and variance using permutational analysis of variance (PERMANOVA) with 999 iterations. Other statistical analyses were conducted with JMP Pro Statistics, version 15.1.0 (SAS Institute, Cary, NC, USA).
The microbiomeSeq R package using the Adonis function was used to calculate the overall significance level and visualization of soil environmental variables in canonical correlation plots. BIOENV analysis in QIIME2 was conducted to select the best possible subsets of soil chemical variables that correlate with the community pattern using Spearman's correlation value [52].
Datasets of QIIME2 output, which include metadata of the microbiome dataset with taxonomic classifications of sequences, were further visualized using Explicit software package [53], microbiomeSeq [54] and microbiomeAnalyst [55].

Microbial Network Analysis
The co-occurrence network was inferred based on the Spearman's rank correlation. The ASVs with a sum relative abundance of at least 20% of all samples were subjected to a pairwise correlation [56]. A correlation between 2 ASVs was considered statistically significant at values of >0.6 for the Spearman's correlation coefficient (r s ) and at <0.01 for the p-value [57,58]. Multiple testing correction using Benjamini-Hochberg standard false discovery rate correction were used to adjust the p-values to reduce the chances of obtaining false-positive results [59]. The co-occurrence network is comprised of nodes and edges which represent the genus and all significant correlations identified from pairwise comparison of genera abundance. Topological properties and statistical analyses were calculated with R using the vegan [60], igraph [61], and Hmisc [62] packages. Network visualization was performed using the Gephi platform (http://gephi.github.io/ accessed on 4 July 2019), [63]). All samples were divided into groups by soil amendments.

Soybean Growth and Yield
In the first year, no significant growth difference was observed among treatments (p = 0.94) (Table 3). However, in the second year, the PL and BC/PL treatments showed higher growth than BC-treated and the untreated control plots (p = 0.03) ( Table 3). Nevertheless, there was no significant difference in yield among the treatments in both years, although the BC/PL had the highest yield among treatments over the 2 years.

Microbial Community Profiles and Structure
In terms of alpha diversity, there was no significant difference among the treatments in either of the years with respect to the ASV features species richness and Shannon's diversity based on Kruskal-Wallis non-parametric tests (Table 4 and Figure S2). However, the no-treatment (control) soil from the second year showed a significant difference as compared to the first year in terms of community richness based on phylogenetic lineage (p < 0.04) (Table 4B), whereas the BC-and BC/PL-treated soils in the second year were significantly different with regard to the abundance of observed features as compared the untreated soil in the first year (Table 4A). The parallel analyses conducted on species richness based on abundance and phylogenetic diversity (Faith's pd) showed measurement deviations but consistently supported that the three kinds of soil treatments (BC, PL, and BC/PL) in the second year were significantly different from the no-treatment control in the first year in terms of Shannon's diversity (p < 0.05) (Table 4C). Further analyses to determine beta diversity were performed using metrics including beta group significance (unweighted Unifrac and weighted Unifrac) and the Adonis test ( Table 5). All of the metrics indicated that there was no difference among treatments in the 2 years, but there was a significant difference in each year's overall data between the 2 years (Table 5). After data filtering of uninformative ASV features with low counts and low variance (below 20% for the prevalence filter with a minimum count of 4 and below 10% variance based on the inter-quantile range), the beta diversity patterns of the bacterial communities at the genus level were depicted in Bray-Curtis dissimilarity heatmaps and non-metric multidimensional scaling (NMDS) ordination using weighted Unifrac distance ( Figure 1). Other distance methods using Bray-Curtis, unweighted Unifrac, and weighted Unifrac are presented altogether in the supplemental figure ( Figure S3). We observed the tendency that the BC-treated soil showed the highest level of similarity with regard to the untreated soil in both years ( Figure 1a). The relationship between the bacterial community and the treatments was not significant for the taxonomic-based distance measure using Bray-Curtis for both years. However, the weighted Unifrac distance, a phylogenetic-based distance measure which accounts for the relative abundance of taxa, showed a significant difference among treatments based on permutational multivariate analysis in the first year (PERMANOVA F-value: 2.0225; R-squared: 0.33582; p < 0.02) and a slightly substantial difference in the second year (PERMANOVA F-value: 1.4692; R-squared: 0.26863; p < 0.087). The beta dispersion using "permdisp" homogeneity in group dispersions revealed no significant difference in the spread in any treatments in the 2 years (PERMDISP F-value: 0.64839; p = 0.598 and PERMDISP F-value: 0.76461; p = 0.535, first year and second year, respectively). Thus, the difference in the microbial community composition is not ascribed to variances within the treatment but among treatments, especially in the first year.
The microbial community structures were well represented in the NMDS plots, as shown by the NMDS stress values of 0.08 and 0.11 in the first year and second year, respectively. In the first year, 33.5% variation in the microbial communities was attributed to the different treatments. There were diverging separate clusters between the BC-, PL-, and BC/PL-treated soils as compared to the untreated control in the first year ( Figure 1b). In the second year, 26.8% of the variation among microbial abundances was ascribed to the different treatments, although with slight to nil significant differences (p = 0.087). The NMDS showed the ordination space of microbial communities in the BC-and PLtreated clusters within the ordination distance of BC/PL-treated soil, suggesting that all the microbial communities in the BC-and PL-treatments were present in the BC/PL treatments. The BC-and PL-treatments showed some portions of un-overlap ordination, suggesting that both communities showed unique microbial compositions. On the other hand, the untreated soil converged across all the treatments, which may indicate that the microbial members existing in the untreated soils could also be present in all treated soils. Thus, the results suggest that the composition of microbial community was similar, but the abundance of microbial community members varied under the different treatments. The significant PERMANOVA of the weighted Unifrac distances supported the differences in community structure among treatments because of the variation in the relative abundance of species/taxa among treatments. The correlation of microbial composition with the environmental factors could be observed by other statistical analyses (below section). difference in the second year (PERMANOVA F-value: 1.4692; R-squared: 0.26863; p < 0.087). The beta dispersion using "permdisp" homogeneity in group dispersions revealed no significant difference in the spread in any treatments in the 2 years (PERMDISP Fvalue: 0.64839; p = 0.598 and PERMDISP F-value: 0.76461; p = 0.535, first year and second year, respectively). Thus, the difference in the microbial community composition is not ascribed to variances within the treatment but among treatments, especially in the first year.

Canonical Correspondence Analysis
Treatments influencing the soil chemical properties can change the composition of associated soil microbiomes. Although there was no significant difference in pH, total N, total C, and the C/N ratio among the treatments with one-way ANOVA (Table 6), we attempted to analyze the correlation of the environmental factors with the microbial community through multivariate statistical techniques and Bioenv and Adonis tests with the QIIME2 ( Table 7). The canonical correspondence analysis (CCA) plots showed the correlation of individual soil chemical properties with the structure of the microbial community (Figure 2), and the Bioenv analysis table showed the best set of environmental factors that described the community variation in the 2 years (Table 7). Total nitrogen (R-squared: 0.124, p = 0.039, Spearman's r s = 0.741) and the combination of total nitrogen and pH (r s = 0.559) best explained the variability of microbial composition in the first year ( Table 7). The BC/PL-treated soil had the highest total nitrogen values in the first year, followed by the BC-and untreated soils ( Table 7). The pH (R-squared: 0.197, p = 0.002, r s = 0.358), the combination of pH and total carbon (r s = 0.422), and the combination of pH, total nitrogen, and total carbon (r s = 0.363, p = 0.031) influenced the microbial compositions in the second year ( Figure 2 and Table 7). Total carbon was higher in the BC-and the BC/PL-treated soils throughout the 2 years, especially in the BC-treated soil in the second year (Table 6).     In our findings, the application of soil amendments in soybean fields showed 12% to 19% variation in the microbial communities affected by total N and pH in the first and second year, respectively ( Figure 2). However, total N was strongly associated with the microbial community composition in the first year (r s = 0.74), while the combination of pH and total C was moderately correlated (r s = 0.42) in the second year.

Bacterial Community Composition and Differential Abundance
Proteobacteria, Acidobacteria, Bacteroidetes, Verrucomicrobia, and Gemmatimonadetes were the five dominant phyla across treatments over the 2 years (Table 8 and Figure 3). Phyla that comprise putative beneficial bacteria, including Nitrospirae, Actinobacteria, Planctomycetes, Chloroflexi, and Firmicutes were also found to be major ones (Figure 3a). There was no statistical difference among treatments at the phylum level in both years, except for Actinobacteria (Kruskal-Wallis test: χ 2 = 9.55, p = 0.022) and Chloroflexi (Kruskal-Wallis test: χ 2 = 9.22, p = 0.026), which were more abundant in the untreated control in the first year of this study (Table 8 and Figure 3a). Interestingly, however, the phylum Bacteroidetes increased in the second year compared to the first year in all soil samples, including the untreated control (Figure 3a). The most abundant classes were Sphingobacteria in the phylum Bacteroidetes; Blastocatellia and subgroup 6 in the phylum Acidobacteria; and classes in the phylum Proteobacteria including Alphaproteobacteria, Betaproteobacteria, Deltaproteobacteria, and Gammaproteobacteria (Table 8 and Figure 3b). Other major classes were Gemmatimonadetes, Nitrospira, and the Soil Crenarchaeotic Group (SCG). There was no significant difference observed among the classes in the first year. However, in the second year, Gammaproteobacteria was more abundant in the BC/PL-treated soil than the untreated soil (p = 0.049) (Figure 3b).
The top 20 orders with a relative abundance of >1% were analyzed using the Kruskal-Wallis test (Table 8 and Figure 3c). In the first year, Myxococcales and SCG uncultured bacterium were both abundant in the untreated and BC-treated soil (p < 0.05) (Figure 3c). During the second year, Xanthomonadales (p = 0.008) and Myxococcales (p = 0.043) were significantly more abundant in the BC/PL-and PL-treated soils (Figure 3c). At the family level, unassigned taxa (p = 0.030) were significantly higher in the untreated soil in the first year, while the Chitinophagaceae (p = 0.037) was significantly higher in the BC-treated soil in the second year (Figure 3d). There were significant differences in the microbial composition from 2016 to 2017, which may be attributed to the proliferation of bacterial taxa favored by the addition of soil amendments, and thus a shift in the microbial community over time.
At the genus level, we identified features strongly influenced by different soil treatments, which were enriched or depleted genera based on relative abundance as compared to the untreated soil. There were 166, 214, and 178 enriched genera (primarily the identifiable genera from the phyla Proteobacteria, Bacteroidetes, Actinobacteria, Acidobacteria, and Planctomycetes (Table 9)), and 154, 234, and 197 depleted genera (primarily from the phyla Proteobacteria, Actinobacteria, Bacteroidetes, Chloroflexi and Firmicutes (Table 10)    In the first year, the most significant enriched genera in the BC-treated soils were Nitrospira_uncultured, Acinetobacter, and Rivibacter (Table 9), while depleted genera were Cellumonas, Erwinia, Clostridium, Flavobacterium, and uncultured members of Proteobacteria, Chloroflexi, and Thaumarcheota (Table 10). The PL-treated soil was also enriched with Acinetobacter, Sphingobacterium, Rhodocista, and other uncultured members of the phyla, as well as Proteobacteria, Verrucomicrobia, and Thaumarcheota (Table 9). Meanwhile, Bradyrhizobium, Roseiflexus, Dactylosporangium, Caldithrix, and Cyanobacteria were depleted in the PL-treated soil ( Table 10). The enriched genera in BC/PL-treated soil included the Chthoniobacterales_DA101 soil group, Flavobacteriales_uncultured, Sphingobacterium, Duganella, Prosthecobacter, and other uncultured members of Verrucomicrobia and Bacteroidetes, while Caldithrix, Asanoa, Roseiflexus, Tumebacillus, Bradyrhizobium, Virgisporangium, and uncultured members of Proteobacteria and Nitrospirae were identified as depleted genera (Tables 9 and 10). There were much fewer enriched and depleted genera in the BC-treatment soil as compared to the PL-and BC/PL-treated soils. In the second year, the most significant enriched genera in the treated soils were Pelomonas, Paenarthtobacter, Emticicia, Chitinimonas, Nannocystis, Cellulosimicrobium, Luteimonas, Stigmatella, Lysinimonas, and several members of Proteobacteria (Table 9), while the main depleted genera were Clostridium, Pedobacter, Nitrospira, Sphingomonas, Pseudogulbenkiania, Dechloromonas, Paucimonas, Aquincola, Dechlorobacter, and several members of Proteobacteria and Gemmatimonadetes (Table 10). Notably, enriched genera in the treated soils increased in the second year regardless of the treatment, in which the PL-treated soil contained the highest number of enriched genera (Table 9 and Figure 4), suggesting that the 2 years of treatment caused the enrichment of certain bacterial groups in the treated soil.

Microbial Network Analysis
The co-occurrence network of the microbial community in each treatment and year was constructed to assess the impact of soil treatments on the structure of the microbial community ( Figure 5 and Figure S4). The entire networks across all treatments in the 2 years were composed of 15,328 significant associations (edges) with 373 nodes, with an average clustering coefficient of 0.68 and an average path length of 2.81 (detailed data not shown). Network edges were comprised mainly of strong positive associations. The key taxa that regulate the network structure were identified based on the node degree (>50), closeness centrality (>0.44), and betweenness centrality scores (<0.18) [64]. These taxa include Verrucomicrobia, Gemmatimonas, Bacteroidetes (env_OPS1_uncultured bacterium), and Acidobacteria (Blastocatellaceae_11-24 and subgroup 4).
In the first year, based on top 35 most abundant taxa, the PL-treated and BC-treated soils had the highest co-occurring taxa belonging to Proteobacteria (including Nitrosomon-adaceae_uncultured bacterium, Sphingomonas, Desulfurellaceae_H16, Rhodospirillales_uncultured, Steroidobacter, and Ramlibacter). Meanwhile, the untreated control and BC/PL-treated soils were predominantly composed of co-occurring taxa belonging to Acidobacteria (including uncultured Acidobacteria bacterium, Candidatus Solibacter, and Blastocatellaceae uncultured bacteria). However, the BC/PL-treated soil had the greatest number of co-occurring taxa under the phylum Gemmatimonadetes (Gemmatimonas, Longimicrobia, and uncultured Planctomycetes). In the second year, there was an increase in co-occurring taxa in the phylum Bacteroidetes across treatments, particularly in the BC-treated soil (Sphingobacteria, Chitinophaga, Cytophaga, and other uncultured Bacteroidetes). The BC/PL-treated soil had the greatest number of co-occurring taxa belonging to phylum Proteobacteria (Sphingomonas, Nitrosomonadaceae uncultured, Steroidobacter, Rhodopirillales, Masillia, and Desulfurellales). The PL-treated soil had the greatest number of co-occurring taxa within Verrucomicrobia (Verrucomicrobia uncultured, Spartobacteria, Chtoniobacter, and Opitutus) ( Figure 5). The co-occurring taxa within Gemmatimonadetes, Bacteroidetes, Verrucomicrobia, and some members in the Acidobacteria are involved in degrading complex materials. At the same time, co-occurring taxa within Proteobacteria are involved in nutrient cycling in the soil.

Microbial Network Analysis
The co-occurrence network of the microbial community in each treatment and year was constructed to assess the impact of soil treatments on the structure of the microbial community ( Figures 5 and S4). The entire networks across all treatments in the 2 years were composed of 15,328 significant associations (edges) with 373 nodes, with an average clustering coefficient of 0.68 and an average path length of 2.81 (detailed data not shown). Network edges were comprised mainly of strong positive associations. The key taxa that regulate the network structure were identified based on the node degree (>50), closeness centrality (>0.44), and betweenness centrality scores (<0.18) [64]. These taxa include Verrucomicrobia, Gemmatimonas, Bacteroidetes (env_OPS1_uncultured bacterium), and Acidobacteria (Blastocatellaceae_11-24 and subgroup 4).  In the first year, based on top 35 most abundant taxa, the PL-treated and BC-treated soils had the highest co-occurring taxa belonging to Proteobacteria (including Nitrosomon-adaceae_uncultured bacterium, Sphingomonas, Desulfurellaceae_H16, Rhodospirillales_uncultured, Steroidobacter, and Ramlibacter). Meanwhile, the untreated control and BC/PLtreated soils were predominantly composed of co-occurring taxa belonging to Acidobacteria (including uncultured Acidobacteria bacterium, Candidatus Solibacter, and Blastocatellaceae uncultured bacteria). However, the BC/PL-treated soil had the greatest number of

Discussion
Several studies associate soil health and crop productivity with below-ground bacterial diversity [5,15,[24][25][26][27]30,65]. Following this line of thought, we hypothesized changes in microbial communities due to cultivation management, specifically the application of different soil amendments in the soybean field. The study revealed the gradual impact of soil treatments with BC, PL, and the BC/PL on the bacterial communities and the agronomic traits of soybean plants in terms of growth and yield. Specifically, the growth was significantly highest in the PL-treated soil in the second year and was comparable to that of the BC/PL-treated soil. Biochar characteristics depend upon what it is made from and how it was processed. Most biochar has a small labile fraction of carbon but a much larger recalcitrant fraction that can persist in soils [66]. The biochar used in this study was pinewood biochar (similar to that of a study on Pinus radiata biochar), and tends to be recalcitrant and mineralizes in soil slowly [67]. Treatment with BC alone was comparable to no treatment, despite the higher C content it provided, indicating that the addition of biochar alone did not add appreciable amounts of nutrients to the soil. The PL had high total N, P, K, micronutrients, and organic matter (Table S1). Similar to previous studies, the PL-treated soil showed less crop yield in the first and second years of application [17]. Of note, the BC/PL treatment combination showed the highest crop yield for the 2 years, although it was not statistically different from the other treatments.
Alpha diversity indicates no significant difference within treatments. Nevertheless, the species richness and diversity slightly shifted over the 2 years among treatments, suggesting that the microbiomes in the field change over time. Other research groups have shown that biochar causes a slight shift in bacterial community compositions but has no significant short-term impact on microbial diversity in soil environments including maize fields and pasture soils [67][68][69][70]. It is still unknown whether the long-term application of biochar can increase microbial diversity [71][72][73]. However, biochar treatment was also reported to significantly alter soybean field microbial abundance and community composition at a low pH soil condition in one growing season. This was more obvious in the soybean rhizosphere than the bulk soil [74]. Similarly, in the acidic paddy soil, it was reported that the microbial community in the rice rhizosphere responded to biochar treatment more sensitively than bulk soil [13]. These studies suggest that the microbial community in the bulk soil compartment is likely affected gradually by biochar treatment, while it has more immediate effects on the microbial community in the rhizosphere. In our study, soil pH tended to be neutral in both years regardless of treatments (Table 6), although the untreated and poultry litter decreased in pH in the second year.
Regarding the previous studies mentioned above, the neutral soil pH condition could be a reason for the minor variations within the treatments observed in this study. Interest-ingly, bacterial community structure variation was observed based on phylogenetic-based distance using weighted Unifrac distance, which accounts for the relative abundance of taxa among treatments. Hence, some variations in the lineage of bacterial community composition could be attributed to other factors influenced by soil treatments aside from pH.
pH was found to be an important factor in the community variation of microbial taxa in different agro-ecosystems reported in previous studies [10,[75][76][77]. The availability of N and C was also a main driving factor in the shift of microbial community by different treatments of organic materials [78][79][80]. Biochar was reported to have a lesser effect on the microbial community structure in neutral or alkaline soil [77]. In contrast to other studies, distinct microbial structures were observed between the biochar-treated and untreated acidic paddy soils [81][82][83]. In our study, the BC and BC/PL-treated soil maintained a neutral pH, while the PL and untreated soil had a decreased pH in the second year.
With regard to the bacterial community composition, the dominant phyla across treatments were Proteobacteria, Acidobacteria, Bacteroidetes, Verrucomicrobia, and Gemmatimonadetes, consistent with the previous studies on the most dominant phyla in agricultural soils [65,[84][85][86]. These phyla account for 71-81% of total ASVs, implying that these phyla may play an important role in the soybean bulk soil. Cellulolytic bacteria involved in decomposition under orders Myxococcales and Xanthomonadales had a significantly greater presence in the BC/PL-and PL-treated soils in the second year. Members of Myxococcales are considered micropredators in the soil and play a key role in soil C cycling [79]. Interestingly, the differential abundance at the genus level based on adjusted p-values and log-fold change revealed a significant variation between the BC-, PL-, and BC/PL-treated soils vs. untreated soils. The differentially abundant taxa in the treated soils were related to organic matter decomposition and nutrient cycling. Taxa linked to C cycling were Pelomonas, Cellulosimicrobium, Chitinomonas Acinetobacter, Nannocystis, and Stigmatella. Pelomonas species produce enzymes that degrade hemicellulose, carbon substrates, and aromatic compounds [87]. Chitinomonas and Cellulosimicrobium are chitinolytic bacterium [88,89]. The Cellulosimicrobium genome is composed of genes responsible for protection against salinity stresses and production of volatiles, with a number of genes for the degradation of plant biopolymers [90]. Acinetobacter have metabolic capabilities for the degradation of various long-chain dicarboxylic acids and aromatic and hydroxylated aromatic compounds that are associated with plant degradation products [91]. Nannocystis is a common myxobacteria in soil, and is usually from decaying substrates [92]. Stigmatella is another myxobacteria capable of breaking down peptidoglycans, polysaccharides, protein, and other cellular detritus. Moreover, Stigmatella produces antibiotics toxic to yeast and filamentous fungi but not most bacteria [93]. On the other hand, differentially abundant taxa in the treated soils associated with N and other nutrient cycling were Nitrospira_uncultured, Duganella, Leptospirillum, and other indicators of soil fertility such as Pilimelia, Ramlibacter, and Chthoniobacter_DA101 [94]. Nitrospira are chemolithoautotrophic organisms capable of complete nitrification [95]. Duganella can solubilize phosphorus, potassium, and zinc in soils to promote plant growth [96,97]. Some species of Leptospirillum, such as Leptospirillum ferrooxidans, are iron-oxidizing bacterium that contain all genes necessary for nitrogen fixation such as Mo-Fe nitrogenase, specific regulator (nifA), global regulators (glnB and ntrC), and other sensors and transport systems related to nitrogen assimilation. On the other hand, most denitrifying bacteria such as Pseudogulbenkiania, Dechlorobacter, and Flavobacterium were depleted in the treated soils [98][99][100], while Rivibacter [99] was enriched in the BC-treated soil.
The key taxa in the co-occurrence network analysis were members of the Verrucomicrobia, Gemmatimonadetes, Bacteroidetes, and Acidobacteria subgroup 4. Verrucomicrobia members are known as degraders of recalcitrant organic matter [101]. Gemmatimonas is highly abundant in soils with pyrogenic organic matter, which likely decomposes polyaromatic C [82,102]. In addition, these bacteria are capable of decomposing cellulose, lignocellulose, and chitin in biochar-amended soils [103]. Bacteriodetes are considered as specialists for degrading high molecular weight organic matter and are indicators of nutrient-rich soil [104,105]. Acidobacteria of subgroup 4, especially Blastocatella, were abundant in organic C content [79]. Overall, the co-occurring taxa associated with treated soils and those that regulate the microbiome structure in the treated soils are known to be related to organic matter decomposition and nutrient cycling, suggesting that soils treated with poultry litter and the combination of biochar and poultry litter favor beneficial bacteria for soybean growth.

Conclusions
This study used high-throughput amplicon sequencing to characterize the soil microbial community profiles in soybean field plots. The community structures were correlated with total N, total C, and pH in the 2 years. The co-occurrence networks revealed key taxa involved in organic matter decomposition and nutrient transformation in the treated soils. The PL treatment had the most profound impact on the microbial community among other treatments, followed by the BC/PL and BC treatments. Our study provides insights into the impact of organic fertilizing materials on the improvement of the microbial community to favor beneficial bacteria that promote plant growth and health. The comparable effect of PL and BC/PL treatments on soybean growth indicates a potential use of BC/PL to prevent nutrient loss and enrich the soil for an ecologically sound crop production strategy. Further studies with extended field trials accompanied by more comprehensive studies of the microbial community structure in the different compartments (e.g., bulk soil, soybean rhizosphere, and root endosphere) would substantially enhance our understanding of the mechanistic basis underlying the beneficial effects of the organic fertilizing materials. Furthermore, a parallel experiment such as a mesocosm study using field soil under controlled conditions would minimize the risk of possible confounding variables in the field.