The study was conducted in seven agroecosystems in the Northeast China Transect (NECT), a climatic zone with significant precipitation gradients and a critical component of the International Geosphere-Biosphere Program terrestrial transects (Zhang et al., 1997; Ni and Zhang, 2000). At latitudes ranging from 42° N to 44° N and longitudes from 123° E to 128° E within the NECT, precipitation and aridity demonstrate steep gradients (Ni and Zhang, 2000). Meteorological data that included monthly observations of the precipitation and temperature from 7 meteorological stations in Jilin Province were collected from the National Meteorological Information Center ¹ for the period from 1971 to 2017. The NECT spans a wide large range of precipitation gradients from west to east, with mean annual precipitation (MAP) increasing from 300 mm to 900 mm, which MAP is the main driver of this transect. And the mean annual temperature (MAT) has slight difference among sites shown in Supplementary Table S1. We standardized our site selection by focusing on maize planting areas. These agricultural sites, which have existed for over 50 years, are located in flat areas and are managed conventionally under the same government and according to the same policies. Conventional tillage included fall plowing, spring cultivation, planting. The depth of agricultural tillage (commonly referred to as plow sole) is about 20 cm (Yan et al., 2017). Fertilizer for maize was applied with 200 kg N ha⁻¹, 70 kg P ha⁻¹ and 90 kg K ha⁻¹ as the base fertilizer at planting. The weeds were controlled using broad-spectrum herbicides and supplemented as required with manual hoeing (Guo et al., 2020). The soil types are classified as salt-alkali, light chernozem, chernozem, black soil, and dark brown soil, based on the Chinese soil classification system (Xiong and Li, 1987). The salt-alkali and light chernozem profile is composed of a humus layer and a compact gray-white calcium carbonate deposition layer. The organic matter content of the surface layer is 1.5 to 4.5%. There are local alkalinization and salinization, and the soil is weakly alkaline. or alkaline. The cross-section structure of chernozem has obvious humus accumulation and calcification process, the organic matter content is high (surface content is 3.5 to 12%), and the soil is neutral. The topsoil of dark brown soil and black soil is rich in humus (5 to 10%, even up to 17%), and the profile is slightly acidic (Ni et al., 1999). The basic properties of the study sites are listed in Supplementary Table S1.

The seven sites, Changling (CL), Sijianfang (SJF), Nongan (NOA), Changchun (CC), Dongliao (DL), Longwan (LW) and Baihe (BH), are positioned along a precipitation gradient. Three plots (approximately 50 × 50 m) situated 50–100 m from each other were randomly selected at each site. After removing the litter and rocks, the soil was sampled by horizon to a maximum depth of 50 cm. Samples were taken discretely at five depth intervals (0–10 cm, 10–20 cm, 20–30 cm, 30–40 cm and 40–50 cm). Three cores of soil, each with a diameter of 5.5 cm, were randomly collected and gently mixed to form one composite sample for each layer of each plot. In total, 105 soil samples (7 sites × 5 layers × 3 plots) were collected prior to maize harvest in September 2018 to minimize the effects of agricultural management. And all sampling was done within 2 days of sunny weather to ensure consistent weather and temperature.

Each sample was separated into two portions. The first fresh soil samples were immediately placed inside coolers containing ice packets for transport, returned to the laboratory, and stored at 4°C. Soil moisture (SM) was determined by oven-drying technique using 10 g of fresh soil samples dried at 105°C for 48 h to a constant weight. The amount of water in the sample can be determined and the moisture content calculated and expressed as a percentage of the dry soil weight (Schmugge et al., 1980). The second portion was air-dried for the analysis of edaphic properties. The soil pH was measured in a soil-water suspension (1:5 soil to water ratio) using a glass electrode with a pH meter (Thermo Fisher Scientific Inc., United States). Soil total carbon (TC) and total nitrogen (TN) concentrations were determined by combustion 25.000 mg of soil subsamples through a 0.149 mm mesh at 1100°C using an automatic elemental analyzer (Flash EA 1112 elemental, Italy) (Zhang et al., 2022). Total phosphorus (TP) concentration of soil was measured by first digesting 0.5 g soil sample with H₂SO₄-HClO₄, making up 100 ml of the digested solution, and then the TP concentration was measured using a continuous flow analyzer (SKALAR SAN++, Netherlands; Wang et al., 2022).

