Temperature-driven vapor pressure deficit structures forest bryophyte communities across the landscape

Atmospheric vapor pressure deficit (VPD) controls local plant physiology and global vegetation productivity. However, at ecologically crucial intermediate spatial scales, the role of VPD variability in forest bryophyte community assembly and the processes controlling this variability are little known.

To explore VPD effects on bryophyte community composition and richness and to disentangle processes controlling landscape-scale VPD variability, we recorded bryophyte communities and simultaneously measured forest microclimate air temperature and relative humidity across a topographically diverse landscape representing a bryophyte diversity hotspot in temperate Europe. Based on VPD importance for plant physiology, we hypothesize that VPD can be important also for bryophyte community assembly and that VPD variability will be jointly driven by saturated and actual vapor pressure across the topographically diverse landscape with contrasting forest types and steep microclimatic gradients.

Contrary to our expectation, VPD variability was dictated by temperature-driven differences in saturated vapor pressure, while actual vapor pressure was surprisingly constant across the landscape. Gradients in species composition, species richness and community structure of bryophyte assemblages followed closely the VPD variability. The average daily mean VPD was a much better predictor of species composition than average daily maximum VPD. The mean VPD also explained significantly more variation in species composition and richness than maximum temperature, indicating that time-averaged evaporative stress is more relevant for bryophyte communities than microclimatic extremes. While mesic forest bryophytes occurred along the whole VPD gradient, species occurring near their distributional limits and locally rare species preferred sites with low VPD. Consequently, low VPD sites represent species-rich microclimatic refugia within the landscape, where regionally abundant mesic forest bryophytes coexist with rare species occurring near their distributional range limits.

Our results showed that VPD variability at ecologically crucial landscape scales is controlled by temperature-driven saturated vapor pressure. Future climate warming will thus increase evaporative stress and reshuffle VPD-sensitive forest bryophyte communities even in topographically diverse landscapes, which are traditionally considered as microclimatic refugia buffered against climate change. Bryophyte species occurring near their distributional range limits in microclimatic refugia with low VPD will be especially vulnerable to the future changes in atmospheric VPD.

1 Introduction

Vapor pressure deficit (VPD) expresses atmospheric water demand as the difference between the amount of water vapor the air can hold at a given temperature and the actual amount of water vapor present in the air. Unlike relative air humidity, VPD accurately expresses plant evaporative stress (Campbell and Norman, 1998). Since air capacity to hold water vapor increases exponentially with temperature, the same relative humidity at different temperatures indicates very different atmospheric moisture conditions (Anderson, 1936). An atmosphere with the same relative air humidity may be very “dry” (when the temperature is high) or it may be very “wet” (when the temperature is low). Relative air humidity therefore does not indicate the atmospheric moisture condition in physiologically meaningful way, despite its popularity in ecological studies (Campbell and Norman, 1998). In contrast, VPD directly expresses the atmospheric moisture conditions in terms of plant evaporative stress (Anderson, 1936).

Atmospheric VPD is a key driver of plant functioning in terrestrial ecosystems (Ruehr et al., 2014; Grossiord et al., 2020), because higher VPD leads to reduced photosynthesis in the short term and drought-induced mortality in the long term (McDowell et al., 2008; Fu et al., 2022). Ongoing climate changes exacerbate VPD-driven evaporative stress because higher temperatures lead to an exponential increase in VPD (Lawrence, 2005; Grossiord et al., 2020). Increasing VPD already limits global vegetation productivity (Yuan et al., 2019; López et al., 2021; Lu et al., 2022) and triggers large-scale forest diebacks (Breshears et al., 2013; Eamus et al., 2013; Williams et al., 2013). Yet, in contrast to the widely recognized role of VPD in local plant physiology and global vegetation functioning, VPD effects on plant community assembly are largely unknown (Novick et al., 2024).

The knowledge about VPD effects on plant communities and the processes that control VPD variability over the landscape are crucial for more realistic predictions of climate change impacts on vegetation and the identification of microclimatic refugia (Ashcroft and Gollan, 2013; Davis et al., 2019; Finocchiaro et al., 2024; Ogée et al., 2024). Because VPD is a difference between saturated vapor pressure (P_sat) and actual vapor pressure (P_air), VPD variability reflects the interplay between spatial patterns in saturated and actual vapor pressures. While saturated vapor pressure is solely an exponential function of air temperature, actual vapor pressure is influenced by many processes operating at different spatial scales – ranging from regional atmospheric circulation and precipitation to local evaporation and plant transpiration (Campbell and Norman, 1998). Yet, despite increasingly recognized VPD importance, it is still unknown how these contrasting processes integrate into the VPD variability over the landscape.

A deeper understanding of the mechanisms behind landscape-scale VPD variability is particularly important for climate change biology. Scientists predict a temperature increase of up to 4.4 °C by 2100 (IPCC, 2023), which would lead to a more than 40 % increase in VPD for the same atmospheric water vapor content (Will et al., 2013). These changes can also modify VPD variability over the landscape, potentially shift the distribution of individual species and therefore alter the composition of plant communities. However, VPD effects on plant distribution and community assembly over the landscape are not sufficiently known.

