Co-varying effects of vegetation structure and terrain attributes are responsible for soil respiration spatial patterns in a sandy forest–steppe transition zone

Forest–steppe habitats in central Hungary have contrasting canopy structure with strong influence on the spatiotemporal variability of ecosystem functions. Canopy differences also co-vary with terrain feature effects, hampering the detection of key drivers of carbon cycling in this threatened habitat. We carried out seasonal measurements of ecosystem functions (soil respiration and leaf area index), microclimate and soil variables as well as terrain features along transects for 3 years in poplar groves and the surrounding grasslands. We found that the terrain features and the canopy differences co-varyingly affected the abiotic and biotic factors of this habitat. Topography had an effect on the spatial distribution of soil organic carbon content. Canopy structure had a strong modifying effect through allocation patterns and microclimatic conditions, both affecting soil respiration rates. Due to the vegetation structure difference between the groves and grasslands, spatial functional diversity was observed. We found notably different conditions under the groves with high soil respiration, soil water content and leaf area index; in contrast, on the grasslands (especially in E–SE–S directions from the trees) soil temperature and vapor pressure deficit showed high values. Processes of aridification due to climate change threaten these habitats and may cause reduction in the amount and extent of forest patches and decrease in landscape diversity. Owing to habitat loss, reduction in carbon stock may occur, which in turn has a significant impact on the local and global carbon cycles.

1 Introduction

Due to global climate change, extreme weather events (i.e., heat waves and extreme rainfalls) frequently occur worldwide. They are expected to have a major impact on the structure of forest ecosystems and on the dynamics of environmental factors and functional variables in forest patches and nearby open areas (Chatterjee and Jenerette, 2011; Cunningham et al., 2006; Shi et al., 2011). Transition zones (e.g., edges, forest patches with grasslands) with distinct vegetation physiognomy and plant species compositions could be sensitive habitats (Erdős et al., 2014), as adjacent communities may respond differently to extreme events, depending on their sensitivity and stability. Despite this sensitivity, the ecosystem functions (e.g., soil respiration) of these transition zones are less studied among temperate ecosystems (Erdős et al., 2018).

Soil respiration (R_s) can be considered to be one of the most important functional variables of ecosystems, consisting of respiration by various organisms and showing the overall biological activity of the soil (Jassal and Black, 2006; Shi et al., 2011). The two main components of R_s are autotrophic respiration (R_aut) by plant roots and root-associated microbes and heterotrophic respiration (R_het) by microbes and animals in the soil (Balogh et al., 2016).

Soil factors generally influence the temporal and spatial variance of R_s rate, e.g., soil organic carbon content (SOC), soil physical and chemical properties (texture, bulk density, pH, etc.). The density and activity of plant roots through R_aut as well as the activity of soil microbial biomass through R_het have great influence on R_s rate too (Lellei-Kovács et al., 2016; Michelsen et al., 2004; Mitra et al., 2019; Thomas et al., 2018). The share of R_s components is largely dependent on these variables (Hu et al., 2001; Moyes and Bowling, 2016; Raich and Tufekcioglu, 2000; Tang et al., 2020), responding to the changes of environmental factors differently (Balogh et al., 2016; Shi et al., 2011). Among the soil parameters, R_s shows a positive temporal relationship with soil temperature (T_s) and soil water content (SWC) (Balogh et al., 2015; Lellei-Kovács et al., 2016; Moyes and Bowling, 2016; Shi et al., 2011; Tang et al., 2020). Thus, seasonal changes in soil parameters play a significant role in defining seasonal differences in soil CO₂ emission (Hao et al., 2010; Herbst et al., 2009; Mitra et al., 2019; Raich and Tufekcioglu, 2000; Thomas et al., 2018).