The potential activities of five hydrolytic enzymes, including two C-acquiring enzymes, β-glucosidase (BG) and β-1,4-cellobiosidase (CB), two N-acquiring enzymes, β-N-acetylglucosaminidase (NAG) and Leucine aminopeptidase (LAP), one P-acquiring enzyme, phosphatase (AP), were measured at a temperature of 25°C based on the methods described by Saiya-Cork et al. (2002) and German et al. (2011). The types, functions and substrates of soil hydrolases are detailed in Table 1. The mass of 1.5 g of fresh soil was accurately recorded and 125 ml of 50 mM acetate buffer was added. We used glacial acetic acid to adjust acetate buffer pH to the pH of this study site. The samples were homogenized by blending on highest speed for 2 min to make a slurry with a magnetic stirrer. Blank wells received 250 μl acetate buffer. Standard wells received 50 μl of standard substrate and 200 μl acetate buffer. Substrate wells received 50 μl of 200 μM enzyme substrate solution in 200 μl of acetate buffer. Sample wells received 50 μl of acetate buffer and 200 μl of soil sample suspension. Quench Control wells received 50 μl of the standard substrate (10 μM 4-methylumvelliferyl or 7-amino-4-methylcoumarin) plus 200 μl of soil sample suspension. Assay wells received 50 μl of 200 μM enzyme substrate solution and 200 μl of soil sample suspension. The pre-experimental results exhibited that the activity was the highest in the dark after 4 h of incubation at 25°C. Fluorescence plates were read at 365 nm excitation and 450 nm emission using a microplate reader (BioTek Synergy HTX, Winooski, Vermont, United States). Eight analytical replicates were conducted per sample. Enzymatic activity was enumerated as nanomoles of substrate released per gram of dry soil and incubation time (nmol h⁻¹ g⁻¹).

The soil nutrient ratios (C:N, C:P and N:P ratios) were expressed by soil elements concentration ratio. The stoichiometry of soil enzymes, such as enzyme C:N ratio, enzyme C:P ratio and enzyme N:P ratio were expressed as ln BG: ln (NAG + LAP), ln BG: ln AP and ln (NAG + LAP): ln AP, respectively. Resource limitation was measured by vector analysis of enzymatic stoichiometry (length and angle) (Moorhead et al., 2013) according to the relative activities of enzymes that acquire C, N and P. The rationale for using the relative activities of C-, N-, and P-acquiring enzymes to interpret the relative limitation of C, N, and P for microorganisms is based on ecosystem stoichiometry and metabolic theory (Allison et al., 2010). Relatively longer vector lengths express a greater C limitation, while vector angles <45° and > 45° indicate the relative degrees of N and P limitation, respectively (Moorhead et al., 2016). The vector length (unitless) and angle (°) are calculated as follows:

The Shapiro–Wilk normality test was performed using the R version 3.5.3 ‘stats’ package (R Development Core Team, 2019). We examined the relationships between precipitation and different soil properties, enzyme activities, and C:N, C:P, and N:P ratios using linear regressions. C-acquiring enzymes were calculated as BG + CB, N-acquiring enzymes were calculated as NAG + LAP, and P-acquiring enzyme was represented by AP. Redundancy analysis (RDA) was performed to analyze the relationships between soil enzyme activities and selected soil properties for 105 samples by the packages “vegan” (Oksanen et al., 2019), “adespatial” (Dray et al., 2018). We used the variation inflation factor (VIF) (“vif.cca” function) to assess multicollinearity among soil properties (VIF < 10). The significance of RDA correlations was tested Monte Carlo permutation test.

Similarly, we explored the relationships between precipitation and vector length and angle for each soil profile. The analysis of covariance (ANCOVA) was used to detect the differences in the slope of regression lines between the angle and length of vector and precipitation among soil profiles by SPSS software version 26.0 (IBM Corporation, Armonk, New York, United Stated). The angle and length of vector were dependent variables, precipitation was covariates, and soil profiles were a grouping factor with five levels. A significant soil profiles × covariate interaction would indicate that the slope of the precipitation–vector relationship differs among soil profiles. When the interaction was significant, the homogeneity of slopes assumption of ANCOVA was rejected (Poorter et al., 2008). Differences at the p < 0.05 level were considered statistically significant. To determine the relative influences of environmental variables on the ecoenzymatic vectors in the five soil layers, we used the rdacca.hp. function to calculate the independent adjusted R-squared (R²) and the percentage of each explanatory variable in the RDA, running 999 permutations in the ‘rdacca.hp’ R package (Lai et al., 2022). This function performs variation partitioning and hierarchical partitioning to calculate the unique, shared, and independent contributions of each predictor to the explained variation (R² and adjusted R²) in canonical analysis (RDA) by applying the hierarchy algorithm of Chevan and Sutherland (1991).

