Sampling took place across an elevational gradient on the windward slope of Mauna Loa, Hawaii Island, a model system for studies of ecosystem development (Vitousek et al. ¹⁹⁹⁵; Raich et al. ¹⁹⁹⁷). Mean annual temperature (MAT) declines with elevation on this slope at a lapse rate of 6.4°C/km (Juvik and Nullet ¹⁹⁹⁴). Bryophytes were sampled at six elevations (287–2,239 m) on young lava flows (126–152 years old) and old lava flows (≥3,400 years old; Lockwood et al. ¹⁹⁸⁸) at ten sites (Table 1; Fig. 1). Bryophytes occupied seven primary substrates (bark, humus, leaf litter, rock, soil, rotted wood, and tree fern trunk), and accounted for 10–80% of surface cover.

Eleven moss and one liverwort species were collected at ten sites (Table 1). Species were chosen that were sufficiently common to allow collection of replicates at given sites. The species varied in elevational range; the extreme cases were Campylopusincarvatus, found only at the 700 m old flow site, and Racomitriumlanuginosum, which, like the dominant tree species M.polymorpha, occurred at all sampled elevations. We collected samples of 40–150 cm² in projected area from three colonies per species, where possible, at each site (in 6 of 53 cases only a single colony was found for a given species). The projected area of each sample was traced on clear plastic and scanned (300 dots per inch; Epson 3170 Photo Scanner; Seiko Epson, Nagano, Japan) for area determination (using ImageJ 1.42q software; US National Institutes of Health, Bethesda, MD, USA). The substrate and the estimated height above ground for epiphytic colonies were recorded to the nearest 10 cm; non-epiphytic colonies were collected from rocks, fallen logs, and the ground surface. Samples were brought to the laboratory in plastic bags (Ziploc; SC Johnson, Racine, WI, USA) on ice and refrigerated until processing.

We estimated the climate at each sampling site using a Geographic Information System model based on climate station data (Cao et al. ²⁰⁰⁷) and Geographic Positioning System coordinates (GPSmap 60CSx; Garmin, Olathe, KS, USA). We estimated mean annual temperature (MAT), precipitation (MAP), relative humidity (MARH) and vapor pressure deficit (VPD), calculated from MAT and MARH (Campbell and Norman ¹⁹⁹⁸). VPD is a measure of atmospheric drought, i.e., the driving force for evaporation, and can be calculated as an absolute pressure difference (in kPa of vapor) or after normalizing by atmospheric pressure (mole fraction; i.e., unitless); the mole fraction VPD accounts for the effects of atmospheric pressure and temperature on diffusibility and thus more clearly indicates differences in driving force across elevations. Across sites the mole fraction VPD correlated strongly with absolute VPD (r = 0.92; P < 0.001), and the two estimates correlated similarly with other measured traits; only relationships with absolute VPD are presented. Forest overstory cover (OC, %), the proportion of the sky obscured from view (equivalent to canopy closure in Jennings et al. ¹⁹⁹⁹), was visually assessed, to the nearest 10%. Visual canopy cover estimates involve a level of uncertainty, but correlate with measurements using a densiometer or hemispherical photography (Schott and Pieper ¹⁹⁸⁵; Vora ¹⁹⁸⁸; Bellow and Nair ²⁰⁰³; Korhonen et al. ²⁰⁰⁶), especially given training as was undertaken previously to this study, to achieve a correspondence of ±5% relative to hemispherical photos (Teti and Pike ²⁰⁰⁵; Korhonen et al. ²⁰⁰⁶; Paletto and Tosi ²⁰⁰⁹).

CH was determined with a ruler to the nearest mm of the green portion of the canopy at three random positions in each sample. The upper, green portions of the stems were cut from the older brown tissue using a razor blade, dried at 80°C for at least 48 h before determination of dry mass (AB204S/FACT Analytical Balance; Mettler Toledo, Columbus, OH, USA). Canopy mass per area (CMA) was determined by dividing dry mass by ground area. Dried samples were analyzed for δ¹³C, nitrogen per mass (N_mass) and carbon per mass using high temperature combustion in an elemental analyzer (Costech ECS 4010; Valencia, CA, USA), with effluent passed into a continuous flow isotope ratio mass spectrometer (ThermoFinnigan Delta V Advantage with a Conflo III interface; ThermoFisher Scientific, Waltham, MA, USA; Fry et al. ¹⁹⁹⁶). Samples were dry ashed in glass vials (Miller ¹⁹⁹⁸), dissolved in 1 N HCL and analyzed for phosphorus per mass (P_mass) using inductively coupled plasma-optical emission spectroscopy (Varian Vista MPX Instrument; Varian, Palo Alto, CA, USA; Porder et al. ²⁰⁰⁵). Concentrations of nitrogen and phosphorus per ground area (N_area and P_area) were determined as, respectively, N_mass and P_mass multiplied by CMA. N:P was determined by dividing N_mass by P_mass; C:N and C:P by dividing carbon per mass by, respectively, N_mass and P_mass.

