This experiment was carried out in northern Sweden on ten experimental sites in the vicinity of Arvidsjaur (65°35′–66°07′N, 17°15′–19°26′E). The sites are distributed along a natural fire-induced gradient of successional stages with time since the last fire disturbance ranging from 48 to 368 years (Zackrisson et al., 1996). Based on the time since last fire disturbance, the sites are classified as young (<100 years since the last fire; N = 3), intermediate (100–260 years since the last fire; N = 4), or old (>260 years since the last fire; N = 3) (Jackson et al., 2013; Supplementary Figure 1 and Supplementary Table 1). The most abundant tree species on all sites are B. pubescens, P. abies, and P. sylvestris (Jackson et al., 2013). The understory vegetation across the ten sites is characterized by ericaceous shrubs (E. hermaphroditum, V. myrtillus, and V. vitis-idaea) and the feather moss P. schreberi (Jackson et al., 2013). The relative abundance of the various over- and understory species is dependent on the time since the last fire: younger sites are characterized by shallow feather moss layer (mean: 29.2 mm), dominant biomass of V. myrtillus (mean: 102.3 g/m²) and high basal area of P. sylvestris (mean: 18.4 m²/ha), whereas older sites have a thick feather moss layer (mean: 51.4 mm), high biomass of E. hermaphroditum (mean: 90.5 g/m²) and dominant basal area of P. abies (mean: 8.9 m²/ha) (DeLuca T. et al., 2002; Jackson et al., 2013; De Long et al., 2015; Supplementary Table 1). The soils of all sites of the chronosequence are either Typic or Entic Haplocryods, which are acidic podzols typical of cold and wet environments (DeLuca T. et al., 2002). The mean annual air temperature across the sites is −2°C (average temperature in January: −14°C; average temperature in July: 12°C) and the mean annual precipitation is 600 mm (Jackson et al., 2013).

We established the experiment in 2010 with a full factorial combination of passive warming, shrub removal and moss removal at each of the ten sites. To apply the treatments, we established hexagonally shaped plots (diameter of 165 cm × 190 cm; total area 2.35 m²), which each contained a homogeneous vascular plant community composition and feather moss (Hylocomium splendens and P. schreberi) coverage specific for each site (De Long et al., 2015). Each of these plots was randomly assigned to one of four full factorial combinations of warming (ambient or increased temperature) and shrub removal (shrubs present or removed). These plots were then halved by a north-south middle line to create subplots randomly assigned to one level of the moss removal treatment (moss present or moss removed). This resulted in 8 different treatment combinations per site (Supplementary Figure 1) and 80 experimental subplots across the 10 sites.

To increase soil and air temperature of the plots assigned to warming we used transparent Perspex Open-Top Chambers (OTCs) (diameter: 165 cm × 190 cm, height: 47 cm, central exposed area: 0.95 m², MakroLife, Arla Plast AB, Sweden); a widely used method to investigate warming effects in tundra (Henry and Molau, 1997; Marion et al., 1997) and forested systems (Munier et al., 2010; Kaarlejärvi et al., 2012). To simulate summer warming during the growing season, we placed the OTCs over each plot assigned to warming each year in June and removed them again in October. The OTCs increased the air temperature by 0.4°C and the soil temperature by 0.2°C during summer (De Long et al., 2015) which is comparable to other studies in forest ecosystems (De Frenne et al., 2009). To avoid interference due to continuous shading of the forest canopy with the warming treatment the experimental plots were placed under canopy gaps at each site. In addition, we severed all roots around the outer edge of each plot with a shovel (to a depth of approximately 25 cm or the nearest rock) to exclude interference from tree roots.

For the shrub and moss removal treatments, we manually removed the aboveground parts of the respective functional group at the beginning of the growing season, but left the belowground parts intact to avoid disturbance to the litter layer and other plants (Wardle and Zackrisson, 2005; De Long et al., 2015). The functional group removal treatment was done once per year.