Among plants, bryophytes are exceptionally sensitive to evaporative stress because they lack roots, lignified water-conducting system, water storage tissues, and active stomata and have a large surface area in proportion to biomass (Rice et al., 2001; Goffinet and Shaw, 2009). When exposed to the air with non-zero VPD, bryophytes therefore inevitably lose water (Hinshiri and Proctor, 1971; Busby and Whitfield, 1978). Because bryophytes transport water only passively, mainly through external capillary spaces between tiny parts of their body (Schofield, 1981), their internal water content is a function of the water availability in the surrounding environment (Vanderpoorten and Goffinet, 2009). Once the external water evaporates, bryophyte cells rapidly lose turgor, metabolic activity slows down, and carbon fixation decreases.

To cope with this evaporative stress, bryophytes developed an evolutionary and ecologically unique desiccation strategy, allowing them to survive drought episodes in a desiccated state (Proctor, 2000, 2001). Despite this ability to survive microclimatic extremes, bryophyte assemblages are potentially sensitive to evaporative stress, because desiccation tolerance widely differs among bryophyte species (Hinshiri and Proctor, 1971; Wagner and Titus, 1984; Oliver et al., 2000; Proctor et al., 2007a, b). Therefore, it can be assumed that the atmospheric VPD – an ecologically meaningful variable expressing evaporative stress – will strongly affect composition, richness and structure of bryophyte assemblages. Yet surprisingly little is known about the VPD effect on bryophyte assemblages in temperate forests (Fenton and Frego, 2005).

To provide this missing knowledge, we combined detailed in-situ forest microclimate measurements with simultaneous bryophyte inventories conducted across topographically diverse landscape representing bryophyte diversity hotspot in central Europe. Using these data, we explored how landscape-scale VPD variability affects bryophyte community composition and species richness in temperate forests, quantified VPD variability over the topographically diverse landscape, and identified which processes drive this variability.

2 Material and methods

2.1 Study area

We recorded bryophytes and measured microclimate in the Bohemian Switzerland National Park in the Czech Republic (Fig. 1). The rugged terrain of this sandstone landscape creates a fine-scale mosaic of contrasting habitats with steep microclimatic gradients over short distances (Wild et al., 2013). The elevation within the national park ranges from 125 to 619 m, and the mean elevation is 340 m. According to the data from the Tokáň weather station (Fig. 1), the mean annual air temperature during the 2011–2019 period was 8.3 °C, and the mean annual precipitation was 765 mm.

Figure 1 We measured microclimate and simultaneously recorded bryophyte species composition at 38 permanent research plots within the Bohemian Switzerland National Park in Central Europe. This forested area has rugged terrain creating steep environmental gradients over short distances.

Most of the Bohemian Switzerland is covered with coniferous forests. Norway spruce (Picea abies) planted mostly during the 19th and 20th century dominates in the valleys and on the plateaus, while patches of semi-natural forests are dominated either by Scots pine (Pinus sylvestris) on the upper slopes and rocky ridges or by European beech (Fagus sylvatica) on more mesic sites.

The nutrient-poor and strongly acidic soils result in a relatively low diversity of vascular plants, which contrasts with the exceptionally rich bryophyte flora (Härtel et al., 2007). The Bohemian Switzerland currently hosts more than 300 bryophyte species and therefore represents a hotspot of bryophyte diversity in Central Europe (Marková, 2008).

The bryophyte flora of the Bohemian Switzerland is dominated by forest species like Tetraphis pellucida, Bazzania trilobata, and Dicranum scoparium. These dominant floristic elements are enriched by disjunct occurrences of (sub)alpine or (sub)montane (e.g., Hygrobiella laxifolia, Geocalyx graveolens, Anastrophyllum michauxii), boreal (e.g., Dicranum majus, Rhytidiadelphus subpinnatus), and (sub)oceanic (e.g., Tetrodontium brownianum, Plagiothecium undulatum) species (Härtel et al., 2007; Marková, 2008).

2.2 Field data collection

We recorded bryophyte species composition and measured microclimate on 38 permanent plots within the Bohemian Switzerland National Park (Fig. 1). These plots were selected through stratified-random sampling to capture the main microclimatic gradients within the core zone of the national park. Specifically, using GIS and LiDAR-based digital terrain model, we first divided the study area into geographical strata defined by the terrain (valley bottoms, lower slopes, upper slopes, and ridges) and further separated the slopes with predominantly northern and southern orientation. Within each stratum, we randomly selected an equal number of locations separated by at least 50 m. In the field, we navigated to the selected location with GPS device and placed the center of plot 1.5 m to the north from the nearest tree.

