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

Abstract

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_{\mathrm{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_{\mathrm{s}}$ are autotrophic respiration ($R_{\mathrm{aut}}$) by plant roots and root-associated microbes and heterotrophic respiration ($R_{\mathrm{het}}$) by microbes and animals in the soil (Balogh et al., 2016).

Soil factors generally influence the temporal and spatial variance of $R_{\mathrm{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_{\mathrm{aut}}$ as well as the activity of soil microbial biomass through $R_{\mathrm{het}}$ have great influence on $R_{\mathrm{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_{\mathrm{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_{\mathrm{s}}$ shows a positive temporal relationship with soil temperature ($T_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ emission (Hao et al., 2010; Herbst et al., 2009; Mitra et al., 2019; Raich and Tufekcioglu, 2000; Thomas et al., 2018).

SWC and $T_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ 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_{\mathrm{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 ${}^{\circ }\mathrm{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 omitted. See PDF]

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.

[Figure omitted. See PDF]

Table 1

Measurement campaigns during the study.

Measurement			Study site		Transects
Year	Campaigns	Number of	Location	Average	Quantity of	Length	Cardinal	Max. elevation
		positions		diameter	transect	of each	directions	difference of
				of grove		transect		positions
2018	Spring, Summer, Autumn	89	First grove	15 m	4	44 m	N–S, EN–SW, E–W, SW–NW	4.72 m
2019	Summer, Autumn	61	Second grove	26 m	2	60 m	N–S, E–W	2.63 m
2020	Spring, Summer, Autumn	61	First grove	15 m	2	60 m	N–S, E–W	5.51 m

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_{\mathrm{s}}$ and $T_{\mathrm{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_{\mathrm{het}}$ and $R_{\mathrm{aut}}$ components, we calculated the mean $R_{\mathrm{s}}$ from the whole dataset of each measuring campaign. From the LAI–$R_{\mathrm{s}}$ linear regression, we determined the $y$-axis intercept (where LAI is theoretically 0) for each measurement campaign, considering it as the average $R_{\mathrm{het}}$ rate. Therefore, we calculated $R_{\mathrm{aut}}$ rate as the difference between mean $R_{\mathrm{s}}$ rate and the $R_{\mathrm{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):

1$$\text{VPD}=(\mathrm{100}-\text{RH})\times \mathrm{6.112}\times {\mathrm{e}}^{\mathrm{17.67}\times t/(t+\mathrm{234.5})},$$ 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 $\mathrm{5}\times \mathrm{5}$ box blur kernels as the neighborhood as follows:

2$$\text{SD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n_R-\mathrm{1}}}\sum _{i=\mathrm{1}}(z_i-\stackrel{\mathrm{‾}}{z}{)}^{\mathrm{2}}},$$ 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${}^{\circ }$ from the north (both 0 and 360${}^{\circ }$ meaning N-facing slope, 90${}^{\circ }$ 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 $-\mathrm{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.

[Figure omitted. See PDF]

Table 2

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

Campaigns	Average of air		Sum of precipitation		Average of global
	temperature (${}^{\circ }$C)		(mm)		radiation ($\mathrm{W}{\mathrm{m}}^{\mathrm{2}}$)
	3 d before	Measuring day	3 d before	Measuring day	3 d before	Measuring day
Spring 2018	20.66	19.63	2.7	1.2	264.27	192.75
Summer 2018	21.13	21.94	0.1	0	246.98	301.32
Autumn 2018	15.68	13.84	0.3	0	143.73	151.08
Summer 2019	26.41	21.65	0	0	314.33	293.15
Autumn 2019	17.71	16.2	1.3	0	138.46	181.96
Spring 2020	16.25	18.47	0.1	2.5	267.97	261.99
Summer 2020	22.38	22.27	5.8	11.9	263.09	264.15
Autumn 2020	13.79	13.78	2.4	0.2	67.35	35.28

3.2 Relationships between the measured variables

$R_{\mathrm{s}}$ showed significant positive correlation with SWC ($P<\mathrm{0.001}$, $R^{\mathrm{2}}=\mathrm{0.13}$) (Fig. 4a) and LAI ($P<\mathrm{0.001}$, $R^{\mathrm{2}}=\mathrm{0.28}$) (Fig. 4c) and showed negative correlation ($P=\mathrm{0.059}$, $R^{\mathrm{2}}=\mathrm{0.004}$) with $T_{\mathrm{s}}$ (Fig. 4b). The highest $R_{\mathrm{s}}$ values were observed in summer 2019 (23.99 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). In general, $R_{\mathrm{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 ${}^{\circ }\mathrm{C}$) and low LAI (0.15 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{\mathrm{2}}$). No significant relationship was found between $T_{\mathrm{s}}$ and $R_{\mathrm{s}}$ (Fig. 4b) due to the high temperatures measured in summer 2018 coupled with low water availability. LAI showed significant negative correlation ($P<\mathrm{0.001}$, $R^{\mathrm{2}}=\mathrm{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).

[Figure omitted. See PDF]

3.3

$R_{\mathrm{s}}$ component rates

A characteristic trend can be observed in the $R_{\mathrm{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_{\mathrm{aut}}$ seemed to closely follow LAI.

Table 3

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

Campaigns	Mean soil	Heterotrophic	Autotrophic	Mean leaf	$R_{\mathrm{aut}}$/mean
	respiration	$R_{\mathrm{s}}$ ($R_{\mathrm{het}}$)	$R_{\mathrm{s}}$ ($R_{\mathrm{aut}}$)	area index	$R_{\mathrm{s}}$
	($R_{\mathrm{s}}$)			(LAI)
Spring 2018	3.21	1.62	1.59	1.09	0.50
Summer 2018	2.27	0.64	1.62	1.36	0.72
Autumn 2018	1.52	0.6	0.92	0.98	0.61
Summer 2019	7.06	2.3	4.76	2.05	0.67
Autumn 2019	3.37	1.1	2.27	1.15	0.67
Spring 2020	1.7	0.36	1.33	1.13	0.79
Summer 2020	3.3	1.91	1.39	0.52	0.42
Autumn 2020	1.54	0.46	1.08	0.66	0.70

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_{\mathrm{s}}$ showed a positive relationship with SWC and LAI, while it showed a negative relationship with $T_{\mathrm{s}}$ and VPD, as shown by PCA loadings along the first axis. The second axis separated SWC, $R_{\mathrm{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_{\mathrm{s}}$ were separate from LAI. Mainly, the loadings of $R_{\mathrm{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_{\mathrm{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_{\mathrm{s}}$), soil water content (SWC), leaf area index (LAI), soil temperature ($T_{\mathrm{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).

[Figure omitted. See PDF]

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_{\mathrm{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_{\mathrm{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_{\mathrm{s}}$–VPD values at south-facing directions and also larger $R_{\mathrm{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<\mathrm{0.001}$, ${}^{\ast \ast }$: $P<\mathrm{0.01}$, ${}^{\ast }$: $P<\mathrm{0.05}$).

Campaigns	Axes	SOC	East	North	Sl	SD
Spring 2018	PCA I	${}^{\ast \ast \ast }$ 0.44	${}^{\ast \ast \ast }$ $-\mathrm{0.55}$	$-\mathrm{0.05}$	$-\mathrm{0.02}$	0.02
	PCA II	$-\mathrm{0.02}$	$-\mathrm{0.19}$	0.02	${}^{\ast \ast }$ $-\mathrm{0.29}$	${}^{\ast \ast }$ $-\mathrm{0.27}$
Summer 2018	PCA I	${}^{\ast \ast \ast }$ 0.54	${}^{\ast \ast \ast }$ $-\mathrm{0.51}$	0.08	$-\mathrm{0.01}$	0.02
	PCA II	${}^{\ast }$ $-\mathrm{0.21}$	0.11	$-\mathrm{0.13}$	$-\mathrm{0.2}$	$-\mathrm{0.19}$
Autumn 2018	PCA I	${}^{\ast \ast \ast }$ 0.4	${}^{\ast \ast \ast }$ $-\mathrm{0.55}$	$-\mathrm{0.08}$	0.084	0.12
	PCA II	0.1	0.1	0.04	$-\mathrm{0.15}$	$-\mathrm{0.17}$
Summer 2019	PCA I	${}^{\ast \ast \ast }$ 0.67	$-\mathrm{0.1}$	${}^{\ast \ast }$ $-\mathrm{0.36}$	${}^{\ast \ast \ast }$ $-\mathrm{0.42}$	${}^{\ast }$ $-\mathrm{0.3}$
	PCA II	$-\mathrm{0.18}$	0.03	${}^{\ast \ast }$ 0.4	0.2	0.21
Autumn 2019	PCA I	${}^{\ast \ast }$ 0.36	${}^{\ast }$ $-\mathrm{0.3}$	${}^{\ast }$ $-\mathrm{0.31}$	$-\mathrm{0.23}$	$-\mathrm{0.14}$
	PCA II	${}^{\ast \ast }$ $-\mathrm{0.4}$	$-\mathrm{0.15}$	0.1	${}^{\ast \ast }$ 0.34	${}^{\ast }$ 0.28
Spring 2020	PCA I	${}^{\ast \ast \ast }$ 0.53	${}^{\ast \ast }$ $-\mathrm{0.39}$	$-\mathrm{0.12}$	$-\mathrm{0.16}$	$-\mathrm{0.11}$
	PCA II	$-\mathrm{0.19}$	0.09	0.01	0.01	0.02
Summer 2020	PCA I	${}^{\ast \ast \ast }$ 0.65	$-\mathrm{0.19}$	$-\mathrm{0.08}$	$-\mathrm{0.06}$	0.004
	PCA II	0.04	0.2	0.14	$-\mathrm{0.2}$	$-\mathrm{0.19}$
Autumn 2020	PCA I	${}^{\ast \ast \ast }$ $-\mathrm{0.44}$	${}^{\ast }$ 0.29	0.18	$-\mathrm{0.05}$	$-\mathrm{0.11}$
	PCA II	${}^{\ast \ast }$ 0.4	0.02	0.05	$-\mathrm{0.04}$	$-\mathrm{0.07}$

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_{\mathrm{s}}$

The responses of $R_{\mathrm{s}}$ to soil moisture and $T_{\mathrm{s}}$ were investigated in several studies, which already reported that the correlation between $R_{\mathrm{s}}$ and $T_{\mathrm{s}}$ as well as SWC and $T_{\mathrm{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_{\mathrm{s}}$ showed positive correlation with SWC and LAI. Although it was often found that there is a positive correlation between $R_{\mathrm{s}}$ and $T_{\mathrm{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_{\mathrm{s}}$ showed negative correlation with $T_{\mathrm{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_{\mathrm{s}}$ and LAI in the PCA analysis of spatial data. The below-canopy area had higher $R_{\mathrm{s}}$ and LAI with lower air temperature and $T_{\mathrm{s}}$ due to the microclimate below the complex vegetation structure and higher biomass with trees and bushes. In contrast, the grassland had higher $T_{\mathrm{s}}$ and VPD with lower $R_{\mathrm{s}}$ due to the higher irradiation (Huang et al., 2020; Moyes and Bowling, 2016; Thomas et al., 2018).

4.2

The components of $R_{\mathrm{s}}$

We found seasonal differences in $R_{\mathrm{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_{\mathrm{s}}$ rates in autumn than in the summer. Besides, the correlation between $R_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ 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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{s}}$ and SWC. Opposite orientations of $R_{\mathrm{s}}$–SWC–LAI and $T_{\mathrm{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_{\mathrm{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_{\mathrm{s}}$ and the components of $R_{\mathrm{s}}$ along the investigated seasons. The vegetation structure also influenced the microclimate, which in turn had an effect on $R_{\mathrm{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.

Data availability

The datasets generated and/or analyzed during the current study are available in Figshare: 10.6084/m9.figshare.14269088.v2 (Süle, 2021a), 10.6084/m9.figshare.14269142.v1 (Süle, 2021b), 10.6084/m9.figshare.14269181.v1 (Süle, 2021c) and 10.6084/m9.figshare.14269214 (Süle, 2021d).

Author contributions

GS, LKö, SF, and JB conceived the study and designed the methodology. GS, LKö, SF, DP and JB performed fieldwork. GS, SF, and LKa analyzed the data. GS, LKö, SF, and JB wrote the paper. All authors have read and agreed to the published version of the paper.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We would like to thank the Institute of Ecology and Botany, Centre for Ecological Research, for use of the Fülöpháza research station during the measurement campaigns. We are grateful to György Kröel-Dulay and Gábor Ónodi for the meteorological data.

Financial support

This research was supported by the Ministry of Innovation and Technology within the framework of the Thematic Excellence Program 2020, Institutional Excellence Sub-Program (TKP2020-IKA-12) in the topic of water-related research of MATE.

Review statement

This paper was edited by Daniel Montesinos and reviewed by two anonymous referees.