Tropical montane cloud forests are among the most biodiverse ecosystems on Earth, and are vulnerable to climate change due to reliance on atmospheric moisture. Epiphytic bryophytes (i.e., mosses, liverworts, hornworts) dominate these ecosystems and drive important ecosystem processes, yet their underlying strategies of resource use and functional structure within the canopy are not well understood. Community-level functional trait analyses along environmental gradients are valuable for understanding patterns of plant resource use in ecosystems. Along environmental gradients, intraspecific trait variation may obscure or drive patterns of functional structure but has often been overlooked. We examined bryophyte community functional structure among three vertically stratified zones in a Caribbean slope tropical montane cloud forest near Monteverde, Costa Rica. We tested how morphological and water-related traits associated with bryophyte economic spectra differ among vertical zones within cloud forest trees and determined the relative importance of intraspecific variation in shaping this structure. Functional structure differed significantly among zones and is suggestive of an economic trade-off whereby structural investment toward water holding capacity for species in the canopy comes at the cost of photosynthetic capacity, and vice versa on the trunk and base. Patterns of functional structure were mostly due to species turnover rather than intraspecific trait variation, which was supported by clear shifts in community composition among zones, and by species with high fidelity to specific zones. We found 171 bryophyte species (50 mosses, 120 liverworts, 1 hornwort), including 9 new species for Costa Rica, and 4 new for Central America. Our results suggest that these extremely diverse epiphytic bryophyte communities exhibit acquisition-conservation trade-offs in resource use like those known in vascular plants.

Materials and methods:

Study site and tree selection

This study was conducted in May–July of 2023 at the private ranch “Finca Agroecoturística BARPE” in mature secondary Tropical Montane Cloud Forest on the Caribbean slope of the Cordillera de Tilarán, 13km north-west of the Monteverde Cloud Forest Preserve, in the Cordillera de Tilarán, Costa Rica (1400 m a.s.l, 10° 22' 5'' N, 84° 49' 30'' W) (Fig. 1A). Research was conducted with permissions from SINAC and ACAT (M-P-SINAC-PNI-ACAT-0022-2023), and the private land owner. Mean annual rainfall, temperature, and relative humidity at the study site were 3720mm, 18.3 ºC, and 96.9% respectively, as determined by two years of continuous monitoring (2022–2024) (Gotsch et al. 2024 unpublished data). The vascular vegetation is characterised by a tall evergreen angiosperm-dominated canopy (to 30 m tall) with abundant lianas and epiphytes in the Piperaceae, Ericaceae, Clusiaceae, Pteridaceae, Hymenophyllaceae, Orchidaceae, and Bromeliaceae.

At Monteverde, like other tropical cloud forests, there is little variation in epiphyte host preferences, and rather tree size, branching architecture, bark texture, and microsite conditions are more important determinants of epiphyte composition (Pócs, 1982; K. Wagner et al., 2015). Given that 75–80% of the total bryophyte diversity for a given tropical forest stand can be adequately recorded from just 3-5 trees (Gradstein et al., 2003), six study trees were selected from an area of approximately 2 hectares that represented the oldest age class available (DBH 66–119 cm), had similar branching structure, and were determined safe to climb by single rope arborist techniques. The trees were located along a north-west oriented ridge, between 1362 and 1421 m above sea level.

Field data collection

In this study, we sampled bryophytes from three vertical strata: 1. Trunk and Base: the lower half of the trunk including the basal trunk, 2. Inner canopy: the first one-third of the crown branches, and 3. Mid-canopy: the middle third of the crown and branches (Fig. 1B). A laser rangefinder was used to determine the height to the crown and the top of the canopy, which was divided into thirds for sampling. Within all zones, three epiphyte mats (n = 54: 6 trees, 3 zones, 3 replicates per zone; Table 1; Fig. S1) were harvested in quadrat plots sized 400 cm² (20 cm x 20 cm), with rectangular quadrats (5 cm x 80 cm) used for small diameter branches (Shen et al., 2022). Vascular plants occurring in the plots were removed in the field, and arboreal soils were carefully removed in the lab to retain only the bryophyte layer and <1cm of soil. In this study, microclimatic differences between zones were assumed (e.g., Shen et al., 2022) and not measured directly, however several epiphytic mat and host-tree characteristics were recorded at each plot location and included branch (or trunk) inclination, diameter, soil depth, bark texture, vascular plant cover, bryophyte cover, average epiphyte layer height, and maximum epiphyte height. Branch diameter measurements were missed for the lower two layers of the first sampled tree. Missing data were estimated using the regression equation of height above ground and the diameter from the nearest tree of similar size and architecture (Fig. S2). Total bryophyte cover was estimated in the lab using Image J software (Schindelin et al., 2015) from images taken of the intact bryophyte community. Bark texture varied only slightly among our study trees and was often obscured by variable depths of arboreal soil; therefore, it was not considered to be important for shaping bryophyte composition or contributing to host-tree specificity.

