This study was carried out in Marojejy National Park in the Sava Region of north-eastern Madagascar and centred on the Marojejy massif. The mountain spans a wide elevational gradient from 250 to 2132 m, making it an ideal site for examining trait filtering across elevations. The rugged topography of the massif creates diverse habitats that transition quickly with changes in altitude. The moist evergreen forests occurring at lower elevations gradually transition at about 1200 m into montane cloud forest, marked by its lower canopy and more dense understory and, finally, a forest line occurs around 1800 m and ericoid montane thickets occur to the peak at 2132 m. Rainfall is abundant, and almost no dry season has been recorded on the eastern slope of the massif. Across this slope, average annual temperatures range from 14.9 to 21.8°C, with relative humidity consistently high, averaging between 94.8% and 98.9% throughout the year (Marline et al. 2023).

Bryophytes were collected along an elevational transect at 200 m intervals from 250 to 2050 m and following the hierarchical sampling method described in Ah-Peng et al. (2007) and Marline (2018). At each elevation, two plots of 10 × 10 m were established. Within each plot, three quadrats of 2 × 2 m were selected randomly. In each quadrat, three trees were randomly chosen, and 5 × 10 cm epiphytic bryophyte samples (microplots) were collected at three heights: 0–50 cm, 50–100 cm, and 100–200 cm above ground level.

Temperature and relative humidity were monitored using MadgeTech RHTemp1000 data loggers, installed at five elevations (450, 850, 1250, 1650, and 2050 m) in December 2013. These loggers recorded temperature and humidity continuously at hourly intervals from December 2013 to December 2014. To estimate temperature and relative humidity at locations without data loggers, a calibration curve based on elevation was applied. For each elevation, vapor pressure deficit (VPD) was calculated as the difference between saturated vapor pressure (SVP) and actual vapor pressure (VP) following Monteith and Unsworth (2013). Minimum and maximum canopy height were also recorded at each site. A summary of the environmental variables, including temperature, relative humidity, vapor pressure deficit, and canopy height across elevations, is presented in Supplementary material 1.

For the species found on the transect, we compiled data on 12 morphological traits (Table 1) that are likely to be related to resource use, life history, species defense, resistance to desiccation, and photosynthetic activity, and could be obtained for at least 95% of the collected taxa. Trait values for each taxon were collected from the literature or by direct measurements on herbarium specimens.

Two multidimensional indices, functional evenness (FEve) and functional dispersion (FDis), were used to characterize functional diversity, following the framework proposed by (Botta-Dukát 2005; Mason et al. 2005; Villéger et al. 2008; Laliberté and Legendre 2010). FEve, a measure of the evenness of species abundance distributions in functional trait space (Mason et al. 2005; Villéger et al. 2008) reflects the degree to which a community can effectively utilize the entire range of resources available to it. FDis measures the weighted (by relative abundance) mean distance of individual species to their weighted community centroid and describes the functional similarity of species in a community in trait space. Its variation across species indicates the degree of functional redundancy, typically measured across multiple traits, among species within a given spatial scale (Laliberté and Legendre 2010).

We used linear regressions to evaluate the relationship between functional diversity indices and species richness across altitudes.

To investigate potential functional responses underlying the process of community assembly along the elevational gradient, we measured the community-weighted means (CWM) and community-weighted variances (CWV) for each trait along the elevational transect. The CWM for each trait is the mean across species present in the community weighted by their relative abundance (Violle et al. 2007; Bernard-Verdier et al. 2012), in this case the frequency of species occurrence across the smallest sampling units (microhabitat). Similarly, CWV_i is calculated following Bernard-Verdier et al. (2012) and Sonnier et al. (2010). For binary traits, CWM_i estimates the total frequency of species sharing the attribute, and CWV_i estimates the variance in this frequency.

A null model approach was used to test whether the observed trait metrics differ from random expectation. Species abundances were randomized across elevations in order to divide existing relationships between trait values and species abundance. Random assemblages were constrained to observed species richness per elevation. Weighted community traits were computed at each iteration, and the 2.5^th and 97.5^th percentiles were obtained from the null distribution to produce null confidence intervals across the elevation gradient.

