Response of Soil Bacterial Community to Application of Organic and Inorganic Phosphate Based Fertilizers under Vicia faba L. Cultivation at Two Di ﬀ erent Phenological Stages

: It is essential to investigate to which extent and how speciﬁcally soil–plant–microbe interactions can be conditioned by di ﬀ erent agricultural practices. Legumes such as Vicia faba is one of the essential functional group in intercropping and crop rotations due to its higher N ﬁxing capacity. Hence, it is important to study the living microbial community of this legume. Further, it is also expected that ﬂuctuations in soil microbial diversity and composition could be complemented by plant phenological stages and di ﬀ erent fertilizer amendments. Thus, community composition in soil treated with phosphate-based inorganic and organic fertilizers, in the presence of Vicia faba plants at ﬂowering and fruiting time using NGS 16S rRNA gene amplicon sequencing. Further, the evaluation of plant biomass parameters under di ﬀ erent fertilizer treatments was also carried out. The presence of the Vicia faba plant increased the abundance of N ﬁxing bacterial such as Bardyrhizobium , Microvirga (Rhizobiales), Arthrobacter and Psuedoarthrobacter (Actinomycetales) in soil. Fluctuation in composition and diversity of bacterial community was further by plant phenological stages. These alterations could be due to changes that occurred in the plant nutrient requirement and varied root exudation patterns at a speciﬁc phenological stage. Further, fertilizer treatments also have a profound e ﬀ ect on the diversity and structure of the bacterial community. Organic fertilizers, especially vegetable tanned leather waste (VTLW), have a stronger e ﬀ ect on the composition and diversity of bacterial community compared to inorganic fertilizer (PT—triple superphosphate). Alpha-diversity was signiﬁcantly decreased by both organic and inorganic amendments, especially a species evenness because each fertilizer tends to stimulate the growth of distinctive microbes that dominated the community of amended soil. Microvirga a most genus that contributed most in co-occurrence pattern, suggests that these generalists are adapted to a variety of time, ecting M.T.C. and M.N.; Formal analysis, S.I.P., C.G., S.S. and M.N.; Methodology: S.I.P., C.G., M.N. and M.T.C.; Writing—original draft: S.I.P.; Writing—review & editing: S.I.P., C.G., M.N., M.T.C., S.S., G.P. and S.O.; supervision, M.T.C.; project administration, M.T.C. and S.O.; funding acquisition, M.T.C. O.L.


Introduction
Legumes are one of the cheaper and higher sources of proteins, minerals, and fibers for animals and humans [1], thus, it is a staple food in many developing countries. Further, legumes are a key functional group in intercropping and crop rotation and are highly valued by providing multiple services in the context of sustainable management [2][3][4]. Along with being served as basic high protein food and forage to humans and animals, legumes further contribute to improving soil health by nitrogen (N) fixation, phosphorus mobilization, water retention, and reducing greenhouse gas emission [1,5]. Thus, they have a high potential for sustainable agriculture being functional either as growing crops or as crop residue and promoting ecosystem efficiency [2,6]. Due to the higher nutritional values, such as protein, carbohydrates, B group vitamins, and minerals [7], the Faba bean (Vicia faba L.) is one of the important legumes around the world [8]. Moreover, Faba bean has high N fixing capacity and can fix approximately 50-330 kg N ha −2 [9,10]. Despite legume's N fixing capability through a symbiotic relationship with Rhizobium leguminosarum, by some means, they also require fertilization, particularly during the seedling time. Further, low nutrient availability in soil could have a negative effect on legume production [11]. Previous studies demonstrated that nodulated legumes have a strong direct positive relationship with soil phosphorus (P) and organic matter content [12][13][14][15]. Moreover, much of the nitrogen benefit of legumes to upcoming crop comes from the plant residue-shoots and roots-which demonstrated the need for enough soil nutrients to improve legume biomass that can be fulfilled by fertilizer amendments.
Soil microbial community is a key driver of energy flow and biogeochemical cycling in terrestrial ecosystems [16,17]. Hence, soil ecological processes preliminary depend on the composition, diversity, and stability of soil microbial communities [18,19]. Furthermore, soil microbial community is a significant soil quality indicator and is highly sensitive to different biotic and abiotic factors, such as soil nutrient, pH, plant species, and soil management practices [20][21][22][23]. Many studies showed that the soil microbial community is profoundly influenced by the use of legume in intercropping or crop rotation. These influence directly [24][25][26][27] or indirectly, through changes in rhizosphere pH, enhanced availability of nutrients, or the effect on the quantity, the quality, and the distribution of soil organic matter in upper-soil horizons [28,29]. Along with plant species, plant growth stages also have a strong influence on soil microbial community, which could be due to different nutrient requirements and varied root exudation patterns during different phenological stages of the plant [30][31][32]. Despite a substantial body of literature, indicating that legumes play a vital role in stimulating microbial community and their activities [1,33], none of them actually looked into the microbial community of living legumes and the time in which these interaction manifest.
Phosphorus (P) is the second most important macronutrient in agroecosystem [34]. Unlike nitrogen (N) and sulfur (S), P does not occur as gas or disappear into the atmosphere but can be easily fixed and bind to iron and aluminium or strongly absorbed on soil particles and become unavailable to plants. Soil microbial community plays a vital role in P solubilization and mineralization [35] and as a group these microbes are named P solubilizing microbes (PSB). Vice versa, organic and inorganic P has a profound impact on the abundance, diversity, and structure of soil microbial communities [23,36]. Recently, many meta-analysis studies demonstrated that organic and inorganic amendment with matched macronutrients might have a similar effect on plant production, while microbial community respond differently [37,38]. These demonstrate that nutrients are a limiting factor for plant production, though the soil microbial community is more correlated with their chemical and biochemical properties and other abiotic factors.
Hence, the present study aims to assess the response of the bacterial community to plant presence and time after the addition of phosphorus-based organic and inorganic amendment to the soil.
Soil bacterial community was analysed with and without Vicia faba cultivation at flowering and fruiting time using NGS 16S rRNA gene amplicon sequencing. We hypothesized that (i) diversity and structure of bacterial community would differ among different fertilizer amendments due to various biochemical characteristic; (ii) changes in the bacterial community would be supplemented by plant presence and sampling time due to differences in plant nutrient requirement and root exudation pattern during flowering and fruiting time; (iii) changes in a microbial community under the organic amendment would be due to the supplied organic material, not by fertilizer-born-microbes.