Statistical procedures were applied using R 2.6.1 (R Development Core Team ²⁰⁰⁷). We tested for differences in δ¹³C across elevations, soil ages and species using three-way analysis of variance (ANOVA) and across substrate types using one-way ANOVA after log-transformation to improve normality and homoscedasticity (Zar ¹⁹⁹⁹). Multiple pairwise camparisons of δ¹³C among substrate types was calculated with the Tukey honestly significant difference test (HSD; Zar ¹⁹⁹⁹). Subsequent correlation analyses used species mean values for each site. Pearson correlations were determined for linear and power law relationships (i.e., linear relationships after log-transformation) for δ¹³C against environmental variables and bryophyte composition and structure variables. Relationships were tested for species mean values across all sites, and separately across only young and old soil ages independently, and for each of the six taxa that were sampled at six or more sites (five species: A. fuscoflacum, Bazzania cf. trilobata, D. speirophyllum, M. microstomum, R. lanuginosum; and one genus, Campylopus). The δ¹³C relationship with OC was tested across substrate types using raw data rather than species means because a single species could be found on several substrates within a single site.

For relationships that were significant for all species across all sites and also for individual soil ages and/or for more than one individual taxon, we tested for shifts in the relationship. We tested for similarity of the slope and intercept between the soil ages or among taxa, using untransformed or log-transformed data according to which fitted best the relationship for all species across all sites (using SMATR; http://www.bio.mq.edu.au/ecology/SMATR; Warton et al. ²⁰⁰⁶). Ordinary least squares regression was used for testing relationships between δ¹³C and climate variables or OC, and standard major axes were used in testing relationships between δ¹³C and structure and composition traits, as appropriate given co-dependence of variables (Sokal and Rohlf ¹⁹⁹⁵).

When bryophyte traits were correlated with climate variables that were themselves inter-correlated, we calculated partial correlations to resolve underlying relationships (using the corpcor package in R; Schaefer et al. ²⁰⁰⁷).

The landscape-scale increase of δ¹³C with elevation in bryophytes on Mauna Loa was similar to that of M. polymorpha trees at these sites (Vitousek et al. ¹⁹⁹⁰; Cordell et al. ¹⁹⁹⁸, ¹⁹⁹⁹) and that observed generally for vascular plants along elevational gradients worldwide (Körner et al. ¹⁹⁸⁸, ¹⁹⁹¹, ²⁰⁰³; Morecroft et al. ¹⁹⁹²; Hikosaka et al. ²⁰⁰²). This finding extends the elevational trends reported for Sphagnum (Menot and Burns ²⁰⁰¹; Skrzypek et al. ²⁰⁰⁷) to a broader set of genera. In vascular plants, the trend was attributed to reduced CO₂ partial pressure, associated with the lower atmospheric pressure at higher elevations, which would reduce c_c and increase δ¹³C (Körner et al. ¹⁹⁹¹), but later work concluded that such an effect would be too small to explain observed variation (Terashima et al. ¹⁹⁹⁵; Körner ²⁰⁰³, ²⁰⁰⁷). The trend in vascular plants has also been attributed to increased leaf mass per area (LMA) with elevation (Körner ¹⁹⁸⁹, ²⁰⁰³; Vitousek et al. ¹⁹⁹⁰; Cordell et al. ¹⁹⁹⁸; Hikosaka et al. ²⁰⁰²; Takahashi and Miyajima ²⁰⁰⁸; Li et al. ²⁰⁰⁹). A higher LMA in vascular plants may lead to lower mesophyll conductance, lower ratio of chloroplast to atmospheric CO₂ partial pressure (c_c/c_a) and higher δ¹³C (Terashima et al. ¹⁹⁹⁵) because a greater leaf thickness increases the diffusion path length from the stomata to the mesophyll cell walls (Parkhurst ¹⁹⁹⁴), and/or because thicker mesophyll cell walls reduce CO₂ diffusion from intercellular air space to the chloroplasts (Niinemets ¹⁹⁹⁹; Niinemets et al. ²⁰⁰⁹). We found that the δ¹³C decreased with decreasing bryophyte canopy mass per area (CMA), analogously to the decrease of δ¹³C with decreasing LMA in M. polymorpha (Vitousek et al. ¹⁹⁸⁸), but whereas in M. polymorpha and in many other vascular plants LMA increases with elevation (Vitousek et al. ¹⁹⁸⁸; Cordell et al. ¹⁹⁹⁸; Hikosaka et al. ²⁰⁰²; Körner ²⁰⁰³), bryophyte CMA did not correlate with elevation but rather with OC (as discussed further below). These results indicate that the altitudinal trend in δ¹³C was likely due to lower CO₂ pressure and/or the “temperature effect” as shown for three Sphagnum species (Menot and Burns ²⁰⁰¹; Skrzypek et al. ²⁰⁰⁷). Such an effect contrasts with tracheophytes, for which there is no evidence for a general direct effect of temperature on δ¹³C (Troughton and Card ¹⁹⁷⁵; Christeller et al. ¹⁹⁷⁶). There is also the possibility that the declining cloud cover and thus increased irradiance, lower MAP and higher VPD at higher elevation might have contributed to the increase of δ¹³C. However, we found no correlation of δ¹³C with MAP or VPD, and thus such an effect was apparently not important at the landscape scale.

We observed a tendency for lower bryophyte δ¹³C on older soils on Mauna Loa. This finding is similar to that reported for leaves of M. polymorpha which have lower LMA and lower δ¹³C on the young flows (Vitousek et al. ¹⁹⁹⁰). However, for the bryophytes, the effect was weak and apparently linked with the lower OC on young flows.

The strong correlation of bryophyte δ¹³C with OC across the landscape would likely have arisen from several factors. A first candidate is the “source air effect” if incomplete mixing resulted in air enriched in ¹²C due to soil-respired CO₂ at bryophyte-level at sites with higher OC; such an effect would be stronger when canopies are dense, hindering airflow (Schleser and Jayasekera ¹⁹⁸⁵; Broadmeadow et al. ¹⁹⁹²). However, the bryophytes in this study were typically collected <0.5 m above ground, and up to only 5.5 m above soil, and tissue δ¹³C did not correlate with vertical height as expected if soil respiration were a major driver of the δ¹³C–OC relationship. The δ¹³C was significantly higher in colonies growing on rock than on organic substrates, as would be expected if substrate respiration had an important effect. However, colonies growing on humus, where substrate decomposition presumably leads to high CO₂ production, did not differ significantly from those growing on bark substrates, and, most importantly, the δ¹³C–OC relationship held on individual substrates, including on rock. Thus, a vertical gradient in source air composition would likely contribute as a subtle effect, rather than as key driver of the relationship. The strong gradient in irradiance was likely the primary driver of the relationship, as reported for leaves of tracheophytes within forests (Buchmann et al. ^1997b; Sternberg et al. ¹⁹⁹⁷) and within the canopies of single trees (Le Roux et al. ²⁰⁰¹). The δ¹³C–OC relationship may arise from the “irradiance-A effect” effect, i.e. a higher δ¹³C due to increased photosynthetic rate. The relationship might also arise from an indirect effect of water stress in more exposed, drier microhabitats leading to the development of thicker cell walls with greater water storage capacity, which during periods of high water availability would reduce CO₂ diffusion (Dilks and Proctor ¹⁹⁷⁹; Waite and Sack ²⁰¹⁰). Such an effect may hold for individual species, but the lack of a significant relationship of δ¹³C with MAP or with VPD despite their varying by two- to three-fold across sites did not support the importance of such an effect at the landscape scale. We note that variation across taxa in cell size might also be expected to play a role; across three species of microalgae, larger cell size was related to reduced CO₂ conductance and higher δ¹³C (Popp et al. ¹⁹⁹⁸). That scenario is unlikely to be important across bryophyte species, in which larger cells tend to be associated with deeper shade (Waite and Sack ²⁰¹⁰) while higher δ¹³C was related to lower OC and hence greater irradiance.

Thus, the results best support a direct effect of irradiance on photosynthetic rate as a key driver of the bryophyte δ¹³C–OC relationship. Such an effect was found experimentally for two liverwort species grown at four irradiance levels; δ¹³C correlated positively with growth irradiance (Fletcher et al. ²⁰⁰⁶). Our findings extend this pattern across multiple species and show its potential importance at the landscape scale, across sites varying in forest successional stages. The effect is quantitatively strong; when the range of OC was divided into five bins, the average δ¹³C shifted from −27.1 to −30.6‰. The effect was also strong on bark, leaf litter and rock substrates individually; the lack of the correlation on humus may be due to the lower range of OC values compared to the other substrate types. Notably, the correlation of OC and δ¹³C observed across all species and sites was also found within two of the six tested taxa, i.e., for Campylopus and in R. lanuginosum, which had the greatest range in OC values. Whereas these within-species relationships emerged from their wide distribution across OC levels, the across-species relationship reflected the assembly of diverse species along a light gradient (cf. Ackerly and Cornwell ²⁰⁰⁷).

We found that bryophyte δ¹³C correlated with the morphological variable CMA, analogously to LMA for tracheophyte leaves. This finding points to integrated adjustment of canopy form and physiology with irradiance for bryophytes across forest successional stages. However, by contrast with tracheophytes, in which the correlation of δ¹³C and LMA may arise in part from their relationship with mesophyll conductance, our partial correlation analysis indicated that in bryophytes, the relationship of δ¹³C with CMA arose indirectly due to relationships with OC. The high scatter in the correlation of δ¹³C with CMA, and their joint relationship with OC indicates that additional data would be needed to support a direct control of δ¹³C by CMA for any given species.

The bryophytes of this study exhibited great similarities with vascular plants in landscape-scale patterns, showing independent and strong trends for lower carbon isotope fractionation with younger soil age, greater canopy openness and increasing elevation. The similarity with tracheophytes is especially remarkable given the centrality of stomata in explaining carbon isotope fractionation in tracheophytes and the lack of these structures in bryophytes. Our findings point to direct environmental impacts of light and temperature on gas exchange rates as primary mechanisms for these landscape-level trends. We inferred the importance of given potential effects from correlation analyses at landscape scale; individual species may differ in the factors that determine δ¹³C. These findings highlight a necessity for additional mechanistic studies to clarify the underlying basis for these landscape-level trends. Indeed, for bryophytes, we lack a physiological predictor of δ¹³C signals, in the way that stomatal water use efficiency can be used for tracheophytes (reviewed in Seibt et al. ²⁰⁰⁸). This study, showing strong δ¹³C signals across the landscape scale, points to the importance of elucidating the determinants of carbon isotope discrimination in bryophyte structure and biochemistry (see Meyer et al. ²⁰⁰⁸). Such studies will need to determine isotope fractionation in real time for given species under ranges of controlled conditions to deepen the understanding of the landscape scale patterns demonstrated here.

We found the variation in δ¹³C in bryophytes at the landscape scale—including diverse taxa across diverse systems—to be information-rich, reflecting the signal of multiple environmental factors. This finding indicates the possibility of extending the usefulness of preserved bryophyte tissue to estimate past environmental conditions. Previous work emphasized that the lack of stomata can link bryophyte δ¹³C more closely with atmospheric CO₂ concentration than for vascular plants (White et al. ¹⁹⁹⁴; Fletcher et al. ²⁰⁰⁵, ²⁰⁰⁶, ²⁰⁰⁸), and pointed to the δ¹³C of fossilized Sphagnum as a proxy for moisture conditions, given that the “liquid film effect” would lead to higher δ¹³C in wetter climates (Loisel et al. ²⁰⁰⁹; Zhu et al. ²⁰⁰⁹). However, the lack of a relationship in our study between δ¹³C and MAP in the Mauna Loa bryophytes suggests that not all species would be good proxies for moisture conditions, though they may be proxies for other environmental factors. An avenue for future work is to determine how to model the combined effects of elevation, irradiance, soil age and temperature on given bryophyte floras. Alternatively, it may be possible to tease apart specific environmental signals, perhaps by using additional isotopes to carbon. Thus, bryophytes, by occupying, and responding to, a very wide range of environments could provide in their tissues recoverable information not only of physiological responses across gradients but also information of their climates and microhabitats.

Below is the link to the electronic supplementary material.

Species nomenclature follows Staples et al. (²⁰⁰⁴); family nomenclature follows Tropicos, Missouri Botanical Garden [http://www.tropicos.org/]. Lava age is designated as ‘young’ (Y) for 1855–1881 flows and ‘old’ (O) for the 3,400-year-old lava flow (and the ca. 400-year-old lava flow at the 1,750 m site). R. lanuginosum was present but not sampled at the 2,200 m site because of difficulties distinguishing live from dead material

Pearson correlation coefficients, with italicized values derived from log-transformed data, providing better fit and stronger significance. Replication is provided in parentheses in the first row, with exceptions denoted in the table

CAMP Campylopus spp., RALA Racomitriumlanuginosum

Significant correlations in bold: * P < 0.05; ** 0.01 ≥ P > 0.001; *** P ≤ 0.001. No within-species correlations were found for Acroporiumfuscoflavum (n = 7) Bazzania cf. trilobata (n = 7) Dicranumspeirophyllum (n = 6) or Macromitriummicrostomum (n = 9)

Each row represents a separate partial correlation analysis of δ¹³C with the factors in the columns, partialling out the effect of all the other factor(s) in that row. For example, the first value in row (1) presents the partial correlation of δ¹³C and OC, controlling for the effect of elevation

Values in italics were derived from log-transformed data (see Materials and methods)

Bold type indicates correlations that remained significant after partialling out the other factor(s)

CMA, canopy mass per area; MAP, mean annual precipitation; MAT, mean annual temperature; N_area, nitrogen per area; OC, overstory cover

Significance: * P < 0.05, ** 0.01 ≥ P > 0.001, *** P ≤ 0.001. The correlations with elevation and MAT remained significant after partialling out the effect of OC for log-transformed data but not untransformed data. Correlations shown without significance were non-significant using log-transformed or untransformed data