Distribution of Genes and Microbial Taxa Related to Soil Phosphorus Cycling across Soil Depths in Subtropical Forests

: Although many studies have focused on the roles of soil microbes in phosphorus (P) cycling, little is known about the distribution of microbial P cycling genes across soil depths. In this study, metagenomic sequencing was adopted to examine the differences in the abundance of genes and microbial taxa associated with soil P cycling between organic and mineral soil in subtropical forests. The total relative abundance of inorganic P solubilizing genes was the highest, that of P starvation response regulating genes was second, and organic P mineralizing genes was the lowest. The soil organic carbon concentration, N:P ratio, and available P concentration were higher in the organic soil than the mineral soil, resulting in abundances of organic P mineralizing genes ( appA and 3-phytase ), and inorganic P cycling genes ( ppa ), whereas those of the inorganic P cycling genes ( gcd and pqqC ) and the P starvation response regulating gene ( phoR ) were higher in mineral soil. The four bacteria phyla that related to P cycling, Proteobacteria , Actinobacteria , Bacteroidetes , and Candidatus _ Eremiobacteraeota were higher in organic soil; conversely, the three bacteria phyla ( Acidobacteria , Verrucomicrobia , and Chloroﬂexi ) and archaea taxa were more abundant in mineral soil. Therefore, we concluded that the distribution of genes and microbial taxa involved in soil P cycling differed among soil depths, providing a depth-resolved scale insight into the underlying mechanisms of P cycling by soil microorganisms in subtropical forests


Introduction
Phosphorus (P) is the key limitation nutrient for tree growth in subtropical forests [1]. Although total soil P content may be sufficient [2], only small quantities of inorganic P-namely orthophosphate (H 2 PO 4 − and HPO 4 2− ) ions-can be directly absorbed by plant roots [3]. The dominant soil P pools, including low-solubility inorganic and highly complex organic P forms, can be transformed into orthophosphates by physical chemical reactions (i.e., dissolution and desorption) and biological activities (i.e., mineralization). Soil microbes play an essential role in regulating P cycling and available P concentration [4]. Furthermore, natural forest ecosystems are more severely P-limited than cropland ecosystems [5], highlighting the importance of the soil microbial regulation of P cycling in forest ecosystems. Three microbial gene groups, including (1) inorganic

Site Description
This study was conducted in Dashanchong Forest Park (28 • 23 58 -28 • 24 58 N, 113 • 17 46 -113 • 19 08 E) in Changsha County, Hunan Province, China. A hilly topography with an altitude of 55−217 m a.s.l. has been characterized in the park. The area has a humid mid-subtropical monsoonal climate with a mean annual precipitation of 1416 mm and mean monthly temperatures ranging from −10.3 to 39.8 • C. The soil is well-drained clay loam that developed from slate and shale parent rock and is classified as Ferralsols [16], which is very shallow [17].
This study operated in two typical secondary forests: (1) Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest (PLF), and (2) L. glaber-Quercus glauca evergreen broadleaved forest (LGF). These two forests are approximately 70 years old. There are almost no herbs and very sparse shrubs and saplings due to the closed canopy. A 200 m wide valley separates the PLF (28 • 24 36 N, 113 • 18 12 E) and LGF (28 • 24 32 N, 113 • 18 18 E). The environmental factors (e.g., elevation, soil moisture, pH, and nutrient availability) differed between PLF and LGF [18]. Four plots with 30 m × 30 m area were established in each forest. The species, diameter at breast height, and height of all trees were identified and recorded. The stand characteristics were detailed in the work of Ouyang et al. [19] and Table S1.

Soil Samples Collection
In June 2022, we randomly sampled seven soil cores (10 cm in internal diameter and 20 cm in depth) in the central 10 m × 10 m area of each plot. Depending on the soil color, every soil core was separated into two soil depths: a surface layer (brownish-black color; organic soil) and a deeper layer (brownish-yellow color; mineral soil). The soils at the same layer were mixed thoroughly, and about 500 g of mixed soil was collected for each sample. A total of 16 soil samples were obtained. The stones and roots of each sample were removed by a 2 mm mesh and then divided into three subsamples: one for measuring soil moisture; one for determining soil pH and the concentrations of soil organic carbon (SOC), total P (TP) and nitrogen (TN), inorganic nitrogen (NH 4 + -N and NO 3 − -N), and available soil P (AP); and one for metagenomic sequencing.