Comparing the estimated CWM and CWV values to expectations from randomized data allows detection of significant functional structuring in observed patterns. Higher values of CWM compared to random simulations indicate higher trait values at the community level, or higher frequency of species with a given attribute for binary traits. Conversely, values lower than expected provide evidence of ecological filtering towards high or low trait values (Hulshof and Swenson 2010; Bernard-Verdier et al. 2012), analogous to directional selection in evolution. Compared to random expectations, higher CWV values indicate significantly high variability in trait values or, for binary traits, high variability in frequency among species sharing similar attributes. Conversely, lower CWV values suggest reduced variability. Such patterns reflect ecological convergence/divergence for the considered traits and, indirectly, are evidence of competition or ecological filtering (Hulshof and Swenson 2010; Kang et al. 2017).

In order to understand the possible drivers of community functional traits, we built statistical models that relate either elevation only, or a combination of elevation and other environmental variables to each functional diversity measure and CWM. A multiple general linear model was used to select the best combinations of environmental variables associated with community functional composition. For the sake of simplicity, we restrict the analyses to metrics related to average functional composition (CWM) and not variability (CWV). AIC scores were compared among three model types: linear, cubic, and quadratic polynomial, to determine the best fit. The model with the lowest AIC score and highest R² value was taken as best explaining the relationship between the observed CWM trait values and environmental variables. CWM for binary traits was logit-transformed prior to analysis, allowing the consideration of models in a classical Gaussian framework for all response variables.

All analyses were conducted with the statistical software R v.2024.12. (R Core Team 2024) and the packages ade4 v.1.7-23 (Dray 2016), permute v.0.9-7 (Simpson 2016), vegan v.2.6-10 (Oksanen et al. 2015), and FD v.1.0-12.3 (Laliberté and Legendre 2010; Laliberté et al. 2015).

The epiphytic leafy liverwort flora along the altitudinal transect comprises 149 species, distributed among 44 genera and 18 families (Supplementary material 2). Lejeuneaceae are the most species-rich family and comprise nearly half of the total species diversity (70 species in 19 genera). Species richness peaks at 1250 m (69 species) and lowest species richness was found at 650 m (19 species) and 2050 m (19 species).

FEve shows a significant increase with altitude (R² = 0.405, p < 0.05), as does FDis (R² = 0.632, p < 0.01) (Fig. 1). Neither of the two functional diversity indices shows a significant relationship with species richness (Fig. 2).

Figure 1. 

Pattern of functional diversity metrics with elevation. A. Functional evenness. B. Functional dispersion. R² value from a linear regression. * p value < 0.05, ** p value < 0.01. Note that vertical axes start above 0.

Figure 2. 

Linear regression between functional diversity measures and species richness. A. Functional evenness. B. Functional dispersion. R² value from a linear regression and 95% confidence level.

Community weighted means (CWM) for size-related traits (i.e. leaf length, leaf width, gametophyte length, gametophyte width, and stem diameter), as well as for papillae and trigone presence increase significantly with increasing elevation. CWMs for lobule, ocelli, and oil body presence decrease significantly with elevation (Fig. 3, Table 2). CWM for underleaf presence did not show a significant relationship with elevation. A hump-shaped relationship with elevation was found for the elongation index (Fig. 3, Table 2).

Figure 3. 

Distribution of community-weighted means (CWM) of traits (black dots) along the elevational gradient and values estimated by linear models from climatic variables (grey squares). Vertical lines indicate the 95% confidence interval for modeled values.

Considering the intraspecific variability did not improve these relationships except for lobule presence and trigone presence, where the trigone presence decreases with increasing elevation and lobule presence increases with increasing elevation (Fig. 3). The traits showing the strongest relationships with environmental variables are leaf length (R² = 0.9, p < 0.001), gametophyte length (R² = 0.92, p < 0.001) and trigone presence (R² = 0.87, p < 0.001) (Table 2).