Soil Properties
Sandy soil was taken from an experimental field (0-10 cm depth) in Florence, Italy. After air-drying, the soil was sieved at 5 mm. After sieving, the soil was used for the pot experiment, and the main characteristic of soil is given in Table 1. is an organo-mineral fertilizer from vegetable tanned leather, certified as organic products. The waste material is hydrolyzed by sulphuric acid and then converted into a liquid form; natural soft rock phosphate is added to solidify and create an intermediate product. The transition from liquid to solid state occurs in about three months, a period called maturation. Then, the product is granulated for commercialization and contains 3% N mainly as amino acids and short-chain peptides, 11% phosphoric anhydride (P 2 O 5 ), 10% organic carbon (OC), 15% calcium oxide (CaO), and 12% sulphuric anhydride (SO 3 ).

Inorganic Fertilizers
Triple superphosphate (TP) and urea (U) (Fertilizzanti Certaldo Srl, Certaldo, Tuscany, Italy). TP, manufactured by letting the phosphoric acid react with phosphate rock, contained 46% of P 2 O 5 . Urea contained 46% of N. The main characteristics of used organic and inorganic fertilizers are given in Table 2.

Experimental Design
The experiment was three-factorial and the experimental plan included four fertilizer treatments, two soil sampling time during the growing season (flowering and fruiting), and one soil with and without plant presence. The treatments were performed in triplicates. To test the effect of the fertilizer on the soil chemical and biological properties, a total of 48 pots (30 cm × 30 cm × 30 cm) were filled with 18 kg of soil. Fertilizer treatments (Supplementary Table S1) comprised control (C) containing only U, VTLW plus U (P), PCM plus NSRP (PP), and TP plus U (PT). To evaluate the effect of the plant presence on soil properties, the test was performed in pots with plants (WP) and no plants (NP) in the open air. Treatments were set up by distributing the respective fertilizers (Supplementary Table S1) in each pot and then homogenously incorporating the fertilizer by manual mixing. Twenty-four pots were sown with three Vicia faba L. seeds per pot in October 2017, while the other 24 pots remained unseeded. All pots were arranged in a randomized pattern and were irrigated with a dripping system, depending on soil moisture. At flowering (FL) and at fruiting (FR) time, the soil was sampled for each of the pots.
Since the root system of the three plants grown in each pot completely occupied the available space, we assimilated the soil as ectorhizospheric. The control was the soil in a pot without plants. All the soil samples were stored at −20 • C prior to further analysis.

Plant Biomass Analysis
At maturity, the plants were separately harvested for analysing the length of stem and roots, root collar diameters, and the number of seeds per plant. For each pot, roots, aboveground biomass (stem and leaves), and seeds inside the pods were placed separately inside paper bags and kept in a drying oven at 60 • C for 3 days until a constant weight was achieved. Then, the same plant parts were separately weighed on an electronic balance to measure the dry weight. Roots were recovered from each pot and maintained separately. The soil was removed from the pot and disaggregated by hand carefully to keep the roots intact and recover as much as possible. The roots were then gently washed with distilled water to remove residual soil particles.

Soil Chemical Analysis
Inorganic nitrogen, ammonium-nitrogen (NH 4 + -N), nitrate-nitrogen (NO 3 -N), and nitrite-nitrogen (NO 2 -N), and urea-nitrogen (U-N) were analysed by the photometric method. Briefly, 40 g of air-dried soil from each sample was added to the 160 mL of a 0.01 M CaCl 2 solution (Sigma-Aldrich, 97%) and then was shaken for 2 h at 150 rev min −1 . Later, the suspensions were filtered through a nitrate-free filter and immediately analysed for N compounds or stored at −40 • C in the dark.

