DOI: 10.5061/dryad.3bk3j9kxt
title: Data from: Community-level trait variation of epiphytic bryophytes supports trade-off aligned with leaf-economic spectrum in vertically stratified tropical montane cloud forest canopies

description: 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.
description: Materials and methods: <em>Study site
 and tree selection</em> 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.
 <em>Field data collection</em> In
 this study, we sampled bryophytes from three vertical strata: 1.
 <em>Trunk and Base: </em>the lower half of the trunk including
 the basal trunk, 2. <em>Inner canopy:</em> the first one-third
 of the crown branches, and 3. <em>Mid-canopy</em>: 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<sup>2</sup> (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 &lt;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).
 <em>Functional trait sampling and
 selection</em> 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). <em>Water-related
 traits</em> 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). <em>Morphological
 traits</em> <em> </em>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<sup>2</sup>/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.  
 <em>Data analysis </em>
 <em>Community composition, epiphytic mat and host-tree
 characteristics</em> 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<sub>pb</sub>)
 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.
 <em>Community-level </em><em>functional
 structure</em> 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 &lt; 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.
description: # Data from: Community-level trait variation of epiphytic bryophytes
 supports trade-off aligned with leaf-economic spectrum in vertically
 stratified tropical montane cloud forest canopies Dataset DOI:
 [10.5061/dryad.3bk3j9kxt](10.5061/dryad.3bk3j9kxt) ## Description of the
 data and file structure ### Files and variables #### File:
 Tucker_et_al_2025_Trait_Data.xlsx **Description:** This dataset is the raw
 trait data. Each trait is a new sheet in the excel file.  ##### Traits: *
 SM: Shoot mass (mg) * SW: Shoot width (mm) * SL: Shoot length (mm) * SLW:
 Shoot length to width ratio (mmL/mmW) * T50: Shoot half-dessication rate
 (mins to 50% saturated weight) * WHC: Shoot water holding capacity ( %
 drymass) * SSA: specific shoot area (mm2/mg) * SA: Shoot area (mm2) #####
 Variables: * Note: For WHC and SSA which are derrived from calculations,
 it was not possible to calculate values for each individual independantly,
 therefore the values represent the mean value across all shoots for that
 species-zone combination. * Note: For T50: 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. * Column Variables:  * Zone: * 1 =
 Trunk and Base * 2 = Inner-canopy * 3 = Mid- canopy * Individuals:
 Individuals (shoots) of species that traits were measured on.
 "NA" indicates for that particular species/zone/trait
 combination it was not possible to measure traits on a full set of five
 individuals (i.e., &lt; less than five individuals were suitable for
 sampling). #### File: Tucker_et_al_2025_ENV_Data.csv
 **Description:** Epiphyte mat and host-tree characteristics #####
 Variables: * Plot: 400 cm2 plots - unique identifiers * Tree: Tree
 identity * Zone: Sampling zone on tree, Zone 1 = trunk and base, Zone 2 =
 inner canopy, Zone 3 = mid-canopy * Branch/trunk diameter (cm) *
 Branch/Trunk Inclination (degrees) * Bark texture * vascular epiphyte
 cover (%) * Bryophyte cover (%) * Average epiphyte height (cm) * Max
 epiphyte height (cm) * Soil Depth (cm) * Date Note: NA = value for average
 epiphyte height was missed for plot 31 during field collection. #### File:
 Tucker_et_al._2025_community_comp.csv **Description:** Raw community
 composition data. Abundance values are visually assessed percent cover.
 ##### Variables * Columns:  * Plot = unique identifiers for the total
 number of samples (400 cm2 qudrats of bryophytes removed from the trees).
 Organization is in ascending order by zone and tree. For example, plots
 1-3 = Tree 1, Zone 1 (trunk and base), plots 4-6 = Tree 1 Zone 2 (inner
 canopy), plots 7-9 = Tree 1, Zone 3 (Mid-canopy). * All remaining columns
 are species * Values: Values are the raw visually estimated percent cover.
 Species with trace abundance in a plot are given a value of 0.5%. ####
 File: Tucker_et_al._2025_Dominant_Species_Relative_Abundance.csv
 **Description:** Abundance matrix for the subset of dominant species
 (those that contribute up to 90% relative abundance in each plot) used to
 calculate trait indices re-normalized to relative abundances.  #####
 Variables * Columns:  * Plot = unique identifiers for the total number of
 samples (400cm2 qudrats of bryophytes removed from the trees).
 Organization is in ascending order by tree and zone. For example, plots
 1-3 = Tree 1, Zone 1 (trunk and base), plots 4-6 = Tree 1 Zone 2 (inner
 canopy), plots 7-9 = Tree 1, Zone 3 (Mid-canopy).  * All remaining columns
 are species * Values: Values are the relative abundance proportions of
 dominant species in each plot ## Access information Other publicly
 accessible locations of the data: * iNaturalist was used as a repository
 to share interesting species observations. No data were derived from
 iNaturalist. Species observations can been found in the
 project [https://www.inaturalist.org/projects/bryophytes-of-monteverde](https://www.inaturalist.org/projects/bryophytes-of-monteverde)