Within each permanent plot, we installed HOBO U23 ProV2 (Onset, USA) microclimatic datalogger at 1.5 m height on the north side of a tree nearest to the plot center. Each HOBO datalogger was protected by a white radiation shield with good ventilation and measured air temperature (resolution 0.02 °C, accuracy ± 0.21 °C) and relative humidity (resolution 0.05 %, accuracy ± 2.5 %) every 30 min from 1 June to 31 August 2022.

Simultaneously with microclimate measurements, we recorded the presence of all bryophyte species in each research plot following the nomenclature of the Czech national checklist (Kučera et al., 2012). We deliberately sampled bryophytes in a relatively small circular plot with 1 m radius (3.14 m²) without any exposed rocks or big stones to reduce the possible effects of within-plot environmental heterogeneity (Rambo and Muir, 1998; Vanderpoorten and Engels, 2002; Schmalholz and Hylander, 2011).

2.3 Microclimate data processing

First, we checked the microclimatic time series visually and then with standard automated procedures implemented in the myClim R package (Man et al., 2023). Air humidity measurement with microclimatic loggers is sensitive to water condensation, resulting in unrealistically high measurements for prolong periods of time (Ashcroft and Gollan, 2013; Feld et al., 2013). We therefore carefully checked microclimatic time series and found no signs of the condensation effect.

Using checked air temperature and relative humidity data, we first calculated the saturated vapor pressure (P_sat) following the updated Buck formula (Buck, 1981, 1996):

$$\begin{array}{rl}{P}_{\mathrm{sat}}& =(\mathrm{1.003}+\mathrm{4.18}\times {\mathrm{10}}^{-\mathrm{6}}\times \mathrm{101}\mathrm{kPa})\times \mathrm{0.61115}\\ & \times e^{\left(\right(\mathrm{23.036}-t/\mathrm{333.7})\times (t/(\mathrm{279.82}+t)\left)\right)},\end{array}$$

where t is air temperature [°C].

Then, we calculated the actual vapor pressure (P_air) using the Tetens's formula (Tetens, 1930):

$$P_{\mathrm{air}}=P_{\mathrm{sat}}\times \left({\displaystyle \frac{\mathrm{rh}}{\mathrm{100}}}\right),$$

where rh is relative humidity [%].

Finally, we calculated atmospheric VPD as the difference between P_sat and P_air (Jones, 2014). Using the resulting microclimatic time series, we calculated three variables representing evaporative stress (Table 1). First, we calculated the average daily maximum temperature (T_max). While T_max is ecologically less meaningful proxy for evaporative stress than atmospheric VPD (Campbell and Norman, 1998; Eamus et al., 2013), several previous studies identified T_max as highly relevant microclimatic variable linked to evaporative stress and affecting species composition and richness of forest vascular plants and bryophytes within the central Europe (Macek et al., 2019; Man et al., 2022). Then, we calculated two variables capturing different aspects of VPD driven evaporative stress. First, we calculated the average daily maximum VPD (VPD_max), which represents site-specific microclimatic extremes (Ashcroft and Gollan, 2013). Second, we calculated the average daily mean VPD, which represents time-aggregated evaporative demand experienced by bryophytes on each site.

To disentangle the drivers of spatio-temporal VPD variability over the landscape, we calculated also plot-specific daily average values of P_sat and P_air (Table 1).

Table 1 Overview of microclimatic variables representing evaporative stress (T_max, VPD_max, VPD_mean) and its components (P_sat, P_air). For each variable, we provide the overall mean and range of plot-specific averaged daily values measured continually during summer 2022 on 38 forest research plots in the Bohemian Switzerland National Park, Czech Republic.

2.4 Data analysis

2.4.1 Bryophyte community composition, richness and structure

In our analysis, we focused on the relationship between microclimatic variables representing evaporative stress and bryophyte community composition, structure, and richness. First, we identified the main gradients in community composition and explored their relationship with variables representing evaporative stress. Then, to explore which variable representing evaporative stress is more closely associated with bryophyte community composition and richness, we calculated the variability in species composition and richness explained by the mean and maximum atmospheric VPD and maximum air temperature. Further, to disentangle the effects of atmospheric VPD from the effects of the maximum temperature, we partition the explained variability into independent and shared fractions. Finally, we tested the link between VPD and bryophyte community structure through nestedness analysis.

To explore the main gradients in the bryophyte community composition, we used non-metric multidimensional scaling (NMDS) to extract the main patterns in bryophyte community composition expressed with the Sørensen dissimilarity index. We calculated two-dimensional NMDS with the weak treatment of ties, a maximum of 500 random starts, and 999 iterations in each NMDS run using metaMDS function from the vegan R package version 2.6-4 (Oksanen et al., 2022). To maximize variance along the first ordination axis, we centered and rotated the resulting two-dimensional configuration with principal component analysis.