Soil DNA Extraction and 16S rRNA Gene Sequencing
Total soil and chicken manure (PP) DNA was extracted using the Fast DNA™SPIN Kit for soil (MP Biomedicals, Santa Ana, California, USA) according to Panico et al. [42]. Quality control and DNA yield were checked by 2% agarose gel electrophoresis and by spectrophotometer (Picodrop limited, Hinxton, UK).

Sequencing Data Processing and Data Analysis
Paired reads were assembled, quality-filtered, and analysed using the pipeline SEED 2.0.3 [44] with the inclusion criteria of mean quality score ≥32 and length ≥250 bp. Briefly, chimeric sequences were detected using the de novo VSEARCH algorithm [45] and removed from the dataset. Sequences were then clustered into operational taxonomic units (OTUs) at a 97% sequence identity threshold using the VSEARCH algorithm; consensus sequences were constructed for all clusters [44]. Low abundant sequences (≤5 of total count) were excluded from further analysis. Identification and the taxonomic assignment were done using representative sequences retrieved from RDP database [46] and the NCBI (https://www.ncbi.nlm.nih.gov/) using a 10-4 E value threshold. Sequences identified other than bacteria were discarded. The remaining sequences were used to create OTU table and then normalized by dividing sequences of individual OTU by dividing the sum of all sequences for each sample. Phylogenetic assignment to bacterial phyla and class level was based on best hits, by dividing the number of sequences belonging to each phylogenetic group by the total number of sequences in the given sample. All Illumina datasets were submitted to the European Nucleotide Archive under the study accession number PRJEB40215.
PAST 3.03 [47] and R Statistical Environment (R development Core Team, 2008) were used for all statistical analysis. Venn diagram was constructed to identify shared and unique OTUs among different fertilizer treatment, time, and plant presence/absence. Rarefaction and alpha diversity of OTUs was performed on resampled datasets with the same number of sequences randomly selected from all samples (20.000 sequences) using the SEED 2.0.3 software [44]. Three-way analysis of variance (ANOVA, fertilizer treatment × plant presence × time) was performed to check any significant effect of fertilizers, plant presence, and their interaction on the variability of alpha diversity data. F values resulting from ANOVA were used to find out if the average between two sets of data were significantly different or not. Further, multiple pairwise comparisons of means was done by Tukey's honestly significant difference (HSD) test at p < 0.05 level of significance to analyse the individual effects of each factor. Correlations between alpha diversity and physicochemical variables or plant production data were examined by means of a Spearman's rank correlation test. OTUs with >0.5% abundance were used to evaluate differences in beta diversity. Principal Coordinate Analysis (PCoA) and PERMANOVA test were conducted based on Bray-Curtis similarity distance to determine the distribution of diversity and statistical significance of beta diversity, respectively. Heatmap was constructed for the top 50 OTUs based on abundance at different taxonomic levels (Phyla to Genus) to identify which OTUs contributed most to discriminate the soil bacterial community within a different treatment, time, and plant presence. Correlations between beta diversity and physicochemical variables or plant production data were determined by canonical corresponding analysis (CCA). The forward selection was based on Monte Carlo permutation tests (permutations D 999) [48]. The ordination in the x-and y-axes and the length of the corresponding arrows indicated the importance of environmental variables in explaining the taxon distribution across communities. The co-occurrence of the 50 most abundant OTUs in microbial communities among all samples was analysed by calculating Spearman's rank coefficients (p) using R package Hmisc [49]. Subsequently, those significant (FDR-adjusted p value < 0.01) and robust (p ≥ 0.6) correlations between OTUs were exported as a GML format network file using R package igraph [50]. Network visualization was conducted using the Fruchtermann-Feingold layout of the interactive platform Gephi version 0.9.2. Possible keystone genera were those that demonstrated high betweenness centrality values [51]. The modular structure of the community was evaluated via the modularity index [52]. Table 2 shows the considered plant biomass data per pot according to soil treatments. No significant differences were observed between treatments in terms of stem length, stem and leaves dry weight, and roots collar diameter ( Table 3). The P treatment significantly reduced (p < 0.05) the roots length as compared to plants grown in C and PT treatments, while no significant difference was observed with respect to plants grown in PP. There was no difference in the roots length among the C, PP, and PT treatments. The heavier roots among those tested were measured in the PP treatment, followed, in decreasing order, by PT, C, and P. The roots dry weight measured in the P treatment resulted significantly lower (p < 0.05) than those measured in the other treatments. Both for the number of seeds per plant and seeds dry weight the highest values were measured in the PT treatment, followed by the PP, C, and P treatments. The number of seeds per plant and the seeds' dry weight measured in the P treatment were significantly lower (p < 0.05) than those measured in the PT treatment, while no significant difference was observed with respect to the C and PP treatment. Quantities are expressed per pot. Significant differences are shown in the column and are indicated in different lowercase letters (One-way ANOVA followed by Tukey post hoc test, p < 0.05). Means followed by a different letter indicate significant differences at p < 0.05.

Soil Inorganic N Compounds
Three-way ANOVA results showed that all three factors-fertilizers treatment, plant presence and sampling time-and their interactions significantly affected the different nitrogen forms (ammonia, nitrite, nitrate, urea, and total nitrogen); time and plant presence did not have a significant effect on urea and nitrate, respectively (Supplementary Table S2). Further, ammonia and nitrate contents were significantly higher in soil treated with VTLW (P) during flowering than the other fertilizer treatments (Supplementary Table S3), while fertilizer treatments did not have significant effect on Urea content (Supplementary Table S3).

Bacterial Alpha Diversity
After quality filtering, chimera cleaning, and removal of low abundant sequences (≤5 total count), 2,655,968 16S rRNA sequences, and 30,198 OTUs were obtained from a total of 51 samples (48 soil samples and 3 fertilizer samples-PCM). The rarefaction curve was reached to saturation for all samples, indicating the sequencing depth sufficient to cover detectable species in all samples (Supplementary Figure S1). Venn diagram revealed that 9744 OTUs (68.7%) were exclusively shared by soil with and without plant presence, while 2619 OTUs (18.3%) were only enriched in soil with plant, and 1949 OTUs (13.6%) in soil without plant ( Figure 1). Assessing time effect, results showed that 10,068 OTUs (71%) were shared by soil during flowering and fruiting time, though 1008 (7.1%) and 3105 (21.9%) OTUs were uniquely enriched in the soil during flowering and fruiting time, respectively. Further, a comparison among fertilizer showed that 8645 OTUs (48.7%) were shared by different Sustainability 2020, 12, 9706 7 of 20 fertilizer treatments. Whilst somehow, each treatment had exclusive OTUs in the range of 1241 to 1005 with the highest number enriched by C treatment (7.0%) and the lowest enriched in PP treatment (5.7%).
samples, indicating the sequencing depth sufficient to cover detectable species in all samples (Supplementary Figure S1). Venn diagram revealed that 9744 OTUs (68.7%) were exclusively shared by soil with and without plant presence, while 2619 OTUs (18.3%) were only enriched in soil with plant, and 1949 OTUs (13.6%) in soil without plant (Figure 1). Assessing time effect, results showed that 10,068 OTUs (71%) were shared by soil during flowering and fruiting time, though 1008 (7.1%) and 3105 (21.9%) OTUs were uniquely enriched in the soil during flowering and fruiting time, respectively. Further, a comparison among fertilizer showed that 8645 OTUs (48.7%) were shared by different fertilizer treatments. Whilst somehow, each treatment had exclusive OTUs in the range of 1241 to 1005 with the highest number enriched by C treatment (7.0 %) and the lowest enriched in PP treatment (5.7 %).  Table 4. Time had a consistent significant effect (p < 0.05), while fertilizer treatment had a significant effect on Shannon index and evenness not on Choa-1. Further, plant presence did not have any significant effect on alpha diversities (Table 4 and Supplementary  Table S4). Interaction of time with plant presence and did not have any significant effect on alpha diversities, while, the interaction between time and fertilizers had only significant effect on evenness. Interaction between plant presence and fertilizers had a significant influence on all three diversity indices (Table 4 and Supplementary Table S4).
Moreover, the interaction between all factors significantly affected the Shannon index and evenness, while this interaction did not influence Chao1 richness. Evenness was decreased considerably in fertilizer treated soils during both flowering and fruiting time compared to control soil, especially in P treatment. While Choa-1 richness was increased in treated soils, but the differences were not significant.   Table 4. Time had a consistent significant effect (p < 0.05), while fertilizer treatment had a significant effect on Shannon index and evenness not on Choa-1. Further, plant presence did not have any significant effect on alpha diversities (Table 4 and Supplementary Table S4). Interaction of time with plant presence and did not have any significant effect on alpha diversities, while, the interaction between time and fertilizers had only significant effect on evenness. Interaction between plant presence and fertilizers had a significant influence on all three diversity indices (Table 4 and  Supplementary Table S4). Moreover, the interaction between all factors significantly affected the Shannon index and evenness, while this interaction did not influence Chao1 richness. Evenness was decreased considerably in fertilizer treated soils during both flowering and fruiting time compared to control soil, especially in P treatment. While Choa-1 richness was increased in treated soils, but the differences were not significant.

Correlation of Alpha Diversity with Soil Edaphic Properties and Plant Biomass
The correlation of soil bacterial alpha diversity with soil chemical properties and plant biomass data are shown in Table 5a,b, respectively. In plant presence, only evenness was significantly and positively correlated with ammonia content during flowering time. In soil without plant, all soil edaphic factors are significantly but negatively correlated with evenness except nitrite during flowering time, although Shannon index and Chao-1 were also negatively (significant) correlated only with ammonia during flowering time. During fruiting time, only urea was negatively correlated with Shannon index and evenness.

Correlation of Alpha Diversity with Soil Edaphic Properties and Plant Biomass
The correlation of soil bacterial alpha diversity with soil chemical properties and plant biomass data are shown in Table 5a,b, respectively. In plant presence, only evenness was significantly and positively correlated with ammonia content during flowering time. In soil without plant, all soil edaphic factors are significantly but negatively correlated with evenness except nitrite during flowering time, although Shannon index and Chao-1 were also negatively (significant) correlated only with ammonia during flowering time. During fruiting time, only urea was negatively correlated with Shannon index and evenness.
Plant biomass parameters, root and stem length, and root dry weight were significantly and positively correlated only with evenness during flowering time. Further, Chao-1 was negatively (significant) correlated with root dry weight during flowering time. Any plant biomass data did not have a significantly correlation with Shannon index during flowering time, while Shannon index was significantly and positively correlated only with root length during the fruiting time. Further, stem and root length were also significantly correlated with evenness during fruiting time. Plant biomass parameters, root and stem length, and root dry weight were significantly and positively correlated only with evenness during flowering time. Further, Chao-1 was negatively (significant) correlated with root dry weight during flowering time. Any plant biomass data did not have a significantly correlation with Shannon index during flowering time, while Shannon index was significantly and positively correlated only with root length during the fruiting time. Further, stem and root length were also significantly correlated with evenness during fruiting time.

Changes in Bacterial Community Structure (Beta Diversity) in Response to Time, Plant Presence and Fertilizer Treatment
Beta diversity allows studying how different the population structure is in various environments. PCoA of Bray-Curtis distance was used to analyse the variation in the bacterial community as affected by plant presence, time, and fertilizers (Figure 2a-c). The significance level of variation was checked by PERMANOVA. The first two principal coordinators explain a high percentage of variance (~60%, coordinate 1:42.24% and coordinate 2:17.50%) with distinction in community structure associated with all three factors. Plots revealed that community clustered differently under all three factors, especially under plant presence and absence. PERMANOVA results confirmed a significant global Regarding the fertilizer effect, PCoA plot and pairwise PERMANOVA showed that the community treated with VTLW (P) was significantly different (p < 0.05) from control soil and the other fertilizer treatments. Moreover, significant differences were also observed between the community structure of C (control) and PP (PCM plus NSRP) treatment.

Correlation between Bacterial Community Structure and Soil Edaphic Properties and Plant Biomass Data
CCA analysis and Mantel test were performed to correlate soil and plant variables to bacterial community structure (Figure 3a,b, Table 6). Output of Mantel test indicated that different soil nitrogen forms (ammonia; p = 0.028, nitrite; p = 0.039, nitrate; p = 0.0002, urea; p = 0.032) have significant correlation with bacterial community. CCA plots showed that the first two axes explain 48.89% (p = 0.013) and 26.77% (p= 0.001) of total variance. Soil bacterial community formed distinct clusters by plant presence and fertilizer treatment, while weakly grouped into a cluster by sampling time (Figure 3a). Ammonia and nitrite had a stronger correlation with the bacterial community under the treatment of organic fertilizers (P and PP) especially, with plant presence, while urea and nitrate had a stronger correlation with the bacterial community under treatment of triple phosphate (PT) and control soil (C). Soil bacterial community without plant presence did not have a stronger correlation with N forms as shown by the community under plant presence. 0.0003), and root dry weight (p = 0.0002) are significantly correlated with bacterial community. CCA plot showed that first two axes explain 58.89% (p = 0.001) and 22.50% (p = 0.043) of total variance. Soil bacterial community are distinctly clustered by fertilizer treatments, while weakly separated by sampling time. Treatment of organic fertilizers has a stronger effect as they are clustering separately from each other and also from C and PT treatments. Stem length has a stronger correlation with a bacterial community of C and PT treatment, while stem dry weight has a stronger correlation with the community treated with PP.  Rt collar dia: roots collar diameters, Roots L: roots length, Rt dry wt: roots dry weight, gr nº: seeds nº, gr dry wt: seeds dry weight. To assess the link with plant variables, only soil bacterial community with plant presence were analysed (Figure 3b, Table 6). Mantel test showed that only stem length (p = 0.001), root length (p = 0.0003), and root dry weight (p = 0.0002) are significantly correlated with bacterial community. CCA plot showed that first two axes explain 58.89% (p = 0.001) and 22.50% (p = 0.043) of total variance. Soil bacterial community are distinctly clustered by fertilizer treatments, while weakly separated by sampling time. Treatment of organic fertilizers has a stronger effect as they are clustering separately from each other and also from C and PT treatments. Stem length has a stronger correlation with a bacterial community of C and PT treatment, while stem dry weight has a stronger correlation with the community treated with PP.

