Two potential equilibrium states in long-term soil respiration activity of dry grasslands are maintained by local topographic features

Soil respiration of grasslands is spatio-temporally variable reflecting the changing biological activities of the soil. In our study we analysed how the long-term soil respiration activities of dry grasslands would perform in terms of resistance and resilience. We also investigated how terrain features are responsible for response stability. We conducted a 7-year-long spatial study in a Hungarian dry grassland, measuring soil respiration (Rs), soil temperature (Ts) and soil water content (SWC) along 15 measuring campaigns in 80 × 60 m grids and soil organic carbon content in 6 of the occasions. Two proxy variables were introduced to grasp the overall Rs activity, as well as its temporal stability: average rankRs, the temporal average Rs rank of a measuring position from the campaigns revealed the persistent spatial pattern of Rs, while rangeRs, the range of ranks of the positions from the campaigns described the amplitude of the Rs response in time, referring to the response stability in terms of resistance or resilience. We formulated a hypothetic concept of a two-state equilibrium to describe the performance of the long-term Rs activity: Rs activity with smaller rangeRs, that is both the lower elevation positions with larger rankRs (“state I”) and the higher elevation positions with smaller rankRs (“state II”) correspond to an equilibrium state with several terrain attributes being responsible for the equilibrium responses. Majority of the measuring positions was belonging to none of these equilibrium states. These positions showed higher rangeRs for medium rankRs, suggesting resilience (not resistance) as a major strategy for this ecosystem.

The key concepts of ecological stability, such as persistence, resistance and resilience are properties hard to quantify and are always context-dependent 20,21 . Following the concepts found in literature 20-23 , we will use those terms as defined in Fig. 1 (see proxy variables in Methods later on). An ecological system is always exposed to a certain level of disturbance, e.g., related to global changes 22,24,25 . In general, an equilibrium system would respond with different amplitude and response time than a perturbed system 20,22 , whether in nutrient cycling or in community dynamics 26 .
Digital elevation models (DEM) are frequently used not only in landform classification [27][28][29][30][31] or soil mapping 32 but in ecological studies 18,33 as well. Detailed spatial analysis of DEM help to capture relevant information about the terrain surface elements, which can have important ecological effects. The slopes and altitudinal differences can be closely related to surface runoff, water accumulation, snow movement or subsurface biophysical processes 19,31,34 , which influence e.g., vegetation patterns or plant species abundance, diversity and distribution 30,33,34 . The aspect as a measure of slope orientation captures different physical and subsequent biological effects related to predominant wind direction and solar radiation (north-and south-facing slopes differ in the duration of shade, snow cover, vegetation period 34 , or a west-facing slope would be warmer than an eastfacing slope late in the afternoon), both affecting landscape formation and microclimate characteristics. Other terrain attributes like local mean elevation, standard deviation of elevation within a specific area or topographic position index (TPI) used in our study revealed no co-varying aspects 35 of the surface.
Compared to the use of DEM in the above-mentioned studies, the effects of the terrain features on the spatial patterns of grassland soil respiration were scarcely studied. Cultivated areas, grazed or restored grassland vegetation types with different aspects and slope positions have mostly been analysed for soil organic carbon, total nitrogen or other nutrient distribution/accumulation/erosion as well as patterns of above and belowground biomass [36][37][38][39][40] . These studies provided evidence for the effects of these terrain features on the differences in the spatial patterns of the soil nutrients or plant biomass, both influencing R s spatio-temporally.
The complexity of the terrain and the study scale have important consequences on the terrain attributes. An ecological phenomenon and an underlying mechanism can have different spatial scales, as in some cases a neighbouring effect can rather act than a single factor at one particular position 34,41 . Matching scales has to be explorative 31 , since it has to be taken into account that some characteristics disappear at broader scales 29 and that the relative importance of an attribute may change across scales 35 . The picture becomes even more complex if www.nature.com/scientificreports/ we consider the fact that an ecological function can be influenced by different terrain characteristics 42 and that an attribute may be involved in different biophysical effects 31 . In our study we conducted a long-term (7-year-long) spatial investigation in a piece of semi-natural grazed grassland in Hungary. The finely undulating (no more than 1.5 m elevation differences within the study site) surface in our study site is formed through the combined effects of wind, water erosion and drought, resulting in uneven soil nutrient (soil organic carbon, SOC) and water distributions together with different above-ground biomass covers between crests and depressions 18 . We measured R s , T s and SWC along 15 measuring campaigns in 80 × 60 m grids and SOC in 6 of these occasions. Some of these datasets have already been used for detailed geostatistical analysis with a different focus ( 18 : effects of grassland management on CO 2 and N 2 O flux spatial patterns). Our current research question formulated on the basis of spatial samplings (15 occasions) was how soil respiration activity of the grasslands would respond to a range of environmental constraints in terms of resistance and resilience in the longterm, and whether the terrain features were responsible for differences in response stability. We hypothesized on the basis of our previous work 18 that lower micro-elevation levels (surface depressions), rather than the crests, could be responsible for more stable R s activities in general through the effect of more persistent water availability even under drought.