Soil Chemical and Physical Analysis
Soil TP concentration was measured using the Mo-Sb colorimetric method [20]. Soil AP concentration was quantified using a solution of 0.05 M HCl and 0.025 M H 2 SO 4 [21]. The SOC concentration was measured using the K 2 Cr 2 O 7 -H 2 SO 4 oxidation method. Soil TN concentration was measured using the Kjeldahl method. The inorganic N was extracted using 0.5 M K 2 SO 4 solution, and then measured the NH 4 + -N and NO 3 − -N concentrations of the filtered extract by a flow injection analyzer (FIAstar 5000, FOSS, Höganäs, Sweden). Soil moisture was quantified by drying the samples at 105 • C to constant weight. Soil pH was determined at a soil-to-water ratio of 1:2.5 by using an FE20 pH meter (Mettler Toledo Instrument Co., Ltd., Shanghai, China). Soil properties are presented in Table 1. Table 1. Physico-chemical properties (mean±SD) of organic and mineral soil in the two subtropic forests. The letters a and b indicate significant differences between organic and mineral soil in PLF, as do x and y in LGF.

Metagenomics Analysis
The adaptors of sequenced reads were trimmed. Reads with lengths shorter than 50 bp or quality value lower than 20 or containing N bases were removed using fastp software (v0.20.0) [22]. Metagenomics data were assembled by MEGAHIT software (v1.1.2) [23] using succinct de Bruijn graphs. The final assembled contigs with lengths ≥300 bp were selected, and their open reading frames were predicted by Meta-Gene [24]. The predicted open reading frames with lengths of 100 bp or more were retrieved and translated into amino acid sequences through NCBI. A non-redundant gene catalogue was constructed using CD-HIT [25] with 90% sequence identity and 90% coverage. After quality control, reads were mapped to the non-redundant gene catalog with 95% identity using SOAP aligner [26], and gene abundance was evaluated for each sample. Overall, a total of~0.77 billion high-quality reads were produced, ranging from 41.24 to 55.99 million reads per sample.
Representative sequences of the non-redundant gene catalogue were aligned to the NCBI NR database with an e-value cutoff of 1 × 10 -5 using Diamond [27] for taxonomic annotations. The KEGG annotation was conducted using Diamond against the Kyoto Encyclopedia of Genes and Genomes database with an e-value cutoff of 1 × 10 -5 . The abundance of a taxonomic group was calculated by summing the abundance of genes annotated to a feature. Relative gene abundances (%) were normalized to the annotated read number across all samples for subsequent analysis. In this study, 27 functional genes (Table S2) related to organic P mineralization, inorganic P solubilization, and P starvation response regulation were selected for further investigation.