In the summer of 2013, when the disturbance impact of the initial vegetation cover removal was considered to have largely ceased (Wardle and Zackrisson, 2005), we planted two cohorts of tree seedlings in each subplot, one in June and one in August (hereafter referred to as “first” and “second” seedlings, respectively). Prior to transplanting, we grew all seedlings from commercial seeds for northern Swedish and Finnish provenances (Svenska Skogsplantor AB): B. pubescens (Patana I SV-421), P. sylvestris (ZON 6800 200), and P. abies (Björkebo FP-19). Before transplantation to the field in the beginning of June, the “first” seedlings were grown in trays with sand; B. pubescens seedlings in warmed (initially 20°C, later 15°C) and lit growth chambers for two and a half weeks, and P. sylvestris and P. abies seedlings in a greenhouse for one and a half weeks. The “second” seedlings were all grown in individual small pots in a 1:1 (vol.) sand:potting-soil mixture in the greenhouse for approximately 2 months before transplantation to the field at the end of August 2013. The evening before transplantation, we gently washed all seedlings free of sand and soil and kept them between moist tissues at 5°C overnight. At the time of transplantation, all first seedlings had only developed their first cotyledons and a limited root system (tap root with beginning lateral roots), while the second seedlings had developed several leaves and a viable root system.

In June 2013, we planted the first cohort: 50 B. pubescens seedlings, 20 P. sylvestris seedlings, and 10 P. abies seedlings per subplot. To account for heterogeneity in soil and vegetation, individuals were planted in one of two locations for each species within each subplot, resulting in two groups of 25 B. pubescens seedlings, two groups of ten P. sylvestris seedlings, and two groups of five P. abies seedlings per subplot. The number of transplanted seedlings varied between species, due to interspecific differences in germination and pre-transplantation survival, but the distance between transplanted seedlings was kept constant across species (approximately 1.5 cm). As there was very little within-species size variation in this cohort of pre-grown first seedlings, we randomly distributed all these seedlings across sites and subplots.

In late August 2013, we planted a second cohort of larger seedlings than in the first cohort to make sure that there would be enough many seedlings survive throughout the years for the final biomass analysis. This was done because the seedlings of the first cohort were very small and thus had a high risk of dying before the end of the study. The second cohort consisted of five B. pubescens seedlings and four of each P. sylvestris and P. abies seedlings per subplot. Each of these individual seedlings was planted at its own location within each subplot, to account for within-plot heterogeneity, and these locations were different than the locations at which the first cohort seedlings were planted. Further, to account for within-species size variation in the pre-grown “second” seedlings, all these seedlings were grouped into 3–4 size classes and seedlings within each size class were then randomly distributed across all sites and subplots. In addition, we oven-dried (48 h at 60°C) and weighed multiple seedlings of each cohort and species to obtain an estimate of the initial species- and cohort-specific total dry biomass (see Supplementary Table 2).

We planted all seedlings at least 30 cm from the outer edge of the hexagonal plots and at least 10 cm from the border between the subplots, as in Lett et al. (2020). The seedlings were planted at a depth that is comparable to that at which naturally germinated seedlings are positioned, i.e., in the feather moss layer with their stem base a few cm below the tips of the mosses, or in the organic soil with their stem base a few mm below the soil surface (Lett et al., 2017). After transplantation, all living seedlings from the “first” cohort were counted twice per year, in early (June/July) and late (August/September) summer, through 2013 until 2016, to assess their survival.

To obtain the biomass growth dataset in spite of a high seedling mortality rate, we combined biomass data from both cohorts of seedlings (see below). In August 2016, we thus harvested the above- and belowground biomass of the largest of the surviving seedlings of each species per subplot for both the first and second seedling cohorts, in all cases in which at least one seedling was still remaining. We used the largest seedlings because the size of the largest seedling more accurately reflects treatment-induced plot conditions on seedling biomass growth than the average biomass of all surviving seedlings per plot, which is biased by seedling mortality (Wardle et al., 2008). All harvested seedlings were washed in water to remove soil particles and any attached plant material. The seedlings were oven-dried for 48 h at 60°C and weighed. To combine the biomass data for the two cohorts, which differed due to age differences and planting date, we fitted linear regression models for each species describing the relation between first and second seedling biomass (Supplementary Figure 2). We then used these relations to convert the biomass values of second seedlings into that of first seedlings. For subplots with both second and first seedlings, we averaged the biomass of the first and converted second seedling to be consistent with the approach used for subplots with seedlings of only one cohort. In total, we thus harvested 276 seedlings for biomass analysis and obtained growth data for 94, 90, and 92 individuals of B. pubescens, P. sylvestris, and P. abies, respectively. We then calculated total relative biomass increase for each individual seedling (hereafter called “growth”) as: (individual seedling total dry mass at harvest − initial species-specific × cohort-specific total dry mass) / initial species-specific × cohort-specific total dry mass.