Community composition was determined in the lab. Percent cover for detectable bryophyte morphospecies was visually estimated under a stereomicroscope and a 2 cm x 2 cm mesh grid overlayed on the sample to aid in percent cover estimates. Morphospecies were separated into specimen packets and identified to the species level. Small shoots of additional species growing intermixed with, but not initially noticed before close microscopic inspection of the field-identified morphospecies were considered to have trace abundance in the plots. Taxonomy follows Gradstein (2021) for liverworts and hornworts and Allen (1994-2018) for mosses. Specimens were deposited at the Costa Rica National Museum (CR), University of Costa Rica Herbarium (USJ), University of Cincinnati (CINC), and the University of British Columbia (UBC).

Functional trait sampling and selection

We selected a combination of ‘soft’ easily measurable morphological and water-related traits that are important proxies for physiological properties aligned with economic spectra in bryophytes (Grau-Andrés et al., 2022; Wang et al., 2017; Wang & Bader, 2018). Morphological traits measured included shoot length (SL), shoot width (SW), shoot length-width ratio (SLW), shoot area (SA), specific shoot area (SSA), and shoot dry mass (SM). Water-related traits included shoot water holding capacity (WHC) and shoot half-desiccation rate (T50).

Water-related traits

Water holding capacity (WHC) was calculated as the percentage dry-weight difference between dry and saturated mass ((saturated mass (mg) – dry mass (mg) / dry mass (mg)) x 100), using a protocol adapted from (Elumeeva et al., 2011). First, air-dried shoots were weighed on a microbalance (Mettler and Toledo MX5) to obtain their dried mass. To account for daily variation in lab humidity at the time of dry-weight measurement (43-63% RH), three species (typically 15 individuals) were selected to be oven dried at 60 ºC for 24 hours, and the percent difference between their air-dry and oven-dried weight was averaged and applied as a conversion factor for all individuals measured on that date (Roos et al., 2019). For saturated mass, dry shoots were submerged in distilled water for 30 minutes and then placed on moist filter paper and sealed in petri dishes for 24 hours. Shoots were then lightly blotted to remove excess external water and weighed to obtain saturated mass. Following saturated mass, morphological shoot traits were measured, and then re-saturated in distilled water for 15 minutes and lightly blotted to remove excess before measurement of the half-desiccation rate. Measurement of half-desiccation rate was modified from the protocol used by (Elumeeva et al., 2011; Michel et al., 2012). Shoots were pooled on plastic dishes by each species-zone combination, and the half-desiccation rate (minutes to 50% water loss (mg)) was interpolated from the linear desiccation rate of the pooled samples. The pooled saturated shoots were air dried at room temperature and humidity in 10-minute intervals for the first 40 minutes, 20-minute intervals for the next 60 minutes, and then 30-minute intervals to a maximum of 210 minutes. Weighing of the pooled shoots ceased on the next time interval after the sample reached 50% saturated weight. To account for mass differences among species, the half-desiccation rate was divided by the pooled dry mass of the shoots for each sample to get the shoot equivalent T50 (minutes to 50% water loss/mg drymass).

Morphological traits

Morphological traits were measured on saturated shoots projected in 2 dimensions at 1200dpi on a flatbed scanner (Epson V600) in WinFOLIA software (Regent Instruments Canada Inc.). Shoot length (mm) corresponded to the furthest distance between the base of the shoot and apex (or branch apex for multi-branched shoots). Shoot width (mm) was measured by the furthest two points perpendicular to the axis of the shoot length. Specific shoot area (mm²/mg) was calculated as the average shoot area of the individuals, divided by the average dry mass of the individuals for each sample.

