Moss phyllid morphology varies systematically with substrate slope

Background and aims – Tracheophyte leaf morphology is well studied but it is unclear if the findings generalize to poikilohydric plants. We tested combinations of hypotheses to determine if microhabitat characteristics, including light exposure, moisture availability, and substrate slope, controlled for morphological differences between upright and prostrate growth forms, affect phyllid surface area and costa length of mosses. Material and methods – We quantified mean phyllid surface-area and costa lengths for four replicates of 38 moss species from Alabama. Phylogenetic comparative methods that model adaptation were used to evaluate the relative evidence for each hypothesis using information criteria. To further explore mechanistic explanations involving substrate slope, we tested whether mosses on vertical substrates differed from those on horizontal substrates in the average amount of water-retaining, nutrient-rich litter they accumulated. Key results – Substrate slope and growth form combined were the best predictors of phyllid surface area. Mosses growing on vertical substrates exhibited smaller phyllid surface area for both growth forms. Although growth form and phyllid length best explained costa length variation, a more complex model including substrate slope performed nearly as well. Within the prostrate growth forms, species growing on vertical substrates exhibit longer relative costa than those on horizontal substrates. We also estimated rapid rates of adaptation for both traits. Conclusion – The smaller phyllid surface area of both upright and prostrate growth forms is possibly an adaptive response to reduced habitat moisture-retention or nutrient quality that vertical substrates offer. The longer costa lengths of prostrate mosses growing on vertical surfaces relative to prostrate mosses on horizontal surfaces, possibly make up for the decreased ability of smaller phyllids to rapidly reabsorb water when it is available. Further work is required to determine if it is truly substrate slope itself that matters or other variables associated with the differences in slope, and to determine how general this phenomenon is.


INTRODUCTION
Bryophytes diverged from the ancestor of extant tracheophytes between 430 and 451 million years ago (Morris et al. 2018). Whereas tracheophytes developed a pronounced vascular system allowing them to attain large sizes, bryophytes remained comparatively avascular, and consequently much smaller. Mosses (phylum Bryophyta, the most speciose non-tracheophyte plant phylum ;Geffert et al. 2013), exhibit significantly smaller leaf area, leaf mass per unit area, and respiration per unit dry mass of leaf than tracheophytes (Proctor 2000;Glime 2007). The miniaturized leaf-like structures of the moss gametophyte, known as phyllids, are usually a single-cell layer thick and many possess a midrib-like structure known as a costa (Schofield 1981). Whereas tracheophyte leaves regulate water loss through the opening and closing of stomata, moss phyllids lack structural equivalents and thus gain or lose water rapidly to equilibrate with environment water conditions, i.e. they are poikilohydric (Proctor & Tuba 2002). Instead, mosses are known to employ characteristics such as leaf surface roughness and canopy density, that reduce edge exposure to minimize evaporative loss (Rice et al. 2001). The phyllids typically lie within the laminar boundary layer of the bryophyte carpet or cushion, or of the substrate on which it grows (Proctor 2000). Evaporative conditions at the boundary-layer between bryophytes tissues and the environment are largely determined by their growth forms (upright or prostrate sensu La Farge-England 1996), substrate type and exposure to moisture, light, and wind in their habitats (Heijmans et al. 2004).
Among tracheophytes, leaf surface-area directly affects important physiological processes such as photosynthetic efficiency, gas exchange, and evapotranspiration (Wright et al. 2004), and variation in this trait is clearly shaped by niche-mediated natural selection (Donovan et al. 2011). Only a handful of studies, however, have examined how surfacearea of poikilohydric plant leaf equivalents, such as moss phyllids, evolve in response to correlated morphological and/or niche variation. Similarly to tracheophytes, moss phyllid surface-area, costa length, and canopy density in some species are correlated with microhabitat irradiance (Waite & Sack 2010). Bond-Lamberty & Gowers (2007) further demonstrated higher leaf area index (LAI = leaf area/ground area) for boreal bryophytes, indicating an ecophysiological need to maximize light interception for photosynthesis (Whitehead & Gower 2001). The veins of vascular plant leaves, comprising xylem and phloem bundles are well-studied and their roles are clearly defined: xylem conducts water and minerals from the roots, and phloem conducts photosynthates such as sucrose from the leaf. The moss costa, as with tracheophyte veins, is a relatively rigid structure compared to the overall phyllid lamina, consisting of a thick-walled hypodermis made up of enlarged guide and hydroid cells (Thomas et al. 1996). Costa are known to be involved in limited water conductance throughout the phyllid lamina (Scheirer 1983;Thomas et al. 1996) and to provide mechanical support (Vitt & Glime 1984) but other than that, their functional significance is not as clear as that of leaf veins in tracheophytes, and very little is known regarding their evolution.
The present study aims to test evolutionary hypotheses, influenced by the tracheophyte literature, for variation in moss phyllid surface-area and costa length across 38 moss species. The hypotheses centre around the notion that variation in 1) growth form (prostrate versus upright-stem growth forms); 2) light exposure; 3) ambient moisture; and 4) substrate slope affect adaptive optima for the traits. Inclusion of the fourth characteristic is motivated by a hypothesis put forward here that moss species specifically growing on steep or vertical surfaces such as cliff faces or tree trunks, likely do not accumulate as much water-or nutrient-retaining litter. Litter can improve moisture retention, and consequently plant growth, by reducing runoff (Dyksterhuis & Schmutz 1947;Rauzi 1960;Meeuwig 1970), increasing water infiltration (Weaver & Rowland 1952;Dormaar & Carefoot 1996), and reducing evaporative losses through shading (in the case of larger pieces of litter blocking sunlight) and the lowering of soil temperatures (Holland & Coleman 1987;Deutsch et al. 2010). Mosses, unlike tracheophytes, absorb their nutrients from the environments over the surface of the whole plant rather than just the roots (Streeter 1970). Using phylogenetic comparative methods based on models of evolution consistent with adaptation to environmental variables (Hansen 1997;Butler & King 2004;Hansen et al. 2008;Kopperud et al. 2020) and a likelihood-based information theoretic framework (Burnham & Anderson 2002), we test hypotheses based on all combinations of these growth form and microhabitat characteristics.