To explore how main compositional gradients correlate with microclimate variables representing evaporative stress, we passively projected vectors of maximum and mean VPD, and maximum temperature into the NMDS ordination space and tested the significance of the fit with 999 random permutations using the envfit function from vegan R package (Oksanen et al., 2022). Finally, we projected bryophyte species richness gradients into the NMDS ordination space using a generalized additive model (GAM) fitted through ordisurf function from vegan R package (Oksanen et al., 2022).To quantify the relationship between the microclimatic variables representing evaporative stress and species richness expressed as a number of bryophyte species recorded in the plot, we used GAM fitted with the R package mgcv 1.9.1 (Wood, 2011). We used GAM with Poisson distribution, log link function, and smooth terms fitted by thin plate regression splines without null space penalization and smoothing parameter estimation using restricted maximum likelihood. To assess the statistical significance, we used a χ² test comparing the fitted model to the only intercept null model.

To calculate the proportion of variability in bryophyte community composition explained by microclimatic variables representing evaporative stress, we used distance-based redundancy analysis (db-RDA) calculated on the same Sørensen dissimilarity matrix as used for NMDS (McArdle and Anderson, 2001). We calculated the db-RDA with dbrda function from vegan R package (Oksanen et al., 2022) and assess the statistical significance using 999 random permutations of the raw data (Legendre et al., 2011).

As all three microclimatic variables representing evaporative stress were correlated (Appendix A), we explored their shared and independent effects on bryophyte community composition and richness through variation partitioning (Legendre, 2008). Because VPD_max and T_max were almost identical (Pearson r=0.98), we disentangled shared and independent effects of substantially less correlated VPD_mean and T_max (Pearson r=0.78). To quantify their independent and shared effects, we partitioned the variation in bryophyte community composition explained by atmospheric VPD_mean and T_max using adjusted R² (Peres-Neto et al., 2006) calculated with the varpart function from the vegan R package (Oksanen et al., 2022).

To quantify the shared and independent effects of atmospheric VPD_mean and T_max on species richness, we partitioned the deviance explained in GAM models. First, we related species richness to atmospheric VPD_mean and T_max in the full GAM, when both variables were used simultaneously as predictors. Then, we fitted two partial GAMs (first with VPD_mean, second with T_max as explanatory variables). To prevent different smoothing parameters in the partial models, we extracted smoothing parameters from the full GAM and used them in both partial GAMs (Hjort et al., 2012). To assess the statistical significance, we compared each model against the null model with only intercept using a χ² test. To assesses the significance of the independent effects of atmospheric VPD and T_max, we compared partial GAMs with the full GAM using χ² test.

Finally, we used nestedness analyses (Ulrich et al., 2009) to test the VPD effects on bryophyte community structure. To directly test the two hypotheses about the bryophyte community structure along the VPD gradient, we first order the community matrix along the gradient of increasing plot-specific VPD_mean. To test the first hypothesis that the bryophyte communities from sites with high VPD are nested subsets of bryophyte communities from sites with low VPD, we used NODF_sites metric (Almeida-Neto et al., 2008). To test the second hypothesis that more frequent bryophyte species occur along the whole VPD gradient, but less frequent species are concentrated on sites with low VPD, we used NODF_species metric (Almeida-Neto et al., 2008). To calculate both NODF metrics, we used nestednodf function from the vegan R package (Oksanen et al., 2022).

We used a null model approach to assess the statistical significance of nestedness patterns (Ulrich et al., 2009). Specifically, we compared the observed NODF values to the distribution of 999 NODF values calculated through the conservative R1 null model, which maintains species richness of the site and uses species frequencies as probabilities of selecting species (Wright et al., 1997). To quantify the difference between the observed NODF values and the NODF values generated by the R1 null model, we calculated the standardized effect size (SES) expressing the number of standard deviations that the observed NODF value differs from the mean NODF value of the simulated assemblages (Ulrich et al., 2009). To construct the null models and to calculate SES, we used the oecosimu function from the vegan R package (Oksanen et al., 2022).

We used R version 4.4.0 (R Core Team, 2024) for complete data analysis and figure preparation. For the color-blind safe visualisations we used the R package scico 1.5.0 (Pedersen and Crameri, 2023).

2.4.2 VPD variability across the landscape

Using the time-series of both VPD components (P_sat and P_air), we explored their spatio-temporal variability and quantify their influence on the VPD variability over the landscape. First, we explored how variable was VPD and both its components over the landscape in a daily timesteps. Then, we averaged this daily variability into the overall measure of spatial variability in VPD, P_sat, and P_air during the whole study period. Finally, we used variation partitioning to quantify how much was VPD variability controlled by P_sat and P_air.

To quantify spatial variability in daily VPD and both its components (P_sat and P_air) over the landscape, we calculated the standard deviation (SD) of the plot-specific daily mean VPD, P_sat and P_air values among all study plots. In this first step, we calculated SD of these microclimatic variables for every day within the study period separately. Then we averaged these daily inter-plot SD values separately for VPD, P_sat and P_air into an overall measure of spatial variability for each microclimatic variable during the whole study period.

Finally, to disentangle the contribution of P_sat and P_air to the VPD variability over the landscape, we performed variation partitioning based on a multiple linear regression model and adjusted R² (Legendre, 2008) with the plot-specific mean VPD as the response variable and the mean P_sat and P_air as the predictors.