SWC and T_s strongly influence the growth and distribution of plant species as well (Cuena-Lombraña et al., 2018; Latif and Blackburn, 2010; Potter, 2014). In addition, the vegetation composition and structure are crucial determinant factors of R_s rate by influencing canopy microenvironment (Chatterjee and Jenerette, 2011; Petrone et al., 2008), soil microclimate and soil organic matter through plant residue input (Stoyan et al., 2000). Consequently, transitions in the vegetation types cause sharp changes in the spatial and temporal variability of the canopy's carbon assimilation, governing allocation patterns and litter production, thus resulting in different levels of stored terrestrial carbon belowground (Chen et al., 2011). Several previous studies (Huang et al., 2020; Michelsen et al., 2004; Shi et al., 2011; Thomas et al., 2018) showed R_s dissimilarities between different vegetation types (e.g., woody vegetation vs. grassland). Grasses allocate greater amounts of carbon for growth and typically produce a larger amount of readily decomposable litter, while trees generally allocate more carbon to structural tissues and produce smaller amounts of longer-lasting, nutrient-poor litter (Moyes and Bowling, 2016; Raich and Tufekcioglu, 2000). Related to relative carbon allocation, R_s rates were usually found to be lower in grasslands than in forests (Moyes and Bowling, 2016; Thomas et al., 2018). Woody vegetation has the largest carbon stocks among vegetation types; thus the transition of woody vegetation to grasslands via natural (i.e., forest–steppe mosaics) or unnatural (deforestation) mechanisms is of great importance. This type of land cover change was considered to be one of the most important drivers of global R_s changes (Huang et al., 2020).

Forest–steppe habitats in central Hungary are naturally fragmented transition zones, where a mosaic of forests with grasslands can develop, as neither of them has a decisive advantage over the other (Erdős et al., 2018). This type of ecosystem has a very sharply separated structure between the forest patches/groups of trees and the open grasslands, which are in close contact with each other through biotic and abiotic factors (Süle et al., 2020). The strong abiotic (i.e., soil water content, SWC) differences in the edges prevent the forest's expansion and development by constraining its growth (Erdős et al., 2014). The sandy forest–steppe vegetation type is in danger of total disappearance in the near future due to the aridification of this region (Pongrácz et al., 2014). Hence, improving our knowledge about ecosystem functioning in this habitat could provide new insights into the effects of global change and a key to preserving these habitats and maintaining landscape diversity (Erdős et al., 2014, 2020; Morecroft et al., 1998).

The vegetation patterns in open landscapes such as sandy grasslands are often influenced by elevation and environmental factors determined by micro-topography. Surface altitudinal heterogeneity not only is a key factor in determining not only the spatial patterns of SWC and T_s but also has a strong connection with the microclimate above the soil surface. Actually, the abiotic factors of the landscape/habitat related to terrain attributes (e.g., slope, altitudinal differences, aspect) and the subsequent distribution of plant species (Cuena-Lombraña et al., 2018; Hohnwald et al., 2020; Latif and Blackburn, 2010; Potter, 2014) mutually influence each other (Alexander et al., 2016). The slopes and altitudinal differences can be closely related to water accumulation and surface runoff, which were found to influence functional spatiotemporal patterns of the vegetation (Fóti et al., 2020; Lassueur et al., 2006). Slopes with different orientation are exposed to different physical effects (e.g., predominant wind direction and irradiation). The main difference between north- and south-facing slopes resides in their shade duration, due to the different amounts of incoming solar radiation (Alexander et al., 2016; Fóti et al., 2020), and the west side of a grove could be colder than the east side during a diurnal period (Süle et al., 2020). This variability in shading causes significant variation in plant composition (Erdős et al., 2014), with the subsequent differences in microclimate giving evidence of a mutual spatiotemporal effect. Although the abovementioned publications show evidence of terrain attributes' effects, little is known about how a tree group acts on their covariation with shading effects.

The most important microclimatic components usually examined with R_s variability include air temperature and air humidity (Matlack, 1993; Potter, 2014) and vapor pressure deficit (VPD) (Bolton, 1980; Chen et al., 1995; Süle et al., 2020). VPD can be an important limiting factor in plant growth because above a certain threshold photosynthetic CO₂ uptake is strongly limited due to stomatal closure to prevent water loss (Novick et al., 2016; Shamshiri et al., 2018; Young and Mitchell, 1994). In our previous study (Süle et al., 2020) we found that at a 3 kPa threshold value of VPD species composition differences between the grove and the adjacent grassland areas were accompanied by differences in their VPD distribution (Süle et al., 2020).

The objectives of this study were (1) to examine the spatial variability of the functional responses of a grove with the surrounding grassland in terms of R_s and canopy characteristics in different phenological stages as well as (2) to link these responses to the background factors given by the spatially co-varying effects of the habitat's terrain features and the grove. Our main question was as follows: which are the main influencing factors of soil respiration in an ecosystem with such a large vegetation structure heterogeneity?