Moss morphology and habitat variables as selective regimes
Previous studies have established strong effects of plant growth forms (e.g. Horn 1971) and microclimate mediating selection on tracheophyte leaves, as well as similarities in some trends between tracheophyte and moss leaf morphology and how they relate to ecosystem function (Rice et al. 2008(Rice et al. , 2011(Rice et al. , 2018Waite & Sack 2010;Nicotra et al. 2011;Niinemets & Tobias 2019). Motivated by these studies, we hypothesized that morphological and microhabitat characteristics of mosses, including growth form, light exposure (Bell 1982;Guerra et al. 1992;Glime 2007), ambient moisture availability (Niinemets & Tobias 2019), and substrate slope exert differential selective pressures on moss phyllid morphology. To begin with, we a priori chose a list of taxa that: 1) are present in the state of Alabama where the authors are currently based (as reported on the Consortium of North American Bryophyte Herbaria Bryophyte Portal -CNABH 2020); and 2) are also included in the molecular phylogeny of Rose et al. (2016) (see below). We then combed through various sources describing microhabitat characteristics for many of these taxa (CNABH 2020;Shaw & Goffinet 2000;Glime 2007;McKnight et al. 2013;Pope 2016;eFloras 2020). A major challenge with obtaining continuous variables from the literature is homogenizing the data collection procedures of various researchers, often leading to lower sample sizes. Discretization offers a solution to this problem by eliminating the need to work with specific values, and instead using categorized ranges.
Based on the following "rules of thumb" for keywords, the microhabitat characteristics were divided into subcategories as follows: light exposure categories included fully exposed or filtered. Fully exposed taxa were those described as found in wide open spaces, such as meadows or fields, or growing exposed on top of boulders or described as growing in full sunlight. Filtered taxa were those described as being found directly under tree foliage, rocky outcrops, growing on tree trunks, shaded cliff sides, crevices, or on the sides of, or next to any man-made structure capable of blocking sunlight. We identified three subcategories for ambient moisture: moist, moist-mesic, or mesic. Mesic taxa were identified as those growing on substrates with no known water source other than ambient moisture of the environment and that are either exposed to wind and sunlight (e.g. on hardwood tree trunks in sparsely forested areas, or on smooth, dry rock surfaces).
Moist-mesic taxa were those found in shaded, drainage basins (e.g. ditches, temporary streams, flood plains, tree hollows) that occasionally collect standing water. Moist taxa were those found growing on substrates that are typically permanently submerged in water (rivers, streams, lakes, swamps, marshes). Substrate slope was derived from the substrate categories: soil; rock; rock face; gravel; tree stump, tree base; tree trunk and tree branches or fallen tree branches along with our observations of the sample to classify our species into one of two substrate slope categories. Vertical substrates include upright tree trunks, rock faces, walls and branches or tree bases with steep angles (above ~45 degrees). Horizontal substrates include soil, gravel, tree stumps, and rocks, branches and tree bases with shallow angles (below ~45 degrees). The species for analysis were chosen a priori from the above-mentioned shortlist to maximize balance between the number of species in each growth type and across the different microhabitat qualities.

