Planting Year-and Climate-Controlled Soil Aggregate Stability and Soil Fertility in the Karst Region of Southwest China

: The effects of long-term monocropping systems combined with climate on soil water aggregate stability (WSA) and soil fertility in the karst region of Southwest China (KRSWC) are unclear. Our research was conducted in the KRSWC, wherein tobacco ( Nicotiana tabacum ) production is characterized by heavy fertilization and continuous monocropping. The tobacco ﬁelds in the study area have similar soil types and fertilization and tillage practices and are spread over an area of 11,500 km 2 . A total of 568 topsoil samples were collected in 2021. Soil fertility was reﬂected using the soil fertility index (SFI), which was calculated using the minimum data set method with six soil fertility-related factors, namely, soil pH, soil organic matter, cation exchange capacity, available nitrogen, available phosphorus, and available potassium. Results showed that long-term planting generally promoted soil fertility levels and WSA content. WSA and SFI had inconsistent spatial distribution patterns likely due to different climate-driven effects. WSA variability was greatly controlled by precipitation (Spearman correlation coefﬁcient [r] = − 0.49, p < 0.01), whereas SFI variability was mostly dominated by temperature (r = − 0.36, p < 0.01). The levels of SFI and WSA were optimal under conditions of low temperature and precipitation and poor under conditions of high temperature and precipitation. Moreover, long-term planting could alleviate the negative effects of climate on SFI and WSA in the KRSWC. The results of this study could provide valuable information on fertilization and climate-adapted strategies for tobacco ﬁelds in the KRSWC.


