Soil CO2 emission and soil attributes associated with the microbiota of a sugarcane area in southern Brazil

The spatial structure of soil CO2 emission (FCO2) and soil attributes are affected by different factors in a highly complex way. In this context, this study aimed to characterize the spatial variability patterns of FCO2 and soil physical, chemical, and microbiological attributes in a sugarcane field area after reform activities. The study was conducted in an Oxisol with the measurement of FCO2, soil temperature (Ts), and soil moisture (Ms) in a regular 90 × 90-m grid with 100 sampling points. Soil samples were collected at each sampling point at a depth of 0–0.20 m to determine soil physical (density, macroporosity, and microporosity), particle size (sand, silt, and clay), and chemical attributes (soil organic matter, pH, P, K, Ca, Mg, Al, H + Al, cation exchange capacity, and base saturation). Geostatistical analyses were performed to assess the spatial variability and map soil attributes. Two regions (R1 and R2) with contrasting emission values were identified after mapping FCO2. The abundance of bacterial 16S rRNA, pmoA, and nifH genes, determined by real-time quantitative PCR (qPCR), enzymatic activity (dehydrogenase, urease, cellulase, and amylase), and microbial biomass carbon were determined in R1 and R2. The mean values of FCO2 (2.91 µmol m−2 s−1), Ts (22.6 °C), and Ms (16.9%) over the 28-day period were similar to those observed in studies also conducted under Oxisols in sugarcane areas and conventional soil tillage. The spatial pattern of FCO2 was similar to that of macropores, air-filled pore space, silt content, soil organic matter, and soil carbon decay constant. No significant difference was observed between R1 and R2 for the copy number of bacterial 16S rRNA and nifH genes, but the results of qPCR for the pmoA gene presented differences (p < 0.01) between regions. The region R1, with the highest FCO2 (2.9 to 4.2 µmol m−2 s−1), showed higher enzymatic activity of dehydrogenase (33.02 µg TPF g−1 dry soil 24 h−1), urease (41.15 µg NH4–N g−1 dry soil 3 h−1), amylase (73.84 µg glucose g−1 dry soil 24 h−1), and microbial biomass carbon (41.35 µg C g−1 soil) than R2, which had the lowest emission (1.9 to 2.7 µmol m−2 s−1). In addition, the soil C/N ratio was higher in R2 (15.43) than in R1 (12.18). The spatial pattern of FCO2 in R1 and R2 may not be directly related to the total amount of the microbial community (bacterial 16S rRNA) in the soil but to the specific function that these microorganisms play regarding soil carbon degradation (pmoA).