Results
Spatial patterns of stability proxies and background variables. Figure 2 a, b show the spatial distribution of our two proxy variables, the average rank of R s per position (rankR s ) and of the range of the ranks per position (rangeR s ) in kriged maps. The middle to southern areas were found to have the largest, whilst the north-eastern areas the smallest rankR s values, whereas a slightly different pattern was characteristic for rangeR s with some additional north-western large values. Similarly, larger average soil organic carbon content (mean-SOC) and average soil water content (meanSWC) (Fig. 2 c, d) were detected at the western-middle-southern regions and smaller at the north-eastern part of the study site. On the basis of the correlation analysis we found an important difference in terrain attribute features between DEM 5 and 6, especially in SD, Sl, North and East. All subsequent results are then based on DEMs 1-5, which were found to be more similar to each other and to the original DEM1. The maps of terrain attributes with the box blur kernel from DEM1-5 can be found in the Supplementary Information (SI) together with the descriptions and calculations. As we couldn't find any of the blur kernels superior to the other when considering correlations, the results hereafter are only presented for the box blur kernel calculations for simplicity.
When we considered the entire dataset (named hereafter: "A" dataset), we could only find significant correlation between rangeR s and TPI at less smoothed DEMs but the correlation was very weak (black symbols and line in Fig. 3).
Any other correlation between the proxies and the terrain attributes could only be deduced indirectly from the positive correlations between rankR s and meanSOC, meanSWC (cf. Table 1b). These correlations were scaleindependent, i.e., we detected them at every DEMs. In general, the larger the soil carbon content and soil moisture . Direct correlation between TPI and stability proxy, rangeR s at less smoothed DEMs, DEM1-2 for datasets A (black symbols and line) and S (blue symbols and line, see the information later on). The correlations were significant at p = 0.0076 and p = 0.0172 levels, although they were weak, r 2 = 0.09, r 2 = 0.42 for A and S (see the information later on), respectively. Table 1. (a) Statistically significant (p < 0.05) linear correlation between terrain attributes, ALT, mALT, TPI, SD, Sl, North, East and background factors, meanSOC, meanSWC for A dataset. (b) Statistically significant (p < 0.05) linear correlation between terrain attributes, background factors and stability proxies, rankR s , rangeR s in A dataset and in the subgroups (see codes in the text). Regular letters mean scale-independent correlations (valid for DEM1-5), italic underlined letters mean correlations valid for less smooth DEMs (DEM1-2 or 1-3), "pos" and "neg" indicate the sign of the correlation. www.nature.com/scientificreports/ at a position (cf. Figure 2c,d, showing quite similar patterns to the proxy patterns in the figure upper row), the larger the R s activity detected and the opposite was true for lower carbon content and soil moisture positions. As we investigated the background of these correlations more thoroughly in dataset A in terrain attributes (cf. the maps in SI, Table 1a), we found that meanSOC showed negative correlation with ALT, mALT, SD and Sl, while positive with North and East at DEMs 1-5. Similarly, meanSWC correlated negatively to SD, Sl, except for DEM5. Several terrain attributes were then responsible for the meanSOC patterns, i.e., higher absolute elevation and neighbouring surface heterogeneity, as well as steeper slope positions facing more South-West could be characterized with lower meanSOC. The opposite features were characteristic for higher meanSOC level positions on lower elevations with lesser neighbouring heterogeneity and gentle slopes facing mostly North and East. Similarly, meanSWC was higher at smaller surface heterogeneity with more gentle slopes, while lower at more heterogeneous surfaces with steeper slopes. rankR s followed these patterns with higher R s activity in the middle-southern part of the study area, while, in contrast, lower R s activities were characteristic at lower meanSOC and meanSWC at the north, north-east facing locations in the study site, mainly on local ridges as found on the basis of the direct TPI correlations.
Correlations between stability proxies and background variables along DEMs: subgroups. We also checked the correlations within different data subgroups corresponding to specific rankR s or meanSOC categories because we hypothesized that these kinds of groupings could enable us to grasp some important characteristics of the stability.
Direct correlation between terrain attributes and proxies showed considerable variation depending on the subgroups and variables (Table 1b collects  It seems that the meanSOC pattern related negatively to ALT, mALT detected in dataset A could have acted as a driver for the negative rankR s and ALT, mALT correlations in the groups M and C1. It was in subgroup C1 that TPI, SD and Sl acted negatively on rankR s as well, most probably more directly through the patterns generated by terrain attributes in meanSWC. Further negative correlations were found between rangeR s and TPI in S data (see also blue symbols and line in Fig. 2), like in dataset A (cf. Fig. 3 black symbols and line), as well as between rangeR s and North in C4. Accordingly in the long run, local valleys but mostly constant slope positions (with TPI close to zero, cf. blue symbols in Fig. 3) with lower neighbouring surface heterogeneity and gentle slopes with more elevated meanSOC and meanSWC could be characterized with larger R s activity with higher variability (through the negative TPI-rangeR s correlation) in these subgroups per se, similarly to dataset A. The opposite was likely to be the case for local ridges.
The subgroups mentioned here were restricted groups of measuring positions, where carbon content in the soil was the lowest of all (C1) or, as in subgroup S, coincided with low meanSOC levels (Fig. 4), and low meanSWC levels as well.
Furthermore, measuring positions, grouped by either on the basis of rankR s or on the basis of meanSOC, occupied more or less well delimited spatial areas within the sampling grid (Fig. 5, positions coloured according to C1-C5, where e.g., C5 category, indicated with asterisk occupied the lowest altitudinal positions, C4 was found mostly around C5, while C1 category could be found along the edges of the study area on the crests), which would also be characteristically different in respect of the terrain attributes, especially in SD, Sl, as found by the correlations (cf. Table 1).
Finally, rangeR s was fitted to rankR s using the following equation: where µ is the mean rankR s (42.61), σ is the standard deviation of rankR s (25.38), and a (4,291.73) is a model parameter. The correlation, approaching a bell-shaped curve (Fig. 6, both the curve and the model parameters are statistically highly significant, p < 0.0001), visualized together with the subgroups showed that both low and high rankR s could be associated with small rangeR s with larger stability and a typically resistant response, while middle values of rankR s corresponded to larger rangeR s with a more flexible, resilient response of R s activity. Furthermore, it was also showed that rankR s categories more or less fitted to C1-C5 meanSOC categories (cf. square symbols of C1 in S-M-group regions, C2 in M, while asterisks mostly in the upper half of the rankR s range), giving strong evidence of SOC as a controlling factor in R s stability. The smallest and the largest rankR s values could correspond to the largest potential stability (in terms of resistance) in the activities, rankR s being either low in general (cf. Fig. 2a north and north-east regions) due to low meanSOC and meanSWC (Fig. 2c,d) or high (cf. Figure 2a more the middle and southern regions within the study plot) in the opposite cases. On the other hand, medium rankR s with larger rangeR s overlapped spatially with C2-C3 and M groups with medium meanSOC levels, and these positions showed a more resilient response. (1)

Discussion
In our analysis we tried to grasp the persistence, as well as resistance/resilience in soil respiration activity, which is regarded as an important ecophysiological functioning of ecosystems. We used two proxies derived from spatial replicate measurements from several years to analyse stability in R s dynamics. R s has an inherent temporal variability due to its environmental control through spatio-temporally varying and co-varying factors 4,10,43 . R s is also exposed to disturbances, which can be either defined as a "sudden shock" 20 or as a constant disturbance regime like shifts in climatic conditions due to global changes (e.g., nighttime or daytime warming, change in amount or timing of precipitation etc. 22,24 ). Our study site is characteristically exposed to the latter, experiencing frequent droughts and heatwaves in summer demonstrated mostly by NDVI data (Fig. 8a), coinciding with earlier predictions to our East-Central European region 44 . We revealed the spatial pattern of long-term, persistent R s functioning by mapping the average rank of R s in space from a series of measurements. We could conclude that rankR s calculated from 15 measurements, irrespective of the actual environmental conditions and plant biomass levels, enabled us to describe the long-term average functioning, while rangeR s allowed us to have an insight into the dynamics of this process. Spatial patterns have already been found to be similar to a certain degree 15,16,18 but the pattern similarities are assumed to be more qualitative than quantitative 45 in the long run. rankR s showed a general picture of the spatial R s activity but did not provide information about the stability of the process. rangeR s could be more informative in terms of stability. As it corresponds to the amplitude (if there is an alteration) of R s response in time or to the efficiency component of disturbance response 20 it can indeed be defined as a measure of the stability in R s response. Smaller rangeR s means more resistant R s activities against the environmental conditions, mostly drought in our region. Larger rangeR s could correspond to a more flexible response, reflecting resilient R s activities.
However, our results based on the correlations of our proxies and the terrain attributes as well as on the bellshaped curve relation between rankR s and rangeR s suggest the existence of two types of resistant responses with smaller rangeR s : one at lower elevations with high meanSOC, meanSWC at gentle slopes with large rankR s , and another at higher elevations, at local ridges with lower levels of meanSOC and meanSWC, at steeper slopes and www.nature.com/scientificreports/ increased surface rugosity with lower rankR s . The first type of response coincided with our hypothesis based on our previous work 18 that R s activity would be more stable at lower elevations where water availability is more constant. However, the other type was a new finding in our system but theoretically equally meaningful because complex systems can have several equilibrium states 19,22 .
In order to characterize the complexity of our system we can rely on the terrain attribute analysis. It was finally restricted to DEMs 1-5 as DEM smoothness was found to have significant effects on terrain attributes 29,46 ,  www.nature.com/scientificreports/ slope, aspect, curvature, size of catchment area, which were all found to be affected by DEM resolution but also by the vertical precision. We found significant variability of the terrain attributes within our study area (cf. SI Figs. 1-3) which implies substantial spatial differences both in the environmental conditions and in community structure 30,33,34 . Thus, this spatial variability in the conditions may simultaneously cover the effects the system generally faces and those experienced due to the disturbances. Slope differences were found to be responsible for the water regime in general 31,34 , while aspect for incoming solar radiation and for the surface formation by wind 34 . Finally, we came to the conclusion that the equilibrium state in our system was dynamic 22 . Todman et al 22 stated that several "smaller-scale domains of attraction" could exist in complex systems. We concluded that both lower elevation positions with larger rankR s ("state I") and higher elevation positions with smaller rankR s ("state II") but both characterized with smaller rangeR s could correspond to an equilibrium state. This theory can also rely on the observations that wet and dry soil moisture patterns with transitory phases between them characteristically occur 45,47 . Wet state is a result of non-local forces, acting on excess water supply, while dry state is locally driven by soil properties, incoming radiation and vegetation 47 . Although semiarid regions are in the dry state most of the time 45 , our C5 positions could correspond to a generally wetter or at least more transitory state which has a generally larger soil moisture variability with some intervention of the above-mentioned non-local forces. C1 and S are typically more locally controlled especially if we consider the importance of TPI in these places: local ridges are the most exposed to sunlight and evapotranspiration is strong. These differences in water availability between the two equilibrium states together with meanSOC differences are well demonstrated in Fig. 4 showing that generally C1-C5 and S-M-L subgroups experience different levels of meanSWC and meanSOC.
On the basis of these observations we attempted to formulate a concept (Fig. 7) concerning the stability of R s activity by trying to grasp the terrain features and background factor characteristics as surrogates 31 .
Several terrain attributes were responsible for the meanSOC, meanSWC patterns which were found to be the direct spatial drivers of R s activity. Higher absolute elevation (mALT) and neighbouring surface heterogeneity (SD), as well as steeper slope (Sl) positions ( Fig. 7 right part: "State II") could be characterized with lower meanSOC, meanSWC, mostly on local ridges together with a more resistant response, while the opposite features ( Fig. 7 left part, "State I") were characteristic of higher meanSOC, meanSWC level positions with lower elevations (mALT) with lesser degree of neighbouring heterogeneity (SD) and gentle slopes (Sl), but also with a more resistant response through a generally larger R s activity. Intermediate levels of the background factors and no specific terrain features characterized the positions (Fig. 7, middle part), where a more resilient R s response was detected.

Methods
Site description. The study site is in the Kiskunság National Park, near Bugac (46.69° N, 19.6° E, 114 m a.s.l.), according to a long-term research permission (Grant No: 60960-1-11/2015). The vegetation, which is a semi-arid sandy grassland, is dominated by Festuca pseudovina Hack. ex Wiesb., Carex stenophylla Wahlbg. and Cynodon dactylon L. Pers. The prevailing wind direction is North-West and the surface is slightly undulating. The mean annual precipitation and temperature was 585 mm and 10.6 °C, respectively, for 15 years following the establishment of the eddy-covariance station in 2002. According to the FAO classification 48 the soil type is Chernozem with a relatively high organic carbon content, the soil texture is a sandy loam with a sand:silt:clay ratio of 81:11:8% in the topsoil layer 49 . The study plot had been under extensive grazing for decades. Grazing intensity was 0.66 ± 0.18 Hungarian Grey cattle animal ha −1 year −1 for the previous few years. The grazing period usually lasted from May to June and from late August to the end of November. The grassland may potentially turn into a source of carbon in drought years 50  Soil respiration (R s , µmol CO 2 m −2 s −1 ) was measured by means of closed chamber systems (Licor 6,400, LiCor, Inc. Lincoln, NE, USA and EGM-4 PPSystems, Amesbury, USA) at 78 sampling locations (arranged as a 80 × 60 m grid 18 ) in each measurement campaign. Target CO 2 concentration was set by placing the soil chamber on its side on the ground to monitor the CO 2 concentration over the surface. Collars were not used with the soil gas exchange chambers (area of 78.54 cm 2 ) to minimize disturbance 53,54 since both measuring systems performed well without collars 55 . Although the sampling positions remained relatively constant for the duration of the experiment, a shift of a few centimetres was applied when selecting the actual patch for each measurement. The standing biomass was removed 1.5 h before starting the soil respiration measurements. To minimize the effects of diurnal temporal patterns the measurements were started at noon and lasted ~ 1.5 h for the grid.
Soil water content (SWC, %) was measured at the same spots as the CO 2 gas fluxes by time domain reflectometry (ML2, Delta-T Devices Co., Cambridge, UK; FieldScout TDR300 Soil Moisture Meter, Spectrum Technologies, IL-USA) in the top 0-5-cm layer of the soil. The measurements were performed usually after the R s measurements in all positions in one run. Soil temperature (T s , °C) was determined at a depth of 5 cm by a digital soil thermometer near the R s chambers parallel with the R s measurements. The soil organic carbon content (SOC, %) of the bulked soil samples from the upper 10 cm was determined by sulfochromic oxidation/loss on ignition. environmental conditions during the study period. Meteorological data were available from the eddy covariance system functioning at Bugac continuously from 2002. The yearly average air temperature, sum of precipitation and NEE data for the study period are shown in Table 2.  www.nature.com/scientificreports/ Annual precipitation sum was lower by 27% in the driest (2012) and higher by 39% in the wettest (2014) year of the study period than the fifteen-year average (585 mm). This variability in water availability resulted in a source activity of + 37.9 g C m −2 year −1 in 2012, while in a general sink activity between −38.4 and −79.8 g C m −2 year −1 for the period 2013-2018.
Daily average temperature, precipitation sum and broad-band normalized difference vegetation index (NDVI) 50 are presented in Fig. 8a with the measuring campaigns for the entire study period between 2012 and 2018. The first measuring campaign was usually scheduled for the spring-early summer active periods, the second was in the period of summer drought, while the third was in autumn, often along a re-greening period. Several very intensive precipitation events could be distinguished. The actual R s of the measuring campaigns covered wide ranges of the potential activity with different variability from time to time due to a corresponding variability in SWC and T s (Fig. 8b).
Stability proxies for resistance/resilience interpretation. In order to gain an insight into the stability of the R s patterns we first calculated the average rank of soil respiration (rankR s ) for the measuring positions on the basis of the full dataset. This way we could grasp the long-term, persistent distribution of the spatial positions with low or high R s activities. The calculation relied on ranking the measured values from the 1st (the smallest value) to the 78th (the largest value), then averaging the 15 campaigns for each position. Small rankR s value corresponds to a generally low R s activity in the position, while large rankR s value means a generally more significant R s activity in that particular position. Another proxy was rangeR s which we calculated as the difference between the maximum and minimum rank (rankR s ) for a position, giving low values with more stable patterns, because the R s rank was found to be more constant for the duration of the measuring campaigns, i.e. these positions showed considerable resistance against environmental constraints. Large rangeR s values referred to more variability in R s activity as well as more resilient functioning.
DEM processing and terrain attribute calculations. Digital elevation models (DEM) contain bareearth elevation data in raster grids. Raster data can be processed by different image processing operations, which serve to extract information from the objects, which is a reasonable approach to capture relevant structures when performed at different scales. Pyramid image processing is a multi-scale approach when the image is processed by smoothing and subsampling steps in several runs. We used a specific case of pyramid image processing along the terrain attribute calculations, called the "mixed scaling" 41 , when DEM is smoothed and subsampled www.nature.com/scientificreports/ and the calculated terrain attributes are upscaled. This method was found to result in fewer artifacts and more straightforward patterns than other techniques 41 .
In our study, we implemented 56 the method on a 0.2 by 0.2 m input DEM raster (originating from laser scanning) as follows: 1. DEM was progressively smoothed, i.e., aggregated by a factor of two, resulting in six different resolution DEM rasters along a geometric series between 0.2 and 6.4 m (DEM1: 0.2 m, DEM2: 0.4 m, DEM3: 0.8 m, DEM4: 1.6 m, DEM5: 3.2 m, DEM6: 6.4 m), and another scale was also calculated to meet the resolution of the measuring campaigns (10 m, DEM7) 2. Terrain analysis was performed on each of the DEMs, giving firstly a series of terrain attribute rasters with different raster cell sizes/resolutions 3. Terrain attributes were then disaggregated, or upscaled to the original resolution of DEM1, resulting in differently smoothed terrain features with the same raster cell sizes Six terrain attributes were calculated following the guidelines 35 for the best characterization of the surface with the least potential co-variance between the selected attributes and the ones which were found to be applicable for a range of terrain complexities: Details about the calculations can be found in SI.
Spatial data processing. In order to visualize stability proxies we performed variography and kriging.
The steps of the spatial data processing, detailed description of variography, kriging and leave-one-out crossvalidation can be found in the SI (as well as some results, which are only background information for the main focus of the present article). In brief, we performed variogram analysis first 57-62 on rankR s , rangeR s , SWC and SOC data. The criterions for variogram model selection were the residual sum of squares (model with maximum SSErr from exponential, Gaussian and spherical), the Nash-Sutcliffe model efficiency coefficient (E > 0.5), and the range of autocorrelation, a (the best fit with a within the spatial scale of the study site, i.e., a < maximum distance of the diagonal of the rectangle that spans the data locations). Two kinds of kriging methods were used for mapping the variables, ordinary kriging (OK), and kriging with external drift (KED). Kriging results were evaluated by means of the leave-one-out cross-validation method 58 , and as the error estimates for OK and KED didn't show important differences (in terms of cross-validation errors, i.e., normalized root mean squared error (nRMSE), mean error (meanErr) and mean squared deviation ratio (MSDR)) but we lack OK map for rangeR s because we lack valid variograms, the presented maps are KED maps.