Hypotheses testing with adaptation-inertia phylogenetic comparative analysis
Modern phylogenetic comparative methods, based on appropriate evolutionary process models, provide powerful statistical tools for testing adaptive hypotheses for trait evolution across related groups of species whilst simultaneously controlling for potential pseudoreplication through common inheritance. Here, we use a suite of comparative methodologies (Hansen et al. 2008;Bartoszek et al. 2012;Kopperud et al. 2020) based on an Ornstein-Uhlenbeck (OU) process that uses a likelihood-based information criteria framework (Burnham & Anderson 2002;Butler & King 2004) to compare user-defined adaptive and non-adaptive models of trait evolution. The models underlying these methods assume that species trait values are close to adaptive fitness peaks and then model the movement of these peaks as a function of various hypothesized selective regimes. Traits can be univariate or multivariate continuous morphological or gene expression measures. The methods then test whether groups of species living in distinct niches show systematic differences in the position of their adaptive peaks by estimating a "primary optimum", defined as the average optimum reached by large numbers of species evolving independently for a long time in a specific niche so that the effects of all secondary factors average out. Differences in the primary optima reflect systematic influences of different niches on the positions of adaptive peaks (Hansen 2012). The methods provide a rich modelling framework and are capable of discerning hypotheses of evolution by drift, fluctuating, directional, stabilizing and reciprocal selection and parallel evolution from interspecific data using various information criteria.
The following hypotheses were set up as competing evolutionary explanations for variation in phyllid surface area and costa length: i) the traits accumulate changes randomly over time with respect to the adaptive hypotheses such that their evolution can be approximated by a Brownian motion process; ii) a single global adaptive optimum exists such that changes in trait values are constrained by stabilizing selection around this optimum; or the traits evolve towards and are stabilized around different adaptive optima that are affected by: (iii) the two different growth types (supplementary file 2A); (iv) the two light exposure niches (supplementary file 2B); (v) the three ambient moisture niches (supplementary file 2C); vi) the two substrate slope niches (supplementary file 2D); (vii) all pairwise combinations of niche characteristics; (viii) all three-way combinations of the niche characteristics; and ix) a four-way combination of growth form, light exposure, ambient moisture and substrate slope. For costa length we also include phyllid length as a covariate that affects adaptive optima for costa length.