To assess the main and interactive effects of warming, moss removal, shrub removal and successional stage on seedling survival from 2013 to 2016, we applied a Cox mixed effects model for each tree species separately, using only the survival data from the first cohort of seedlings. Cox models evaluate the survival probability for each individual subject (here tree seedling) at a given time and allow for integration of random factors. Here we used subplot nested within plot nested within site as a random factor to account for the nested experimental design. We assessed the significance of the effects of the treatments and their interactions using Chi-square tests (Therneau, 2020).

To analyze the main and interactive effects of warming, moss removal, shrub removal, and successional stage on seedling growth, we applied linear mixed effect models separately for each specie, using the combined data from both seedling cohorts. We used subplot nested within plot nested within site as random factor to account for the nested experimental design. Due to 100% mortality in some of the treatment combinations (i.e., no seedling was available for biomass harvest) we excluded four-way interactions from the models. The significance of treatment effects was assessed with F-tests. We used model diagnostic plots to check the homogeneity of variances and the normality of errors (Crawley, 2007) and the data were log- or square-root-transformed whenever necessary (Tables 1, 2).

We used the “coxme” package (Therneau, 2020) to run Cox mix effect models, the “lme4” package (Bates et al., 2015) to run linear mixed effect models, the “tidyverse” and the “reshape2” package (Wickham, 2007; Wickham et al., 2019) to arrange data and the “ggplot2” package (Wickham, 2016) for plotting in R statistical software (R Core Team, 2020 version 4.0.2).

At the last survival counting occasion after 3 years 361 (9%) out of 4,000 planted B. pubescens, 357 (22.3%) out of 1,600 planted P. sylvestris and 105 (13.1%) out of 800 planted P. abies seedlings had survived. Throughout the duration of the experiment, warming had no significant effects on survival of any of the three tree seedling species either by itself or in combination with any other factor (Table 1), although there was a trend (p < 0.1) for a weak positive warming effect in the younger successional stages changing into a negative effect in the older stages for B. pubescens (Figure 1A, warming × successional stage interaction, Table 1).

Moss removal, shrub removal and post-fire successional stage all exhibited some significant species-specific main and interactive effects on seedling survival. As such, moss removal increased survival of all three species, with B. pubescens receiving the greatest benefit (Table 1 and Figures 1A–C). For this species, moss removal moreover increased survival more when shrubs were present compared to when shrubs were removed (moss removal × shrub removal interaction, Table 1) and this effect was more pronounced in the late post-fire successional stage (moss removal × shrub removal × successional stage interaction, Figure 1A and Table 1). For P. sylvestris seedlings, the positive effect of moss removal on survival was also highest in old successional stages (moss removal × successional stage interaction, Figure 1B and Table 1) but not dependent on shrub removal. The positive effect of moss removal on P. abies seedlings was independent of shrub removal and post fire successional stage (Figure 1C and Table 1).

Shrub removal reduced survival overall for B. pubescens seedlings, increased survival overall for P. sylvestris seedlings and had no effect on P. abies seedlings (Figures 1A–C and Table 1). The effect of shrub removal on B. pubescens survival changed from slightly positive to negative with increasing time since last fire (shrub removal × successional stage interaction, Figure 1A and Table 1), which was primarily due to a strong negative effect of shrub removal on B. pubescens seedling survival when combined with moss removal in the older successional stages (shrub removal × moss removal × successional stage interaction, Figure 1A and Table 1). There was no interactive effect of shrub removal with any other factor on seedling survival of either P. sylvestris or P. abies.

