Environmental drivers of Sphagnum growth in peatlands across the Holarctic region

The relative importance of global versus local environmental factors for growth and thus carbon uptake of the bryophyte genus Sphagnum—the main peat‐former and ecosystem engineer in northern peatlands—remains unclear. We measured length growth and net primary production (NPP) of two abundant Sphagnum species across 99 Holarctic peatlands. We tested the importance of previously proposed abiotic and biotic drivers for peatland carbon uptake (climate, N deposition, water table depth and vascular plant cover) on these two responses. Employing structural equation models (SEMs), we explored both indirect and direct effects of drivers on Sphagnum growth. Variation in growth was large, but similar within and between peatlands. Length growth showed a stronger response to predictors than NPP. Moreover, the smaller and denser Sphagnum fuscum growing on hummocks had weaker responses to climatic variation than the larger and looser Sphagnum magellanicum growing in the wetter conditions. Growth decreased with increasing vascular plant cover within a site. Between sites, precipitation and temperature increased growth for S. magellanicum. The SEMs indicate that indirect effects are important. For example, vascular plant cover increased with a deeper water table, increased nitrogen deposition, precipitation and temperature. These factors also influenced Sphagnum growth indirectly by affecting moss shoot density. Synthesis. Our results imply that in a warmer climate, S. magellanicum will increase length growth as long as precipitation is not reduced, while S. fuscum is more resistant to decreased precipitation, but also less able to take advantage of increased precipitation and temperature. Such species‐specific sensitivity to climate may affect competitive outcomes in a changing environment, and potentially the future carbon sink function of peatlands.