Moss samples
We collected four specimens for each of the 38 different moss species, chosen as described above. Most samples were collected by the authors on publicly-accessible lands in the state of Alabama from June 2018 until May 2020. The climate of Alabama is classified under the Köppen climate classification system as humid subtropical (Belda et al. 2014). Alabama experiences yearly temperature averages of 17.4°C. It has long, hot, humid summers (June to September), averaging 26.2°C, and mild to cold winters (December to March), averaging 8.1°C (Herbert & Caudill 2012) and an average relative humidity of 71.6%. To validate the general ecological niche categories for our specimens, we cross referenced our own sampling observations with the categorical descriptions derived from field guides and text books. Upon preliminary identification in the field, the samples were brought back to the lab for definitive identification based on microscopic features using dichotomous identification keys (McKnight et al. 2013;Vitt & Buck 2019). Some specimens for analysis were obtained from the Herbarium at the University of Alabama (UNA), either because we could not find any in the field, or to bump up replicate numbers to four for each species (see supplementary file 1 for list of herbarium species). The established protocol for bryophyte herbarium collections uses air drying with light pressure, and is designed to preserve the distinctive leaf shape of bryophytes (Flowers et al. 1945). In addition, insects rarely attack bryophyte herbarium specimens. Together, these facts suggest that moss herbarium specimens do not significantly differ morphologically from their fresh counterparts (see also Espinosa & Pinedo Castro 2018). To validate this assumption, we performed t-tests on the mean area and phyllid length for species that include both fresh and herbarium specimens. These tests did not reveal significant differences (see supplementary file 3).

Litter accumulation between vertical and horizontal substrates
The rationale behind our substrate slope hypothesis is that mosses growing on horizontal substrates, irrespective of growth form, may accumulate more moisture-retaining and/ or nutrient-rich litter than mosses on vertical substrates. To validate this assumption, we collected three fresh moss specimens from three randomly chosen species per category under consideration (i.e. horizontal substrate, upright growth forms: Rosulabyrum capillare, Fissdens bryoides, Dicranella heteromalla; vertical substrate, upright growth forms: Barbula orizenbensis, Bryum pseuodoquetrium, Schistidium apocarpum; horizontal substrate prostrate growth forms: Anomodon attenuatus, Campyliadelphus chrysophyllus, Schwetschkeopsis fabronia, and vertical substrate, upright growth forms: Amblystegium varium, Brachythecium oxycladon, Pylaisia selwynii) and measured out 4.3 cm 3 of material from each moss specimen. The moss material was placed on filter paper and allowed to dry in an incubator at 65°C for 1 hour. Afterward, the material was gently washed with 30 ml of DI water by placing the filter paper in a glass funnel. Washed moss material was removed from the filter paper using forceps, and only organic debris was allowed to remain. The filter paper was returned to the incubator to dry out at 65°C for 1 hour then placed on microbalance to record the dry mass of organic debris. We tested for differences in mulch accumulation across the four different niches using one-way ANOVA, and used Tukey's post hoc tests for all pairwise comparisons.

Imaging and measurement protocols
Fresh moss samples collected from the field were allowed to air dry for 48-72 hours prior to measurement. In general, bryophytes are known to equilibrate rapidly with their surrounding environmental moisture and have the ability to survive desiccation for prolonged periods of time (Proctor et al. 2007). Both fresh samples and herbarium samples were subjected to the same rehydration protocol (dropping 1 ml of DI H 2 O on top of moss sample and allowing to equilibrate for 5 min. prior to photography). We visually inspected branches for signs of heteroblasty and took care not to sample any perichaetal leaves if observed. Other than that, ten hydrated moss phyllids randomly selected from different branches from each of four individual plants for each species were imaged with an Olympus SC180 18MP camera attached to an Olympus SZX16 dissecting microscope at 10× to 20× magnification. ANOVA revealed no systematic branch or specimen related difference in leaf shape or size within a species (see supplementary file 4). Phyllid area (mm 2 ) and costa length (mm) were quantified using ImageJ software (Schneider et al. 2012). Phyllid length was measured from lamina apex to the point of intersection of the lamina with the stem axis as per established protocols (Cho et al. 2007;Rouphael et al. 2010) (see fig. 1 for examples of the different leaf morphologies measured). Rose et al. (2016) provides a maximum likelihood, time calibrated phylogeny based on a supermatrix approach for 3121 moss species constructed from sequenced data for 21 loci including 14 plastid markers, 5 mitochondrial markers, and 2 ribosomal markers. Time calibrations were performed with penalized likelihood based on 12 secondary date ranges obtained from recent estimates in the literature as well as fossil data (Rose et al. 2016). We pruned the tree down to the 38 species found in our dataset using the  Pl. Ecol. Evol. 154 (3), 2021 "keep.tips" function in the ape package (Paradis et al. 2004) for the R statistical environment (R Core Team 2013), which maintains the branch lengths estimated by Rose et al. (2016). The hypothesized selective regimes were reconstructed on the phylogeny using maximum likelihood ( fig. 2 and supplementary file 2) so as to provide estimates, contingent on the Rose et al. (2016) phylogeny, of how long each lineage spent evolving in a given niche. We used the continuous-time Markov model for discrete character evolution (Pagel 1999) as implemented in the ape package for this purpose. For each regime, we contrasted the simplest possible transition rate model (i.e. equal rates of transition for all states) with alternatives encompassing various degrees of complexity depending on the number of niches and types of change. Models were favoured based on information criteria (Burnham & Anderson 2002) and their ability to fully resolve all nodes. Scaled likelihoods at each node were used to specify a given niche in the preceding branch, selecting the ancestral state with highest probability.