The physicochemical properties and soil nutrient ratios displayed various patterns in topsoil (0–20 cm) and subsoil (20–50 cm) along the precipitation gradient (Figure 1). The soil pH decreased significantly from 9.36 to 4.46 with increasing precipitation. In contrast, soil moisture and nutrients (TC, TN, and TP) were positively correlated with precipitation, but this strong correlation was not observed in the subsoil for TC. Interestingly, for these components of soil nutrient stoichiometry components (TC:TN, TC:TP, and TN:TP), significant positive relationships with precipitation were observed in the subsoil but not in the topsoil. The soil properties and the variation of soil properties with soil depth at seven sites were shown in Table 2 and Supplementary Figure S1, respectively.

The soil enzyme patterns were dissimilar to soil nutrient patterns along the precipitation gradient (Figure 2). C-acquiring enzyme activities were unaffected by precipitation. The N-acquiring enzyme activities decreased across the precipitation gradient in all soil profiles, while P-acquiring enzyme activities were positively correlated with precipitation at a depth of 0–30 cm but not at that of 30–50 cm. The C:P enzyme stoichiometry pattern was the opposite of that of the P-acquiring enzyme. With increasing precipitation, the enzymatic C:N and N:P ratios demonstrated increasing and decreasing trends, respectively; however, both stoichiometries responded to precipitation more drastically in the subsoil than in the topsoil. The variation of soil enzyme activities and stoichiometry with soil depth at seven sites was shown in Supplementary Figure S2. The major environmental variables shaping the soil enzyme pattern were determined through RDA (Figure 3). According to the RDA results, approximately 33.58% of the variance in total soil enzyme activity can be explained by environmental variables. The first two RDA axes explained 30.04% of the relationships between enzymes and the environment (24.79 and 5.25% for the first and second axes, respectively). According to the model selection results, soil pH and moisture were the most important factors affecting the enzyme activity.

The nutrient limitation of soil microbes was calculated using enzyme stoichiometry and quantified by vector angle and length (Figure 4). Vector angles ranging from 36.1° to 80.5° were positively related to the precipitation, and they displayed a more intense correlation with precipitation in the subsoil than in the topsoil (Figure 4A; Supplementary Table S2). We also subdivided the driving factors for nutrient limitation according to the soil layer (Figure 5). We observed that pH had the largest independent effects on vector angle in topsoil, explaining 69.41 and 44.78% of the variance in the depth ranges of 0–10 cm and 10–20 cm, respectively (Figure 5A). From the topsoil to the subsoil, the soil physicochemical properties (soil pH and moisture) variance explanation ranged from 85.64 to 16.79% (Figure 5A). In contrast, the interpretations of soil nutrients and their stoichiometry ranged from 14.36 to 83.21% (Figure 5A). Vector length was positively correlated with precipitation only in the subsoil (Figure 4B). The interpretation of soil nutrient stoichiometry with respect to the vector length grew stronger with depth (Figure 5B). Furthermore, it was not found that the temperature would have an effect on our results by redundancy analysis and hierarchical partitioning analysis.

Soil moisture, total carbon, total nitrogen, total phosphorus increased as the MAP increased along the NECT, but the soil pH decreased as the MAP increased. These results are consistent with the results at the regional scale (Wang et al., 2002; Cui et al., 2019a). Due to the steep moisture gradient, the virgin vegetation changes gradually from mixed deciduous broad-leaved forest to meadow steppe along the transect (Prentice et al., 2011). The forest can improve the availability of soil nutrients by nitrogen fixation and phosphorus strategies (Nasto et al., 2014, 2017). In addition, soil properties are related to the interaction of precipitation and soil parent materials (Huston, 2012). The high pH of the soil in west of the transect is due to the fact that the soil is salt-alkali and the evaporation is three to four times higher than the precipitation (Wang and Ba, 2008; Yang et al., 2021). Therefore, less precipitation, high surface evaporation and parent salinity are the major reasons for high soil pH. Understanding soil properties at the regional scale is fundamental to predicting climate changes at the global scale.

