Elucidating the association between the phoD-harboring bacterial community and soil ecosystem multifunctionality, which is crucial for the comprehension of the phoD-harboring bacterial role and contribution in agro-ecosystems, is an essential but rarely investigated subject. Here, we explored the phoD-harboring bacterial community in long-term fertilized soils using amplicon sequencing and multiple analysis methods including the null, neutral, and niche breadth models. We found distance-decay relationships of community similarities against geographical distance on a large spatial scale. Community dissimilarity was significantly lower in the organic fertilization treatment (M) than that in the no (CK) and mineral (NPK) fertilizer treatments. Dispersal limitation governed community assembly in CK, M, NPK, and whole samples, with corresponding relative contributions of 58.2%, 58.3%, 52.8%, and 54.4%, respectively. Electrical conductivity, total carbon, total nitrogen, total phosphorus, organic phosphorus, and available phosphorus were responsible for the community assembly of phoD-harboring bacteria. Multiple model analysis revealed that the phoD-harboring bacterial community was less constrained by the environment and presented flexible metabolism in soils with the M fertilization treatment. phoD-harboring bacteria presented more conflicting interaction and exhibited significantly higher ecosystem multifunctionality in soils with the M fertilization treatment than that in the CK and NPK fertilization treatments. To our knowledge, this is the first study to report a less environment-constrained phoD-harboring bacterial community might lead to a larger difference in ecosystem multifunctionality in fertilized soils. Therefore, we suggest phoD-harboring bacterial community assembly could be a biotic indicator for evaluating soil ecosystem multifunctionality.