Bacterial Community Composition of PP and P Fertilizers
Extracting DNA from the pure PP fertilizer and sequencing the 16S rDNA, results revealed that the most represented phyla were Firmicutes (68%) and Actinobacteria (27.7%); the remaining were unknown (4.3%) (Supplementary Figure S2). Due to hydrolysis of sulphuric acid, the microbial community was not able to survive in VTLW (P) and, thus, we were not able to extract or sequence the bacterial community of VTLW (P) fertilizer.

Influence of Plant Presence, Time, and Fertilizer Treatment on Taxonomic Composition
To analyse the effect of all three factors on soil bacterial composition, we assessed the bacterial relative abundance at two different taxonomic levels, Phylum, and Class. More than 50 phyla were identified in the present study, but we only showed those present >1% (Figure 4). At the class level, only the ten most abundant classes are showed (Supplementary Figure S3). Overall, phylum Proteobacteria (~30%) with classes alpha, beta, delta, and gamma was the most abundant followed by Actinobacteria (~20%), Acidobacteria (~20%), Planctomycetes (Planctomycetia), Bacteriodetes (Chitinophagia and Cytophagia), Chloroflexi, Firmicutes, Verrucomicrobia, etc. All three factors had various influences on soil bacterial taxonomic compositions.   Actinobacteria, Bacteroidetes, Proteobacteria with classes alpha, beta, and gamma, and Verrucomicrobia significantly increased in soil with plant present, while Acidobacteria, Chloroflexi, Cyanobacteria Gemmatimonadetes were significantly enriched in soil without plant (Figure 4 and Supplementary Figure S3). Further, bacilli class was significantly enriched in soil with plant (Supplementary Figure S3). Time did not have a strong effect as plant at phylum level but had a stronger effect at class level taxonomy (Supplementary Figure S3). At class level only Bacteroidetes and Verrucomicrobia were enriched during flowering time, whereas relative abundance of Candidatus, Cyanobacteria Gemmatimonadetes increased during the fruiting time. Beta and gamma Proteobacteria were significantly increased during flowering time, while alpha and delta Proteobacteria were strongly enriched during the fruiting time.
Relative abundances of Acidobacteria and Planctomycetes and unknown bacteria were depleted in all fertilizer treatments (Figure 4 and Supplementary Figure S3). Actinobacteria were significantly enriched in soil treated with VTLW (P). The relative abundance of Verrucomicrobia was strongly increased in soil treated with chicken manure (PP), while the amendment of triple phosphate (PT) significantly increased the abundance of Cyanobacteria. In general, the plant presence, time, and treatments, especially the organic fertilizers (P and PP), induced the largest variations within classes. Delta-proteobacteria was significantly affected by treatments (PP > P > PT > C) while these did not have a significant influence on other Proteobacteria (Supplementary Figure S3).
Differences in the relative abundance of the 50 most abundant OTUs were shown by heatmap analysis (Supplementary Figures S4-S6). Evidently, these OTUs differed among plant presence, time, and fertilization but with different strengths. OTUs belonging to Actinomycetales (such as OTUs 1,2,21) and Rhizobiales (such as OTUs 24, 25) were significantly higher in plant presence, while OTUs identified as cyanobacteria (such as OTUs 9, 10, 14) or unknown (OTUs 13, 38, 59) were highly abundant without plant.
Time did not have a strong effect as plant presence. OTUs identified for cyanobacteria phylum were more abundant during fruiting time, while the abundance of OTU 54 (Paeniglutamicibacter) and OTU60 (Pseudomonas) were significantly higher during flowering time.
Among fertilizers, OTUs belonging to Cyanobacteria were significantly more abundant in PT treatments; Arthrobacter and Pseudoarthrobacter were significantly more represented in P samples. The remaining OTUs were more or less represented at the same amount in PP and C.