Traits were sampled for species that comprised up to 90% relative abundance (cover) in plots. Five shoots of each species (when possible) were carefully removed from the plot where they had the highest relative abundance in each zone (Carmona et al., 2012) and trimmed to remove dead basal tissue. Traits were measured on species by ‘environmental levels’ where species trait values are fixed within each level (i.e., the trait value for a species is the same for all plots in that zone) (de Bello et al., 2021). This resulted in trait data for 58 species found in one or more zones (101 unique species x zone combinations) (Table S1). Traits were measured in batches of 10–12 species with the order randomized. Implicit in this approach is the assumption that intraspecific trait variation within a zone is lower than intraspecific variation among the zones.

Data analysis

Community composition, epiphytic mat and host-tree characteristics

We used linear mixed models in the R package lme4 (Bates et al., 2025) to analyze differences in species richness, and epiphytic mat and host-tree characteristics among zones. Zone was used as our fixed effect predictor, with species richness, or the epiphytic mat and host-tree characteristics as the response variables. Tree identity was included in each model as a random factor (intercept). We tested the significance of the fixed effect using Type III ANOVA with Satterthwaite’s method (Satterthwaite, 1946). When a significant effect was detected, we performed post-hoc pairwise comparisons using the estimated marginal means with Tukey adjustment for multiple tests. Average epiphyte height, and maximum epiphyte height were log-transformed, vascular epiphyte cover, bryophyte cover, and inclination (sine component) were logit transformed, and soil depth was inverse hyperbolic-sine transformed to normalize the distribution of the residuals and meet the model assumptions.

To examine if community composition differed between zones, we used non-metric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarities, followed by Permutational Multivariate Analysis of Variance (PERMANOVA) (Anderson, 2001) and post-hoc pairwise permutational t-tests (with p-values Bonferroni adjusted for multiple comparisons) with 999 iterations to test for differences between zones using the R packages vegan (Oksanen et al., 2024) and RVAideMemoire (Hervé, 2023). Prior to ordination, extremely rare species (occurring in less than two plots) were removed to reduce noise introduced by rare species and small sample size in the analysis (Cao et al., 2001), and absolute species abundances were standardized by plot to relative abundances. A 2-dimensional solution was found (stress = 0.164). Species fidelity towards zones was assessed on the reduced data set using the point-biserial correlation coefficient (r_pb) which provides a correlation value (between -1 and 1) for habitats based on species abundances (Cáceres & Legendre, 2009). Significance of species associations were tested permutationally with 999 iterations. To account for potential non-independence within the trees, plot data were averaged by zone for each tree.

Community-level functional structure

Traits were log-transformed prior to analyses to normalize the distribution of the residuals and meet model assumptions. To examine how functional structure varied among different zones, the community-weighted means (CWMs) were calculated for each trait and every plot (n = 54). To quantify the influence of intraspecific trait variability (ITV) on the functional structure, we used the procedure described by Lepš et al. (2011) where two different CWMs are calculated, one where species trait averages are calculated at each plot (specific average), and the other using the averages of species trait values across all plots (fixed average). The difference between these averages should isolate the intraspecific effect. Three ANOVAs were run in parallel using the trait.flex.anova function in R on the separate averages (specific, fixed, and their difference) to test for significance differences among zones for each effect (Lepš et al., 2011). When the ANOVAs of specific averages were significant (p < 0.05), pairwise Tukey HSD tests were used to further examine the differences between zones. The procedure outlined by Lepš et al. (2011) does not support multiple levels of error, therefore we plotted the residuals of each ANOVA model against each study tree to assess if there was structure in the data due to trees (Figs. S3, S4). As well, all models were also fit with the response as averaged values for each tree to remove the potential non-independence, but this did not change our interpretation of the results, so we proceeded with the original models.

To determine the variation explained by the environment for each component (i.e., turnover, ITV, covariation) sum of squares of each ANOVA (specific and fixed averages and their difference), were decomposed into the variability explained by the environmental factor and the unexplained residual variance (Lepš et al. 2011). Two plots were found to be extreme outliers for all traits, classified here by having CWM values 2.5 times below or above the first and third quartiles and were excluded prior to analyses.