2 Materials and methods

2.1 Study area

Our study sites were situated in central Hungary (Fülöpháza region of the Kiskunság National Park). The 10-year (2010–2020) temperature and precipitation averages of this region were 11.5 ^∘C and 622.17 mm, respectively (National Meteorological Service, Fülöpháza Meteorological Station). Humus-poor sandy soils are characteristic in this region, which have extremely low water-holding capacity and soil organic matter content (Erdős et al., 2014). Due to these poor soil features, precipitation drains quickly, so the amount of water in the soil available for plants is very small. Two groves of poplar (Populus alba) and their surrounding grasslands were selected for this study (Fig. 1) to represent the sharp vegetation structure differences of the forest–steppe habitat. Relief maps as well as sampling positions of the first grove ($\mathrm{46}{}^{\circ }{\mathrm{53}}^{\prime }{\mathrm{28.18}}^{\prime \prime }$ N, $\mathrm{19}{}^{\circ }{\mathrm{24}}^{\prime }{\mathrm{46.91}}^{\prime \prime }$ E, 107 m a.s.l) and second grove ($\mathrm{46}{}^{\circ }{\mathrm{53}}^{\prime }{\mathrm{06.11}}^{\prime \prime }$ N, $\mathrm{19}{}^{\circ }{\mathrm{24}}^{\prime }{\mathrm{29.17}}^{\prime \prime }$ E, 106 m a.s.l.) can be seen in Fig. 2. Measurements of environmental and functional variables were performed in intersecting transects with different cardinal directions, forming together a star and cross-shape sampling arrangement with the group of trees in the middle (Fig. 2, Table 1).

Figure 1 Study sites with the first grove (a) and the second grove (b).

Figure 2 Digital elevation models and sampling arrangements of the sites (a: first grove in 2018; b: second grove in 2019; c: first grove in 2020; white dashed line: visual tree edges; black dots: sampling positions). Coordinates refer to the Hungarian Unified National Projection system and are given in meters.

Table 1 Measurement campaigns during the study.

2.2 Measurements of soil parameters and functional variables

Measurements of soil parameters and functional variables were performed along transects with 2 m intervals between the measurement positions. R_s and T_s were measured by an EGM-4 infrared gas analyzer (PP Systems, UK). SWC was measured by a FieldScout TDR 300 (FieldScout, USA) with a 3.0 in. (7.62 cm) rod. Leaf area index (LAI) was measured by an ACCUPAR LP-80 ceptometer (METER Group, USA). Manual measurements were taken near noon and lasted about 1.5 h (Fóti et al., 2018). All of the statistical evaluation was performed in R (R Core Team, 2018). We used standardization (Hmisc package: Harrell and Dupont, 2020) on the variables mentioned above for principal component analysis (PCA biplot; vegan package: Oksanen et al., 2019).

To separate R_het and R_aut components, we calculated the mean R_s from the whole dataset of each measuring campaign. From the LAI–R_s linear regression, we determined the y-axis intercept (where LAI is theoretically 0) for each measurement campaign, considering it as the average R_het rate. Therefore, we calculated R_aut rate as the difference between mean R_s rate and the R_het rate.

In each measurement year, bulk soil samples were taken from the upper 10 cm at the measurement positions along the transects. SOC (%) of sieved soil was determined by sulfochromic oxidation and loss on ignition (Ivezić et al., 2016).

2.3 Microclimate measurements and vapor pressure deficit