Material and methods
Location and characterization of the study area. This study was conducted in an agricultural area destined for sugarcane production located in Barrinha, São Paulo State, Southern Brazil. The geographical coordinates of the area are 21°13′ S and 48°07′ W, with a mean altitude of 555 m above sea level. The soil is classified as an Oxisol (Eutrorthox, USDA Soil Taxonomy), with a slope lower than 0.5%. The area occupies a geomorphological province named basaltic cuestas. The regional climate is defined as B 1 rB′ 4 a′ according to the Thornthwaite classification, that is, a humid mesothermal climate with a small water deficit and summer evapotranspiration lower than 70% of the annual evapotranspiration 32 . The mean annual precipitation registered during the last 30 years was 1560 mm, concentrated from October to March, and a mean annual temperature of 22.2 °C.
Sugarcane has been cultivated in the area for more than 10 years under the mechanized harvesting system, with the last field reform carried out in 2008. Ratoons were mechanically eliminated before tillage operations using an implement consisting of rotary hoes, which cuts the soil and ratoons at a high rotation, and throws them against an impact plate of the implement at a high speed, with the breaking of clods and separation of ratoons from the soil. Limestone and gypsum were applied in the area after ratoon elimination, followed by a leveling harrow operation.
Soil tillage operation consisted of the use of an offset disc harrow with 28-notched discs of 28″, half of them in the front section and the other half in the rear section, a working width of 3.5 m, and a working depth of 0.25 m. Two operations were carried out with this implement at a mean speed of approximately 7 km h −1 , the second operation immediately after the first one to simulate the effect of the disc harrow 11 . A regular grid of 90 × 90 m containing 100 points with minimum distances between sampling points of 10 m was installed after these operations.
Measurements of soil CO 2 emission, soil temperature, and soil moisture. Ten measurements were performed over a period of 28 days (September 24, 26, and 30 and October 3, 7, 10, 14, 16, 18, and 21, 2013) from 8:00 to 10:50 h. A portable LI-8100 automated soil CO 2 flux system (LI-COR, Lincoln, NE, United States) was used to measure the soil CO 2 emission (FCO 2 ) in the experimental area. This system monitors changes in CO 2 concentration inside a closed chamber using optical absorption spectroscopy in the infrared spectrum. The chamber is a closed system with an 854.2-cm 3 internal volume and an 83.7-cm 2 circular soil contact area coupled to PVC collars previously inserted at each sampling points to a depth of 3 cm. In the measurement mode, FCO 2 was determined inside the chamber every 2.5 s and 1.5 min were required to record it at each point 11,16,17,23,33 .
Soil temperature (Ts) was measured using a temperature sensor from the LI-8100 system, which consists of a 20-cm probe that was inserted into the soil near the PVC collars. Similarly, soil moisture (Ms) was measured using a Time Domain Reflectometry (TDR) system (Hydrosense, Campbell Scientific Inc., Logan, UT, United www.nature.com/scientificreports/ States), which consists of two 12-cm probes that were also inserted into the soil near the PVC collars. The measurements of Ts and Ms were carried out concomitantly with FCO 2 assessments.
Determination of soil physical and chemical attributes. Soil samplings were performed at the end of FCO 2 assessments at a depth of 0-20 cm and each grid point. Therefore, 100 disturbed soil samples were collected for chemical analysis, and 100 undisturbed soil samples were taken for physical analysis. For soil chemical analysis, samples were collected with a Dutch auger, being then dried, decloded, and sieved through a 2-mm mesh sieve. The analyses included soil organic matter content (SOM), phosphorus availability (P), K, Ca, Mg, and H + Al contents 34 , sum of bases (SB), and cation exchange capacity (CEC). The total organic carbon (TOC) was estimated by dividing SOM by 1.724. The total soil nitrogen (N) was quantified by the Kjeldahl method 35 .
The particle size distribution of sand, silt, and clay were determined by the pipette method using 0.1 mol L −1 NaOH as a chemical dispersant and low rotation mechanical stirring for 16 h 36 . Fractions containing particles larger than 0.1 mm were separated by sieving (0.105 mm sieve) and those of smaller size by sedimentation, according to the Stokes law; silt was determined by difference. Soil bulk density (Ds) was determined from undisturbed soil samples collected with a sampler adapted to cylinders with an internal diameter of 5 cm and height of 4 cm. The total pore volume (TPV) was calculated based on Ds, with the pore size distribution determined by a porous plate funnel under a 60-cm water column tension in previously saturated samples. The pore volume retained in the sample corresponded to the micropores, and the difference calculated between TPV and micropores corresponded to the macropores 36 . Air-filled pore space (AFPS) was calculated as the difference between the porosity fraction filled with water (Ms) and TPV.
Soil carbon stock (Cstock) was calculated based on the methodology described by Veldkamp 37 , using Eq. (1).
where Cstock is the soil carbon stock (Mg ha −1 ), OC is the organic carbon content (g kg −1 ), Ds is the soil bulk density (Mg m −3 ), and E is the soil layer depth (20 cm). Soil carbon stability in the sugarcane field reform area was obtained by Eq. (2).
where k is the soil carbon decay constant and C-CO 2 is the labile carbon emitted into the atmosphere as CO 2 .
Microbiological analyses. Two regions with distinct emission potentials were identified after FCO 2 interpolation, being named as R1 (Region 1) and R2 (Region 2), from which soil samples were collected for the microbiological characterization.
Quantitative real-time PCR. DNA extraction from environmental samples. Soil samples were collected from nine points at each experimental area (R1 and R2). Each of these points came from composite sampling, which means that 15 samples were collected in the area around a central point to obtain its results. Pre-sterilized PVC tubes with dimensions of 50 cm in length and 5 cm in diameter were used for the soil sampling. These PVC tubes were inserted vertically into the soil to collect soil samples. Subsequently, all the PVC tubes were sealed to avoid losses, stored in iceboxes, and sent to the Laboratory of Biochemistry of Microorganisms and Plants of the School of Agricultural and Veterinarian Sciences of the São Paulo State University (FCAV-UNESP), Jaboticabal campus, São Paulo State, Brazil, where they were stored in an ultra-freezer at − 80 °C to be used for genomic DNA extraction. The genomic DNA of the soil collected at each region was extracted by the FastDNA SPIN Kit for Soil (MP Biomedicals), following the manufacturer's instructions, and stored at − 20 °C. A 5-µL aliquot was submitted to electrophoresis at 1% (w/v) agarose gel stained in GelRed (Uniscience) in SB buffer to confirm the DNA extraction quality 38 . A 2-µL Low DNA Mass Ladder (Invitrogen) was used as a molecular standard. This gel was submitted to an electric field of 85 V for approximately 30 min. Subsequently, the DNA was quantified in a Nanodrop 2000c spectrophotometer (Thermo Scientific) with an optical density ratio of 1.0 at 260 nm (OD 260 ) equal to 50 ng of DNA μ −139 .
Detection and quantification of bacterial 16S rRNA, pmoA, and nifH genes by quantitative real-time PCR. Quantitative real-time PCR reactions were performed on the StepOnePlus Real-Time PCR System (Applied Biosystems) for bacterial 16S rRNA, pmoA, and nifH genes in the SYBER Green I system. The established conditions for each gene are shown in Table 1. A melting curve was added at the end of each reaction under the following conditions: 95 °C for 15 s and primer annealing temperature for 1 min, gradually increased with the reading of the data at every 0.7 °C until 95 °C.
Standard curves were obtained by performing amplification with the copy number of template DNA of Pseudomonas fluorescens for the bacterial 16S rRNA gene, environmental sample DNA for the pmoA gene, and Bradyrhizobium liaoningense for the nifH gene. The data obtained by amplification of the DNA extracted from soil were interpolated to determine the copy number of the genes under study. The real-time PCR reaction for each gene was prepared at a final volume of 10 µL, containing 5 µL SYBR Green ROX qPCR Master Mix (Thermo Scientific), 5 pmol of each forward and reverse primers, 1 µL DNA of test sample, and sterile ultrapure water. The data of qPCR were obtained by using the StepOne Software 2.2.2 (Applied Biosystems), being subsequently exported to Excel (Microsoft) to calculate the number of gene copies per gram of dry soil.  46 , respectively. Moreover, the enzymatic activity of amylase was determined according to Barbosa 44 , in which substrate extraction was carried out by the Cole (1977) method, followed by the determination of reducing sugars by the Somogyi 47 method. The enzymatic activity of cellulase was determined following the method proposed by Kanazawa and Miyashita 48 and the reducing sugar method of Somogyi 47 .
Data processing and statistical analysis. The data were initially analyzed using descriptive statistics (mean, standard error of the mean, standard deviation, minimum, maximum, coefficient of variation, skewness, and kurtosis). The analysis of variance was performed using the software R 49 . The spatial dependence was analyzed using geostatistical techniques, with the estimate of the experimental variograms and permissible model adjustments. The variogram was estimated under the assumption of the intrinsic hypothesis by Eq. (3).
where N(h) is the number of pairs of points of the experimental observations Z(x i ) and Z(x i + h) separated by the h distance. The parameters nugget effect (C 0 ), sill (C 0 + C 1 ), and range (a) were estimated for each variable in the adjustment of mathematical models to the experimental variograms. The spatial dependence index (SDI) was used to classify the spatial dependence as weak, moderate, or strong, according to Seidel and Oliveira 50 . This index considers the C 0 , a, and contribution (C 1 ) values, the variogram model, and the maximum distance between sampling points. The choice of the best fit to the experimental variogram was based on the lowest sum of squared residuals, the highest coefficient of determination (R 2 ), and the cross-validation through the calculation of the root mean square error (RMSE) 51,52 .
The estimate of ordinary kriging (OK) at the non-sampled point x 0 is given by Eq. (4).
in which ẑ(x 0 ) is the OK estimate at the point x 0 , z(x i ) is the neighboring value at the site x i , i = 1, 2, ..., n , and i is the weight of observations associated with the neighboring values, estimated based on the adjusted variogram. The construction of maps of spatial patterns (interpolation by OK) was performed using the software Surfer version 9.11.947 (Golden Software Inc., Golden, CO, United States). Soil CO 2 emission and the data analyzed in the regions R1 and R2 were submitted to the multivariate exploratory analyses of hierarchical clustering and principal components. The hierarchical clustering analysis is an exploratory multivariate technique to assemble the sample units into groups, allowing the homogeneity within the group and the heterogeneity between groups. The structure of groups in the data is found in a dendrogram constructed with the similarity matrix between samples 53 using the Euclidean distance, and the linkage of groups was conducted by the Ward method. Hotelling's T 2 test was performed to test a possible significant difference between groups observed in the dendrogram. Heatmaps were generated using the Euclidean distance as the distance method, being processed using the package gplots in the software R 49 .
Principal component analysis (PCA) is also an exploratory multivariate technique that condenses the information contained in a set of original variables into a smaller set consisting of new latent variables, preserving a relevant amount of original information. The new variables are the eigenvectors (principal components) generated by linear combinations of the original variables, constructed with the eigenvalues of the covariance matrix 54 . Principal components whose eigenvalues were higher than the unity were considered in this analysis, according to the criterion established by Kaiser 55 . The coefficients of linear functions, which define the principal components, were used to interpret their meaning using the signal and relative size of coefficients as an indication of the weight to be assigned to each variable. Only coefficients with high values, usually those higher than or equal to 0.70 in absolute value, were considered in the interpretation. The multivariate analyses were processed using the software Statistica 7.0 (StatSoft Inc., Tulsa, OK, United States).