There was a significant main effect of post-fire successional stage on P. sylvestris seedling survival, as more seedlings survived in old compared to young and intermediate sites (Figure 1C and Table 1). There was no significant main effect of successional stage on seedling survival for P. abies, but a weak trend (p < 0.1) of higher survival in young compared to intermediate and old sites for B. pubescens (Figure 1A and Table 1).

Warming did not significantly affect total seedling growth for any of the three tree seedling species, although there was a positive trend (p < 0.1) for B. pubescens (Table 2). Interestingly, despite its strong effect on survival, moss removal did not have a significant effect on growth of the surviving seedlings in any of the species (Table 2).

Shrub removal increased growth for P. abies and P. sylvestris and there was a weak positive trend (p < 0.1) for B. pubescens as well (Figures 2A–C and Table 2). Furthermore, growth of P. sylvestris seedlings exhibited weaker growth in intermediate and old successional sites compared to young sites, which effect seemed (p < 0.1) primarily to occur in combination with moss removal and warming (Figure 2B and Table 2). However, growth of B. pubescens and P. abies seedlings was not affected by post-fire successional stage (Table 2). There were also no significant interactive effects among factors for any of the three species (Table 2).

In contrast to the prediction of our first hypothesis, experimental warming of air and soil temperature did not significantly increase seedling survival or growth. This is consistent with the results of Marty et al. (2020) who found no effect of experimental warming on seedling growth, but is in contrast to Nissinen et al. (2020) and Okano et al. (2021) who observed a positive effect on growth, and Wright et al. (2018) who found a negative effect on survival. These divergent responses can origin from different magnitudes of warming in different experiments, ranging from 0.4°C in our experiment, over 2°C in the experiment of Marty et al. (2020) to 3.4°C in the experiment of Wright et al. (2018). Boreal forests are expected to experience a 1.5–4°C increase of temperature (Gauthier et al., 2015; IPCC, 2021) suggesting that the magnitude of warming in our study was comparatively small and therefore may not have triggered a seedling response. It is also possible that other essential factors such as nutrient, water and light availability could act as limiting factor for seedling responses to warming (Luo et al., 2019). We observed weak trends of positive warming responses on B. pubescens survival and growth and P. sylvestris growth that were limited to the young sites. This may be related to higher nutrient availability, because an earlier study conducted in the same experiment measured higher ammonium amounts in warmed and young successional stage plots (De Long et al., 2015). Besides nutrient availability, light availability, soil moisture, and/or surrounding vegetation might be other factors more strongly affecting seedling performance than warming (Lett and Dorrepaal, 2018) especially when they are a limiting factor and not alleviated by warming.

We did not find support for the second part of our first hypothesis that warming responses in seedlings would depend on the presence of mosses and shrubs. Still, the main and some interactive effects of moss and/or shrub removals strongly impacted seedling performance. Interestingly, survival and growth were strongest affected by different understory functional groups: seedling survival of all three species responded positively to moss removal and seedling growth of P. sylvestris and P. abies responded positively to shrub removal. This may indicate different mechanisms for seedling survival and growth, as well as different modes of action for mosses and shrubs, such as facilitation or competition. Earlier findings have shown that the presence of mosses can be a great barrier for seedling germination and survival in the boreal forest (Valkonen and Maguire, 2005; Hypponen et al., 2013; Kokkonen et al., 2018) and relate that to changes in moisture (Lindo and Gonzalez, 2010; Soudzilovskaia et al., 2013) and competition for space as mosses can also suffocate tree seedlings (Stuiver et al., 2014). Shrubs have also been shown to suppress seedling growth (Wardle et al., 2008), e.g., E. hermaphroditum is known to produce poor quality litter with allelopathic compounds hampering seedling establishment (Nilsson and Zackrisson, 1992). This assumption is supported by positive effects of shrub removal on P. sylvestris root and shoot biomass and P. abies shoot biomass (Supplementary Table 3). Furthermore, interspecific competition, for nutrients and light, could be another mechanism, explaining why shrubs have a negative effect on seedling growth because established shrubs have a better developed root system and higher photosynthetic biomass than young tree seedlings (Kuuluvainen and Rouvinen, 2000). Seedling survival might thus be governed by competition for water and space while seedling growth might be more affected by competition for light and nutrients, overall highlighting the importance of different understory functional groups and their influences on resource availability as drivers of tree seedling performance.