Introduction
Soil water aggregate stability (WSA) is an important indicator of soil erosion risk and soil sustainability [1][2][3][4].Soil fertility, which is defined as the ability of soil to sustain plant growth, is generally assessed using a scheme proposed in the 1990s [5][6][7].The level of soil fertility determines the performance of sustainable agriculture [7][8][9].According to Jenny's factorial model of soil formation, soil properties are influenced by numerous factors, including climate and human activities [10].Climate determines the moisture and temperature of soil formation [10].Planting years, an important human factor for arable soils, are related to the duration of material input and anthropogenic disturbances [11].A clear understanding of the effect of planting years and climate on soil fertility and structure is important for the sustainable management of regional land use and early warning of unfavorable trends.
Climatic conditions, including temperature and precipitation, are regarded as the key factors affecting soil properties on a regional scale [2,[12][13][14].WSA is governed by four well-known mechanisms (i.e., raindrop impact, slaking, differential swelling, and physicochemical dispersion) [15,16], all of which are closely related to precipitation.Our previous studies confirmed that precipitation in the karst region of Southwest China (KRSWC) exerts a strong negative control on WSA, whereas temperature has a weaker effect than precipitation on WSA [2].The relative importance of temperature and precipitation for soil fertility indicators may depend on regional biological constraints [14].For example, compared with other factors, temperature has a greater influence on soil fertility indicators in China [13] and Europe [12], whereas precipitation contributes more highly in eastern Australia, which is under strong water limitation [17].Past studies revealed positive correlations between WSA and soil fertility with field experiments [18][19][20].However, given the different effects of climate on soil fertility and structure, their positive correlation may break down at regional scales.Therefore, exploring how climate affects the spatial distribution of WSA and soil fertility may help the proposal of climate-adapted agricultural strategies and is thus necessary.
Continuous soil management might alter soil properties, and such change varies across planting years [11,[21][22][23][24][25].For example, the storage of soil carbon, nitrogen, and phosphorus in tea gardens [11] and ginseng fields [25] increases with an increase in planting years.Under conventional farming practices, the soil structure in Brazilian tobacco (Nicotiana tabacum) fields continuously declines with increasing planting years [21].By contrast, longterm no-tillage was found to be an effective strategy for improving soil fertility in a Brazilian oxisol [24].Overall, planting years affect soil fertility and aggregate stability differently under various land use types and soil management strategies.The KRSWC is one of the least underdeveloped regions in China, with a per capita rural income of 86% of the national average in 2016 [26].Tobacco cultivation, as an important income source for farmers in the KRSWC [27], is characterized by high fertilization and continuous monocropping [28,29].Previous long-term field experiments revealed the effects of different tobacco rotation systems on soil fertility and aggregate stability [29,30] but paid little attention to the effects of cropping systems that are widely used in actual tobacco production.Moreover, climate and planting year might interactively influence soil fertility and WSA, which has received scant attention.
We assumed that long-term tobacco monocropping with heavy fertilization affects soil fertility and structure and that climate might break the positive correlation between soil fertility and structure at regional scales.This study was implemented on KRSWC tobacco fields with similar soil types, fertilization, and tillage.Specifically, after quantifying soil fertility (as indicated by the soil fertility index (SFI)) using a classical process [5][6][7], we measured the single and joint controls of planting years and climate on SFI and WSA.

Study Area
The study site is located in southwestern China, which is part of the largest contiguous karst zone in the world [26] (Figure 1).The spatial location of the study area is 107 • 80 -109 • 31 E and 28 • 33 -28 • 99 N.The region is predominantly mountainous and hilly with a humid subtropical monsoon climate.In accordance with the Chinese Soil Genetic Classification System, the vast majority of the soils in the tobacco fields in the study area are classified as yellow soil that developed from Permian or Triassic limestones.In the World Reference Base for Soil Resources system, yellow soil is equivalent to haplic alisols or haplic luvisols.Additional information on the climate, topography, geology, and soil types of the region can be found in our previous publication [2].Tobacco was first planted in the study area in the 1980s.Tobacco, an important cash crop that helps local farmers escape poverty, is under the unified, professional, and careful guidance of the local agricultural department.Tobacco nurtured in the greenhouse is transplanted to the field in early May each year and harvested approximately 120 days later.Tobacco fields are fallow during the rest of the year to maintain soil quality.Conventional tillage processes are used to manage tobacco fields, and agricultural machinery is uniformly equipped by the local agricultural department.The tobacco farms in the study area rely solely on rainwater as their source of moisture and lack supplementary irrigation systems.Field surveys in 2017 showed that fertilization practices in the region are similar, with an inorganic nitrogen fertilizer application rate of 112.5 kg per hectare, which is roughly four times the application rate of organic nitrogen fertilizers.Organic fertilizers include oil cake fertilizer, farmyard manure, and amino acid fertilizer.The amount of fertilizer applied has remained almost unchanged in recent years.
manage tobacco fields, and agricultural machinery is uniformly equipped by the local agricultural department.The tobacco farms in the study area rely solely on rainwater as their source of moisture and lack supplementary irrigation systems.Field surveys in 2017 showed that fertilization practices in the region are similar, with an inorganic nitrogen fertilizer application rate of 112.5 kg per hectare, which is roughly four times the application rate of organic nitrogen fertilizers.Organic fertilizers include oil cake fertilizer, farmyard manure, and amino acid fertilizer.The amount of fertilizer applied has remained almost unchanged in recent years.

Data and Pre-Processing
Topsoil (0-20 cm) samples were collected using random sampling prior to tobacco planting and fertilization (March 2021).For ensuring representativeness, 3-5 subsamples were collected within a 5 m radius from the center of a field-determined sampling point.Each sample was geotagged with GPS.Approximately 1 kg of soil was collected from each subsample.Approximately 2 kg of well-mixed soils was obtained using the quartering

Data and Pre-Processing
Topsoil (0-20 cm) samples were collected using random sampling prior to tobacco planting and fertilization (March 2021).For ensuring representativeness, 3-5 subsamples were collected within a 5 m radius from the center of a field-determined sampling point.Each sample was geotagged with GPS.Approximately 1 kg of soil was collected from each subsample.Approximately 2 kg of well-mixed soils was obtained using the quartering method and placed in ziplock bags for subsequent chemical analysis.Meanwhile, undisturbed soils were collected at the same location, and each subsample was carefully stripped off along the natural profile.Ultimately, a total of 1 kg of soil was stored in aluminum containers.Information about planting years was recorded at the site through field observations and farmer interviews (Table 1).Samples with missing information were deleted.Finally, 568 samples collected from yellow soil were used in the current study.Aggregates with different diameters were measured on the basis of the wet sieving method [31] using undisturbed soil, and WSA was defined as the sum of aggregates with particle sizes larger than 0.25 mm.A high WSA generally indicates high soil aggregate stability.Refer to Zhang et al. for additional details [2].Various soil chemical characteristics, including soil pH, soil organic matter (SOM), available nitrogen (AVN), phosphorus (AVP), and potassium (AVK), as well as cation exchange capacity (CEC), were measured.The following analytical methods were utilized to measure the above characteristics: SOM, potassium dichromate oxidation; AVN, NaOH hydrolysis; AVP, HCl-H 2 SO 4 extraction and colorimetry; AVK, ammonium acetate extraction and colorimetry; soil pH, glass electrode method (soil/water 1:2.5); and CEC, titration [32].
Two climate raster maps of mean annual temperature (MAT) and mean annual total precipitation (MAP) with a spatial resolution of 1 km were sourced from our previous study [2].The raster maps of five commonly used topographic factors, including aspect (Asp), slope (Slp), topographical wetness index (TWI), plan curvature (PLC), and profile curvature (PRC), derived from the 30 m digital elevation model were used for data analysis [2].Altitude was not considered because the strong co-linearity between altitude and MAT may disturb analyses.

Data Analytical Tools
Descriptive statistics were used to examine the degree of concentration and variability across the original data set.Spearman correlation coefficients were calculated between response variables (WSA and SFI) and environmental factors because not all data were normally distributed.Given that the assumptions of homogeneity of variances and normality always were violated, the statistical comparisons of WSA and SFI were performed using the Kruskal-Wallis test with Games-Howell post-hoc analysis.These statistical analyses were completed with SPSS v22.
Principal component analysis (PCA) is a popular technique for data dimensionality reduction that aims to replace a large number of original variables with a few variables (i.e., principal component (PC)) and retain the maximum amount of information [33].The PCs with eigenvalues greater than 1 were used to determine how many PCs should be retained [6,34].Variance partitioning analysis (VPA) was applied to quantify the interpreta-tion ratio of different environmental factors to WSA and SFI.If covariance existed among the factors in the same group, the data were down-dimensioned using PCA before VPA was performed.PCA and VPA were executed with the vegan package of R 4.0.5.
Local indicators of spatial association (LISA) analysis, which was developed by Anselin [35], was used to analyze the local spatial autocorrelation and existence of clusters in the spatial arrangement of a given variable.LISA analysis provided the ability to detect the local spatial associations and clustering of WSA or SFI [36].Five clustering patterns were identified: (i) no significant spatial dependence; (ii) high-high, high value in a highvalue neighborhood; (iii) high-low, high value in a low-value neighborhood; (iv) low-high, low value in a high-value neighborhood; and (v) low-low, low value in a low-value neighborhood.LISA analysis was conducted with ArcGIS v.10.3.Furthermore, we obtained the spatial distributions of SFI and WSA with ordinary kriging in ArcGIS v.10.3.

Soil Fertility Evaluation
The total data set for soil fertility evaluation consists of 6 soil fertility-related factors, including pH, SOM, CEC, AVN, AVP, and AVK.A widely used scheme (minimum data set (MDS) method) for soil fertility evaluation was adopted [5,6,34] and it consists of three steps.
(1) MDS construction: PCA and Spearman correlation analysis were used to select the most appropriate factors indicating the fertility characteristics.As described by Doran and Parkin [5] and Qi et al. [6], only PCs with eigenvalues greater than 1 were used to represent the local soil fertility characteristics.Factors with weighted loading values within 10% of the highest weighted loading among each PC were selected.Redundant factors with strong correlations were eliminated.(2) Factor scoring and weighting: based on MDS, the scores of soil fertility-related factors were assigned using two standard scoring functions (upper limit function [Equation (1)] and peak limit function [Equation ( 2 )].The weights of each factor were determined according to the communality [6].Commonality is the sum of the squared loading for each variable.It spans from 0 to 1, and a value close to 1 denotes greater variance.(3) Calculation of SFI: SFI was assessed for each sample site according to Equation (3).The higher the SFI, the better the soil fertility.The SFI was divided into five grades using the isometric division method [34]: where f (x) is the score, x is the factor value, and L and U are the lower and upper threshold values, respectively.The upper and lower limits used in this work referred to Qi et al. [6] and Qu et al. [7] (Table S1): where W i is the assigned weight, which was determined in accordance with communality [6], N i is the factor score, and n is the number of factors in MDS.

Descriptive Statistics
The descriptive statistics for soil fertility-related factors and the WSA of topsoil (0-20 cm) are shown in Table 2. WSA accounted for 51.67 ± 4.06% of the soil bulk.Soil fertility-related factors were graded (Tables 3 and 4) in accordance with the MLR [37], CNSSO [38], and tobacco growing experience.The samples had average AVN, AVP, and AVK contents of 155.73, 64.63, and 441.37 mg/kg, respectively (Table 2), with more than half of them assigned to Grade I (Table 3).Most of the samples had Grade II CEC contents with a mean ± standard deviation of 17.44 ± 3.13 cmol/kg.The SOM content was low, and 54% of the samples had Grade III SOM contents.Approximately 63.2% of the samples had soil pH values at sub-suitable and unsuitable levels for tobacco growth (Table 4).The study area has a complex topography but little climatic variability (Table 2).WSA, pH, MAT, and MAP exhibited low variability (CV < 15%) [39].SOM, CEC, and AVN showed moderate variability (16-35%).By contrast, AVK, AVP, Asp, Slp, and TWI presented high variability (CV > 35%).

Soil Fertility Assessment
The first three PCs had eigenvalues >1 and explained 82.18% of the variation in all soil fertility-related factors (Table 5).For each PC, the variables within 10% of the highest factor loading were reserved as key indicators for SFI calculation.AVN was removed due to its strong correlation with SOM (r > 0.9, p < 0.01).The other five factors were all chosen.Their weights were quantified using their communality values.SOM was assigned the highest weight (0.24), followed by pH (0.20), CEC (0.20), AVP (0.18), and AVK (0.17).Approximately 64.85% of the soil samples were in Grade I (Table S2).This result implies that local tobacco fields have excellent soil fertility.Note: PTV: percentage of the total variation; CPTV: cumulative PTV; (a) : the Spearman's correlation coefficient between AVN and SOM was greater than 0.9 (p < 0.01), so AVN was removed.Variables in bold type have a high factor loading within 10% of the highest loading in a principle component.The abbreviations and descriptions of soil factors are listed in Table 1.

WSA Content and Soil Fertility in Different Planting Years
The differences in WSA, SFI, and key soil fertility-related factors in different groups of planting years (p < 0.05) were investigated (Figure 2).The numbers of soil samples in the short, middle, and long planting year groups were 46, 327, and 195, respectively.No significant differences in AVP and AVK were found between different planting years.WSA, SFI, pH, SOM, and CEC significantly differed between the long and short planting year groups but were not always significantly different between the middle and short planting year groups or middle and long planting year groups.The general trend of WSA, SFI, and key soil fertility-related factors showed that the contents of WSA, SFI, SOM, CEC, AVP, and AVK increased and pH declined with increasing planting years.

Climate Affects SFI and WSA Differently
The spatial distribution of WSA differed from that of SFI (Figure 3).WSA values were low and high in the east and west, respectively, whereas low SFI values were concentrated in the central area.The spatial structure derived from LISA analysis further confirmed the spatial differences between WSA and SFI (Figure 3).For example, in the central area, SFI exhibited the significant spatial distribution of high-high relationships, whereas WSA generally displayed the spatial distribution of low-low relationships.This result implies that the soils of the tobacco fields in the study region are mainly characterized by low WSA values and high SFI values.In general, the spatial distributions of WSA and SFI were inconsistent at the regional scale, likely due to the different influences exerted by climate on SFI and WSA.The VPA method was used to verify the pure and joint effects of precipitation, temperature, and topography on WSA and SFI.MAP, MAT, and terrain factors accounted for 18.26%, 0.05%, and 1.15% of the total variation in WSA, respectively, and for 0.76%, 11.18%, and 0.55% of the total variation in SFI, respectively (Figure 4a,b).Precipitation dominantly controlled WSA variability, whereas temperature governed SFI variability.Moreover, linear and nonlinear fits revealed monotonically negative correlations (p < 0.01) between MAP and WSA and between MAT and SFI with Spearman correlation coefficients of −0.49 and −0.36 (p < 0.01), respectively (Figure 4c,d 1.

Climate Affects SFI and WSA Differently
The spatial distribution of WSA differed from that of SFI (Figure 3).WSA values were low and high in the east and west, respectively, whereas low SFI values were concentrated in the central area.The spatial structure derived from LISA analysis further confirmed the spatial differences between WSA and SFI (Figure 3).For example, in the central area, SFI exhibited the significant spatial distribution of high-high relationships, whereas WSA generally displayed the spatial distribution of low-low relationships.This result implies that the soils of the tobacco fields in the study region are mainly characterized by low WSA values and high SFI values.In general, the spatial distributions of WSA and SFI were inconsistent at the regional scale, likely due to the different influences exerted by climate on SFI and WSA.The VPA method was used to verify the pure and joint effects of precipitation, temperature, and topography on WSA and SFI.MAP, MAT, and terrain factors accounted for 18.26%, 0.05%, and 1.15% of the total variation in WSA, respectively, and for 0.76%, 11.18%, and 0.55% of the total variation in SFI, respectively (Figure 4a,b).Precipitation dominantly controlled WSA variability, whereas temperature governed SFI variability.Moreover, linear and nonlinear fits revealed monotonically negative correlations (p < 0.01) between MAP and WSA and between MAT and SFI with Spearman correlation coefficients of −0.49 and −0.36 (p < 0.01), respectively (Figure 4c,d).1.

Climate Affects SFI and WSA Differently
The spatial distribution of WSA differed from that of SFI (Figure 3).WSA values were low and high in the east and west, respectively, whereas low SFI values were concentrated in the central area.The spatial structure derived from LISA analysis further confirmed the spatial differences between WSA and SFI (Figure 3).For example, in the central area, SFI exhibited the significant spatial distribution of high-high relationships, whereas WSA generally displayed the spatial distribution of low-low relationships.This result implies that the soils of the tobacco fields in the study region are mainly characterized by low WSA values and high SFI values.In general, the spatial distributions of WSA and SFI were inconsistent at the regional scale, likely due to the different influences exerted by climate on SFI and WSA.The VPA method was used to verify the pure and joint effects of precipitation, temperature, and topography on WSA and SFI.MAP, MAT, and terrain factors accounted for 18.26%, 0.05%, and 1.15% of the total variation in WSA, respectively, and for 0.76%, 11.18%, and 0.55% of the total variation in SFI, respectively (Figure 4a,b).Precipitation dominantly controlled WSA variability, whereas temperature governed SFI variability.Moreover, linear and nonlinear fits revealed monotonically negative correlations (p < 0.01) between MAP and WSA and between MAT and SFI with Spearman correlation coefficients of −0.49 and −0.36 (p < 0.01), respectively (Figure 4c,d).data were categorized into four groups, namely, HT-HP, HT-LP, LT-HP, and LT-LP, in accordance with the median values of MAP and MAT.HT and LT denote high and low temperature groups, respectively, and HP and LP indicate high and low precipitation groups, respectively.WSA differed significantly between the HP and LP groups, and SFI differed significantly between the HT and LT groups (Figure 5).This result is consistent with the finding that WSA and SFI were dominated by precipitation and temperature, respectively (Figure 4).The soil in the LT-LP group exhibited high WSA and SFI contents.Raw data were categorized into four groups, namely, HT-HP, HT-LP, LT-HP, and LT-LP, in accordance with the median values of MAP and MAT.HT and LT denote high and low temperature groups, respectively, and HP and LP indicate high and low precip itation groups, respectively.WSA differed significantly between the HP and LP groups and SFI differed significantly between the HT and LT groups (Figure 5).This result is con sistent with the finding that WSA and SFI were dominated by precipitation and tempera ture, respectively (Figure 4).The soil in the LT-LP group exhibited high WSA and SF contents.Raw data were categorized into four groups, namely, HT-HP, HT-LP, LT-HP, and LT-LP, in accordance with the median values of MAP and MAT.HT and LT denote high and low temperature groups, respectively, and HP and LP indicate high and low precipitation groups, respectively.WSA differed significantly between the HP and LP groups, and SFI differed significantly between the HT and LT groups (Figure 5).This result is consistent with the finding that WSA and SFI were dominated by precipitation and temperature, respectively (Figure 4).The soil in the LT-LP group exhibited high WSA and SFI contents.The differences in WSA and SFI across various climates and planting years were further explored (Figure 6).In each planting year group, WSA significantly differed between different precipitation groups, whereas SFI significantly differed between different temperature groups (similar to the result in Figure 5).In the HT-HP group, the levels of WSA and SFI were significantly higher in the long planting year group than in the middle and short planting year groups.Planting year had no significant effect on WSA and SFI in the other climate groups.
and low precipitation group, LT-HP is the low temperature and high precipitation group, and LT-LP is the low temperature and low precipitation group).
The differences in WSA and SFI across various climates and planting years were further explored (Figure 6).In each planting year group, WSA significantly differed between different precipitation groups, whereas SFI significantly differed between different temperature groups (similar to the result in Figure 5).In the HT-HP group, the levels of WSA and SFI were significantly higher in the long planting year group than in the middle and short planting year groups.Planting year had no significant effect on WSA and SFI in the other climate groups.

Discussion
As indicated by SFI values, the tobacco fields in the study area exhibited excellent soil fertility, with 64.79% and 32.92% of the soil samples classified as Grades I and II, respectively (Table S2).By contrast, the tobacco fields in the study area showed low soil aggregate stability with a mean WSA content of 51.8%, and almost all samples (98.98%) were in an unstable state [2], likely due to their low SOM content.Previous studies have discovered that the effect of SOM on aggregates always exceeds that of other soil chemical properties [4,42,43].Therefore, for the protection of soil aggregates, fertilization strategies should be modified to increase SOM content.For example, the application of organic fer-

Discussion
As indicated by SFI values, the tobacco fields in the study area exhibited excellent soil fertility, with 64.79% and 32.92% of the soil samples classified as Grades I and II, respectively (Table S2).By contrast, the tobacco fields in the study area showed low soil aggregate stability with a mean WSA content of 51.8%, and almost all samples (98.98%) were in an unstable state [2], likely due to their low SOM content.Previous studies have discovered that the effect of SOM on aggregates always exceeds that of other soil chemical properties [4,42,43].Therefore, for the protection of soil aggregates, fertilization strategies should be modified to increase SOM content.For example, the application of organic fertilizer should be increased.The benefits of organic fertilizers for SOM promotion and aggregate protection have been widely confirmed [30,44,45].To illustrate, Zou et al. [29] revealed that manure amendment increases soil macroaggregates and soil organic carbon.Tobacco straws, as valuable organic fertilizer materials, are also effective in enhancing SOM and physical structure [46,47].
WSA, SFI, SOM, and CEC increased significantly with the increase in planting years (Figure 2).Overall, the current cropping patterns are beneficial to tobacco soil.Consistent with previous results [29,48,49], this above finding suggests that the long-term application of organic fertilizers (e.g., manure and oil cake fertilizers) might neutralize the negative effects of continuous monocropping on soil fertility and aggregate stability.For example, soil fertility in tobacco fields significantly improved after 11 years of consecutive tobacco cropping with three different organic fertilizers [48].However, in line with a reported finding [28], soil pH declined with increasing planting years (Figure 3).This phenomenon is an early warning that needs to be taken seriously because declining soil pH potentially threatens soil aggregate stability [3] and fertility [50].Generally, although the current cropping patterns could improve soil fertility and structure in the long term, soil pH decline should be controlled.
Previous works at the field scale usually considered fertility and aggregate stability to be positively correlated [18][19][20].However, our study found that the spatial distribution of soil fertility and WSA was inconsistent at the regional scale (Figure 3), presumably due to the different effects of climate on SFI and WSA.Temperature controlled the vast majority of variability in SFI, whereas precipitation mainly governed WSA.Our previous work demonstrated the critical negative control of WSA by precipitation [2].The relative contribution of temperature and precipitation to soil fertility varied across environmental conditions [12,13,17,51], which depend on the main limiting factor of biological activity [14].For example, Ballabio et al. [12] and Song et al. [13], respectively, revealed that temperature has great effects on soil fertility factors in Europe and China.The role of precipitation is greater in areas under strong water limitation (e.g., eastern Australia) [17].Given that our study area is not under water limitation (annual precipitation of 1000-1500 mm), temperature dominated SFI variability.The negative correlation between temperature and soil fertility found in our research is consistent with that discovered with studies in southwestern China [50,52,53] and may be closely related to soil microbial activity [14].Increased planting years significantly and positively affected WSA and SFI in environments with high temperature and rainfall (Figure 6).This result shows that long-term tobacco planting could alleviate the negative effects of climate on soil fertility and structure.

Conclusions
Our regional-scale study was implemented on tobacco fields in the KRSWC, wherein tobacco production is characterized by high fertilization and continuous monocropping.SFI values showed that in the study area, the soils under tobacco production are of good fertility but have low SOM levels, which may result in unstable aggregates.Long-term tobacco planting promotes soil fertility and WSA.The spatial variabilities of WSA and SFI are largely negatively controlled by precipitation and temperature, respectively.Climate controls WSA and SFI differently, leading to inconsistencies in the spatial distribution of these factors.WSA and SFI are highest under the LT-LP condition and lowest under the HT-HP condition.Interestingly, considering the combination of climate and planting years, long-term tobacco planting exerts a significant positive effect only under the HT-HP condition.In general, the above results reveal that long-term tobacco planting could promote soil fertility and structure and alleviate the negative effects of climate on soil fertility and structure.However, fertilization strategies should be modified, for example, by increasing the application of organic fertilizers to improve SOM content, and thus protect aggregates.Evidently, long-term planting caused temporal variability of tobacco soil fertility and structure.Therefore, continuous investigations are necessary in the future.

Figure 1 .
Figure 1.(a) Soil sample points, digital elevation model, and spatial distribution of tobacco fields on Permian and Triassic limestone; (b) geographic location of the study area.

Figure 1 .
Figure 1.(a) Soil sample points, digital elevation model, and spatial distribution of tobacco fields on Permian and Triassic limestone; (b) geographic location of the study area. ).

Figure 2 .
Figure 2. Boxplots of WSA, SFI, and key soil fertility-related factors in different planting years.Different letters indicate significant differences at p < 0.05 between different planting years.Planting year grouping: short: <10 years; middle: 10-20 years; long: ≥20 years.SFI: soil fertility index.The abbreviations and descriptions of soil factors are listed in Table1.

Figure 2 .
Figure 2. Boxplots of WSA, SFI, and key soil fertility-related factors in different planting years.Different letters indicate significant differences at p < 0.05 between different planting years.Planting year grouping: short: <10 years; middle: 10-20 years; long: ≥20 years.SFI: soil fertility index.The abbreviations and descriptions of soil factors are listed in Table1.

Figure 2 .
Figure 2. Boxplots of WSA, SFI, and key soil fertility-related factors in different planting years.Different letters indicate significant differences at p < 0.05 between different planting years.Planting year grouping: short: <10 years; middle: 10-20 years; long: ≥20 years.SFI: soil fertility index.The abbreviations and descriptions of soil factors are listed in Table1.

Figure 3 .
Figure 3. Ordinary kriging interpolation maps and LISA cluster maps of WSA (a) and SFI (b).

Figure 3 .
Figure 3. Ordinary kriging interpolation maps and LISA cluster maps of WSA (a) and SFI (b).

Figure 4 .
Figure 4. Variation partitioning analysis for WSA (a) and SFI (b).X1-X3: pure effect of each explan atory variable; X4-X7: joint effects of two or three explanatory variables.Linear and non-linear re gressions between MAP and WSA (c), and between MAT and SFI (d).Linear fits are indicated by red solid.Nonlinear relationships are represented by the generalized additive model (GAM), shown by the black dashed line.The shaded area represents 95% confidence interval.**: significant at 0.01 level; r: Spearman correlation coefficient; AIC: Akaike information criterion.

Figure 4 .
Figure 4. Variation partitioning analysis for WSA (a) and SFI (b).X1-X3: pure effect of each explanatory variable; X4-X7: joint effects of two or three explanatory variables.Linear and non-linear regressions between MAP and WSA (c), and between MAT and SFI (d).Linear fits are indicated by red solid.Nonlinear relationships are represented by the generalized additive model (GAM), shown by the black dashed line.The shaded area represents 95% confidence interval.**: significant at 0.01 level; r: Spearman correlation coefficient; AIC: Akaike information criterion.

Figure 4 .
Figure 4. Variation partitioning analysis for WSA (a) and SFI (b).X1-X3: pure effect of each explanatory variable; X4-X7: joint effects of two or three explanatory variables.Linear and non-linear regressions between MAP and WSA (c), and between MAT and SFI (d).Linear fits are indicated by red solid.Nonlinear relationships are represented by the generalized additive model (GAM), shown by the black dashed line.The shaded area represents 95% confidence interval.**: significant at 0.01 level; r: Spearman correlation coefficient; AIC: Akaike information criterion.

Figure 5 .
Figure 5. Boxplots of (a) WSA and (b) SFI in different climate groups (different letters indicate significant differences at p < 0.05 in different climate groups.Red bolded numbers indicate the number of soil samples in different climatic groups.WSA is water-stable aggregates, SFI is the soil fertility index, HT-HP is the high temperature and high precipitation group, HT-LP is the high temperature

Figure 5 .
Figure 5. Boxplots of (a) WSA and (b) SFI in different climate groups (different letters indicate significant differences at p < 0.05 in different climate groups.Red bolded numbers indicate the number of soil samples in different climatic groups.WSA is water-stable aggregates, SFI is the soil fertility index, HT-HP is the high temperature and high precipitation group, HT-LP is the high temperature and low precipitation group, LT-HP is the low temperature and high precipitation group, and LT-LP is the low temperature and low precipitation group).

Figure 6 .
Figure 6.Boxplots of WSA and SFI in different climate and plant year groups (different uppercase and lowercase letters indicate significant differences at p < 0.05 across planting years and climate groups, respectively.Red bolded numbers indicate the number of soil samples in different groups, WSA is water-stable aggregates, SFI is the soil fertility index, HT-HP is the high temperature and high precipitation group, HT-LP is the high temperature and low precipitation group, LT-HP is the low temperature and high precipitation group, and LT-LP is the low temperature and low precipitation group).

Figure 6 .
Figure 6.Boxplots of WSA and SFI in different climate and plant year groups (different uppercase and lowercase letters indicate significant differences at p < 0.05 across planting years and climate groups, respectively.Red bolded numbers indicate the number of soil samples in different groups, WSA is water-stable aggregates, SFI is the soil fertility index, HT-HP is the high temperature and high precipitation group, HT-LP is the high temperature and low precipitation group, LT-HP is the low temperature and high precipitation group, and LT-LP is the low temperature and low precipitation group).
Note: SD: standard deviation; CV: coefficient of variation; WSA: water-stable aggregates; /: The mean values for PLC and PRC are negative, which typically means that the CVs are misleading.The abbreviations and descriptions of soil factors, topographic factors, and climatic factors are listed in Table1.

Table 3 .
The grading criteria for SOM, CEC, AVN, AVP, and AVK, and percentage of different grades.The grading criteria are derived from the Specification of Land quality geochemical assessment.b : The grading criteria are derived from the Second National soil census of China.
a :

Table 4 .
The grading criteria for pH and percentage of different grades.

Table 5 .
Principal component analysis for soil fertility-related factors.