Leaf length, leaf width, gametophyte length, gametophyte width, papillae presence, and trigone presence exhibited significant negative correlations with mean canopy height (cmea). Conversely, oil body presence showed a significant positive correlation with cmea. Stem diameter and underleaf presence were negatively correlated with the mean temperature (temp.m), whilst lobule and ocelli presence exhibited significant positive correlation with this variable. Only the elongation index showed a significant correlation with vapor pressure deficit (vpd) and a weak correlation with relative humidity (rh.m) (Table 3). This indicates that a high percentage of variation in these traits is explained by the corresponding environmental variables.

Trait filtering is indicated when the mean value of a trait within a community deviates from null expectations. This occurs when certain environmental conditions at a given elevation selectively favour species with specific traits, causing the community’s trait composition to diverge from what would be expected under random assembly. All traits that we studied, with the exception of the presence/absence of underleaf and ocelli, exhibited departures from the null model at at least one elevational band (Fig. 4). Filtering for lobule absence is evident at 1250–1650 m and oil body presence at 1250–1650 m, and at lower elevation for low values of stem diameter (250–850 m) and gametophyte width (250 m, 560–850 m). Observed values of community-weighted means (CWM) were higher than expected at higher elevation for leaf length (1250–2050 m), leaf width (1650–1850 m), gametophyte length (1250–2050 m), gametophyte width (1650–1850m), stem diameter (1250–1650 m), papillae presence (1450–2050 m), and trigone presence (1250–2050 m). Only the elongation index at mid elevation showed significant deviation above the null model.

Figure 4. 

Distribution of community-weighted means (CWM) of traits (black dots) along the elevational gradient. Envelope (grey) shows the 95% confidence interval under the null model of random community assembly. Diamonds (grey) indicate values estimated by linear models against elevation with 95% confidence interval.

Community weighted variances (CWV) for all traits except gametophyte width and stem diameter also departed significantly from those expected under a null model at least at some elevations. Thus, leaf length (1650–1850 m), leaf width (450 m), papillae presence (1420–2050 m), elongation index (450–1250 m), gametophyte length (1250–1650 m and 2050 m), ocelli presence (250 m), and oil body presence (1250–1650 m) were more variable than expected, indicating that those traits have scattered distributions within communities. This means that abundant species tend to exhibit dissimilar functional trait values in communitites along the elevational gradient. By contrast, trait convergence, indicated by lower-than-random CWV values, is detected at 850–1250 m for underleaf presence, 1250–2050 m for trigone presence, and 250–850 m for lobule presence (Fig. 5).

Figure 5. 

Distribution of community-weighted variances (CWV) of traits (black dots) along the elevational gradient. The envelope (grey) shows the 95% confidence interval under the null model of random community assembly.

Much attention has been focused on the use of functional traits, and their abundance and distribution in communities, in the exploration of relationships between biodiversity and ecosystem functioning (Halpern and Floeter 2008; Cadotte et al. 2009; Mouillot et al. 2011; Roscher et al. 2012; Correia and Lopes 2023). The positive relationship between functional evenness and elevation found here suggests a more regular functional distance among species as elevation increases (Villéger et al. 2008). The fact that more individuals of the common species are recorded as elevation increases suggests that liverwort functional traits are not evenly distributed along the elevational gradient. Functional dispersion quantifies the extent of functional similarity among species in the trait space. The increase in functional dispersion with elevation indicates low redundancy in the liverwort community. This is consistent with the idea that new species added at each elevation are less similar than the existing species and reflects dissimilar functional traits across abundant species (Karadimou et al. 2015; Ferrara et al. 2024).

A main goal of this work was to determine which environmental factor(s) influence the variation and distribution of CWM of liverworts traits across the elevational gradient. This study suggests that canopy height and temperature are the most powerful environmental variables shaping the pattern of trait distributions. Relative humidity and vapor pressure deficit, on the other hand, have a weaker effect on trait variation and distribution among communities.