| INTRODUC TI ON
Net primary production in peatlands is relatively low, but because, in general production exceeds decomposition, peatlands have remained important carbon (C) sinks throughout the Holocene (Yu, 2012). As a result, northern peatlands store c. 500 Gt C (Loisel et al., 2014;Yu, Loisel, Brosseau, Beilman, & Hunt, 2010). A recent analysis suggests that the stock may even be >1,000 Gt (Nichols & Peteet, 2019), which is more than the C in the atmosphere today (829 Gt;IPCC, 2013). In northern peatlands, the growth of peat mosses (genus Sphagnum) is vital for C sequestration and storage, as they effectively engineer a wet and acidic environment that inhibits decomposition (Rydin & Jeglum, 2013). However, despite their crucial role in peatland C cycling, we still lack a clear understanding of the factors controlling Sphagnum growth at global and local scales.
The important role of Sphagnum in C sequestration suggests a need to include Sphagnum growth in terrestrial ecosystem models.
For example, Charman et al. (2013) suggested that biomass input (i.e. NPP) has driven the peat accumulation rate over the last millennium, while variation in decomposition has mattered less. Peatland models have recently included Sphagnum-specific growth functions (Turetsky et al., 2012), but peatland modules are generally lacking in Earth system models (e.g. ORCHIDEE, Qiu et al., 2018), but see Bechtold et al. (2019). A literature survey revealed that attempts to link variation in Sphagnum NPP to important environmental drivers have come to different conclusions (Table 1). Gunnarsson (2005) found that temperature, precipitation, altitude and latitude explained 40% of the variation in productivity. Moore (1989) identified annual mean temperature as a key driver of Sphagnum production, while Krebs, Gaudig, and Joosten (2016) reported that growth of S. papillosum was primarily influenced by water availability (precipitation frequency and precipitation/temperature quotient). Loisel, Gallego-Sala, and Yu (2012) found that cumulative photosynthetically active radiation (PAR) for days over 0°C was the most important driver of growth.
Globally, terrestrial NPP is largely determined by precipitation and water availability, and in northern regions (above 50°N) temperature and solar radiation become increasingly important (Gallego-Sala et al., 2018;Schloss, Kicklighter, Kaduk, Wittenberg, & The Participants of the Potsdam NPP Model Intercomparison, 2001).
The contrasting findings of these previous studies may result from inconsistencies in methodology. First, meta-analyses use data collected using different methods, in different habitats (e.g. pH and degree of ground-water influence) and over different time periods (Gunnarsson, 2005;Krebs et al., 2016;Limpens et al., 2011;Loisel et al., 2012). Second, many of the previous studies are geographically limited (Asada, Warner, & Banner, 2003;Moore, 1989;Walker et al., 2017) and consequently cover a narrow range of climatic conditions. Third, many studies include different species for different areas and species can covary with environmental conditions, which may add bias and uncertainty to the results. To overcome these problems we performed a Holarctic, coordinated sampling effort, focusing on common but ecologically divergent species.
In addition to climatic drivers, local conditions and air pollutants such as nitrogen (N) and sulphur (S) deposition affect Sphagnum growth. Limpens et al. (2011) examined the effects of increased N deposition in a meta-analysis comprising 29 studies across Europe and North America. They showed that small N additions to sites with low background N deposition stimulated growth, while adding higher levels of N depressed growth. Interestingly, increased vascular plant cover caused the shift from a positive to a negative growth response to N to occur at lower N addition levels (Limpens et al., 2011). Experiments have shown that N addition leads to increased cover of vascular plants that will reduce light availability, increase litter and lower surface temperature, thereby reducing Sphagnum productivity (Berendse et al., 2001;Chong, Humphreys, & Moore, 2012). Moreover, Hayward and Clymo (1983) reported reduced moss mass growth when PAR was lowered. In addition to vascular plant cover, the position above the water table, which is a proxy for water availability, has been correlated with Sphagnum NPP.
Production was lower for Sphagnum growing higher above the water table (e.g. Gunnarsson, 2005;Krebs et al., 2016). These local factors and N availability are known to modify Sphagnum NPP, but to what extent relative to climatic drivers has not been investigated.
Here, we set out to test the importance of global and local environmental drivers of length growth and C sequestration potential (i.e. NPP) of two common peat-forming Sphagnum species in northern peatlands. Length (height) increment is essential for the shoot to keep up with vascular plants and to compete for light with neighbouring Sphagnum shoots. NPP is more directly related to peatland C sequestration, and depends on the combination of length increment (LI) and the bulk density of the material produced. The two TA B L E 1 Previous studies relating Sphagnum growth to climatic variables in temperate and boreal peatlands Sphagnum species, S. fuscum and S. magellanicum were selected to represent different niche preferences within peatlands, i.e. lawns and hummocks. Our aim was to explain and predict variation in Sphagnum productivity within sites and across the Holarctic using statistical modelling-including structural equation models (SEMs, Shipley, 2009). By sampling in a consistent and standardized manner across sites and years, we tested the relative contributions of previously identified environmental drivers (e.g. water table position, temperature, precipitation, PAR) to growth. We did this by exploring both direct and indirect effects (e.g. impact on vascular plant cover, shoot density) of these drivers. This allows us to test to what extent climate-change effects on Sphagnum will be modified and/or mitigated by local conditions and if the response is speciesspecific. More broadly, improving our understanding of the drivers of Sphagnum NPP will support forecasting of changes in peatland carbon dynamics under ongoing climate change.

| Study species and sampling design
Two Sphagnum species belonging to different subgenera (sections) were selected for this study: Sphagnum fuscum (Schimp.) H. Klinggr. and S. magellanicum Brid. Sphagnum fuscum in the subgenus Acutifolia is characterized by its brown colour and its small and densely packed capitula ( Figure 1). It typically forms hummocks in bogs, but also in minerotrophic peatlands where its hummocks often are high enough to be ombrotrophic. Recently, some individuals of S. fuscum in Europe have been identified as conspecific to the North American S. beothuk R. E. Andrus (Kyrkjeeide et al., 2015). Genetic analyses performed on one sample from each of 26 of our 102 sites indicated a small proportion (two samples) of S. beothuk (N. Yousefi, K. Hassel, & H.K. Stenøien, unpubl.). Sphagnum magellanicum belongs to the subgenus Sphagnum, which is a group of large and stout species with cucullate leaves. Sphagnum magellanicum stands out in its subgenus due to its colour that varies from deep red to almost completely green ( Figure 1). Also, recently, S. magellanicum in northern peatlands has been split into two species: S. medium Limpr. that commonly occupies the bog expanse, where it grows close to the water table and at relatively intense light conditions, and S. divinum Flatberg & Hassel that occupies bog margins and poor fens, growing further from the water table and at lower light intensity (Hassel et al., 2018). Genetic analysis (N. Yousefi, K. Hassel, & H.K. Stenøien, unpubl.) on one sample each from 62% of our sites indicates that both species were included in approximately equal proportions. Given that S. beothuk/ S. fuscum and S. medium/S. divinum are sister species and were treated collectively in the past, and that they had not been differentiated at the time of data collection, we here follow Granath et al. (2018) and treat our study species as S. fuscum coll. and S. magellanicum coll.
In 2013 and 2014 we sampled throughout the species' Holarctic distributions to cover the climatic variation of sites where the species co-occur. In total, we sampled at 102 sites, of which three sites were excluded from the analysis due to inconsistent sampling (Hani and Mangui, China) or lack of sample tissue N concentration data  Table S1). The sample sites covered a large proportion of the boreal and temperate biomes where northern peatlands are present ( Figure 3).
Our aim was to sample from similar habitats across the sites, and we therefore sampled species in their typical open, nutrient-poor habitats, ranging from predominantly bogs to poor fens. At each site, four mono-specific, homogeneous non-treed patches (1 m 2 ) were chosen for each species, with each patch at least 10 m apart.