Bacterial Co-occurrence Network Analysis
Across all the considered variables, the correlation network analysis ( Figure 5) showed 71 strong positive correlations among 31 genera (R > 0.6, p < 0.01) with 24 internal type and 47 external type concurrence. The network of positive correlations formed three distinct major modules (8 ≥ OTUs). The average path length (APL) between two nodes was 2.67 edges with a diameter of 6 edges. The clustering coefficient (CC) was 0.68 and the modularity index (MD) was 0.67, where MD > 0.4 suggests that the network has a modular structure. All genera in the network were assigned to a different bacterial taxonomic level. Among these, Proteobacteria, Actinobacteria, and Cyanobacteria made up the largest proportions, accounting for 32.28, 26.01, and 19.25% of all nodes, respectively. Based on betweenness centrality scores, the top genera identified were Chelatococcus (OTU 23), Arthrobacter (OTU 1), Cyanobacteria (OTU 16), and Microvirga (OTUs 44, 25), which indicated the most influential bacteria and the critical roles they play in the co-occurrence network. 2020, 12, x FOR PEER REVIEW 15 of 23

5.
Network of co-occurring OTUs in different fertilizer treatments, based on correlation analysis. Each edge stands for a strong (Spearman's ρ > 0. cant (p < 0.01) correlation. The size of each node is proportional to the number of connections (i.e., degree). The nodes are also coloured by degree prop ickness of the edges is proportional to the robustness of a given Spearman's ρ, ranging from 0.6 to 0.9. Each edge stands for a strong (Spearman's ρ > 0.6) and significant (p < 0.01) correlation. The size of each node is proportional to the number of connections (i.e., degree). The nodes are also coloured by degree proportion. The thickness of the edges is proportional to the robustness of a given Spearman's ρ, ranging from 0.6 to 0.9.

Discussion
The variation in bacterial community composition in terms of beta diversity and taxonomy was manifested by plant presence. Driving force-underlying variances in the bacterial community can be directly due to easily available energy sources in the rhizosphere zone or indirectly due to changes in soil chemical and biochemical properties through root exudations [53,54]. Increased relative abundance of Actinobacteria, Proteobacteria, and Verrucomicrobia was also previously identified in the soil of Faba beans [18,55]. Increased abundance of N fixing bacteria such as Bardyrhizobium, Microvirga (Rhizobiales), Arthrobacter, and Psuedoarthrobacter (Actinomycetales) is anticipated in a leguminous plant. However, alpha diversity was not influenced by plant presence, which is supported by previous studies, which demonstrated that plant have less control over bacterial alpha diversity [54,56].
Fluctuation in composition and diversity of bacterial community was further supplemented by plant phenological stages. These alterations could be due to changes that occurred in the plant nutrient requirement and varied root exudation patterns at specific phenological stages [37,57]. Henceforth, the availability of varied organic substrates during flowering and fruiting time could selectively attract specific microorganisms and, thus, stimulate changes in the bacterial community. Results suggest the need of further analysis on the identification of root exudation patterns of Vicia faba during different phenological stages and consequences of a changed pattern on soil microbial community and their functions. The genus Pseudomonas is one of the more diverse genera; because of their catabolic versatility, excellent root-colonizing abilities, and capacity to produce a wide range of antifungal metabolites, the soil-borne Pseudomonas had received particular attention among PGPRs [58]. Its presence in the soil at flowering could be positive for the flower's formation of the plant. Cyanobacteria are one of the major components of the nitrogen fixing biomass in soil and provide a potential source of nitrogen fixation at no cost. Due to the important characteristic of nitrogen fixation, Cyanobacteria have a unique potential: to contribute to productivity in a variety of agricultural and ecological situations [59,60]. Temperature is determinant to enhance the growth rate of Cyanobacteria, as they surpass other taxonomical groups under high temperature [61][62][63]. Therefore, in addition to the plant presence effect, an increase in temperature between flowering and fruiting could also have been a factor that promoted the presence of Cyanobacteria.
Concerning the fertilizer effect, there is evidence that organic and inorganic amendments with equivalent total nutrient content have comparable fertilizer effect not only on crop yield [64,65] but also on soil microbial community [38]. Alpha-diversity was significantly decreased by both organic and inorganic amendments, especially species evenness because each fertilizer tends to stimulate the growth of distinctive microbes that dominated the community of amended soil [38]. Inorganic fertilizer (PT) contains easily available micro-and macronutrients, while organic fertilizers P and PP, comprise of both C and nutrients. Further, in VTLW, 90% of nitrogen is organic N such as amino acids and peptides, while chicken manure contains up to 70% of nitrogen in the form of urea and uric acid. Slow nutrient release fertilizer P had a larger effect on community structure than PP and PT. Previous findings also showed that the addition of nutrients or labile carbon with different C:N:P ratio to soil has frequently been shown to have a divergent and selective effect on the diversity and structure of microbial community [38,66]. Considering that, varying C compound and nutrient components tend to influence the microbial community [67,68], further fluctuations in community compositions are due to the contrasting chemical composition of the used fertilizers amendments. Recent findings by Tao et al. [69] support increase of Pseudomonas in VTLW treatment, as authors reported that bioorganic fertilizers tend to increase indigenous Pseudomonas. Cyanobacteria abundance tends to increase in the optimal combination of biogenic elements in soil, such as reduced nitrogen content under a sufficient amount of phosphorus [60], which could underpin the higher abundance of these bacteria in PT treatment, which contain higher phosphorus.
We have further looked into the effect of organic-fertilizer-borne microorganism on soil microbial community because, in some way, these fertilizers contain harmful microbes, antibiotics, etc. Due to the hydrolysis of sulphuric acid, the microbial community might be not able to survive in VTLW and, thus, we were not able to extract or sequence the bacterial community of VTLW fertilizer. While we successfully sequenced and identified poultry litter, bacterial community and results clearly showed that native soil bacterial community was not adopting fertilizer-borne microorganisms. Thus, this finding suggests that poultry litter-borne microorganisms entering the soil be outcompeted by native soil-borne bacteria, thus, indicating the suppressiveness of the native microbial community to harmful pathogenic microbes.
Network analysis of OTUs based on correlation analysis provides insight into co-occurrence patterns among the 50 most abundant taxa that were influenced interactively by all three studied factors; plant, time, and fertilizers. Seventy-one positive correlations among 32 different genera indicate metabolic cooperation may play an important role in modelling species co-occurrence [70]. Proteobacteria, Actinobacteria, and Cyanobacteria were the most abundant phyla that contributed most to co-occurrence patterns, which stipulate that these generalists are adapted to a variety of environments [71,72]. Chelatococcus, Cyanobacteria, Sphingomonas, and Microvirga had the top betweenness centrality (degree of the nodes stand together), indicating their high degree of the importance of these nodes in the co-occurrence network. Cyanobacteria has already been mentioned above. Chelatococcus is an aerobic denitrifying bacteria genus and attracted a lot of attention for its denitrification rate, comparable, if not higher, to that of anaerobic denitrification in biotrickling filters [73]. Moreover, Chelatococcus can reduce atmospheric dinitrogen when living in symbiosis with leguminous plants [74]. Sphingomonas is a bacterial genus isolated from a variety of anthropogeneously contaminated environments, including soil, sediments, and aquatic habitats. Species belonging to this genus were shown to possess abilities to degrade a variety of pollutants including polycyclic aromatic hydrocarbons and herbicides [75]. This metabolic function within agricultural fields has important implications for the longevity of polluting compounds and their mobility to the ground and surface water [76,77]. Microvirga is an Alphaproteobacterial root nodule bacteria that fix nitrogen and have been isolated from nitrogen-fixing nodules of legumes. Because Microvirga is not a close relative of Rhizobium, the close affinity of nif genes to those of Rhizobium suggests that these genes were acquired through horizontal transfer [78]. These results, clearly point out that plant, time, and fertilizers not only had singular but also shared effects on bacterial community as the top five genera found in co-concurrence network were somehow also affected by each factors individually. Further, there are few genera; OTUs 6,22,24,26,30,45, and 45 were not affected by any of the three factors, which specifies not only adoptive but resistance responsiveness of soil microbial community under different environmental circumstances.
The bacterial community was further affected by soil-plant properties, with plant roots (length and dry weight), stem length, and soil inorganic Nx as key drivers. While other plant variables did not have a significant correlation with the bacterial community. Microbial diversity was significantly correlated with soil inorganic N forms, especially the community treated with VTLW. N in VTLW was mainly available in amino acid and small peptides forms, which allow constant stimulation to the N cycling community due to slow-release tendency of given fertilizer. Thus, future work is required to study the molecular ecology of bacteria, encoding different N cycling-genes to identify responsible microbes that are performing these functions under Faba bean plantation. The plant roots and the microbial community influence each other [79]. In fact, they are able to micro-engineer their habitats by changing the porosity and clustering properties of the soil pores. Biota act to significantly alter their habitat toward a more porous, ordered, and aggregated structure that has important consequences for functional properties [80]. Our results showed that the bacterial community was significantly affected by two plant properties: root length, and root dry weight. These observations support the hypothesis that the soil-plant-microbe system is self-organized.
Moreover, the plant root system is involved in the acquisition of nutrients and water, synthesis of hormones as well as anchoring to the plant; and crop growth and yield formation is significantly affected by roots. So, we considered some root and seed parameters for a comparison between organic and inorganic treatments. Plant biomass data of PP samples (root length and dry weight, seed number, and dry weight) were significantly higher than C and P samples. On the other side, the same parameters (except root dry weight) in PP were significantly lower than in PT samples. However, considering the number and the dry weight of seeds in PP, the values were only 18.8% and 28.3% lower than in the chemical PT samples. This gap can be considered in line with the one calculated using a new large meta-dataset [81]. A question still discussed among the proponents of sustainability in agriculture concerns crop yield. There are, in fact, those who argue that sustainability should be as productive as conventional agriculture and others who argue that sustainable could be less productive than conventional in exchange for better protection of the soil, the environment, and, consequently, of human health.

Conclusions
In this study, we determined the effect of phosphorus-based fertilizers with equivalent total nutrients on soil bacterial community of Vicia faba. Our results suggest that variation in diversity and structure of bacterial community among different fertilizer amendment was due to fertilizer biochemical characteristics and further changes were supplemented by plant and sampling time due to differences in nutrient requirement and root exudation patent during the flowering and fruiting time. Further, changes in the bacterial community under organic amendment were due to supplied organic material, not by fertilizer-born-microbes. With this study, we offer a contribution also to reflect on the importance of improving the fertility of agricultural soil by stimulating soil microbial community, especially different plant growth-promoting bacteria, even if this involves a slightly lower production yield. Maintaining soil fertility for future generations is well worth a small renunciation. However, future research should investigate the soil functional diversity of microbial groups involved in different nutrient cycling under given treatments to the eventual selection of specific biochemical mechanisms.