Air temperature and air humidity were measured in the herb layer with a sensor network along the transects at 2 m intervals for 48 h (1 min resolution) during each measurement campaign. The data loggers were placed 20 cm above the soil surface, at the average height of the herbaceous vegetation. Crossbow MICA XM2110CA mote (Crossbow Technology Inc., Milpitas, CA, USA), UNI-T UT330B Mini USB temperature humidity logger (UNI-TREND Technology CO Ltd., Guangdong, China) and Voltcraft DL-120TH USB temperature humidity logger (Voltcraft, Hirschau, Bavaria, Germany) were used. The sensors were shielded with a white plastic plate to avoid direct solar radiation. The sensors were calibrated before the measurements. We selected precipitation-free measurement periods, but it rained during the measurements in October 2020. The main changes in the weather conditions were described by observation (e.g., clouds' shading and movement) and recorded manually. The location of the visual edges of the grove, the positions of bushes and trees in the surrounding area, and the shadow of the groves were described by observation and were also recorded manually.

VPD was computed (Hmisc package: Harrell and Dupont, 2020) from the relative air humidity (RH) and air temperature (t) according to the formula developed by Bolton (1980):

$$\begin{array}{}\text{(1)}& \text{VPD}=(\mathrm{100}-\text{RH})\times \mathrm{6.112}\times {\mathrm{e}}^{\mathrm{17.67}\times t/(t+\mathrm{234.5})},\end{array}$$

with t in degrees Celsius, RH in percent, and VPD in pascals.

We used the average of per minute VPD in the 11:00–13:00 LT period for this study to examine the same period covered by the other measurements during the day. We also used standardization (Hmisc package: Harrell and Dupont, 2020) on VPD for principal component analysis (PCA biplot; vegan package: Oksanen et al., 2019).

2.4 Measurement with GPS, digital elevation model (DEM) processing and terrain attribute calculations

Coordinates and altitude of measuring positions along the transects were recorded by high-precision STONEX S8 PLUS GPS (STONEX Srl., Paderno Dugnano, Italy).

We used spline interpolation (akima: Akima et al., 2016, and field packages: Nychka et al., 2017) of the measured altitude data to generate DEMs of 0.4 m resolution for the groves, which could serve as input raster data for terrain attribute calculations. Although DEMs in general are built on the basis of data collected by remote sensing techniques, they may equally be built from land surveying by GPS receivers. These interpolation maps can be seen in Fig. 2 together with the measuring positions. The maps lie on 183 (first grove) and 106 (second grove) randomly measured data, respectively.

Four terrain attributes were calculated 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. For a more detailed description of the calculations please refer to Fóti et al. (2020). We used RSAGA (Brenning et al., 2018) and raster packages (Hijmans, 2018) for the terrain attribute analyses mentioned below. We calculated Spearman correlation coefficients (Hmisc package: Harrell and Dupont, 2020) between PCA scores of the biplots mentioned earlier, SOC and these terrain attributes.

Standard deviation of elevation (SD) was calculated by “mixed scaling” (see Behrens et al., 2018) for a smoothed DEM of 1.5 m (to approximately meet the sampling resolution) by 5×5 box blur kernels as the neighborhood as follows:

$$\begin{array}{}\text{(2)}& \text{SD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n_R-\mathrm{1}}}\sum _{i=\mathrm{1}}(z_i-\stackrel{\mathrm{‾}}{z}{)}^{\mathrm{2}}},\end{array}$$

where z_i values are the elevations of the correspondent R radius while $\stackrel{\mathrm{‾}}{z}$ is the mean elevation within R. SD describes the heterogeneity and local surface roughness within the raster.

For slope- and aspect-derived easterness and northness, we used eight neighbors, as suggested by Lecours et al. (2017):

Slope (Sl) expresses the rate of change in elevation between positions; it is the tangent (vertical “rise” divided by horizontal “run”) of a surface angle to the horizontal in degrees.

Easterness and northness (east, north) aspect is the compass direction that a slope faces, derived from the maximum downslope rate of change in a value from a raster cell to its neighbors. It is a circular variable ranging clockwise from 0 to 360^∘ from the north (both 0 and 360^∘ meaning N-facing slope, 90^∘ meaning E-facing slope; Ritter, 1987). In ecology, use of sine (easterness) and cosine (northness) of aspect is more frequent because they provide a continuous gradient of east–west and north–south directions, respectively. Northness and easterness with the values close to +1 mean that the slope is northward and eastward in general, while values close to −1 mean a generally south- and west-facing slope.

3 Results

3.1 Measuring circumstances (National Meteorological Service, Fülöpháza meteorological station)

Daily average temperatures and daily precipitation sums are presented in Fig. 3 for the study period between May 2018 and October 2020. The first measuring campaign was usually scheduled for spring, the second was in mid-summer and the third was in autumn. We present the summary of the preceding meteorological conditions of the measurement campaigns in Table 2. During spring 2018 and all of the 2020 measurement campaigns it had rained before the measuring day. The amount of daily precipitation observed on the day of the measurement in summer 2020 fell during the previous night, but no precipitation occurred during the measurements. During the autumn 2020 measurement campaign, the global radiation was very low compared to the other measurement periods.