We did not observe a significant relationship between C-acquiring enzymes and precipitation (Figure 2), and surprisingly, the vector length increased with precipitation, indicating that microbial C limitation was gradually exacerbated along the precipitation gradient in 0–50 cm. Moreover, N-acquiring enzyme activity decreased and P-acquiring enzyme activity increased with increasing precipitation (Figure 2), which alters the stoichiometry of enzymes leading to C limitation. We attributed this pattern of microbial C limitation to changes in soil nutrients and their stoichiometry across precipitation regions (Figure 5B). Microbial resource limitation is relative (Chen et al., 2018); thus, microbial C limitation can depend on the combined effect of other nutrients.

We detected a decrease in the enzymatic C:P and N:P acquisition ratios with increasing precipitation (Figure 2). Microorganisms can shift from secreting more C-acquiring enzymes and N-acquiring enzymes to secreting more P-acquiring enzymes, which implies an increased demand for P. Moreover, the vector angle became greater than 45° with increasing precipitation (Figure 4A). Our results notably indicated that microbial P limitation was gradually exacerbated along the precipitation gradient. A perpendicular projection of P-acquiring enzyme onto the line overlaying the environmental variable arrow revealed that P-acquiring enzyme was influenced by soil pH and moisture (Figure 3). P is predominately derived from mineral aerosol deposition and weathering, and it is highly dependent on pH (Thingstad et al., 2005; Mahowald et al., 2008). High soil moisture and low soil pH increase P-acquiring enzyme activity, which facilitates P solubilization (Collavino et al., 2010; Xu Z. et al., 2020). In addition, soil pH and moisture were the main factor determining the pattern of soil enzymes along the precipitation gradient (Figure 3). Therefore, soil pH and moisture play an essential role in microbial P limitation throughout soil profiles.

Furthermore, soil nutrients and their stoichiometry and soil physicochemical properties are equally important for microbial phosphorus limitation in 0–50 cm (Figure 5A). More C and N might lead to microbial P limitation (Cui et al., 2020). This observation may be explained for two reasons. On the one hand, if the C and N accumulated in the soil reduces the relative concentration of soil P, the microorganisms may be limited by P due to the elemental stoichiometric balance of the microorganisms (Sinsabaugh et al., 2008, 2009; Chen et al., 2017). Thus, despite the increase in soil P, microbial P limitation remains because of the slow growth rate of P relative to that of C and N. On the other hand, phosphatases have relatively high C and N concentrations, and microorganisms require sufficient C and N to construct extracellular phosphatases (Treseder and Vitousek, 2001; Houlton et al., 2008). Phosphatase may be highly responsive to changes in soil elements, a critical strategy for organism adaption to human modifications of the global biogeochemical cycle that can delay the onset of P limitation (Güsewell and Freeman, 2005; Marklein and Houlton, 2012). These findings demonstrate that precipitation regulates soil physicochemical properties and nutrient status (Smith, 2011; Ma et al., 2015), which consequently greatly enhances microbial C and P limitation. Our current results take a new step and are valuable for understanding microbial nutrient limitation along the precipitation gradient in agroecosystems.

Our findings are in support of our hypothesis that microbial nutrient limitations were stronger in the subsoil than in the topsoil along the precipitation gradient. Also, we show that the driving factors of microbial nutrient limitation are different in topsoil and subsoil. Soil properties and enzyme activities can change with soil profiles (Stone et al., 2014; Li et al., 2020; Luan et al., 2020), and we observed wide differences in soil and enzyme stoichiometry between the topsoil (0–20 cm) and subsoil (20–50 cm). Our study found that enzymatic C:P acquisition ratios were negatively correlated with precipitation only in the topsoil, which is attributable to the phosphatase distribution pattern. Enzymatic N:P acquisition ratios decreased, but enzymatic C:N acquisition ratio increased along the precipitation gradient, and changes substantially in the subsoil. The vector angles visualized that microbial P limitation increased and responded more strongly to precipitation in the subsoil than in the topsoil. In topsoil, soil pH and moisture are the major drivers of microbial P limitation (Figure 5A). Soil pH governs phosphatase spatial variations and is positively correlated with enzymatic C:P and N:P acquisition ratios (Xu et al., 2017). It is interesting to note that, compared to the topsoil, soil nutrient availability and their stoichiometry are important factors in determining microbial P limitation in subsoil (Figure 5A). Soil substrate availability greatly influences enzyme activity patterns (Stock et al., 2019), and low resources (C and N) and O₂ concentrations can limit microbial metabolism (Fierer, 2017; Gu et al., 2017). Zhou et al. (2020) also found that the soil C:N ratio explains a large proportion of the variation in microbial nutrient acquisition. Emerging evidence suggests that subsoils, exhibit more varied organic matter sources, microbial communities, and substrate availability compared to those of the topsoil (Rumpel et al., 2012; Chen et al., 2019). Differences in soil properties across soil depths can lead to greater environmental filtering or competition for nutrients among microorganisms than in the topsoil (Bahram et al., 2015; Schlatter et al., 2020). Especially in agroecosystems, this difference between topsoil and subsoil may be attributed to the agricultural tillage depth, which is typically 0–20 cm. Cropland activities contribute to the homogenization of microbial communities (Rodrigues et al., 2013), and alter microbial resource limitations (Cui et al., 2021). On the contrary, the original soil conditions were preserved in the subsoil, resulting in microbes secreting soil enzymes that are less responsive to precipitation in the topsoil than that in the subsoil and attenuate the reaction of microbe nutrients limitation to precipitation in the topsoil.

