Effects of nitrogen fertilization and bioenergy crop species on central tendency and spatial heterogeneity of soil glycosidase activities

Extracellular glycosidases in soil, produced by microorganisms, act as major agents for decomposing labile soil organic carbon (e.g., cellulose). Soil extracellular glycosidases are significantly affected by nitrogen (N) fertilization but fertilization effects on spatial distributions of soil glycosidases have not been well addressed. Whether the effects of N fertilization vary with bioenergy crop species also remains unclear. Based on a 3-year fertilization experiment in Middle Tennessee, USA, a total of 288 soil samples in topsoil (0–15 cm) were collected from two 15 m2 plots under three fertilization treatments in switchgrass (SG: Panicum virgatum L.) and gamagrass (GG: Tripsacum dactyloides L.) using a spatially explicit design. Four glycosidases, α-glucosidase (AG), β-glucosidase (BG), β-xylosidase (BX), cellobiohydrolase (CBH), and their sum associated with C acquisition (Cacq) were quantified. The three fertilization treatments were no N input (NN), low N input (LN: 84 kg N ha−1 year−1 in urea) and high N input (HN: 168 kg N ha−1 year−1 in urea). The descriptive and geostatistical approaches were used to evaluate their central tendency and spatial heterogeneity. Results showed significant interactive effects of N fertilization and crop type on BX such that LN and HN significantly enhanced BX by 14% and 44% in SG, respectively. The significant effect of crop type was identified and glycosidase activities were 15–39% higher in GG than those in SG except AG. Within-plot variances of glycosidases appeared higher in SG than GG but little differed with N fertilization due to large plot-plot variation. Spatial patterns were generally more evident in LN or HN plots than NN plots for BG in SG and CBH in GG. This study suggested that N fertilization elevated central tendency and spatial heterogeneity of glycosidase activities in surficial soil horizons and these effects however varied with crop and enzyme types. Future studies need to focus on specific enzyme in certain bioenergy cropland soil when N fertilization effect is evaluated.

Bioenergy crops have the potential to reduce fossil fuel consumption 1 and the energy crops such as switchgrass (SG: Panicum virgatum L.) and gamagrass (GG: Tripsacum dactyloides L.) are key for supplying biofuel plant biomass 2,3 . Bioenergy crop yields are enhanced by nitrogen (N) fertilizers 4 , and previous studies frequently focused on aboveground crop yield and less so on belowground features. Nevertheless, N fertilization substantially alters microbial community composition and structure in soil 5 and consequently impacts soil extracellular enzymes that microbes produced and excreted to the environment 6 . As important proxies to soil health and management 7,8 , soil extracellular enzyme activities mirror soil community's metabolic requirements and available nutrients 9 . Given the N fertilizer overuse worldwide, investigation of spatial pattern of soil microbial functions such as extracellular enzymes is imperative. Knowledge of spatiotemporal variations of soil extracellular Scientific Reports | (2020) 10:19681 | https://doi.org/10.1038/s41598-020-76837-1 www.nature.com/scientificreports/ soil microbial biomass and the potentially positive relationship of microbial biomass and extracellular enzymes, our second hypothesis was that in soils that have never been fertilized for years, N fertilization would restructure spatial heterogeneity of AG, BG, BX, CBH and C acq . Last, we hypothesized that N fertilization effects on central tendency and spatial heterogeneity varied with enzyme type (e.g. AG, BG, BX and CBH) due to the unique characteristics of each individual enzyme.