We also found species-specific effects of shrub removal on seedling survival. B. pubescens seedlings generally survived better with shrubs, while P. sylvestris seedlings survived better without them. This difference in response between these two species could be due to a trade-off between avoiding competition and being eaten. Notably, B. pubescens seedlings are more palatable for small mammalian herbivores than conifer seedlings (Lyly et al., 2014) and potentially benefit from being surrounded by shrubs of low palatability as this will lower their accessibility to herbivores and thus reduce their mortality (Callaway et al., 2005; Barbosa et al., 2009; Jensen et al., 2012). This can also explain why the positive effect of moss removal on B. pubescens survival was more pronounced when shrubs were present than when they were absent. In addition, the presence of shrubs could have increased moisture in understory vegetation (Angell and Kielland, 2009), promoting slug herbivory particularly targeting P. sylvestris seedlings (Nystrand and Granstrom, 1997). Overall, this further implies the importance of understory functional groups for forest regeneration in terms of seedling survival and illustrates potential facilitating effects which are, however, species-specific.

We also found no support for our second hypothesis that the effect of warming on seedling survival and growth, and how this is governed by understory functional groups, is mediated by the post-fire successional stage. However, we found stronger positive effects of moss removal on seedling survival in late compared to early post-fire successional stages for B. pubescens and P. sylvestris. Old post-fire sites in our chronosequence have higher moss biomass than do intermediate and young sites (DeLuca T. et al., 2002; Jackson et al., 2013), therefore the removal of mosses likely relaxed the competition from mosses for space and resources more in the older sites (Nilsson and Wardle, 2005). Furthermore, B. pubescens survival also benefitted from the combined effect of moss removal and shrub presence in old post-fire successional stages. These older successional stages support a higher biomass of E. hermaphroditum relative to that of V. myrtillus (DeLuca T. et al., 2002; Jackson et al., 2013; De Long et al., 2015), and it could therefore be expected that seedling survival in these older forests would be impaired because of the ability of E. hermaphroditum to produce phytotoxic allelochemicals (Nilsson et al., 2000). However, as the effects of shrub presence on B. pubescens seedlings was positive and not negative, a different mechanism must have been involved. For example, a high abundance of E. hermaphroditum, which is a species of low nutritional quality, might have protected B. pubescens seedlings from consumption by small mammalian herbivores like hare, vole and/or lemmings (Callaway et al., 2005; Barbosa et al., 2009; Jensen et al., 2012; Lyly et al., 2014).

In general, our results show that understory functional groups are stronger drivers of tree seedling performance than direct effects of experimentally generated warmer temperatures ranging at the lower end of other experimental warming studies (Rustad et al., 2001). However, warmer temperatures at the level expected in boreal forests (1.5–4°C) and an increase in precipitation (Gauthier et al., 2015; IPCC, 2021) can also affect the cover and abundance of understory functional groups (Hedwall et al., 2015, 2019) and thus indirectly influence seedling responses. Overall, climate change is expected to decrease the cover of shrubs (Hedwall et al., 2015; Jonsson et al., 2021), while climate change effects on moss cover are less clear; potentially leading to an increase (Jonsson et al., 2021), but also a reduction (Alatalo et al., 2020). A decrease in shrub cover in future forests may in particular favor P. abies and P. sylvestris regeneration, because these species profited from shrub removal in terms of growth and P. sylvestris additionally benefitted from shrub removal in terms of survival.

In addition, climate change is also expected to promote fire events due to warmer temperatures (de Groot et al., 2013; Flannigan et al., 2013), resulting in post fire understory vegetation change following successional trajectories (Nilsson and Wardle, 2005). However, fire suppression (Pinto et al., 2020) and likely increasing precipitation in the boreal forests of northern Europe might also dampen the risk of forest fires (Krikken et al., 2021), leading to a higher number of old forest sites (>260 years since fire disturbance). In our study, time since last fire disturbance modified B. pubescens and P. sylvestris responses to moss and shrub removal such that performance was always worse at old sites. However, this effect was clearly weaker than the effect of functional group removal leading to the conclusion potential climate change effects on future forest regeneration may primarily act through changes in understory functional groups.