3 Results

3.1 Bryophyte community composition, richness and structure

In total, we recorded 39 bryophyte species: 14 liverworts and 25 mosses (Appendix C, Table C1). The species richness was highly variable among the plots – while the average number of species per plot was 8, the minimum was 1 and the maximum 21. The most frequent species were Dicranum scoparium (n=32), Leucobryum juniperoideum (n=26) and Hypnum cupressiforme (n=24).

Main patterns in community composition and species richness reflected the gradient of evaporative stress (Fig. 2). Gradient in T_max was highly correlated to the gradients in VPD (Fig. 2), but main patterns in community composition were less related to T_max than to VPD (vegan::envfit – T_max: R² = 0.32, p = 0.003; vegan::envfit − VPD_mean: R²=0.52, p = 0.001; VPD_max: R²=0.37, p=0.001).

Figure 2 Nonmetric multidimensional scaling (NMDS) of the bryophyte community composition showing main gradients in bryophyte assemblages sampled at 38 temperate forest plots. Points show the positions of the individual plots within the NMDS ordination space, and the vectors show the gradients in the maximum air temperature (T_max), maximum VPD (VPD_max) and mean VPD (VPD_mean). The smooth surface and associated contours fitted into the NMDS ordination space with a generalized additive model show the pattern of decreasing species richness with increasing evaporative stress.

The number of bryophyte species was higher in plots with low VPD and declined with an increasing VPD (Fig. 2). Both atmospheric VPD and maximum temperature were significantly associated with species richness, but maximum temperature explained substantially less deviance (Table 2). The mean VPD explained slightly more deviance than the maximum VPD (Table 2).

Table 2 Variation in community composition and species richness explained by three microclimatic variables representing evaporative stress. To quantify variation explained by each variable, we used distance-based redundancy analysis (db-RDA) for community composition and generalized additive models (GAM) for species richness.

The mean VPD explained substantially more variation in species composition than the maximum VPD and the maximum temperature (Table 2). When used independently, both VPD_mean and T_max were significant predictors of bryophyte community composition (Table 2). However, the effect of T_max almost completely overlaps with VPD_mean (Fig. 3). When we controlled for the effect of mean VPD, maximum temperature did not explain significant part of variation in community composition (vegan::dbrda – adj. R² = 0 %, p = 0.764) or in species richness (mgcv::gam – D²=3.72 %, p = 0.174). In contrast, the mean VPD explained a significant part of variation in species composition and richness even after the controlling for maximum temperature (vegan::dbrda − adj. R² = 6 %, p = 0.003; mgcv::gam – D²=22 %, p = 0.001). Therefore, the mean VPD explained substantially more variation in bryophyte community composition and richness than maximum temperature and maximum temperature did not have any significant effects independent from the mean atmospheric VPD (Fig. 3).

Figure 3 Variation partitioning showing independent and shared effect of mean VPD (VPD_mean) and maximum air temperature (T_max) on bryophytes species composition and richness in 38 forest plots. Values represent adjusted R² from distance-based redundancy analysis for species composition and explained deviance from a generalized additive model for species richness. While VPD_mean has significant effects even after the controlling for T_max both for species composition (p = 0.003) and richness (p <0.001), the unique effects T_max were non-significant both for species composition (p = 0.764) and richness (p = 0.174).

Bryophyte community structure was closely related to the gradient of mean atmospheric VPD (Fig. 4). Bryophyte communities from plots with higher VPD were generally impoverished and compositionally nested subset of the communities from sites with lower VPD (vegan::oecosimu – NODF_sites=39.17, SES = 4.26, p = 0.001). Moreover, while frequent species occurred along the whole VPD gradient, rare species occurred preferably on sites with low VPD (vegan::oecosimu – NODF_species=29.97, SES = 3.34, p = 0.003).

At the species level, small liverworts (e.g. Riccardia multifida, Lophozia ventricosa) and hygrophilous bryophytes (e.g. Polytrichum commune, Bazzania trilobata), as well as species with boreal (e.g. Dicranum majus) and (sub)oceanic (e.g. Mylia taylorii, Plagiothecium undulatum) distribution preferred plots with low atmospheric VPD (Fig. 4). In contrast, regionally frequent species like Hypnum cupressiforme, Polytrichum formosum or Dicranum scoparium occurred also in plots with higher atmospheric VPD (Fig. 4).

Figure 4 Occurrences of all recorded bryophyte species along the gradient of the mean VPD measured at 38 forest plots. Plots are sorted from the lowest to highest mean VPD and each filled square shows the presence of the focal species within the plot. While rare and species near their distributional range limits prefer sites with low VPD, mesic forest species occur along the whole VPD gradient.

3.2 VPD variability across the landscape

VPD in the forest understory was highly variable across the landscape (Fig. 5). While the variability in saturated vapor pressure was comparable to the variability in VPD, actual vapor pressure was much less variable among the sites (Fig. 5). In average, the landscape-scale spatial variability of P_sat (average daily SD = 0.20 kPa) was almost three times higher than the spatial variability of P_air (SD = 0.07 kPa).