Materials and methods
Study site description and experimental design. In 2011, a bioenergy crop field fertilization experiment was established located at the Tennessee State University (TSU) Main Campus Agriculture Research and Education Center (AREC) in Nashville, TN, USA. Prior to the croplands, the land was mowed grassland for several decades with no amendment of fertilizers. The experimental site marks a warm humid temperate climate with an average annual temperature of 15.1 °C, and total annual precipitation of 1200 mm 42 . The crop type and N fertilization treatments were included in a randomized block design 37,39,43,44 . The two crop types were Alamo SG (Panicum virgatum L.) and GG (Tripsacum dactyloides L.). The three N levels included no N fertilizer input (NN), low N fertilizer input (LN: 84 kg N ha −1 year −1 as urea), and high N fertilizer input (HN: 168 kg N ha −1 year −1 as urea), and each treatment had four replicated plots with a dimension of 3 m × 6 m. The low N fertilization rate was determined as the optimum N rate to maximize cellulosic ethanol production in established northern latitude grasslands 45 . The high N rate doubled the low rate in order to create appreciable gap and detectable effect between the two levels. The fertilizer was manually applied in June or July each year after cutting the grass. The soil series for the plots is Armour silt loam soil (fine-silty, mixed, thermic Ultic Hapludalfs) with acidic soil pH (i.e., 5.97) and intermediate organic matter content of 2.4% 39,46 .
Soil collection and laboratory assay. In this study, soil subsamples were adopted from soil collections based on our former study 39 . Here a brief introduction was presented regarding the former soil collection. On June 6th, 2015, soil samples (0-15 cm) were collected from 12 plots (2 crop × 3 N × 2 replicates). In each plot, soil sampling location was determined in a spatially explicit way accounting for randomization in both sampling direction and distance ( Fig. 1). Given this sampling design, the unique x, y coordinates were assigned to each sample. Twenty-four cores were collected from each plot yielding 288 soil cores in 12 plots. Soil samples stored in coolers filled with ice packs were immediately transported to TSU lab and subsequently stored at 4 °C until chemical analysis; and subsamples were stored in − 20 °C for enzymatic assay. The visible roots and rocks were removed from soil cores by passing through a 2 mm soil sieve prior to chemical analysis and enzymatic assay. For each soil sample, soil gravimetric moisture content was determined by oven drying subsamples for 24 h at 105 °C. Water extractable soil pH was measured given soil: water = 1:5. Four glycosidase activities were quantified by soil fluorimetric enzymatic assay methods in each core. Briefly, soil samples for each plot were assayed for α-glucosidase (AG), β-glucosidase (BG), β-xylosidase (BX), and cellobiohydrolase (CBH) using 4-methylumbelliferyl (MUB)-α-d-glucopyranoside, MUB-β-d-glucopyranoside, MUBβ-d-xylopyranoside, and MUB-β-d-cellobioside with concentrations of 200 mmol/L as substrates, respectively, following published protocols 47,48 . Sample suspensions were prepared by placing 1.0 g soil in a 125 ml Nalgene bottle. Acetate buffer (50 mM, pH 5) was added to the bottle and the resulting suspension was homogenized using a Brinkmann Polytron for approximately 1 min. Additional buffer was added to the bottle to bring the final suspension volume to 125 ml.
The plates were placed in an Echotherm incubator at 20 °C, for 18-24 h given enzyme type. The assay of glycosidase activities was conducted on black 96-well microtiter plates. The assay design included reference standards (eight wells) and quench controls (eight wells per sample) added to each plate. The 10 µM MUB was used as the reference standard for AG, BG, BX and CBH. Quench control wells contained 200 µl of sample suspension and 50 µl of the reference standard. The assay was incubated at 20 °C. The reactions were terminated by adding 10 µl of 1.0 M NaOH to each well. Fluorescence was measured using a Molecular Devices (Multi-Mode Microplate Reader, FilterMaxF5) with excitation set to 365 nm and emission set to 460 nm. All enzyme activities were calculated as µmol activity h −1 g soil −1 . A total of 1152 enzymatic activity data were collected for 288 soil samples and 4 enzymes. Laboratory tests were conducted and specific protocols were optimized to secure sufficient soil mixing. As a result, the variation of each measurement (i.e., coefficient of variation) in multiple tests ranged from 2 to 8% based on our protocol. Statistical analysis. We use both descriptive and geospatial analytical methods to illustrate the central tendency and spatial heterogeneity of enzymes assayed. Mean, frequency distribution, plot-level variance and with-plot coefficient of variation (CV) were estimated to describe central tendencies and variations for enzyme activities in each plot. The two-way ANOVA was used to test whether N fertilization, crop species and their interaction significantly affected each enzyme. To avoid the pseudo-replication impacts, the plot means were used in the two-way ANOVA test. The statistically significant level was set at P < 0.05.
Cochran's C test was performed to test the assumption of variance homogeneity. The test statistic is a ratio that relates the largest empirical variance of a particular treatment to the sum of the variances of the remaining treatments. The theoretical distribution with the corresponding critical values can be specified. Soil properties that exhibited non-normal distributions were log-transformed to better conform to the normality assumption of the Cochran's C test 49,50 .
The sample size required in a research plot can be determined quantitatively under given desired sampling error 51 . That is, under a desired sampling error, the sample sizes derived can be used to evaluate the plot-level variations between different research plots. In this study, the sample size requirement ( N ) in each plot was derived given specified relative error ( γ ), which was defined as the ratio of error term ( t 0.975 × s www.nature.com/scientificreports/ with a range of 0-100% (Eqs. 1-3). To evaluate how sample size requirement varied with N fertilization or crop types at certain relative error, the average of sample size ( N ) in two plots was derived and plotted. Under a relative error of 10%, the sample sizes were also derived from each plot and compared between different plots. For comparison, the higher sample size the greater plot-level variation under the same relative error.
where CI, X, s, n, N, CV, and γ denote confidence interval, plot means, plot standard deviation, sample number (n = 24), coefficient of variation, sample size requirement and relative error, respectively. t 0.975 = 1.96. The log-transformed sample size requirement ( N) has a negative linear relationship (i.e. slope = 2) with the logtransformed relative error (γ ). Figure 1. Illustration of an efficient clustered random sampling design within a plot (2.75 m × 5.5 m). The plot was divided into eight subplots (grey zone) and there was a centroid (dark solid circle) in each subplot (1.375 m × 1.375 m), where three soil sampling points (+) were determined from random directions and distances from a centroid in each sampling region (grey area). The extent of an interpolation map was thus determined by the minimum and maximum values at horizontal and vertical axes, and each map can attain its extent less than or equivalent to a plot area. www.nature.com/scientificreports/ Geostatistical analysis. Three different geostatistical tools were applied to describe the spatial structure of soil exoenzyme activities within and among plots. The methods were briefly described below and more details could be found in Li 52 . First, trend surface analysis (TSA) is the most common regionalized model in which all sample points fit a model that accounts for the linear and non-linear variation of an attribute. Relationships between soil properties and x and y coordinates of their measurement location within the sampling plots are estimated with the trend surface model (Eq. 4): The presence of a trend in the data was determined by the significance of any of the parameters β 1 to β 5 , while the β 0 was the intercept 53,54 . Linear gradients in x or y directions were indicated by the significance of β 1 or β 2 . A significant β 3 indicated a significant diagonal trend across a plot. Significant β 4 and β 5 parameters indicated a more complex, nonlinear spatial structure such as substantial humps or depressions. Trend surface regressions were estimated using R program 55 . Model parameters were determined to be significant at a level of P < 0.05.
Second, residuals from the trend surface regressions were saved for subsequent spatial analysis using a Moran's I index 56 . The Moran's I analysis [57][58][59] was used to quantify the degree of spatial autocorrelation that were present in each plot. The resulting local Moran's I statistics is in the range from -1 to 1 With a positive Moran's I values indicating similar values (either high or low) are spatially clustered, and a negative Moran's I values indicating neighboring values are dissimilar. No spatial autocorrelation or spatial randomness was reached with a Moran's I value of 0. Given that the observed Moran's I value is beyond the projected 95% confidence interval at a certain distance, this is identified as a significant autocorrelation. In this study, correlograms were produced for soil variables in all plots given a range of 0-5.5 m with 0.25 m incremental interval.
Third, an ordinary kriging method was usually used to produce maps which offered direct and visual assessments from which to compare the spatial distributions of the soil properties among the plots 60 . The ordinary kriging method required a large sample size (i.e., a few hundred or more) in order to achieve reliable interpretation maps 60 . Due to the fine-scale sampling region (1.375 × 1.375 m) and a relatively small sample size per plot (n = 24), inverse distance weighting (IDW) interpolation was used in this study. The IDW maps were formerly used to distinguish the effects of different land uses on spatial distributions of soil biogeochemical features in South Carolina, USA 52 . Briefly, the weights for each observation were inversely proportional to the power of its distance from the location being estimated. Exponents between 1 and 3 we typically used for IDW. Tests with different IDW exponents indicated that 2 was optimal with data collected in this study, as the exponent of 2.0 showed the best fit between estimated values and actual data in cross-validation tests 61 . ArcGIS 10.6 (ESRI, USA) was used to generate the IDW maps and perform cross-validations.

Results
Central tendencies and within-plot variances. There were significant main and interactive effects of N fertilization and crop species on BX activity (Table 1), and post hoc tests showed that relative to unfertilized treatment (NN), N fertilization treatments significantly escalated BX activity by 14% (LN) and 44% (HN) in SG ( Table 2). There were also significant effects of crop species on the activities of BG and BX (Table 1) and the effect of crop species on CBH was marginally significant (P = 0.056; Table 1). Relative to SG, GG was higher by 15%, 31%, 32%, and 39% in BX, BG, C acq and CBH, respectively (Table 2). There were no significant effects of N fertilization or interaction of N fertilization and crop species on activities of AG, BG, CBH or C acq ( Table 1).
The frequency diagrams of four glycosidase activities and C acq showed nearly normal distributions under all treatments in two croplands except BG under HN in SG (Fig. 2). The frequency distributions contrasted substantially among different N fertilization treatments for SG, whereas, the frequency distributions of different N fertilization treatments showed relatively similar ranges for GG (Fig. 2). Remarkably, for BX in SG, NN appeared to have the highest frequency in lower values, whereas, HN had the highest frequency in higher values. This frequency distributions between NN and HN was consistent with a significantly escalated BX activity in HN compared to NN (Table 2). On the other hand, the Cochran's C tests showed that N fertilization little changed plot-level variation for most glycosidase activities in both bioenergy croplands, except that HN induced the highest plot-level variance for AG in SG ( Table 3).
The within-plot CVs of four glycosidases and C acq ranged from 15 to 47% in all treatments (Fig. 3). The CVs of AG, CBH and C acq were higher in SG than those in GG (Fig. 3). In 12 plots, the number of plots with CVs larger than 40% for AG, BG, BX, CHB and C acq , were 2, 0, 0, 4 and 2 in SG, and 0, 1, 0, 1 and 0 in GG, respectively; Table 1. P-values of two-way ANOVA tests for the main and interactive effects of N fertilization and crop species on AG, BG, BX, CBH, and C acq (µmol g −1 soil h −1 ). Bold numbers denote significant treatment effects at P < 0.05. AG: α-glucosidase; BG: β-glucosidase; BX: β-xylosidase; CBH: cellobiohydrolase; C acq : the sum of AG, BG, BX and CBH. www.nature.com/scientificreports/ Accordingly, the number of plots with CVs less than 20%, were 0, 1, 0, 0 and 0 in SG, and 0, 4, 3, 1 and 0 in GG, respectively. The sample size requirement (SSR) for all enzymes was generally higher for NN than LN or HN in both croplands, except CBH in GG (Fig. 4). The plotted lines of SSR against relative sampling error departed from one to another in SG in a much wider extent than those in GG for all glycosidases except BX (Fig. 4). In general, a larger number of samples were required in SG than that in GG for all glycosidases under the same desired relative error. Given the same desired sampling error of 10%, the sample size required were always larger in SG than that in GG; A total of 123 samples were required for CBH under NN in SG and only 14 samples were required for BG under HN in GG to achieve the desired error of 10% in both plots (Table 4).
Surface trend, autocorrelation and spatial map. Trend surface analysis results showed only a few significant linear or nonlinear trends in each plot, and about half of plots showed no significant linear or nonlinear trends (Table 5; Table S1). Given the detected significant surface trends, there was a contrasting pattern of surface trend with N fertilization between SG and GG (Table 5). In SG, there were no significant linear or nonlinear trends in any of HN plots for all enzymes. Relative to NN plots, there were more significant linear or nonlinear surface trends of AG and BG in LN plots, whereas there were comparable number of surface trends of BX, CBH, and C acq between NN and LN plots (Table 5). In GG, there were no significant linear or nonlinear trends in any of NN plots for all enzymes except AG, and there were no any significant surface trends in any plot for BX or in any of NN and LN plots for CBH. Relative to NN plots, there were more significant linear or nonlinear surface trends of BG, CBH, and C acq in LN or HN plots; whereas there were comparable or larger number of surface trends of AG in NN than those in LN or HN plots (Table 5). Under the same treatment, the number of significant linear or nonlinear trends varied between the two replicated plots, and in all cases except AG in LN plot in GG, there was significant surface trends in one plot, but none in another plot (Table S1).
The number and distance of significant spatial autocorrelations varied with N fertilization treatments, bioenergy croplands, and among variables ( Table 6). The number of significant spatial autocorrelations in SG was identified more frequently than that in GG for AG, BX, CBH, and C acq across three N fertilization treatments. Compared to NN, fertilized treatments (LN and HN) possessed a higher number of significant spatial autocorrelations for BG, BX, CBH in SG, and for AG, BX, CBH in GG. The distance of significant spatial autocorrelation appeared to be positive or negative in any plot for any enzyme studied. The distances in which the significant spatial autocorrelations appeared ranged from − 5.25 to 5 m across all enzymes. Relative to other enzymes, AG showed significant spatial autocorrelations in more diverse distances in almost all plots except P1 in NN and LN plots in GG (Fig. 5).
With the same scale for each enzyme in two crops, the IDW maps of all enzymes exhibited higher activities (e.g., darker color) in GG than those in SG, and this was true in unfertilized and fertilized plots (i.e., NN, LN and HH) (Figs. 6 and 7). In SG, all IDW maps exhibited low to high activities (e.g., shallower and gradually darker colors) from NN plots, through LN, to HN plots (Fig. 6). Also, the contrast between dark and shallow color regimes in a plot were more pronounced in LN or HN relative to NN, and this was particularly evident for BX (Fig. 6). In GG, the IDW maps exhibited different patterns from those in SG. Large variations of color regime were identified among different enzymes with darker color for BG, CBH and C acq (Fig. 7). The color regimes were comparable and evident among all plots for each enzyme (Fig. 7).

Discussion
N fertilization elevated BX activity in switchgrass cropland soils. Our current study identified that only BX activity was significantly affected by N fertilization. This only partially supported our first hypothesis that all studied glycosidases would increase with N fertilization. The responsiveness of BX may lie in several possible mechanisms. First, N fertilization could increase production and excretion of BX given the elevated relative abundance of Gram-negative bacteria in cropland soils 62,63 because of the close association of BX with Gram-negative bacteria 64 . Second, BX was found out to be significantly correlated with soil C contents but not significantly correlated with microbial biomass 65 . In our plots, N fertilization significantly increased SOC by up to 16% but little changed microbial biomass 39 . The positive response of BX thus was likely associated with the elevated SOC stock under N fertilization. Third, N fertilization stimulated plant growth and subsequent Table 2. Means (± SE) of AG, BG, BX, CBH, and C acq (µmol g −1 soil h −1 ) under three N fertilization treatments (NN, LN and HN) in two bioenergy croplands (SG and GG). SG: switchgrass; GG: gammagrass; NN: No nitrogen fertilizer input; LN: low nitrogen (84 kg N ha −1 year −1 in urea); HN: High nitrogen (168 kg N ha −1 year −1 in urea). In each column, different lowercase letters denote significant difference between fertilization treatments at P < 0.05 (N = 2).

Crop
Fertilization  www.nature.com/scientificreports/ stimulation of microbial activity in our plots 37 , this likely contributed to the elevated BX but it remained unclear why other glycosidases studied showed no positive response to N fertilization. One possible explanation for the lack of response in CBH may lie in the N fertilization effect on the Gram-positive bacteria, fungi and the actinomycetes 63,66,67 , which correlated with CBH 68 . Interestingly, the detectable response of BX suggested that the fundamental attributes of enzymatic reactions could be extrapolated from molecule through community to ecosystem scales 69 . In addition, BX and BG were significantly and CBH was marginally significantly higher (by > 30%) in GG than that in SG. This rejected part of our first hypothesis that there was no significant difference of glycosidase activities between SG and GG. Despite similar characteristics of both SG and GG roots, i.e., massive root volume, the contrasting root chemistry of the two plants may induce different strategy to compete with soil microbes for nutrients 70 . Given the more structurally complex nature of GG root than SG root 37 , this may slow GG root to acquire readily available N (e.g., ammonium or nitrate) and thus result in high nutrient availability to microbes. The strategy of microbial nutrient acquisition may thus shift from the control by nutrient deficiency to nutrient abundance, resulting in less BX production and expression under N fertilization in GG. Besides the crop species, the effect of N fertilization on glycosidases were also reported to co-vary with other factors, such as soil depth and sampling location (rhizosphere vs. bulk soil) 71,72 and soil and ecosystem types 6 .  www.nature.com/scientificreports/  Table 5. The number of significant regression coefficients of trend-surface analysis for AG, BG, BX, CBH, and C acq (µmol g −1 soil h −1 ) under three N fertilization treatments (NN, LN and HN) in two bioenergy croplands (SG and GG). Values represent the sum of significant regression coefficients in two replicated plots under each treatment. The regression coefficients denote parameters β 1 to β 5 in Eq. (4). The significant coefficients of trend-surface analysis for each plot was presented in Table S1. The abbreviations are referred to Tables 1 and 2.  www.nature.com/scientificreports/ Based on the post hoc tests of the interaction term, BX activity was significantly escalated under N fertilization in SG. That is, the positive response of BX to N fertilizer appeared much less in GG than SG. This may lie in several different mechanisms. First, relative to SG, the percentage increases of SOC under N fertilization were much smaller in GG 44 , which could limit the response of BX given the relationship of BX and soil C contents 65 . Second, soil microbial biomass showed no significant change in SG but significantly increased under N fertilization in GG 39 . We speculate that despite the increased microbial biomass, the substrate availability to microbes may have been limited due to no change of SOC under N fertilization 44 , which could result in microbial C limitation and subsequently a constraint on enzyme productions. N fertilization restructured spatial heterogeneity of glycosidases. In support of our second hypothesis, N fertilization resulted in more pronounced spatial heterogeneity of glycosidases in both bioenergy croplands. Despite tremendous efforts committed to avoid nutrient hotspots due to uneven fertilizer application, the elevated spatial heterogeneity were still likely induced by fertilizer application in bioenergy crops because manual spread of N fertilizers would create significant irregularity of nutrient deposit and clusters 39,44 , and consequently favor the formation of hotspots of soil microbial communities 73 . In these hotspots, microbes grew faster and SOC appeared high so that these conditions may eventually create hotspots of microbial functions, such as greater extracellular enzyme activities. Consistent with this speculation, the positive effect of N fertilization on spatial heterogeneity of SOC were also found in the same field experiment 44 . This intricate associations of bulk soil C and glycosidases were also supported by their significant correlation coefficients (P < 0.05; Table 7). Meanwhile, there was no significant correlations between microbial biomass C and any of these glycosidases studied ( Table 7). The correlation between SOC and glycosidase and no correlation between MBC and glycosidase are consistent with the findings revealed in Waldrop, Balser 65 . These suggested that the N fertilization induced heterogeneity of glycosidases may be largely moderated by the microbial substrate availability over space (e.g., SOC), not necessarily by the microbial abundance (e.g., MBC). This strategy was consistent with the recent discovery of the negative correlation between maximum microbial growth rate and soil extracellular hydrolytic enzymes under high resource conditions 74 . With readily available nutrients (e.g., nitrate and ammonium) under N fertilization, the resource acquisition strategists invest heavily in extracellular enzyme production while other microbial groups (e.g., growth strategist and maintenance strategist) either compete for fast growth or limit investment in both enzymes and growth 74 . Noted that a positive effect of N fertilization on spatial heterogeneity of MBC was also found in the same experiment 39 , suggesting a decoupling of spatial distributions of MBC and glycosidases under N fertilizations. Collectively, these results suggested that the change and distribution of glycosidases may primarily be driven by the substrate availability for microbes, not the abundance of microbes themselves. N fertilization effects on central tendency and spatial heterogeneity of glycosidases varied with enzyme type. In support of our third hypothesis, this study showed that the N fertilization effects on central tendency and spatial heterogeneity of glycosidases varied with enzyme type. This result is not surprising due to the generally high indigenous soil heterogeneity. N fertilizer applied in each point in a plot impact glycosidase activities in the specific location via a suite of biogeochemical reactions involving the N-containing molecules, plant root, soil, and microbes; Then, the glycosidase activity could increase or decrease in each location affected by N fertilizer input, leading to re-distribution of enzymes and other soil features (e.g., SOC and MBC), i.e., restructuring spatial heterogeneity in the plot level; As a consequence of spatial restructuring, the plot level mean or central tendency and variation were likely changed as well. The N fertilization-induced changes in spatial heterogeneity of soil enzymes will also likely depend upon the indigenous site condition and the legacy effect that it exerted for years to decades 75 .
In response to N fertilization, both central tendency and spatial heterogeneity of glycosidases contrasted among different enzymes. Besides the unique chemical characteristics of each enzyme, the indigenous variation of each enzyme in a plot may regulate the effect of N fertilization. Glycosidases showed medium CVs (i.e., Table 7. Pearson correlation coefficients between SOC, TN, C/N, MBC, MBN, MBC/MBN, and glycosidase activities (i.e., AG, BG, BX, CBH) under three N fertilization treatments (NN, LN and HN) in two bioenergy croplands (SG and GG). Bold values denote significant correlation coefficients at P < 0.05. www.nature.com/scientificreports/ 15-47%) as they appeared to be higher than that for soil moisture, total pore space, pH, SOC, TN, δ 13 C and δ 15 N 30,44,76 and substantially lower than that for extractable soil Fe and Mn 52 . However, substantially different CVs were evident between enzymes in this study. For instance, BG showed the most extreme CVs of 15-47% while AG showed narrow CVs of 23-29% in GG. The contrasting plot-level variations led to a difference of sample size required up to an order of magnitude (123 vs. 14; Table 4). The same number of soil sampling design in each plot unavoidably induced differential error terms for different enzymes that could have contributed to the insignificant treatment effects for most enzymes. However, BX responded significantly to both N fertilization and crop type suggesting an essentially strong and a suite of high order interaction of this enzyme with plant and soil, mediated by microbial community strategy. Due to the complexity of interactions, the fertilization effects on spatial heterogeneity of BG and CBH were identified. This suggested that the N fertilization effects on central tendency and spatial distribution of enzymes decoupled for the same enzyme.

Conclusions
Our study demonstrated that N fertilization significantly enhanced central tendency and resulted in more pronounced spatial heterogeneity of glycosidase activities in top soil horizons in bioenergy croplands, but the N fertilization effects varied with crop species and enzyme type. N fertilization significantly enhanced BX by 14-44% in SG, and consistently restructured the spatial heterogeneity of BG in SG and CBH in GG. Two bioenergy crops showed contrasting glycosidases' activities (GG > SG) and their plot-level variations (GG < SG). The unique enzyme characteristics and their interactions with plant root and soil, mediated by soil microbial community, possibly explained the enzyme-specific responses under N fertilization. Though N fertilization elevated mean activities and spatial heterogeneity of glycosidases, these effects varied with crop species and enzyme type. Future studies should focus on specific enzyme when evaluating N fertilization effect in certain bioenergy cropland soil.