| Length growth and production
Sphagnum growth is easy to measure as they have no roots or rhizoids (Shaw et al., 2010) and new growth is confined to the top of the shoot. We assessed Sphagnum growth using two different approaches; LI (mm), and biomass production (NPP; g/m 2 ). To measure LI we inserted a minimum of three 'brush wires' in each patch at the start of the growing season ( Figure S1). The brush wires resemble bottlebrushes and were constructed by attaching paintbrush bristles to c. 10-15 cm long wires. This is an improved version of the cranked wire method described in Clymo (1970), and is now the most commonly applied method (Rydin & Jeglum, 2013). The brushes were inserted c. 5-10 cm into the Sphagnum carpet without disturbing the moss surface by initially pressing the bristles into a narrow tube, which was then pulled out while the brush was left in the Sphagnum carpet allowing the bristles to spread and 'anchor' in the vegetation. The height of the part of the brush wire protruding from the Sphagnum carpet was measured each year at the beginning of the season, and then remeasured at the end of the season. The difference between spring and autumn measurements gave us the LI (mm). Measurements with large negative values (<-5 mm growth) indicated physical disturbance and were removed from analyses.
To estimate stem bulk density and shoot density, we collected one moss core (5-10 cm diameter) from each sampling patch at the end of the season at the same time as the final wire length measurements. The capitula (i.e. the moss top with a height of c. 0.5 cm in S. fuscum and 1 cm in S. magellanicum) were removed and dried, and the 3-cm stem sections immediately below the capitula were dried (24 hr at 60-65°C) and weighed to obtain stem bulk density (BD; kg/m 3 ). The area-based biomass growth (NPP; g/m 2 ) was calculated as BD × LI. In addition to measurements of LI and production for the growing season, we calculated LI and NPP per day during the sampling periods. We also counted the number of capitula in the core to obtain numerical shoot density (cm -2 ).

| Local predictors: Water table and vascular plant cover
For each patch and year, we recorded two variables that we expected would relate to growth: height of the moss surface above the water table (HWT; cm) and vascular plant cover (vertical projection on 1 m 2 surrounding the patch, %) at the start and end of the season. For a few sites, these measurements were only taken once per year at either the start or end of the season. For sites where these metrics were made both early and late in the season, the correlation between start and end values of HWT was r ≈ 0.80 (df = 83 in 2013, and 85 in 2014) and between start and end values of vascular plant F I G U R E 2 The sampling design included 102 sites distributed across the Holarctic region where either Sphagnum fuscum or S. magellanicum, or both were sampled. At each site and for each year, we aimed to sample four patches (red circles) of each species if they were both present at the site, and at each patch we took three measurements of length increment, and one measurement of bulk density to calculate biomass accumulation. Illustration: F. Bengtsson Holt Lodge, UK. For total C and N, around 10 mg (6 decimal balance) of milled and dried (70°C) sample was analysed using the combustion method with a Carlo Erba CN analyser (Flash1112 series). For cations and trace metals, 0.6 grams of plant material was digested with a microwave system (Anton Paar Multiwave 3000) using 7.5 ml of HNO 3 at 185°C for 30 min (20 min ramp from ambient to 185°C).
The extract was then adjusted to 25 ml and analysed with an ICP-OES dual view (Thermo ICap 6500). Most samples were from 2013 and in general two patches per species and site were analysed. Total element concentrations were expressed as percentages or mg/kg. In this study we use mean N concentration per site, but include all data in our publication for completeness. W/m 2 ), time during an hour with no precipitation (HOURNORAIN; s) as predictors. HOURNORAIN was recalculated as the mean rain free period (average time between rain events; days). These data are derived from satellite and meteorological station observations and are gridded with a grid cell resolution of longitude 0.667° and latitude 0.5°. We extracted data for the periods of growth measurements at each site, and calculated the average temperature (°C), total precipitation (kg/m 2 ), evaporation (kg/m 2 ) and PAR (PARDF+PARDR, W/m 2 ) and the average number of consecutive days without rain (discarding days with <1 mm precipitation; d; Table S2). Data on nitrogen deposition were extracted from the model synthesis of Lamarque et al. (2013). They modelled global total N deposition for the time period 1995-2005 (model output given as an estimate

| Global predictors: Weather and nitrogen deposition
for the year 2000) with a resolution of 0.5° longitude by 0.5° latitude. Elevation and coordinates were extracted from Google Maps (Table S1).