Results and discussion
The mean value of FCO 2 (2.91 µmol m −2 s −1 ) over the 28-day period ( Table 2) was similar to that observed by Iamaguti et al. 11 in a study also conducted on an Oxisol in sugarcane areas under the most intense soil tillage (conventional soil tillage). Moreover, the values observed for Ts (22.6 °C) and Ms (16.95%) were in accordance with the values reported in other studies carried out after sugarcane field reform 8,9 . These values are consistent with the characteristics of areas under reform, as soil tillage provides low soil moisture values and high soil temperatures. Tillage operations fractionate soil and incorporate air into it, affecting its temperature regime and accelerating the drying process 8 . Soil tillage also mixes the surface layer, more fertile and rich in organic matter, with deeper soil layers, which is associated with better aeration conditions, adequate moisture, high temperatures, increases in FCO 2 due to favorable conditions to the aerobic microbial activity in the decomposition of soil organic matter 56 , and consequent release of CO 2 from soil to the atmosphere 7,57,58 .
In addition, soil aggregate breakdown due to mechanical tillage activities provides high emissions due to the release of the CO 2 stored in deeper soil layers, an effect usually observed in the short term, which comprises intervals from five 59 to 24 h after soil tillage 4 . Soil physical attributes and crop growth are also affected by management systems 2 . Soil density and total porosity reflect the impact to the soil caused by tillage systems and the machinery traffic in the area 4 . The soil under study had a high density (1.33 g cm −3 ) and predominance of micropores (39.55%) compared to macropores (11.08%). Agricultural areas under conventional soil tillage tend to have a high macroporosity and low microporosity 2 , justifying the results found in our study. Moreover, soil management does not change soil texture, and soils with a clayey texture (61.58%) ( Table 2) present high microporosity 60 .
The sugarcane field reform area had the following mean values for chemical attributes: SOM of 32.  Table 2). The variation of soil attributes can be classified according to its coefficient of variation (CV) ( Table 2). The CV value of 19.33% observed for FCO 2 , although considered Table 2. Descriptive statistics of soil physical and chemical attributes assessed in the sugarcane field reform area. N = 100, except for *, in which N = 900. SE standard error of the mean, SD standard deviation, CV coefficient of variation (%), FCO 2 soil CO 2 emission, Ts soil temperature, Ms soil moisture, Ds soil bulk density, Macro macroporosity, Micro microporosity, TPV total pore volume, AFPS air-filled pore space, Clay clay content, Silt silt content, Sand sand content, SOM soil organic matter, Cstock soil carbon stock, k soil carbon decay constant, P phosphorus content, K exchangeable potassium content, Ca exchangeable calcium content, Mg exchangeable magnesium content, Al aluminum content, H + Al potential acidity, SB sum of bases, CEC cation exchange capacity, V base saturation.  11 , also conducted after sugarcane field reform. In addition, the CV values found for microporosity (12.18%), AFPS (15.10%), silt content (16.62%), and Cstock (12.99%) are also considered moderate (12% < CV < 24%). The variables Ts (2.62%), Ms (10.29%), Ds (6.91%), TPV (8.94%), clay content (7.85%), sand content (9.11%), SOM (10.56%), and pH (5.23%) presented CV values considered low (CV < 12%). Moreover, the CV values for P (114.05%), K (31.18%), Ca (78.54%), Mg (47.07%), Al (162.16%), H + Al (24.78%), SB (66.26%), CEC (32.46%), and V (24.13%) are considered high (CV > 24%) ( Table 2). Most of these attributes listed by their high CV are affected by soil chemical management. The high CV values of some attributes suggest a high heterogeneity around the mean. This heterogeneity may have several causes, such as the accumulation and distribution of soil particles as a function of the relief shape and water flow and the lack of homogeneity during the chemical management with soil corrective application in the total area. Considering that the study area presents a smooth relief with a slope lower than 0.5%, the chemical management may assist in understanding the high variability of some soil chemical attributes.
In general, classical statistical methods use both measures of central (mean) tendency and measures of dispersion (variance) to describe a given soil attribute, considering that soil variability occurs entirely random and assuming that its attributes have a normal frequency distribution 62 . However, several studies have shown that soil attributes have spatial dependence 16,17,23,63 , and geostatistical techniques are required to capture the spatial pattern of these attributes 16 . In this sense, the geostatistical analysis showed that the models adjusted to the variograms were mainly spherical, except for clay and silt contents, which presented an exponential model (Table 3). Exponential models are best adjusted to erratic phenomena on the small scale, while spherical models describe properties with a high spatial continuity or less erratic in a short distance 51 .
Spherical models 24,25,28,[63][64][65][66] or the alternation between spherical and exponential models [67][68][69] have been used to describe the spatial variability of  Table 3). A high range value for the spatial variability structure indicates a distribution with high spatial continuity. Table 3. Models and parameters adjusted to the experimental variograms of soil attributes. N = 100; *Logarithmic transformation; C 0 = nugget effect; C 0 + C 1 = sill; a = range (m); SDI spatial dependence index, classified according to Seidel and Oliveira 50 , SSR sum of squared residuals, RMSE root mean square error, obtained through cross-validation, FCO 2 soil CO 2 emission, Ts soil temperature, Ms soil moisture, Ds soil bulk density, Macro macroporosity, Micro microporosity, TPV total pore volume, AFPS air-filled pore space, Clay clay content, Silt silt content, Sand sand content, SOM soil organic matter, Cstock soil carbon stock, k soil carbon decay constant, P phosphorus content, K exchangeable potassium content, Ca exchangeable calcium content, Mg exchangeable magnesium content, Al aluminum content, H + Al potential acidity, SB sum of bases, CEC cation exchange capacity, V base saturation. www.nature.com/scientificreports/ Changes in the ranges of spatial variability models of FCO 2 have been observed between seasons 30,68 and months 70 , after precipitation events 63,69,71 , or even according to the size of the sampling grid 28 . In a similar study conducted in the same area as the present study, Silva et al. 69 observed that the space-time variation of soil CO 2 emission, soil temperature, soil moisture, and soil aeration were affected by three periods related to the same precipitation event. Moreover, the authors incorporated a correlation analysis in the spatial component and identified sites where soil moisture and air-filled pore space were the only factors that influenced soil respiration.
The spatial patterns (maps) of FCO 2 and other soil attributes (Fig. 1) indicated a subdivision (top and bottom) of the studied area regarding FCO 2 , macroporosity, AFPS, silt content, SOM, K, Ca, Mg, and k. This behavior was already expected for k, as this constant was determined from FCO 2 (C-CO 2 ). The similarity in the spatial patterns of these attributes indicates a relationship between them. The spatial structure of soil attributes can be affected by numerous factors in a highly complex way due to the spatial and temporal covariation between influencing factors 72 . Also, soil management practices contribute to increasing the variability of soil attributes because their characteristics are affected 33 . Different agricultural practices induce spatial heterogeneity mainly by affecting the ability to retain carbon, water, and nutrients 25,72 . In this context, as observed in the maps, the direct relationship between FCO 2 and AFPS, macroporosity, SOM, Ca, and Mg is probably due to the activities in the sugarcane field reform, as soil tillage increases the total porosity 4 .
Moreover, soil correction practices increase Ca and Mg contents 56 , promoting favorable conditions for SOM decomposition and, consequently, an increased FCO 2 . Liming carried out before soil tillage operations contributes to increasing FCO 2 due to the chemical dissolution of carbonate, releasing bicarbonate (2HCO 3 − ), which turns into CO 2 , besides to changes in biological processes that increase SOM mineralization in response to an increased pH 73,74 . The direct spatial relationship observed between SOM and FCO 2 occurs because the organic matter represents the main energy reservoir for microorganisms 75,76 . Thus, the increased supply of energy to the microbial metabolism results in an increased release of CO 2 into the atmosphere. According to 25 , soils without vegetation cover present texture properties that contribute to the spatial and temporal variation of FCO 2 since they condition the spatial variability of soil water contents. As in our study, these authors also observed a relationship between FCO 2 and silt content.
Regarding Ms, a division in the east and west was observed in the area. On the other hand, Ts presented a certain homogeneity throughout the area, except for the central region. The maps did not show a relationship between FCO 2 and Ts and Ms. Several authors have observed that the contribution of these factors is not so great when we analyze the spatial variability of FCO 2 67,77 . Soil temperature is characterized as one of the main drivers of the temporal variability of soil respiration 18,78 . However, its covariation with soil moisture masks the spatial correlation between soil respiration and soil temperature 25,72,79 .
As observed in the study area, soil respiration is exclusively due to microbial activity without the presence of plants 25 . It increases the variability of the flow of soil greenhouse gases, such as CO 2 , because they are produced or consumed by a great variety of organisms 80 . In this context, most of the studies on FCO 2 variation have observed relationships with microbial activity 28 , microbial biomass and organic matter content 20 , and variations in the biomass of living and dead organisms associated with the total soil porosity 15 . Although the spatial pattern of soil microorganism diversity is often associated with the vegetation present in the ecosystem 26 , changes occurring in the environment, especially those associated with soil physicochemical management practices, significantly affect its microbiological community 27,81 . This effect occurs because physicochemical management influences the abundance and selection of specific communities of soil microorganisms. These microorganisms are the main agents responsible for GHG production, being highly responsive to any change in the environment. Thus, processes that result in GHG emissions depend directly on the functional diversity and structure of soil microbial communities 82,83 .
The abundance of the bacterial 16S rRNA, pmoA, and nifH genes ( Fig. 2b-d), the enzymatic activity of dehydrogenase, amylase, urease, and cellulase, microbial biomass carbon, and soil C/N ratio were determined after mapping FCO 2 and identifying the regions R1 and R2 (Fig. 2a), with different potential emissions (Fig. 3). No significant difference was observed in the copy number of the bacterial 16S rRNA gene when the regions R1 (4.3 × 10 9 g −1 soil) and R2 (3.1 × 10 9 g −1 soil) were compared (Fig. 2b). Therefore, the abundance of microorganisms present in the soil of the respective regions was similar. The same was verified regarding the copy number of the nifH gene associated with the nitrogen cycle in R1 (823.33 g −1 soil) and R2 (3,541.67 g −1 soil) (Fig. 2d).
However, a difference (p < 0.01) was observed between R1 (9.5 × 10 4 g −1 soil) and R2 (2.9 × 10 4 g −1 soil) when assessing the abundance of the pmoA gene, related to the carbon cycle (Fig. 2c). It indicates that the spatial pattern of FCO 2 in R1 and R2 may not be directly related to the total amount of the microbial community (16S rRNA) present in the soil, but to the specific function that microorganisms play related to soil carbon degradation (pmoA). Also, the abundance of pmoA genes increases during the dry season 41 , the period when the present study was conducted. Sengupta and Dick 84 assessed the methanotrophic bacterial diversity in two different soils under varying land-use practices as determined by high-throughput sequencing of the pmoA gene and observed that the study of the diversity of specific groups, such as pmoA, can lead us to more accurate interpretations of the correlations between the microbial community and edaphic variables, mainly in studies involving changes in land use aiming at GHG mitigation strategies.
Microorganism populations play an essential role in recycling soil chemical elements by controlling the dynamics of the decomposition and stabilization of carbon 85 and, consequently, spatiotemporal variability patterns of FCO 2 into the atmosphere. Moreover, the functional diversity and structure of soil microbial communities are of particular importance for soil functioning as an ecosystem 82,83 . In this sense, the pmoA gene has been widely used to detect the presence of methanotrophic bacteria in the soil because they have the ability to use methane as the sole source of carbon and energy and hence play an important role in the global carbon cycle, being potential biodegrading agents 86   www.nature.com/scientificreports/ The determination of soil enzymes also indicates the functional diversity of the microbial community 87 because soil enzymes are mostly produced from microorganisms. The enzymatic activity has a great potential to indicate biological transformations of the soil in response to changes in its management, as they are sensitive to soil management and directly related to nutrient transformations 88 . Significant differences (p < 0.01) were observed for the enzymatic activity of dehydrogenase (Fig. 3a), urease (Fig. 3b), and amylase (Fig. 3c) when R1 and R2 were compared. On the other hand, no significant difference was observed for the enzymatic activity of cellulase (Fig. 3d). The enzymatic activity of dehydrogenase was higher in R1 (33.02 µg TPF g −1 dry soil 24 h −1 ) when compared to R2 (22.90 µg TPF g −1 dry soil 24 h −1 ) (Fig. 3a). R1 also presented the highest FCO 2 (Fig. 2a). Dehydrogenase occurs intracellularly in all living cells of microorganisms, not accumulating extracellularly in the soil, that is, it is only present in living and active organisms 89 .
Determining the enzymatic activity of urease is also important for understanding the spatial pattern of FCO 2 since this enzyme is responsible for catalyzing the hydrolysis of urea exoenzymes to form CO 2 and ammonium 90 . The differences (p < 0.01) between R1 (41.15 µg NH 4 -N g −1 dry soil 3 h −1 ) and R2 (31.07 µg NH 4 -N g −1 dry soil 3 h −1 ) (Fig. 3b), with a higher biosynthesis of urease in R1, reinforces the discussion on the FCO 2 variability as a function of the higher and lower emission potential be related to soil microbial characteristics in these respective areas.
The enzymatic activity of amylase was also higher in R1 (73.84 µg glucose g −1 dry soil 24 h −1 ) than in R2 (64.81 µg glucose g −1 dry soil 24 h −1 ) (Fig. 3c). Amylases, among other enzymes synthesized mostly by soil microorganisms, are responsible for SOM mineralization 91 , that is, in this process of organic matter transformation, for instance, the carbon of carbohydrate molecules is released as CO 2 75 , explaining its higher activity in R1, which presented a higher FCO 2 .
Cellulase is an enzyme responsible for catalyzing cellulose hydrolysis, being partially responsible for the litter decomposition rate and process 92 . Its increase in the soil indicates the entry of the cellulose-enriched substrate into the agrosystem 93 . Cellulase did not differ between the regions R1 and R2 (Fig. 3d), which may be an indication that the sugarcane residue (straw and dead roots) was homogeneously incorporated into the  www.nature.com/scientificreports/ soil at the time of tillage (Fig. 1b), not favoring a given region compared to other regarding the incorporation of residues in the soil. Soil microbial biomass carbon represents the active and biodegradable fraction of SOM, being partially composed of several microorganism species, such as fungi, bacteria, protozoa, nematodes, and algae, acting as agents in the organic matter mineralization 94 . A significant difference (p < 0.01) was observed for MBC between R1 (41.35 µg C g −1 soil) and R2 (17.87 µg C g −1 soil) (Fig. 3e). Therefore, R1 presented the highest potential for CO 2 emission (Fig. 2a) and also the largest active fraction of SOM, being 131% higher in R1 than in R2. On the other hand, the soil C/N ratio was lower in R1 (12.18) than in R2 (15.43) (Fig. 3f). The dynamics of microbial activity is regulated, among other factors, by the C/N ratio since the relationship between these two elements in the soil interferes with the degree of humification and stability of SOM 95 . The fact that this ratio is lower in R1 indicates that the organic matter decomposition process is more accelerated in this region.
The analysis of variance indicated differences between the regions R1 and R2 when soil microbiological attributes were univariately assessed (Figs. 2 and 3). In addition, the fact that FCO 2 was also higher in R1 and lower in R2 (Fig. 2a) allowed establishing indirect relationships. However, the characterization of a pattern or behavior depends on several interactions among the assessed factors. In this context, multivariate analyses of the data can be very efficient since they allow visualizing natural correlations and multiple influences on behavior, especially when using multivariate techniques of interdependence, that is, when no variable or group is treated as dependent or independent 50 . Therefore, when the multivariate analysis was carried out by the hierarchical clustering method, the dendrogram constructed from samples of FCO 2 and microbiological attributes reaffirmed the existence of two groups: Group 1 (R2), linked to a shorter Euclidean distance, and Group 2 (R1), linked to a larger Euclidean distance (Fig. 4). In addition, a significant effect (F = 150.76; p = 0.007) was observed when carrying out Hotelling's T 2 test, confirming that R1 and R2 are distinct groups. Multivariate analysis techniques are only efficient when there is a correlation structure between variables 50 . In our study, this structure is shown in Fig. 5 and Table 4. Figure 5 shows the biplot with the first two principal components (PC1 and PC2), while the variance assigned to them is shown in Table 4. PC1 represents 64.11% and PC2 represents 12.23% of the total variance of the original data (Fig. 5). All the assessed attributes were retained in PC1 according to the cut-off criterion (> 0.70 in absolute value) considered for interpreting the principal component (Table 4). In PC1, there is a direct dependence relationship among the enzymatic activity of dehydrogenase (0.73), urease (0.85), and amylase (0.70), microbial biomass carbon (0.75), FCO 2 (0.76), and pmoA functional gene (0.88), all directed to the same direction as the samples from R1 (right side of the two-dimensional plane). On the other hand, there is a relationship of indirect dependence of these attributes with the C/N ratio (− 0.86), located in the opposite direction, where the samples of R2 are also located (left side of the two-dimensional plane) (Fig. 5 and Table 4).
These relationships indicate that the process occurring in PC1 may be related to the soil CO 2 emission potential in R1 and R2. Therefore, for interpretation purposes, PC1 can be termed as FCO 2 potential as a function of soil microbial diversity because there is a close relationship between the spatial pattern of FCO 2 and www.nature.com/scientificreports/ soil microbiota dynamics. This dynamic is regulated by the C/N ratio, that is, the readily available C content, which, as soil pH, has been described in the literature as controlling the taxonomic and functional structure of the microbiota [96][97][98] . In addition, adequate conditions of aeration, nutrient supply, soil moisture, and soil temperature ( Table 2) favored the diversity of soil microorganisms (Fig. 2) and organic matter biodegradation by stimulating the synthesis of soil enzymes (Fig. 3). In this context, all the activities carried out in the reform area, which included soil correction and conventional tillage with residue incorporation, certainly affected the diversity of soil microorganisms and, consequently, FCO 2 values, providing conditions for the establishment of the regions R1 and R2. Thus, region R1 presents a higher FCO 2 potential than R2 (Fig. 2a), which may be associated with a higher soil microbial activity in R1 (Figs. 2c and 3a-e) favored by the lower C/N ratio (Fig. 3f) and better nutritional status provided by the higher content of SOM, Ca, and Mg (Fig. 1), as well as better soil aeration conditions due to the higher percentage of macropores and air-filled pore space in R1 (Fig. 1).