Ecosystems cannot be limited by carbon or nutrients alone (Soong et al., 2020), and most microbial communities are co-limited by energy and key nutrients (Moorhead et al., 2016). Another marked trend is that the vector length increases in the subsoil with precipitation. This result indicates that microbial C limitation increased in the subsoil along the precipitation gradient. The detected relationships might result from that soil C:N and C:P ratios progressively decline across the precipitation region in the subsoil (Figure 1), which can lead to insufficient soil C (Chen et al., 2018). In contrast, there is no relationship between microbial C limitation and precipitation in topsoil. This might occur because the effect of precipitation is masked by cropland activity, causing soil nutrient stoichiometry in topsoil that exhibited no response to precipitation. Moreover, soil nutrients and their stoichiometry increase in their interpretation of C limitation with depth (Figure 4). Therefore, soil C:N:P stoichiometry tremendously affected microbial groups involved in major biogeochemical processes (Luo et al., 2020). Additionally, balanced nutrient stoichiometry is essential for maintaining microbial and elemental homeostasis (Cui et al., 2019b).

Further, P limitation changes rapidly in the subsoil, which contributes to C limitation. Nutrients are interdependent and they have a complex coupling, so changing in one nutrient cycle type can alter the availability of another nutrient (Marklein and Houlton, 2012). The role of P limitation in the storage capacity of the C ecosystem has previously been substantially underestimated (Peng and Thomas, 2010). P affects the C-storage capacity both directly and indirectly by affecting the plant growth capacity and limiting biological activities for N fixation (Peñuelas et al., 2012). A model predicted less efficient microbial growth in P-limited soils, indicating that P availability may affect C cycling (Waring et al., 2014). These nutrient interdependencies confirm that resource limitation coupling is complex, and excess precipitation leads to depletion of P and C and induces imbalances in agroecosystems. Our results highlight the contrasting mechanisms governing microbial nutrient limitation in topsoil and subsoil. Microbial P limitation was controlled by the soil pH and moisture in the topsoil, while microbial C and P limitation in the subsoil were regulated by the soil nutrients and their stoichiometry (Figure 6). The findings suggest that ecosystem functions and microbial nutrient limitation cannot be comprehensively assessed from topsoil observations alone.

Ecological stoichiometry is commonly applied to explore relationships between brown food webs such as soils and microbes (Bell et al., 2014). When nutrient-limited, microorganisms reduce their uptake of C for nutrients and excrete excess C through enzymes and respiration, so that the ability to assimilate C into biomass tends to decrease (Manzoni et al., 2017). Nutrient-deprived foods are stoichiometrically unbalanced, due to differences in elemental concentrations in food and consumers reflecting stoichiometric mismatches that limit consumer growth and development, thereby affecting the metabolism and function of the brown food web (Filipiak, 2016). Further studies on soil profiles, we found that global change-induced stoichiometric mismatch in agricultural subsoil is more severe. Our results are an important step forward in better understanding soil microbial metabolism across the whole profile of agroecosystems under changing precipitation patterns. Future research could aim to study the metabolic capacity of brown food webs under precipitation changes through metagenomic next-generation sequencing to better understand how precipitation affects the function of whole-profile brown food webs in agroecosystems.