| Data analysis
To test the importance of the controlling factors for Sphagnum growth, we selected variables based on previous studies (see Table 1). Temperature and precipitation have been used in most studies and were the main predictors in the meta-analysis by Gunnarsson (2005). He also included the geographical variables latitude, elevation and distance from the sea, but we did not test these variables as they are not mechanistically related to growth, and their direct effects on growth should be included in the climatic data. Loisel et al. (2012) compiled global growth data for S. fuscum and S. magellanicum, and modelled seasonal LI as a response to growing season F I G U R E 3 Mean annual temperature (°C) and annual precipitation (mm) extracted from WorldClim (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) for the study sites (black rings) superimposed on the Whittaker biome chart (Biome chart adopted from Ricklefs (2008) plotted with the R package plotbiomes, https://github.com/valen tinit nelav/ plotb iomes) length, PAR0 and P/E q (index of moisture balance: the quotient between annual precipitation [P] and the annually integrated equilibrium evapotranspiration [E]). We included PAR as a predictor and tested if P/E or P-E would perform as better predictors compared to P and E individually. An additional climatic variable was average length of rain free periods (Krebs et al., 2016). Nitrogen availability has been related to Sphagnum growth (Granath, Strengbom, & Rydin, 2012;Limpens et al., 2011) and our model included N deposition (site level) and tissue N concentration (site level but speciesspecific). Within-site predictors were vascular plant cover (Gavazov et al., 2018), shoot numerical density and height above the water Since we sampled approximately four patches per species and site, we used 'site' as a random factor in our models, and the site effect was allowed to vary between years (random intercept-slope model).
As a first step, we fitted LMMs (classic multivariable regression) including the global (temperature, precipitation, evaporation, PAR, average length of rain free periods, N dep) and local (HWT, vascular plant cover, N tissue) predictors, species and year as fixed factors, including interaction effects between global variables and species.
We ran the models without interaction effects to explore main effects and their R 2 values. p values for interaction and main effects were extracted from type II ANOVAs using the car package (Fox & Weisberg, 2011). R 2 values can quantitatively evaluate the fit of the model and were produced for site (R 2 site , explained between site variation) and within-site level (residual error) by calculating the proportional reduction in these variance components under the fitted model when compared to the null model (following Johnson, 2014;Nakagawa & Schielzeth, 2013 and using code from the r package mumin, Bartoń, 2018). To give a measurement of overall model fit we also estimated the amount of total variance explained by fixed effects (Marginal R 2 ). A parametric bootstrap procedure was employed to produce approximate confidence intervals (CIs) of R 2 s and we present the 90% percentile CI.
As a second step we extended the above analyses to a pSEM to also examine indirect effects. A piecewise SEM is flexible as it allows multilevel structures and that each response (endogenous variables; here NPP or LI, shoot density, tissue N concentration, vascular plant cover and HWT) is modelled individually (Shipley, 2009).
We performed SEMs for LI and NPP, and due to important species interactions found in initial multivariable regressions, the SEMs were done for each species. The paths included in the SEM (i.e. model specification) were the same for each species and LI and NPP. In principle, all paths were included, but we removed paths for which there were no clear mechanisms. For example, many of the variables cannot affect HWT, such as tissue N concentration, and tissue N concentration is expected to be modified by vascular plant cover, but vascular plant cover is not affected by tissue N. The complete model specifications can be found in Tables S7 and S8. To simplify the figure illustrations of models we removed links that had minor impact on the response (p ≥ 0.1). To be able to compare direct and indirect effects we standardized both endogenous and exogenous variables by dividing by one standard deviation.
Standard residual analyses were performed to check homogeneity and normality of errors, distribution of random effects and influence of individual data points. Residuals showed increasing variance and Box-Cox transformations were performed using the function powerTransform in the car package (ver 3.0-0, Fox & Weisberg, 2011) to achieve homogeneous errors. We used correlation analysis to examine potential issues with multicollinearity.

| Descriptive statistics
Across species and years, growth measurements showed as much variation within sites as between sites. The proportion of within-site variance to total variation was 47% for LI and 57% for NPP measured on an annual basis (year −1 ; Table S3). The larger species, S. magellanicum, had on average greater LI than S. fuscum, while S. fuscum had somewhat higher NPP than S. magellanicum due to its higher bulk density ( Figure 4; Table S5).

| Multivariable regression models
Our full models with LI as response explained 53% of the growth calculated per day (day −1 ; CI 40%-64%) and 52% on an annual basis (year −1 ; CI 38%-63%) of variance between sites (Table S4). The models explained less variation in NPP between sites than in LI, and models showed similar R 2 site values on a daily basis (31%; CI 11%-46%) and annual basis (32%; CI 12%-47%). The correlations between LI and NPP were r = 0. We detected species-specific responses to precipitation (LI and NPP), temperature (LI), evaporation (NPP) and numerical density (NPP ; Table S4). These interactions suggest different effect sizes and similar directions of the response (e.g. the positive effect of precipitation was weaker and often not distinguished from zero in S. fuscum compared to a clear positive effect on S. magellanicum in all models).
Using precipitation minus evaporation (precipitation surplus, P-E), as a predictor instead of using the two variables as separate predictors gave similar results (Table S6).
For comparison with Loisel et al. (2012) we ran the LI model using PAR, species and year as predictors, which gave an R 2 site of 22% and an R 2 marginal of 12% respectively.

| Structural equation modelling
Using structural equation modelling built from mixed effects models, we investigated direct and indirect effects of studied variables on length (LI; Table 2; Figure 5) and biomass growth (NPP; Table 3;

| Direct and indirect effects of global predictors
Regardless of species, temperature had positive direct effects on growth (Tables 2 and 3), and the effect was strongest for LI in

| Direct and indirect effects of local predictors
Among predictors varying within site, shoot density (i.e. Sphagnum

| Climatic predictors
Precipitation was expected to increase growth (Asada et al., 2003;Backéus, 1988;Bragazza et al., 2016;Gunnarsson, 2005;Moore, 1989) as it helps the poikilohydric Sphagnum maintain a sufficient moisture content to remain active and photosynthesise (Schipperges & Rydin, 1998). In our study the effects of precipitation were speciesspecific and a positive effect on growth was only observed for the wetter-growing species, S. magellanicum. However, more general support for the importance of precipitation was indicated in our data by a consistent negative effect on growth of mean length of rain free periods (cf. Krebs et al., 2016). Interception by vascular plants means that a minimum amount of precipitation per rain event is required to add moisture to Sphagnum; to account for this we counted days F I G U R E 5 Graphical representation of piecewise structural equation models using mixed effects models to describe effects of the environmental variables from the models in Table S4 on annual length increment. Num density = numerical density of Sphagnum shoots, PAR = photosynthetic active radiation, No rain = mean length of rain free periods, HWT = height above the water with precipitation <1 mm as rain free (Farrick & Price, 2009;Price, Rochefort, & Quinty, 1998). Hence, our results imply that extreme dryperiods affect Sphagnum growth but that the amount of precipitation is a better predictor of growth in non-hummock species.
Evaporation is often incorporated with precipitation into a moisture index to reflect water availability. We calculated the amount of available water (P-E), but our results corroborate those of Loisel et al. (2012) who found no effect of this differential on growth.
Keeping precipitation and evaporation as separate predictors resulted in a better fitting model.
Temperature is a major climatic driver of plant production. There is a mechanistic relationship between temperature and photosynthesis (Farquhar, von Caemmerer, & Berry, 1980), exhibiting a unimodal pattern where temperature optimum seems to be rather high in Sphagnum (Harley, Tenhunen, Murray, & Beyers, 1989;Skre & Oechel, 1981). As expected, our results consistently show that warmer temperatures result in increased growth, but this effect was weaker in S. fuscum than in S. magellanicum. In fact, the summed effect of temperature on NPP in S. fuscum was zero due to the indirect effect of reduced density. The strong temperature effect on LI in S. magellanicum illustrates that this species can grow well in length even with little increase in biomass accumulation (Mazziotta, Granath, Rydin, Bengtsson, & Norberg, 2019), the result of which would be patches with low bulk density. For S. fuscum this would not be a viable strategy since maintaining high BD is crucial to maintaining a high water content (McCarter & Price, 2014;Nijp et al., 2014).
PAR is strongly related to temperature and directly controls NPP (Chapin III, Matson, & Vitousek, 2011) because PAR is a key driver of carbon fixation in plants. Using PAR as the only predictor aside from year and species, we found some support for Loisel, Gallego-

Sala, and Yu (2012) results that PAR is an important driver of NPP in
Sphagnum. When controlling for other factors (i.e. our full multivariable models), there was no or little effect of PAR on length growth, but still a positive effect on annual NPP. This indicates that other variables, most notably temperature, with which PAR was correlated and which also control evaporation, were the main drivers of growth.

| Effects related to vascular plant cover
Low levels of N deposition can increase Sphagnum growth in N-limited systems (Limpens et al., 2011), but evidence of large-scale negative effects of greater N deposition on Sphagnum growth are lacking.
Experimental evidence shows that high levels of N deposition will promote vascular plant growth, and thereby intensify the effects of competition for light from vascular plants (Berendse et al., 2001;Bubier, Moore, & Bledzki, 2007;Limpens et al., 2011). Indeed, one of Water table position has frequently been associated with Sphagnum growth, suggesting that a more surficial water table promotes F I G U R E 6 Graphical representation of piecewise structural equation models using mixed effects models to describe effects of the environmental variables from the models in Table S4 on annual net primary productivity. Num density = numerical density of Sphagnum shoots, PAR = photosynthetic active radiation, No rain = mean length of rain free periods, HWT = height above the water

Net primary production
Sphagnum growth (Gunnarsson, 2005;Li, Glime, & Liao, 1992). Our models did indicate a weak negative effect on growth with deep water tables through indirect effects (vascular plant cover and shoot density). Interestingly, for S. fuscum there was a weak positive effect of HWT on NPP due to higher shoot density. Likely, the observed weak effect of water table on growth cannot be easily disentangled from the species effect, as the HWT effect appears stronger between species than within species (Bengtsson, Granath, & Rydin, 2016). These findings suggest that the effects reported in experimental studies (e.g. Li et al., 1992) are weaker in the field. It is possible that strong effects are only apparent during a short period of the season (Walker et al., 2017), or in dry years when the water table is significantly lowered for an extended period (Bragazza, 2008). Unless the drought is severe, hydrological feedbacks will maintain sufficiently high moisture content for photosynthesis across the water table gradient within a species niche (Waddington et al., 2015).

| Comparing species, years and growth responses
The two species show large differences in LI and bulk density, but less so in NPP. In fact, NPP was higher in S. fuscum even though this species had lower LI. This has been reported previously (Lindholm & Vasander, 1990) and is a result of the high BD in S. fuscum, and NPP being the product of LI and BD. The underlying mechanism behind the rather similar NPP is a trade-off between maximizing length growth and maintaining high water content (Laing, Granath, Belyea, Allton, & Rydin, 2014;Mazziotta et al., 2019). Generally, S. fuscum grows in higher hummocks, where its higher bulk density and shoot density increases capillary water uptake and reduces evaporative loss. Sphagnum magellanicum on the other hand, mainly grows in lawn habitats and relies on being closer to the water table to maintain a high water content.
Hence, S. fuscum is less dependent on a continuous wet climate than other peat mosses, and can remain photosynthetically active during drought events and during dry years (Rydin, 1993). This can explain the stronger effect of precipitation and temperature on S. magellanicum and it follows that S. fuscum has a more stable growth and is more resistant to changes in climate, while S. magellanicum can take advantage of higher amounts of rain and grow faster when rain events occur (Mazziotta et al., 2019). This is in line with a paleoecological study showing that wetter climate caused a vegetation dominance shift from S. fuscum to S. magellanicum with subsequent reduction of carbon sequestration (Belyea & Malmer, 2004). The results may be affected by ecological differences between the newly described taxa within S. magellanicum. The two, S. medium and S. divinum, have been described to occur in contrasting habitats (Yousefi et al., 2017). However, we have observed both species in open as well as shaded habitats (Bengtsson, Granath, Cronberg, & Rydin, 2020). The ecological boundaries between these species need to be further studied, since it seems that both could be plastic enough to occur in either environment.
Our models for LI explained more variance than those for NPP.
LI varied less between the two years than NPP, which depended on stem BD. Hence NPP is also somewhat affected by the previous year's growth. Often, in multi-year studies on Sphagnum NPP, a constant BD is assumed, even though it may differ between years.
We collected data on BD for each year, which resulted in a higher between year NPP correlation (0.48-0.58) than often reported (e.g. Lindholm & Vasander, 1990). It is clear that annual NPP can be challenging to measure accurately in Sphagnum and perhaps a time-lag model is required to better capture the environmental effects (e.g. Backéus, 1988 included previous year's weather for predictions).

| Growth models
When Sphagnum moss growth is modelled, HWT is regularly used This is also the approach used in previous attempts to understand broad patterns of Sphagnum productivity (Table 1). Comparisons are difficult as most of these studies report only one main predictor that fitted the data best. If and how other variables were tested is often not mentioned in these papers. Because we sampled Sphagnum growth in the same two species from a broad geographical range, within narrow habitat definitions and with a common methodology, our study is better suited to identify the strongest climatic and environmental predictors for growth. Our dataset and analyses constitute an important step to inform current and future ecosystem models. Indeed, we are able to highlight the important difference between length growth and NPP, and provide direction for what variables to include in ecosystem models. In addition, our study enables comparison and parameterization of local drivers, such as HWT and vascular plant cover, which are two central factors in peatland models (Hayward & Clymo, 1983).
Our models for daily productivity showed the same significant predictors as for the entire growing season, except for the effect of PAR. Using daily growth we could investigate if effects were solely driven by conditions during the growing season, and the observed effect of PAR (summed over the season) supports the previous studies that have highlighted growing season length as an important factor for terrestrial ecosystems in general (Michaletz, Kerkhoff, & Enquist, 2018), and specifically for peatlands (Charman et al., 2013;Koebsch et al., 2020;Loisel et al., 2012). Still, our results suggest that climatic variables do affect growth directly and indirectly, and that growing season length may not be the main driver.
By measuring across the Holarctic we covered a wide range of environmental conditions. However, only 2 years growth limits our ability to study within-site variability and more consecutive years of measurements is needed to understand this variation, and to test if proposed mechanisms hold within a site or region.

| Implications for peatlands
The scale of our study and the consistent design provide a unique dataset to explore the drivers of growth of two common Sphagnum species across Holarctic peatlands. We conclude that variation in growth often is as large within as between peatlands, highlighting the importance of both global climatic factors (precipitation and mean temperature) and local factors (here vascular plant cover).
However, it was challenging to determine which climatic or other environmental variables influenced the growth of Sphagnum, since the most important variables covaried. Surprisingly, water table depth was not a good predictor of growth within a species.
Undoubtedly, the growth of Sphagnum and hence peatland function will be affected by environmental change. Sphagnum magellanicum generally has a wider niche in relation to water table level and nutrients, and showed a stronger direct response in LI to precipitation and temperature, which can give it a competitive advantage in warmer and wetter environments, than S. fuscum, which engineers an environment that is more resistant to desiccation . However, due to indirect effects and differences in bulk density the responses of the two species in productivity (NPP) were more similar, suggesting that the response of peatland moss NPP to climate change may be robust to species turnover (cf. Jassey & Signarbieux, 2019;Robroek et al., 2017).

ACK N OWLED G EM ENTS
We dedicate this work to the memory of Maria Noskova and Richard Payne, collaborators on this project who tragically passed away before the submission of this paper. The project was supported by

AUTH O R S ' CO NTR I B UTI O N S
G.G. and H.R. initiated and designed the study; All authors collected data, and F.B. and G.G. performed the analyses; F.B. wrote the first draft and all authors commented on the manuscript.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and R-script available from the Dryad Digital Repository https://doi.org/10.5061/dryad.1ns1r n8rm .