Phylogenetic comparative analysis
Primary optima sensu Hansen (1997) for log transformed phyllid surface area and costa length were modelled on the hypothesized selective regimes using the R package slouch (Stochastic Linear Ornstein-Uhlenbeck models for Comparative Hypotheses) (Kopperud et al. 2020), that includes methods capable of modelling adaptation towards primary optima that are dependent on changing environmental variables. These methods assume that the traits of interest for groups of species in the same niche evolve towards, or are kept at a primary niche optimum, θ, (the average expected trait value for numerous species evolving in that niche if adaptation were instantaneous). These niche optima can then be estimated using: current species mean trait values, an estimate of how much time different lineages spent in different niches (the likelihood reconstruction of past niches on the phylogeny), and a mathematical model (based on an Ornstein-Uhlenbeck process) that allows for trait values to move towards changing niche optima through time at a rate determined by a parameter α. The α parameter can be interpreted more easily by transforming it to a "phylogenetic half-life" (t 1/2 = ln(2)/α). The phylogenetic half-life has the same units as the phylogeny branch lengths (i.e. is on a linear scale rather than the exponential scale that the α parameter is on, and in this case the units are millions of years) , and measures how long it takes for the influence of half the ancestral trait values to "decay" as they adapt to new niches (i.e. a measure of phylogenetic inertia). Small half-lives relative to the tree height indicate rapid adaptation to new niches whereas longer half-lives indicate lingering phylogenetic inertia (Hansen et al. 2008). The relative support for hypothesized models for traits evolution can be compared to each other using likelihood-based information criteria (Burnham & Anderson 2002;Butler & King 2004).

Mulch accumulation on horizontal and vertical substrates
Mulch accumulation was significantly different amongst vertical and horizontal substrates nested within upright and prostrate growth forms (One-Way ANOVA, F = 40.23, p < 0.05, fig. 3). A post-hoc Tukey's test for multiple comparisons shows that horizontal upright forms accumulate significantly more mulch than vertical upright forms (p < 0.001) and similarly that horizontal prostrate forms accumulate significantly more mulch than vertical prostrate forms (p = 0.002) (see supplementary file 5 for all pairwise comparisons).

Phylogenetic comparative analysis model selection
The AICc scores (table 1) demonstrate that i) the variation in log-transformed phyllid area is best explained by the combined growth form and substrate slope hypothesis (AICc = 18.25) where 42% of the variation in the trait is explained by this combination, and ii) whereas growth form with phyllid length included as a randomly evolving covariate was the best model for costa length (AICc = 554.78, r 2 = 91%) a model of combined growth form and substrate slope along with phyllid length as a random covariate was within one AICc unit of the best model (AICc = 555.32, r 2 = 93%, 2 AICc units difference is indicative of less supported models in this framework).