Figure 5.
Biplot showing the assessed soil attributes and regions. Amylase = enzymatic activity of amylase; Dehydrogenase = enzymatic activity of dehydrogenase; pmoA = pmoA functional gene; MBC = microbial biomass carbon; Urease = enzymatic activity of urease; FCO 2 = soil CO 2 emission. Principal components were processed using the package ggplot2 and factoextra in the software R version 3.6.3 49 (https:// cran.r-proje ct. org/ bin/ windo ws/ base/ old/3. 6.3/). Table 4. Correlation between soil attributes and the first two principal components (PC1 and PC2). *Value referring to the percentage of variation of the original set of data retained by the respective principal components. Correlations in bold (higher than 0.70 in absolute value) were considered in the interpretation of the principal component.

Conclusions
In sugarcane field reform areas, the spatial patterns of soil CO 2 emission are similar to those of the attributes macropores, air-filled pore space, silt content, organic matter, and soil carbon decay constant. In our study, differences in soil CO 2 emission were associated with the dynamics of the microbial activity regulated by the abundance of the pmoA gene (related to the carbon cycle) and the enzymatic activity of dehydrogenase, urease, and amylase, microbial biomass carbon, and soil C/N ratio. Considering that one of the greatest challenges for new directions in world agriculture is to increase the efficiency of conventional practices of soil management, this study reinforces the importance of understanding the spatial pattern of soil attributes and how they influence soil CO 2 emission dynamics in agricultural areas. This aims for future mitigation actions that involve less intense tillage activities and restrained and homogeneous applications of soil inputs, reducing production costs and the contribution of these activities to the additional CO 2 emission during the sugarcane field reform period.