Figure 3 Daily averages of air temperature and daily sums of precipitation during the study period with the measuring days (black arrows), spring 2018–autumn 2020, at the Fülöpháza site.

Table 2 Averages of air temperature, sums of precipitation and averages of global radiation at the measuring campaigns and 3 d before them.

3.2 Relationships between the measured variables

R_s showed significant positive correlation with SWC (P<0.001, R²=0.13) (Fig. 4a) and LAI (P<0.001, R²=0.28) (Fig. 4c) and showed negative correlation (P=0.059, R²=0.004) with T_s (Fig. 4b). The highest R_s values were observed in summer 2019 (23.99 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). In general, R_s was the lowest (0.13 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) in the autumn phenological stages at relatively low temperatures (13.2 ^∘C) and low LAI (0.15 m² m²). No significant relationship was found between T_s and R_s (Fig. 4b) due to the high temperatures measured in summer 2018 coupled with low water availability. LAI showed significant negative correlation (P<0.001, R²=0.08) with VPD (Fig. 4d). No separated data group could be found in these regressions (Fig. 4), except autumn 2020, which was completely separated from the other phenological stages in Fig. 4d. However, characteristic groupings can be seen by phenological stages in Fig. 4b, pointing to the differences in the environmental conditions during measuring campaigns. The points of Fig. 4b and d have a hyperbolic shape.

Figure 4 Relationships between variable pairs of the measurement campaigns. Pairs: soil respiration vs. soil water content (a), soil respiration vs. soil temperature (b), soil respiration vs. leaf area index (c) and leaf area index vs. vapor pressure deficit (d).

3.3 R_s component rates

A characteristic trend can be observed in the R_s component rates of the measurement campaigns (Table 3). Within a year, spring and summer phenological stages were similar to each other, with a slightly higher autotrophic respiration component in summer, while this component was the lowest in autumn. The year 2019 had higher values than the other two surveyed periods. R_aut seemed to closely follow LAI.

Table 3 Rates of soil respiration and its components with leaf area-index during the measurement campaigns.

3.4 Connections between the variables and vegetation physiognomy

PCA biplots (Fig. 5) were used to analyze relationships among variables. Based on this analysis, in every measurement campaign the positions below the groves were notably separated from the positions in the surrounding grasslands, especially the positions lying eastward and southward from the trees. However, in the autumn phenological stage, the positions of the two vegetation types became less separated, especially in the first 2 years, when grassland positions north of the trees overlapped remarkably in PCA space with below-canopy positions. Among the studied variables, R_s showed a positive relationship with SWC and LAI, while it showed a negative relationship with T_s and VPD, as shown by PCA loadings along the first axis. The second axis separated SWC, R_s and LAI in any potential combination, reflecting the same underlying causes with some variation in the actual signs of the loading, although the most frequent case was in which SWC and R_s were separate from LAI. Mainly, the loadings of R_s–SWC–LAI correlated positively with the distribution of the new PCA scores originating from the positions of the groves, while the distribution of the grasslands' positions followed T_s–VPD loadings, especially those that can be considered warmer, including the E, SE and S directions from the trees, except for at the 2020 autumn phenological stage due to loading of SWC.

Figure 5 PCA biplots of the measurement campaigns with separation of the measuring positions below (open square symbols) and surrounding (closed spots with cold–warm colors according to the cardinal directions) the groves. Variables: soil respiration (R_s), soil water content (SWC), leaf area index (LAI), soil temperature (T_s) and vapor pressure deficit (VPD). Cardinal directions: west (W), northwest (NW), north (N), northeast (NE), east (E), southeast (SE), south (S) and southwest (SW).

3.5 Correlations between the variables, soil organic carbon and terrain attributes