Best model parameter estimates
The estimated phylogenetic half-life (t 1/2 ), stationary variance (v y ), and primary optima (θ) for the two best models and the second-best model for costa length (see above for why it is included) are given in table 2. We note that these estimates are contingent on the Rose et al. (2016) phylogeny topology and branch lengths and the ancestral state reconstructions.
Standard errors for primary optima are included in table 2 whereas estimation accuracy for the t 1/2 and v y parameters are given as a joint log-likelihood support surface in fig. 4. The best t 1/2 estimate for the log-transformed phyllid area was essentially zero, indicating rapid adaptation of this trait to new niches. The best t 1/2 estimate for costa length was less than one million years for both the best and second best models, which relative to the tree height of 433 million years can be considered fast adaptation with little phylogenetic inertia. Simulation experience (Hansen et al. 2008) has taught us that when true phylogenetic half-life values are small, variation in tree topology, branch lengths, and internal node states typically do not affect the estimates much, and also that when true half-life values are large, it is highly unlikely that small half-lives will be estimated. The primary optima of the best models for costa length are analogous to a common slopes ANCOVA, where estimated primary optima are regression lines (these are plotted in fig. 5A for the model of growth form + phyllid length and fig. 5B for growth form + substrate slope + phyllid length). The latter shows a clear difference between prostrate growth forms growing on vertical surfaces versus those on horizontal surfaces, but the same distinction is not significant amongst upright growth forms.

DISCUSSION
Mosses are an understudied group of miniaturized plants that can potentially increase understanding of the extensive plant form/function relationships described for tracheophytes. Amongst tracheophytes, particularly the angiosperms, it is generally assumed that the highly variable leaf shapes , and all other models were outside of two AICc units of the lowest AICc score. For costa length, different growth forms was the best model, but the more parameter rich substrate slope combined with growth form model was only ~0.5 AICc units higher than the best model and therefore cannot be excluded. Note that all models for costa length include phyllid length (pl) as a covariate to control for overall phyllid size difference effects.
Table 2 -Parameter estimates for the phylogenetic half-lives (t 1/2 ), stationary variance (v y ), primary optima (θ), and optimal and evolutionary slopes (β) for the best models for each trait. Support for the primary optima and covariate slopes are given as standard errors whereas joint support surfaces for t 1/2 and v y are given in fig. 4.  reflect natural selection operating on their function. Several hypotheses have been proposed to explain leaf shape diversity in tracheophytes, including thermoregulation, hydraulic constraints, patterns of leaf expansion in deciduous species, biomechanical constraints, adaptations to avoid herbivory, and adaptations to optimize light interception (Nicotra et al. 2011).