Figure 5 Spatio-temporal variability of VPD and its components – saturated and actual atmospheric vapor pressures. Each data point shows the standard deviation of the plot-specific daily mean values simultaneously measured at 38 forest plots, and density plots summarize this spatio-temporal variability over the summer season. The individual data points were slightly jittered for better visibility.

The dominant driver of VPD variability across the landscape was temperature-driven saturated vapor pressure (Fig. 6). In a univariate linear regression model, P_sat explained 93 % of VPD variability, while P_air explained 30 %. However, P_sat and P_air were negatively correlated (Pearson $r=-\mathrm{0.31}$) and variation partitioning based on multiple regression model showed that the P_air uniquely explained only 7 % of variability in VPD (Fig. 6). Therefore, temperature-driven P_sat was the dominant driver of VPD variability, while spatial variation in P_air contributed surprisingly little to the overall VPD variability across the landscape.

Figure 6 Atmospheric vapor pressure deficit (VPD) was driven by temperature-dependent saturated vapor pressure, while actual vapor pressure was weakly related to local VPD. Each dot represents the mean VPD, and the mean saturated and actual vapor pressure measured during the summer at 38 forest plots established over topographically diverse landscape. Venn diagram shows variation (adjusted R²) in mean VPD explained solely by mean saturated (P_sat) and mean actual (P_air) vapor pressure and the variation explained jointly by both predictors.

4 Discussion

We found that community composition and richness of forest bryophytes was significantly affected by atmospheric VPD. Our findings have important implications both for theoretical and applied ecology. First, the variation in VPD over the landscape was largely controlled by air temperature. Therefore, air temperature and VPD are tightly coupled at biologically relevant scales, and their effects are hard to disentangle with observational data. Interestingly, this coupling was strongest between maximum VPD, and maximum temperature and maximum temperatures was previously identified as a key driver of bryophyte and vascular plant species distribution in temperate forests (Macek et al., 2019; Man et al., 2022). Unfortunately, these studies did not measure VPD. Considering our results, the importance of maximum temperature does not necessarily stem from its direct effects on plant ecophysiology, but more likely from strong temperature control of VPD variability over the landscape. Nevertheless, this new hypothesis needs further testing.

Interestingly, we also found that mean VPD was a much better predictor of bryophyte community composition and richness than maximum VPD or maximum temperature. At the same time, maximum temperature did not explain any additional variation in species composition and richness not explained by mean VPD. Our results thus provide strong evidence that the mean VPD is more relevant predictor of bryophyte community composition and richness than maximum temperature or maximum VPD. The unique effects of mean VPD, not reflected by the maximum temperature or maximum VPD, suggest that bryophyte communities are more sensitive to the long-term characteristics of site microclimatic conditions, rather than to short-term microclimatic extremes captured by maxima.

Second, our results showing that actual vapor pressure is relatively constant across the landscape imply that it is possible to estimate VPD from local microclimate air temperature measurements combined with non-local measurements of air relative humidity, for example from a nearby weather station. While the general applicability of this approach should be further tested across spatial scales (Dahlberg et al., 2020), in various environmental settings and different vegetation types, our findings suggest that local VPD can be reasonably estimated (Appendix B, Fig. B1). This finding thus opens exciting possibilities for further research as local temperature measurements are increasingly available all over the world (Lembrechts et al., 2020). However, it should be stressed that this approach generates VPD estimates which provide reasonable ranking of the sites along the VPD gradient, but generally overestimate the VPD (Appendix B, Fig. B1), likely because it does not account for locally higher actual vapor pressure, for example near springs, water bodies or on permanently waterlogged soils.

4.1 VPD variability across the landscape

Large spatial variability in atmospheric VPD structured forest bryophyte communities across the landscape. Interestingly, VPD variation was driven by temperature-controlled P_sat, while P_air was relatively constant across the landscape. This finding is important, as the actual vapor pressure should also be variable across the landscape (Ogeé et al., 2024; Johnston et al., 2025). However, our findings suggest that the local and spatially highly heterogeneous processes like evaporation from soil and water surfaces and plant transpiration contribute little to the landscape-scale variation in VPD, even in the topographically diverse landscape with steep microclimatic gradients.

While maximum VPD was solely driven by saturated vapor pressure and therefore maximum temperature, the mean VPD was more affected by actual vapor pressure. However, saturated and actual vapor pressures were negatively correlated and therefore the unique effect of actual vapor pressure on spatial pattern in atmospheric VPD was surprisingly small. The landscape-scale variation in atmospheric VPD was therefore controlled by microclimate temperature variation.