Statistical Analysis
The statistical analyses were conducted by R software (version 4.2.2) and on the Majorbio Cloud Platform. The variations in the gene composition and associated microbial taxa for P transformation between samples were evaluated by ANOSIM with the "vegan" R package [28], and then displayed by nonmetric multidimensional scaling plots (NMDS) via the Bray-Curtis dissimilarity matrix. Significant differences in the relative abundance of P cycling genes between organic and mineral soil were quantified by oneway analysis of variance (ANOVA). Circus has developed to quantify the corresponding relationships between soil microbial taxonomic groups involved in P cycling and samples using Circos-0.67-7 software (http://circos.ca/), and phyla with an abundance of <0.01 were merged with others. Species and functional contribution analysis were performed for the major phyla of archaea (Candidatus_Bathyarchaeota, Candidatus_Thermoplasmatota, Euryarchaeota, Thaumarchaeota), bacteria (Acidobacteria, Actinobacteria, Bacteroidetes, Can-didatus_Eremiobacteraeota, Chloroflexi, Gemmatimonadetes, Planctomycetes, Proteobacteria, Verrucomicrobia), and fungi (Ascomycota, Basidiomycota). The key genes associated with soil P cycling to determine the available soil P were identified using random forest analysis [29].

Differences in Relative Abundance of P Cycling Genes across Soil Depths
The composition of P cycling genes varied significantly between organic and mineral soil (p = 0.005), whereas the differences in forest types were insignificant ( Figure S1a). There were substantial differences among the total relative abundance of genes that function for inorganic P solubilizing, organic P mineralizing, and P starvation response regulating ( Figure 1). With individual genes functioning for organic P mineralizing, only the appA gene coding for phytase was more abundant in surface organic soil than in deeper mineral soil. As for inorganic P solubilizing genes, the ppa gene in PLF and LGF, as well as ppx in LGF, were more abundant in surface organic soil than in deeper mineral soil, whereas the gcd gene in PLF and LGF, as well as the pqqC gene in PLF, were lower in organic soil. The genes functioning for P starvation response regulating, including phoB in LGF and phoR in PLF and LGF, were more abundant in deeper mineral soil than in surface organic soil ( Figure 1).
Forests 2023, 14, x FOR PEER REVIEW 5 of 15 Figure 1. Differences in the relative abundance of microbial genes related to P cycling between organic and mineral soil in PLF (a) and LGF (b). PLF represents the Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest, LGF represents the L. glaber-Quercus glauca evergreen broadleaved forest. The red stars * represent significant differences between different soil depths for a given gene. The black stars * and *** represent significant differences between the total relative abundance of microbial genes involved in P-mineralization, P-solubilization, and P-starvation response regulation.

Differences in Taxonomic Assignments of P Cycling Genes across Soil Depths
The soil microbial taxa associated with P cycling genes varied significantly between organic and mineral soil (p = 0.002), whereas the differences in forest types were insignificant ( Figure S1b). In the similarities analysis via the NCBI-NR database, 99.86% of sequences that were associated with P cycling genes were matched as bacteria, 0.13% as archaea, and 0.01% as fungi in organic soil, whereas in mineral soil, 99.68% were compared as bacteria, 0.31% as archaea, and 0.01% as fungi. Although the total abundance of these three taxa varied insignificantly among organic and mineral soil, archaea abundance was higher in the latter than in the former (Figure 2b). The phyla Proteobacteria, Acidobacteria, and Actinobacteria were the dominant soil bacterial communities (Figure 2a), which were the most crucial taxonomic groups associated with soil P cycling ( Figure 3 and Table  S3). Organic soil has a higher abundance of the phyla Actinobacteria, Bacteroidetes, and Can-didatus_Eremiobacteraeota than mineral soil, whereas Acidobacteria, Verrucomicrobia, and Chloroflexi were higher in mineral soil ( Figure 2b). Differences in the relative abundance of microbial genes related to P cycling between organic and mineral soil in PLF (a) and LGF (b). PLF represents the Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest, LGF represents the L. glaber-Quercus glauca evergreen broadleaved forest. The red stars * represent significant differences between different soil depths for a given gene. The black stars * and *** represent significant differences between the total relative abundance of microbial genes involved in P-mineralization, P-solubilization, and P-starvation response regulation.

Differences in Taxonomic Assignments of P Cycling Genes across Soil Depths
The soil microbial taxa associated with P cycling genes varied significantly between organic and mineral soil (p = 0.002), whereas the differences in forest types were insignificant ( Figure S1b). In the similarities analysis via the NCBI-NR database, 99.86% of sequences that were associated with P cycling genes were matched as bacteria, 0.13% as archaea, and 0.01% as fungi in organic soil, whereas in mineral soil, 99.68% were compared as bacteria, 0.31% as archaea, and 0.01% as fungi. Although the total abundance of these three taxa varied insignificantly among organic and mineral soil, archaea abundance was higher in the latter than in the former (Figure 2b). The phyla Proteobacteria, Acidobacteria, and Actinobacteria were the dominant soil bacterial communities (Figure 2a), which were the most crucial taxonomic groups associated with soil P cycling ( Figure 3 and Table S3). Organic soil has a higher abundance of the phyla Actinobacteria, Bacteroidetes, and Candidatus_Eremiobacteraeota than mineral soil, whereas Acidobacteria, Verrucomicrobia, and Chloroflexi were higher in mineral soil (Figure 2b).
Except for the phoA and phoD genes found in bacterial and fungi communities and phnL located in bacterial and archaea communities, all the organic P mineralization genes were present in bacteria ( Figure 3). The genes functioning for inorganic P solubilizing and P starvation response regulating were present in bacteria and archaea. A significant positive relationship was shown between the beta (β) diversity of genes and the microbial taxa that are associated with soil P cycling (Figure 4), and mineral soil has a higher β diversity of both phyla and genes than organic soil. Except for the phoA and phoD genes found in bacterial and fungi communities and phnL located in bacterial and archaea communities, all the organic P mineralization genes were present in bacteria ( Figure 3). The genes functioning for inorganic P solubilizing and P starvation response regulating were present in bacteria and archaea. A significant positive relationship was shown between the beta (β) diversity of genes and the microbial taxa that are associated with soil P cycling (Figure 4), and mineral soil has a higher β diversity of both phyla and genes than organic soil.   Table S3.  Table S3. Correspondence between species (phyla taxa) and functional (KEGG gene). The contribution of a species indicates the primary species composition for a specific function; the contribution of a function represents the main functions in which a specific species is involved. The contributions are the average of organic and mineral soil. Detailed data are shown in Table S3.

Linkages between P Cycling Genes and Soil P Status
The relative abundance of P cycling genes positively and significantly correlated with soil N:P ratio (p < 0.001; Figure 5a) and available soil P concentration (p < 0.05; Figure 5b).

Linkages between P Cycling Genes and Soil P Status
The relative abundance of P cycling genes positively and significantly correlated with soil N:P ratio (p < 0.001; Figure 5a) and available soil P concentration (p < 0.05; Figure 5b). A random forest analysis revealed that 7 of the 27 P cycling genes contributed to mediating the concentration of available soil P (Figure 5c). A Pearson correlation analysis displayed that the available soil P concentration significantly increased with gene abundances of appA, 3-phytase, and ppa (p < 0.01; Table S4), whereas decreased with that of genes pqqC, phoR, phnP, and gcd (p < 0.05; Table S4). A random forest analysis revealed that 7 of the 27 P cycling genes contributed to mediating the concentration of available soil P (Figure 5c). A Pearson correlation analysis displayed that the available soil P concentration significantly increased with gene abundances of appA, 3-phytase, and ppa (p < 0.01; Table S4), whereas decreased with that of genes pqqC, phoR, phnP, and gcd (p < 0.05; Table S4). Figure 5. The relationships between soil microbial P cycling genes and soil P status. The relationship between soil N:P ratio and the relative abundance of all genes involved in soil microbial P cycling (a), and between the relative abundance of all genes involved in soil microbial P cycling and available soil P concentration (b), are shown. Panel (c) shows the significant (p < 0.05) gene predictors of available soil P, which were identified by random forest analysis. The data of organic and mineral soil dates were analyzed together. Soil N:P ratio represents total soil nitrogen concentration divided by total soil P concentration. PLF represents the Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest, LGF represents the L. glaber-Quercus glauca evergreen broadleaved forest. R 2 represents the coefficient of determination. * p < 0.05, *** p < 0.001.

Key Genes Associated with Soil P Cycling in Organic and Mineral Soil
The inorganic P solubilization genes were the most abundant in both organic and mineral soil in subtropical forests (Figure 1), indicating that inorganic P solubilizing may Figure 5. The relationships between soil microbial P cycling genes and soil P status. The relationship between soil N:P ratio and the relative abundance of all genes involved in soil microbial P cycling (a), and between the relative abundance of all genes involved in soil microbial P cycling and available soil P concentration (b), are shown. Panel (c) shows the significant (p < 0.05) gene predictors of available soil P, which were identified by random forest analysis. The data of organic and mineral soil dates were analyzed together. Soil N:P ratio represents total soil nitrogen concentration divided by total soil P concentration. PLF represents the Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest, LGF represents the L. glaber-Quercus glauca evergreen broadleaved forest. R 2 represents the coefficient of determination. * p < 0.05, *** p < 0.001.

Key Genes Associated with Soil P Cycling in Organic and Mineral Soil
The inorganic P solubilization genes were the most abundant in both organic and mineral soil in subtropical forests (Figure 1), indicating that inorganic P solubilizing may be the dominant resource for the supply of available soil P (soluble orthophosphate), which contrasted with our hypothesis that surface soil has a higher potential to mineralize organic P compared to deeper soil. Owing to subtropical forest soils being highly weathered and acidified [30], P is often adsorbed or precipitated as inorganic forms (Fe-P and Al-P) [31]; as a result, the available soil P pool is tiny [32], and the insoluble inorganic P pool is large [33]. In response to low soil P availability and a sizeable inorganic P pool, the inorganic P solubilization genes will trigger to increase their abundance.
The inorganic P solubilizing genes, pqqC and gcd, were the key predictors of enhanced soil P cycling in mineral soil (Figure 6), which verified the previous results [8,9]. This result can be attributed to the following three aspects. First, the genes pqqC and gcd were more abundant in mineral soil than in organic soil (Figure 1). The gcd gene encodes the glucose dehydrogenase (GCD) and the pqqC gene codes the pyrroloquinoline quinone synthase C (PqqC) that is involved in the synthesis pathway of pyrroloquinoline quinone (PQQ) [34]. The compound forming with GCD and the redox cofactor PQQ is essential to produce gluconic acid by microbial [35]. Gluconic acid is considered the most important organic acid in the solubilization of the recalcitrant inorganic P [36]. Second, there was a negative relationship between the available soil P concentration and the relative abundance of genes pqqC and gcd (Table S4). Consistent with a previous study [37], we observed that the available soil P concentration was higher in surface organic soil than in deeper mineral soil, which indicated that the pqqC and gcd genes were more sensitive to a low soil P availability than other P cycling genes. Third, mineral soil tends to increase inorganic P content in the forest ecosystem. Soil organic carbon has been identified as a key driver of soil TP concentration and the distribution of its forms [38]. In this study, the SOC concentration was significantly lower in deeper mineral soil than in surface organic soil ( Table 1), indicating that more inorganic P existed in the mineral soil. Therefore, the larger inorganic P pool would stimulate the inorganic P solubilization genes, including pqqC and gcd.
The phoR gene (P starvation response regulation gene) may also be a key marker of the soil microbial regulation capacity of P cycling in mineral soil [4,14]. In both PLF and LGF, the phoR gene was more abundant in deeper mineral soil than in surface organic soil (Figure 1), which agreed with previous studies, which showed that low P conditions activate the phoR gene [39], and our hypothesis. The phosphate (Pho) regulon is controlled by the genes phoB, phoR, and phoU which regulates soil inorganic P transformation [40]. Thus, we inferred that the activated phoR gene might be necessary for regulating the pqqC and gcd genes. This assumption was supported by the fact that the genes phoR, pqqC, and gcd were considered as determinants of the concentration of available soil P via the random forest analysis (Figure 5), and that these three gene abundances significantly and negatively correlated with available soil P concentration through the Pearson correlation analysis (Table S3). In addition, the genes phoB and phoR were more abundant in deeper mineral soil than in the surface organic soil of LGF with low soil TP. In contrast, there were no significant differences in the genes phoB and phoU between the organic and mineral soil of PLF with high soil TP (Figure 1 and Table 1). This result is explained by the fact that phoB is activated by phoR under P-low environments, whereas phoB is triggered by phoU in P-rich environments [4]. The phoR gene (P starvation response regulation gene) may also be a key marker of the soil microbial regulation capacity of P cycling in mineral soil [4,14]. In both PLF and LGF, the phoR gene was more abundant in deeper mineral soil than in surface organic soil (Figure 1), which agreed with previous studies, which showed that low P conditions activate the phoR gene [39], and our hypothesis. The phosphate (Pho) regulon is controlled by the genes phoB, phoR, and phoU which regulates soil inorganic P transformation [40]. Thus, we inferred that the activated phoR gene might be necessary for regulating the pqqC and gcd genes. This assumption was supported by the fact that the genes phoR, pqqC, and gcd were considered as determinants of the concentration of available soil P via the random forest analysis (Figure 5), and that these three gene abundances significantly and negatively correlated with available soil P concentration through the Pearson correlation analysis (Table S3). In addition, the genes phoB and phoR were more abundant in deeper The key genes that functioned in P cycling differed between organic and mineral soil (Figures 1 and 6). The ppa gene, which encodes an inorganic pyrophosphatase to hydrolyze inorganic polyphosphate compounds, was the key inorganic P solubilization gene in organic soil [41]. Inorganic polyphosphate, in chains of tens to hundreds of phosphate residues, is environmentally ubiquitous and abundant, such that it is found in every cell in nature [42]. Accordingly, higher levels of inorganic polyphosphate are accumulated in surface organic soil than in deeper mineral soil, as it is derived from forest floor litter decomposition and significant microbial population turnover [43]. Moreover, the ppx gene was more abundant in the surface organic soil than in the deeper mineral soil of LGF, whereas for PLF, that of ppx and ppk1 was slightly higher in the organic soil ( Figure 1). Exopolyphosphatase-which releases orthophosphate anions from inorganic polyphosphate-is encoded by the ppx gene [41], and polyphosphate kinase-which catalyzes the formation of polyphosphate-is coded by the ppk1 gene [44]. Therefore, the high levels of inorganic polyphosphate in organic soil may trigger inorganic P solubilization genes, especially ppa.
Despite the total relative abundance of organic P mineralizing genes being the lowest in subtropical forest soils (Figure 1), the appA and 3-phytase genes were found to have a vital role in determining the available soil P concentration in organic soil ( Figure 5 and Table S4). The appA and 3-phytase genes encode phytases that hydrolyze the phytate [45]. Phytates are the dominant organic P forms in soil [46]. Previous studies suggested that the 3-phytase gene was abundant in P-deficient soil [47], which indicates that mineralizing phytate as a P-source is critical in P-deficient soils [14]. Based on N:P stoichiometry, a high soil N:P ratio implies a high N concentration and/or low P concentration, which may increase the P-starvation of soil microorganisms [48]. Our results revealed that the total relative abundance of P cycling genes significantly increased with the soil N:P ratio ( Figure 5), and the appA gene was more abundant in the surface organic soil (higher N:P ratio) than in deeper mineral soil (lower N:P ratio), which was in line with a previous study [8]. We concluded that the high soil N:P ratio in organic soil stimulated the appA and 3-phytase genes to encode phytase, which improved the mineralization of phytate and then increased the concentration of available soil P.

Phosphorus Cycling Genes Harboring Microbial Taxa Change with Soil Depths
Consistent with previous studies, P-cycling genes harboring microbial taxa distributed among bacteria, fungi, and archaea [14]. As bacteria are better studied and understood in soil P cycling [e.g., primers targeting key P cycling genes have been designed for bacteria, numerous bacterial genomic data], this has resulted in most P-cycling genes being found in soil bacteria [49]. Interestingly, several genes were also presented in fungi and archaea. For example, the phoA and phoD genes presented in fungi agreed with previous studies showing that phoD is widely distributed in different soil microbial taxa [50]. Genes associated with inorganic P solubilizing and P starvation response regulating were also presented in archaea, suggesting that archaea have a critical role in soil inorganic P solubilizing and P starvation responding [14,51]. The species and functional regression analyses also denoted that the soil microbial taxa harbor more distinct P cycling genes in mineral soil (Figure 4b), indicating that more diverse soil microbial taxonomic groups have a higher potential to mineralize and solubilize soil P in subtropical forests.
The bacteria phyla Actinobacteria, Bacteroidetes, and Candidatus_Eremiobacteraeota decreased in abundance with depth ( Figure 2), which may be attributed to their copiotrophic behavior [52,53]. The significant increases in the abundance of Acidobacteria, Verrucomicrobia, and Chloroflexi with depth ( Figure 2) have been previously evidenced by their tolerance for nutrient-poor conditions [54,55]. The copiotrophic hypothesis states that copiotrophic taxa seem to increase in nutrient-rich environments, whereas oligotrophic taxa would likely decrease [56,57]. Furthermore, copiotrophic taxa exhibit fast growth rates with high soil C availability, whereas oligotrophic taxa grow slowly by metabolizing recalcitrant C and nutrient-poor substrates [57]. The variations in copiotrophic and oligotrophic taxa with soil depth may be directly influenced by the decreases in the available soil P concentration, N:P ratio, and SOC, or they are indirectly affected by the increase in soil pH in mineral soil [6,53,58]. For example, the abundance of oligotrophic Acidobacteria significantly decreased with the soil N:P ratio (R 2 = 0.35, p < 0.05) in this study. In contrast, it increased with the soil pH (R 2 = 0.42, p < 0.05), which is consistent with a previous study [8]. Contrary to the copiotrophic hypothesis, the abundance of copiotrophic proteobacteria did not significantly differ between organic (43.45% in PLF, 43.17% in LGF) and mineral (40.24% in PLF, 39.06% in LGF) soil (Figure 2), which suggested that proteobacteria is the predominant contributor in soil P cycling in subtropical forests.
According to previous studies [53,55], the archaea taxa associated with soil P cycling was more abundant in deeper mineral soil than in surface organic soil (Figure 2). Several reasons have been documented for the increasing abundance of archaea taxa with soil depth, including their adaptation to nutrient-limited conditions as slow-growing oligotrophs, adaption to chronic energy stress, preference as methanogens for anaerobic conditions, and function as ammonia oxidizers to stimulate autotropic nitrification in deeper soil [53,54,59,60]. In addition, the total relative abundance of fungi taxa that harbored soil P cycling genes was lower than 0.01% (Figure 2). Compared to bacteria and archaea, the potential of fungi associated with the underlying mechanisms stimulated by changes in soil P availability seems to be more limited [14].

Conclusions
The key genes associated with soil P cycling in organic soil were ppa, appA, and 3-phytase, indicating that soil microorganisms had the potential to both mineralize organic P and solubilize inorganic P in organic soil, whereas those in mineral soil were phoR, gcd, and pqqC, suggesting that soil microbial P starvation response regulating genes were stimulated to improve inorganic P solubilization in mineral soil. The P cycling genes not only contained bacteria, but also harbored archaea and fungi. The bacteria that function in soil P cycling were Proteobacteria, Actinobacteria, Bacteroidetes, and Candidatus_Eremiobacteraeota in organic soil, and Proteobacteria, Acidobacteria, Verrucomicrobia, and Chloroflexi in mineral soil. The relative abundance of archaea associated with soil P cycling was higher in mineral soil than organic soil. These results showed that the distribution of P cycling genes and microbial taxonomic groups significantly differed between the organic and mineral soil. Thus, a depthresolved scale investigation can deeply reveal the functions of the soil microorganisms and the benefits to soil P nutrition management in P-deficient subtropical forests.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/f14081665/s1, Figure S1: Nonmetric multidimensional scaling plots (NMDS) of P cycling genes (a) and soil microbial taxa involved in P cycling (b) in soils. NMDS had a good explanation of variation as stress < 0.1. ANOSIM showed significant differences in variables between organic and mineral soil as p < 0.05. Table S1: Stand characteristics (mean ± SD) in the Pinus massoniana-Lithocarpus glaber coniferous and evergreen broadleaved mixed forest (PLF) and L. glaber-Quercus glauca evergreen broadleaved forest (LGF). DBH represents tree diameter at breast height, H represents tree height. Table S2: Details on the 27 functional genes related to P cycling studied in this study. Table S3: The abundance (RPKM) of species (phyla taxa) for a specific functional (KEGG gene) in organic and mineral soils. Table S4: Pearson correlations between the available soil P concentration and the relative abundance of the P cycling genes selected by Random Forest analysis.

Data Availability Statement:
All data that support the findings of this study are available from the corresponding author (GJW) upon request.