Trait
In this study, substrate slope was found to significantly affect the primary adaptive optima for phyllid surface area for both types of growth form. Although the upright growth form mosses in this study on average exhibit larger phyllid surface areas than the prostrate mosses (similarly to the study of Niinemets & Tobias 2019), species within each growth form were found to exhibit significantly smaller phyllid surface areas when growing on vertical substrates than those that grow on horizontal substrates. Furthermore, within the prostrate growth form mosses, species growing on vertical surfaces exhibit significantly longer relative costa. Growth form was included as a factor primarily to control for the a priori observed differences in phyllid morphologies associated with upright and prostrate growth forms. Below we discuss possible reasons for why our other hypotheses did not pan out as well as some potential means by which substrate slope may select for smaller phyllid surface areas amongst both growth forms, how it might affect relative costa length within the prostrate growth forms, and how these results relate to findings amongst tracheophytes. First however, a note on the potential role of phenotypic plasticity for the traits under consideration as a caveat and how it may impact the selection-based arguments.
Although our trait measurements are averaged over four different individuals for each species, we recognize that all our samples come from the state of Alabama whereas most of the study species have far wider distributions. Seeing as environmental plasticity in mosses is common, and can affect both phyllid size and costa length, there is no guarantee that our measures are the typical or most representative for each species. Our adaptive reasoning below should therefore be seen in terms of hypotheses that need testing with more data points across a wider range before they can be seen as conclusive.
Our light exposure model performed extremely poorly relative to other models for both phyllid surface area and costa length, and in both cases explained the least variance (r 2 = 0.00 in the case of phyllid surface area). We note however that our categories for light exposure, filtered light vs unfiltered is extremely course and at present it is unclear how well it correlates with microhabitat irradiance or whether it correlates with light interception, which would be more relevant for plant growth, but also more complicated to measure and the latter is affected by numerous other factors, including substrate slope. Waite & Sack (2010) quantified microhabitat irradiance as the % of photosynthetically active radiation above the moss colony relative to open habitat for ten Hawaiian moss species using matched quantum sensors (i.e. all their mosses were found in filtered habitats). They found a 15-fold difference across the species in habitat irradiance and that individual phyllid area was negatively correlated with irradiance. Waite & Sack's (2010) results suggest that much of the variation in light use may lay within the filtered vs open light exposure categories rather than between them, suggesting that future work should quantify microhabitat irradiance differences within each category. Waite & Sack (2010) did not, however, find any relationship between costa length and irradiance, but rather that costa length was correlated with phyllid lenght/width ratio, indicating their role in structural support of longer leaves. Given their poikilohydric nature, we further hypothesized that ambient moisture could exert differential selection on leaf surface area and costa lengths, but as with light exposure, this model did not perform well relative to the best models. Water restriction is one of the major limitations for growth and survival of photosynthetic organisms (Fernandez-Marin et al. 2016). Tracheophytes and poikilohydric plants use different strategies to cope with water restrictions. For tracheophytes, their more sophisticated vascular system, waxy cuticles, and stomatal apertures guard against evaporative loss, significantly easing the burden of the tradeoff between photosynthetic area and evaporative loss (Franks & Brodribb 2005). In the case of mosses, evaporative loss places a significant constraint on the size of moss leaves, necessitating a much smaller phyllid area to mass ratio (Waite & Sack 2010). Some mosses in desert environments, for example, demonstrate greater limitations in terms of phyllid surface area compared to mosses in aquatic habitats (Frahm 1978;Priddle 1979).
We further showed that the horizontal substrates accumulate significantly more litter than vertical ones. Haughian & Frego (2017) tested the connection between substrate moisture retention and moss growth. They found that surface humidity was positively associated with organic debris. In addition, moss species richness and growth have been found to be positively correlated with substrate-level moisture and negatively correlated with temporal fluctuations in substrate level moisture (Oberndorfer 2006;Stewart & Mallik 2006). Together, these findings suggest that mosses are highly dependent on the immediate micro-humidity of their substrates, which in turn is modulated by the degree of moisture retaining litter they accumulate. All else being equal, vertical substrates should drain more rapidly than horizontal ones, and given that accumulated litter acts as a moisture-retaining mulch, vertical and horizontal substrates are expected to differ in the rate at which water drains from the microhabitat. Our models that included substrate slope along with growth form outperformed all others tested here for phyllid surface area. We suggest that one possible reason for the smaller surface areas exhibited by mosses growing on vertical substrates reflect the inability of poikilohydric mosses to generate increased photosynthetic surface area without increasing net evaporative water loss, a factor compounded by rapid water drainage from the substrate. However, the opposite problem arises when it is time to reabsorb water from the environment, where a larger surface area would be more advantageous. Below we discuss how the costa may play a role in mitigating this problem.
There is some evidence that leaf costa serve to conduct water in moss phyllids from apex to base and vice versa similar to leaf veins in tracheophytes (Vitt & Glime 1984). Thus far, evidence that relative costa length evolves in response to niche conditions has been lacking (Proctor 1979;Glime 2007). Our best model for costa length was one where adaptive optima were modelled on a combination of growth form and phyllid length. On average, for a given phyllid length, costa length were longer in upright than prostrate growth forms, and all species with no costa were prostrate growth forms in this study ( fig. 5A). The prostrate growth form species studied here all also have much shorter phyllid lengths, suggesting that either water transport, or structural support provided by the costa are less important for this growth form. Water conductance in mosses, especially the prostrate growth forms studied here, that all belong to the Hypnales order within the homocostate pleurocarps, is known to be enhanced by various morphological features including phyllid lamellae, papillae, mammillae, paraphylia, rhizoids, axillary hairs, and differentiated alar cells (Proctor 1979(Proctor , 1982, and future studies should account for these when further examining the hypothesis put forward here. Furthermore, Huttunen et al. (2018) argue that i) since the first major burst of diversification amongst the Hypnales occurred after the evolution of a homogenous costa with undifferentiated cells, ii) there are no known reversals of this character state, and iii) since most Hypnales are found in humid or wet environment, the simplified costa are an adaptation to such environments and are therefore also associated with the evolutionary success of this group. This also raises the interesting, but untested possibility that the phyllids amongst prostrate growth forms are in general smaller than upright growth forms, because they are constrained by a simplified costa. In support of the costa's role in water relations, moss species growing submerged in water, tend to exhibit reduced or absent costa compared to species growing on dryer substrates (Bell 1982;Guerra et al. 1992;Glime 2007). A similarly verifying phenomenon is observed within some aquatic species that have been shown to suppress costa growth as a plastic trait when grown in submerged conditions (Vitt & Glime 1984). The next best model for costa length could not be rejected as it fell well within 2 AICc units of the best model and shows that within prostrate mosses, species growing on vertical substrates exhibit increased relative costa lengths. This suggests the possibility that the increased relative costa length amongst prostrate mosses on vertical surfaces possibly allow moss phyllids to spread residual water across more of its cells thereby buying more active metabolism time to enact cellular desiccation response mechanisms or that it makes up for the smaller phyllid areas by facilitating better absorption of water from the microhabitat when water becomes available. These suggestions are not necessarily mutually exclusive, but do require further experimental evidence for verification. Thus, the results of the present study suggest that, at least within the prostrate stem growth forms, costa length responds to environmentally determined optima where longer relative costa make up for smaller phyllid surface area when it comes to reabsorbing water from the environment. The non-zero, but relatively small (~1 myr) phylogenetic half-life estimate for both phyllid surface area and costa length evolution also suggests some small degree of phylogenetic inertia on shorter time scales, providing some evidence against the notion that we are merely observing a plastic response, which should result in a half-life of zero.
Bryophyte leaves often exist within a canopy, and while this topic remains understudied, some important efforts seek to characterize the interplay between leaf and canopy scales (e.g. Rice et al. 2008Rice et al. , 2011Rice et al. , 2018Waite & Sack 2010;Niinemets & Tobias 2019). The canopy leaf area index (LAI) is a whole plant measurement compromised of leaf size, leaf frequency, leaf area per shoot height, shoot height, shoot leaf area, and shoot number per area, that characterizes the ability of plant canopies to exchange energy and matter with the environment. LAI allows for assessing the contribution of moss to total ecosystem productivity. Niinemets & Tobias (2019) found that even though individual phyllid surface area exhibits high variation across species, it trades off with phyllid number and therefore does not contribute as much to LAI as would be expected. Here, we found a systematic decrease in individual phyllid area with increased substrate slope in both growth forms. We acknowledge that other characteristics of the substrates will likely affect moisture retention ability (soil cover, tree bark vs bare wood, rock type, etc.) and that these could be a fruitful avenue of future research. However, our results also suggest that substrate slopes should be investigated further across a wider range of taxa as potentially very general morphological and ecological characteristics affecting phyllid features and therefore subsequently other traits, such as phyllid number that contribute to moss canopy features. To close with, we echo the argument of Huttunen et al. (2018) that "understanding of the potential function of morphological traits and knowledge of the distribution of morphological characters in various environmental conditions is necessary to understand functional differences among moss species and the communities they belong to".

SUPPLEMENTARY FILES
Supplementary file 1 -An excel file with sheets containing field collection descriptions, raw measurements of all specimens, mean values for raw measurements, and niche descriptions gleaned from observations, online data bases, field guides, and the primary literature.