Soil organic carbon content showed strong positive correlation with the PCA I axes of Fig. 5 of every measurement campaign, except in autumn 2020, where it showed strong negative correlation with the PCA I axis and positive significant correlation with the PCA II axis (Table 4). This difference was in good agreement with the orientation change of the loadings in Fig. 5 (autumn 2020 biplot). In the campaigns of 2018 and 2020 eastern orientation had a significant negative correlation with PCA I axes (and the opposite again in the case of autumn 2020), similarly to the findings in the biplot analysis, with “warmer” grassland positions matching T_s–VPD loadings the most. In addition, Sl and SD showed negative correlations with PCA II axes in the 2018 spring phenological stage, which means that larger R_s and SWC were measured at smaller elevation differences with slight slopes, although this effect was not detectable in the measuring campaign of 2020. This was probably due to the rainy days before the measuring occasion (see Table 2), masking the effects of the surface. At the study site in 2019, north showed a negative correlation with PCA I axes on both occasions, while Sl and SD showed a negative correlation with PCA I axes in spring and a positive correlation with PCA II axes in autumn. These correlations again suggest larger T_s–VPD values at south-facing directions and also larger R_s and LAI values in those positions with less surface heterogeneity.

Table 4 Spearman correlation coefficients between PCA scores of Fig. 4, soil organic carbon (SOC) and terrain attributes (${}^{\ast \ast \ast }$: P<0.001, ${}^{\ast \ast }$: P<0.01, ^∗: P<0.05).

4 Discussion

The vegetation composition and horizontal and vertical structures of ecosystems are crucial determinant factors of the environmental and functional variables influencing, among others, the air and soil microclimate and soil chemical parameters (Huang et al., 2020; Mitra et al., 2019; Thomas et al., 2018). The structural differences in the vegetation types cause spatiotemporal variability in the assimilation and allocation patterns, with subsequent effects on the soil carbon stocks.

4.1 The spatiotemporal explanatory factors of R_s

The responses of R_s to soil moisture and T_s were investigated in several studies, which already reported that the correlation between R_s and T_s as well as SWC and T_s differed in temporal and spatial studies (Allaire et al., 2012; Almagro et al., 2009; Fóti et al., 2020; Herbst et al., 2009; Lellei-Kovács et al., 2016; Moyes and Bowling, 2016; Shi et al., 2011; Tang et al., 2020). Our study is consistent with these findings, because R_s showed positive correlation with SWC and LAI. Although it was often found that there is a positive correlation between R_s and T_s in temporal datasets (Balogh et al., 2019; Savage et al., 2013), based on our entire dataset with sharp vegetation structural differences (Fig. 4b), we could not detect this kind of relationship. However R_s showed negative correlation with T_s and VPD in space by measuring campaigns (Fig. 5) similarly to the findings of other spatial studies (Fóti et al., 2014; Herbst et al., 2009). VPD showed negative correlation with LAI both in time and in space and very often had the opposite loading of SWC, R_s and LAI in the PCA analysis of spatial data. The below-canopy area had higher R_s and LAI with lower air temperature and T_s due to the microclimate below the complex vegetation structure and higher biomass with trees and bushes. In contrast, the grassland had higher T_s and VPD with lower R_s due to the higher irradiation (Huang et al., 2020; Moyes and Bowling, 2016; Thomas et al., 2018).

4.2 The components of R_s

We found seasonal differences in R_s and component rates within a year. In sandy forest–steppe habitat, which is an ecosystem with a semiarid climate and predominantly open vegetation with low biomass (Erdős et al., 2020), the autumn phenological stage may already be a dormant period. It contrasts with other dry sandy grasslands (Balogh et al., 2016) and mesic meadows with trees (Moyes and Bowling, 2016), which have higher biomass and higher organic carbon content. In a closed dry sandy grassland (Balogh et al., 2016) there could be an active re-greening period with higher R_s rates in autumn than in the summer. Besides, the correlation between R_aut and the seasonal variation in LAI reflected the effect of vegetation structure on the ecosystem functioning and revealed the strong modifying effect of plant CO₂ uptake on soil respiration.

4.3 The influence of cardinal directions and topography

Several studies referred to north–south as the coolest–warmest temperature gradient (Erdős et al., 2014; Hohnwald et al., 2020; Latif and Blackburn, 2010; Matlack, 1993) and also highlighted that northern areas of edges and the surrounding open areas showed similarities in microclimate and vegetation composition with the forest interior (Erdős et al., 2014). In our study, the reason for the grouping of grassland positions is the difference in T_s and VPD values, which can be traced back to the temperature gradient. This also caused the scatter of positions of the grassland side north of the trees near the grove.