Microclimate temperature variation over the landscape, crucial for community ecology, is largely dictated by land-surface topography (Dobrowski, 2011). Land-surface topography controls also maximum air temperatures in the forest understory (Vanwalleghem and Meentemeyer, 2009; Macek et al., 2019) and therefore spatial variability in saturation vapor pressure. However, we were surprised that the highly localized processes like evapotranspiration did not contribute much to the spatial variability in absolute air humidity despite our study area with extremely rugged topography and contrasting forest vegetation types. Therefore, spatial variability in absolute air humidity seems to be determined mostly by processes operating at much larger scales like atmospheric circulation and precipitation patterns (Campbell and Norman, 1998). Nevertheless, local topographic depression with waterlogged soils and especially the proximity to flowing water or permanent water bodies can locally elevate actual vapor pressure and therefore decrease atmospheric VPD (Wei et al., 2018; Ogeé et al., 2024) However, our results suggest that the overall pattern in atmospheric VPD will generally follow changes in air temperature and therefore future climate warming will result in non-linear increase in evaporative stress across the landscapes. Given the growing recognition of VPD importance for many ecosystem processes, plant distribution, and community assembly (Grossiord et al., 2020; Kopecký et al., 2024; Novick et al., 2024), the approach we developed here to disentangle the contribution of saturated versus actual vapor pressure can provide new insights into the drivers of VPD variability across spatial and temporal scales. So far, the knowledge of the relative importance of saturated versus actual vapor pressure is limited, therefore it is difficult to compare our results with other studies. Nevertheless, a comparison of the drivers of VPD variability across agricultural fields in Germany supports our conclusion that temperature-driven variability in saturated vapor pressure is a dominant control of VPD variability at finer scales (Wörlen et al., 1999).

4.2 VPD effects on bryophytes

In contrast to vascular plants, bryophytes tolerate desiccation and become metabolically inactive in the absence of water (Proctor, 2000). When conditions improve, bryophytes quickly reactivate physiological processes such as respiration, photosynthesis, cell cycle, or normal cytoskeleton function (Proctor et al., 2007a, b). However, this reactivation requires a lot of energy, for example to produce specific repair proteins (Oliver and Bewley, 1984; Zeng et al., 2002) or to maintain the integrity and normal function of cell organelles and membranes (Platt et al., 1994). Prolonged periods without evaporative stress are therefore key for bryophyte growth and long-term survival (Proctor et al., 2007b; Merinero et al., 2020).

Bryophyte cells at full turgor have osmotic potential rarely more negative than −2 MPa (Proctor, 2000). An osmotic potential of −1.36 MPa is in equilibrium with air at 20 °C and 99 % relative humidity (i.e. VPD <0.03 kPa). If the temperature remains at 20 °C, but the relative humidity drops to 90 %, the water potential outside the bryophyte body decreases to −14 MPa (Proctor, 2000) and bryophytes start to lose water. To maintain full turgor and normal cell function, bryophytes thus need free liquid water close to the cells. However, this external water completely evaporates within 45–50 min if atmospheric VPD reaches 1.22 kPa (León-Vargas et al., 2006). Once the external water evaporates, bryophyte cells rapidly lose turgor, metabolic activity slows down, and carbon fixation decreases.

In our study region, microclimatic conditions with VPD lower than 0.03 kPa and therefore without evaporative stress for bryophytes (Proctor, 2000) occurred on average only 9 % of the measurement time. However, there was large variability among the sites, resulting in fine-scaled landscape mosaic of sites with widely different evaporative stress. We found that this fine-scale VPD variation structured bryophyte communities. Regionally rare species with disjunct distribution in central Europe generally preferred sites with low VPD. These species – otherwise typical for (sub)montane, boreal, or (sub)oceanic regions – are approaching their distributional limits within our study area (Hill and Preston, 1998). For these species, sites with low VPD serve as microclimatic refugia within an otherwise unsuitable landscape matrix. In contrast, regionally widespread bryophytes occurred along the whole VPD gradient. Fine-scale variation in VPD thus functions as an environmental filter for bryophyte community assembly over the landscape. Sites with low atmospheric VPD, hosting simultaneously rare as well as widespread bryophytes, thus represent hotspots of bryophyte diversity in the landscape.

Our findings of dominant temperature control on VPD variability across the landscape suggest that even the sites which can be considered as buffered against climate warming because of locally higher actual vapor pressure will be negatively affected by warming. With climate warming, areas with low VPD will likely shrink, and their bryophyte diversity will become more vulnerable (Pardow and Lakatos, 2013). Moreover, the increasingly frequent and severe canopy disturbances will likely further increase understory temperatures and consequently also VPD (Wolf et al., 2021; Máliš et al., 2023). Our results suggest that such changes will reshuffle bryophyte communities, supporting widespread mesic bryophytes at the expense of regionally rare species near their distributional limits. Such changes will likely decrease not only local and regional bryophyte species richness but also trigger biotic homogenization of bryophyte assemblages across larger spatial scales (Staude et al., 2020).

4.3 Disentangling atmospheric VPD and temperature