The leafy liverwort communities showed a clear functional response along the elevational gradient, as demonstrated by the observed clear pattern in community-weighted means for nine of the 12 traits (Fig. 3). The abundance-distribution of all traits related to plant size significantly increased with elevation, indicating that at high elevation there is more variation in plant size (Henriques et al. 2017). Papillae and trigone presence are thought to be related to water retention and xerophytic adaptation respectively. Papillae may facilitate water uptake by creating capillary spaces (Proctor 1979), whilst elongated cells with trigones enhance water uptake from the surrounding environment (Vitt et al. 2014). This study shows that the abundance of species with those traits also increases significantly with elevation, potentially as a response to the more challenging environmental conditions at elevations higher than 1800 m. Lobules are often referred to as water sacs and, like papillae and trigones, an increasing CWM for lobule presence would have been expected. However, its CWMs were negatively related to increasing elevation. A similar result was found for ocelli presence, although the functional role of this feature is yet to be elucidated. Since the lobules and ocelli are characteristic of Lejeuneaceae, this pattern can be related to the fact that the Lejeuneaceae are most abundant at lower and mid-elevations.

Both habitat filtering and niche differentiation appear to be involved in structuring species abundances in the communities studied here. Our results add to the growing body of evidence for the joint effect of these two processes on community structure (Cornwell and Ackerly 2009; Jung et al. 2010; Mason et al. 2011; Maire et al. 2012).

In the system studied herein, trait filtering appears to occur toward both ends of the elevational gradient. This is in line with hypothesised (Mayfield and Levine 2010) ecological filtering at the extremities of environmental gradients. According to Bernard-Verdier et al. (2012), this is due to the differential filtering of dissimilar traits along a gradient. The detection of ecological filtering is critically dependent on the studied trait. Different traits, related to various functional roles, exhibited filtering at different parts of the elevational gradient. At low elevation, where diversity is lower, the range of traits directly linked to growth and nutrient acquisition (gametophyte and leaf size) exhibits filtering (Herben and Goldberg 2014). This can be related to competition for light, given the fact that at low elevations in a tropical forest the canopy is continuous, and the canopy height is the greatest. This may account for the abundance of small and narrow species at these elevations. At low elevations, for instance, the most abundant species include: Bazzania nitida (F.Weber) Grolle, Ceratolejeunea papuliflora Steph., Cololejeunea appressa (A.Evans) Benedix, Lejeunea abyssinica (Gola) Cufod., and Prionolejeunea grata (Gottsche) Schiffn.

None of the studied traits, except for lobule and oil body presence, were specifically filtered at mid-elevation, suggesting that these provide the most favorable habitats for epiphytic liverworts. Nearly all liverworts, except complex thalloid species, possess oil bodies, their filtering at mid-elevations is unexpected and suggests an ecological or physiological constraint. This pattern may result from environmental factors and selective pressure, leading to a shift in survival strategies. Understanding this anomaly may require further understanding and study of the functional role of oil bodies for liverworts and the selective pressures shaping their distribution along elevational gradients.

Interestingly, at mid-elevations trait variation is either randomly distributed or divergent (elongation index only) due perhaps to the relaxation of environmental constraints allowing for the coexistence of a wide range of functional strategies and therefore a peak of diversity at the center of the gradient. Trait range reduction detected at mid-elevation suggests that trait filtering may also occur in more stable habitats (Bernard-Verdier et al. 2012).

Towards higher elevations, leaf length, leaf width, and gametophyte length tend to be divergent in communities. This suggests that traits related to plant size coexist successfully, and species are less similar to each other (Funk et al. 2008). This trend can be explained by the wide range in size exhibited by the most abundant species. Tall liverwort species such as Herbertus dicranus (Taylor ex Gottsche et al.) Trevis, Mastigophora diclados (Brid. ex F.Weber) Nees, and Conoscyphus trapezioides (Sande Lac.) Mitt. ex Schiffn. are about as abundant as small species such as Drepanolejeunea physaefolia (Gottsche) Steph. at higher elevations.

The divergence in leaf length, leaf width, gametophyte length, papillae presence, elongation index, and oil body presence provide strong evidence of niche differentiation, although mechanisms such as limiting similarity or facilitation cannot be inferred here.