In addition, the warmer areas matched with the easterness terrain attribute, and the T_s–VPD variables aligned with the south-facing slope terrain attribute. Also, the spatial distribution of SOC followed the spatial change of vegetation structure because the grove had higher carbon content. Our results were in good agreement with the findings that woody vegetation has the largest carbon stock among the vegetation types (Huang et al., 2020; Moyes and Bowling, 2016; Tang et al., 2020); thus the transition of woody vegetation to some short vegetation type, for example in forest–steppe mosaics, is of great importance. On the other hand, surface heterogeneity was lower under the grove due to the fact that the seeds of tree species in this habitat can only germinate in ideal depressions, where a sufficient amount of water can accumulate (Erdős et al., 2014). In the context of spatial heterogeneity, deeper positions contain more carbon due to accumulation, resulting in higher R_s activity and SWC (Fóti et al., 2020). These parameters indicated the functional differences within the study sites, mostly between the grove and the surrounding grasslands but also between the grassland areas based on their cardinal directions.

4.4 The effect of the heterogeneous vegetation

We found that the groves were separated from the surrounding grasslands based on the measured variables, except in the autumn phenological stage due to the lower irradiation. In this phenological stage, the grassland heated up less, T_s was lower and soil moisture evaporated at a lower rate. The shade of the group of trees also lasted longer on the grassland in autumn, yet a less marked difference could be observed between the grassland and the grove. This is in good agreement with the results of the study conducted by Thomas et al. (2018), where solar radiation, air and soil microclimate parameters showed the same trend during the phenological stages of a year due to the shade of the trees. A similar seasonal trend was also observed by Moyes and Bowling (2016) in the case of T_s and SWC. Opposite orientations of R_s–SWC–LAI and T_s–VPD loadings along PCA I axes in the 2020 autumn phenological stage can probably be attributed to the low global radiation on the measuring day.

4.5 The main influencing factors of R_s in forest–steppe habitat

Topography sometimes had a very clear modifying effect on abiotic and biotic factors of ecosystems in connection with cardinal directions, exposure, spatial heterogeneity and carbon spatial distribution (Alexander et al., 2016; Fóti et al., 2020; Lassueur et al., 2006). Based on our results, the co-varying effects of the habitat's terrain features and of the grove could be observed. We detected mixed effects of the topography and the shading of the grove on the measured variables; thus topographical effects can be masked by conditions related to the physiognomy of the forest–steppe. We highlight that the extent of the effect of topography may also differ for groves with slightly different locations and vegetation structure, as well as for slightly different sampling resolution and arrangement (Körmöczi et al., 2016). In addition, the distribution of the “warmer” and “colder” grassland positions in the PCA biplots in the context of the environmental and functional variables was a great indicator of the shadowing effect of the grove, which, in addition to the topography, also had a strong effect on these factors (Erdős et al., 2014, 2017; Süle et al., 2020). Topography and shading caused variability in the functional responses of the ecosystem. Overall, the variation in the functional responses was quite small, and the temporal trends were the same, although rainy weather and temporal variability had a strong influencing power.

5 Conclusions

This study provided information on the ecosystem functioning of a small group of trees with the surrounding grassland as a transition zone in a forest–steppe habitat. We described the functional responses and the influential effects of topography and shading.

The terrain features of the habitat and the physiognomy of the grove had a co-varying influence on the abiotic and biotic factors of this habitat. Topography had a clear effect on the ecosystem functioning and carbon spatial distribution, but we also highlighted the importance of locations, vegetation structure, weather and seasonal differences. The responses generally had similar characteristics and temporal trends, such as the vegetation structure showing a strong correlation with R_s and the components of R_s along the investigated seasons. The vegetation structure also influenced the microclimate, which in turn had an effect on R_s. Significant microclimatic differences were observed between the groves and the surrounding grassland. A difference in carbon allocation could also be observed along the transition zone, but these differences were primarily influenced by topography as well as by the effect of the microclimate resulting from differences in foliage.

Our observations are valuable for assessing the dynamics and spatiotemporal patterns of functional and driving variables in this type of ecosystem, which is a natural transition zone in the temperate vegetation. Due to climate change, this habitat will be threatened by total aridification and desertification, which will cause reduction in the amount and extent of forest patches. This process should result, for example, in a reduction in carbon stock, which in turn should strongly affect the carbon cycle.