The close coupling between VPD and temperature clearly shows the need – and simultaneously the difficulty – of disentangling the influences of VPD and temperature on plant communities. While temperature affects basic life functions of bryophytes like photosynthesis, respiration (Dilks and Proctor, 1975), and growth (Furness and Grime, 1982), bryophytes thrive in a wide range of temperatures – from less than −30 °C (Dilks and Proctor, 1975) to over 40 °C in a dry state (Hearnshaw and Proctor, 1982). For most bryophytes, the optimal growth temperature ranges from 12 to 25 °C (Vanderpoorten and Goffinet, 2009). However, many bryophyte species grow even at temperatures around 5 °C (Dilks and Proctor, 1975), and some can even photosynthesize at temperatures below 0 °C (Lösch et al., 1983). Therefore, temperature is hardly a direct limiting factor of bryophyte distribution and community composition in temperate regions.

Our results fully support this conclusion, as we found that mean VPD was much better predictor of bryophyte community composition and richness than maximum temperature or maximum VPD. Bryophytes probably survive the most extreme conditions represented by maximum VPD in desiccated state. However, the time required to recover from desiccation increases and degree of recovery decreases with the length of desiccation (Proctor et al. 2007b). Bryophytes are therefore probably more sensitive to time-averaged characteristics of site microclimatic condition than to short-term extremes captured by maximum VPD. The open question is whether these findings apply also to vascular plants, which cannot survive microclimatic extremes in desiccated state and can be therefore more sensitive to the microclimatic extremes (Schönbeck et al., 2022).

Several studies of vascular plants have attempted to distinguish the independent effect of VPD from other microclimatic factors affecting plant functioning and distribution (Eamus et al., 2013; Denham et al., 2021; Flo et al., 2022; Fu et al., 2022; Kopecký et al., 2024), highlighting the critical importance of VPD (Novick et al., 2016; Schönbeck et al., 2022). Unfortunately, no physiological studies addressed the independent effects of VPD on bryophytes, despite clear indications that VPD plays a key role (Busby et al., 1978; Sonnleitner et al., 2009). So far, studies of bryophyte physiology concentrated on desiccation tolerance (Morales-Sánchez et al., 2022). While desiccation tolerance is an adaptation to cope with the external lack of water, the ultimate driver of desiccation is atmospheric VPD. A deeper focus on atmospheric VPD can therefore bring a new insight into bryophyte ecology and distribution.

5 Conclusions

Atmospheric VPD controls community composition, richness and structure of bryophyte assemblages in temperate forest understory. Even across the landscape with extremely rugged terrain, spatial variability in atmospheric VPD was controlled by temperature-dependent saturated vapor pressure. Maximum air temperature and VPD are thus tightly coupled at biologically relevant scales, and their effects are hard to disentangle. Nevertheless, we found that the time-averaged mean VPD was much better predictor of bryophyte assemblages than maximum temperature (or closely related maximum VPD) representing microclimatic extremes. This points toward the mean atmospheric VPD as the most important variable representing time-averaged evaporative stress and highlights so far overlooked importance of atmospheric VPD for bryophyte community ecology and distribution. With climate warming, the tight coupling between VPD and local air temperature will cause nonlinear increases in VPD-driven evaporative stress, which will subsequently reshuffle bryophyte community composition and decrease species richness. Especially vulnerable will be bryophyte species occurring near their distributional range limits in microclimatic refugia with low VPD.

Appendix A: Correlation of variables representing evaporative stress and its components

Figure A1 Pearson correlation matrix of microclimatic variables representing evaporative stress (maximum temperature, maximum and mean VPD) and its components (mean P_sat and mean P_air).

Appendix B: VPD estimate from in-situ air temperature and regional air humidity

Based on our results, we speculated that local atmospheric VPD can be reasonably estimated using the in-situ air temperature measurements paired with relative air humidity measurements representative for the whole region (and therefore the same for all plots situated within that region).

To explore this idea, we estimated the mean VPD using in-situ measured air temperature (HOBO U23 ProV2 dataloggers in 1.5 m height) and relative air humidity measured in the Tokáň weather station located in the study area (Fig. 1).

While the measured and estimated VPD were closely correlated (Pearson r = 0.97), estimated VPD were consistently higher than in-situ measured VPD (Fig. B1).

Therefore, we conclude that the relative position of the site on the VPD gradient can be reasonably estimated from in-situ microclimate temperature measurements paired with regional relative air humidity measurements. However, it should be stressed that this approach generates VPD estimates which provide reasonable ranking of the sites along the VPD gradient, but generally overestimate the VPD (Fig. B1), likely because it does not account for locally higher actual vapor pressure, for example near springs, water bodies or on permanently waterlogged soils. Therefore, this approach cannot fully replace local air humidity measurements.

Figure B1 Relationship between in-situ measured mean vapor pressure deficit (VPD) and mean VPD estimated from in-situ measured air temperature and relative air humidity measured in regional weather station (June–August 2022). While the measured and estimated VPD are closely correlated (Pearson r=0.97), estimated VPD tends to be higher than in-situ measured VPD, likely because of locally higher air humidity in topographically sheltered sites near valley bottoms.

Appendix C: List of bryophyte species

Table C1 Complete list of bryophyte species recorded at 38 study plots.