Disturbance-effects on treeline larch-stands in the lower Kolyma River area (NE Siberia)

Tree stands in the boreal treeline ecotone are, in addition to climate change, impacted by disturbances such as fire, water-related disturbances and logging. We aim to understand how these disturbances affect growth, age structure, and spatial patterns of larch stands in the north-eastern Siberian treeline ecotone (lower Kolyma River region), an insufficiently researched region. Stand structure of Larix cajanderi Mayr was studied at seven sites impacted by disturbances. Maximum tree age ranged from 44 to 300 years. Young to medium-aged stands had, independent of disturbance type, the highest stand densities with over 4000 larch trees per ha. These sites also had the highest growth rates for tree height and stem diameter. Overall lowest stand densities were found in a polygonal field at the northern end of the study area, with larches growing in distinct “tree islands”. At all sites, saplings are significantly clustered. Differences in fire severity led to contrasting stand structures with respect to tree, recruit, and overall stand densities. While a low severity fire resulted in low-density stands with high proportions of small and young larches, high severity fires resulted in high-density stands with high proportions of big trees. At waterdisturbed sites, stand structure varied between waterlogged and drained sites and latitude. These mixed effects of climate and disturbance make it difficult to predict future stand characteristics and the treeline position.


Introduction
Temperatures in the Northern Hemisphere have increased over the last decades, especially at high northern latitudes (IPCC 2014).While treelines are known to react to climate change by advancing in warmer or retreating in colder times (e.g.Payette and Lavoie 1994;Barnekow 1999;MacDonald et al. 2000), worldwide treelines show varying responses to recent climate warming, with half of them advancing, while two are retreating (Harsch et al. 2009).Treeline advance or densification of the treeline ecotone depends on successful recruitment within or above the treeline ecotone, which on the local scale, is influenced by spatial patterns.These include inter-and intraspecific competition and facilitation (Šrůtek et al. 2002;Lingua et al. 2008) or seed dispersal patterns (Nathan and Muller-Landau 2000), but also microsite conditions after fire (Bonnet et al. 2005) or created by water-logging (Holtmeier and Broll 2012).Various biotic and abiotic influences or past climatic changes may even amplify or mask the response of treelines to climate change (Holtmeier andBroll 2005, 2007).In particular, natural disturbances such as fire or flooding, but also human impact, can completely change the forest stand patterns at a local to regional scale (Bergeron and Archambault 1993;Dale et al. 2000).However, most research is conducted on undisturbed treelines (Cullen et al. 2001).
Fires are a common phenomenon in boreal forests and fire frequency is predicted to increase under climate warming (Kasischke and Turetsky 2006;Kharuk et al. 2008;Balzter et al. 2010); (but see Bergeron and Archambault 1993).Generally, fires have the ability to suppress treeline advance (Payette and Gagnon 1985;Shankman and Daly 1988), although the strength of this effect might depend on the taxa forming the treeline (Landhäusser and Wein 1993) and the actual rate of climate warming (Arseneault and Payette 1997).For example, post-fire recruitment after high-intensity fires depends on taxa-specific traits such as the ability to survive fires as trees and/or seeds or the potential for long-distance seed dispersal (Galipeau et al. 1997;Brown and Vellend 2014).It has been suggested that fire can even promote establishment beyond the current treeline by creating new seedbeds (Landhäusser and Wein 1993).Furthermore, highest recruitment rates in mixed forests of spruce and pine have been found within a few years after a fire (Charron and Greene 2002).Stand patterns initiated in that early phase may continue for many decades (Johnstone et al. (2004).In north-eastern Siberia, fires led to variable tree densities and biomass accumulations in mid-successional stands, ranging from 0 to 4 trees m -2 (Alexander et al. 2012).A simulation study of treelines in the Swiss Alps further suggests that the influence of fires as a shaping factor may almost be equal to that of climate (Schumacher and Bugmann 2006).
A peculiarity of the Siberian treeline eastwards of the Yamal Peninsula is that it lies over continuous permafrost (Abaimov 2010).Thus, besides fire, changing hydrological regimes such as those originating from permafrost degradation (Jorgenson and Osterkamp 2005;Ye et al. 2009) are of relevance for local tree stand development.While studies of North American and north-western Russian treelines suggest that permafrost degradation and subsequent drainage might favour tree establishment (e.g.Lloyd et al. 2003;Bhiry et al. 2007;Wilmking et al. 2012), comparable studies from Siberian treelines are lacking.Permafrost degradation and wetter soils, however, may lead to decreasing tree performance (Iijima et al. 2014).Furthermore, Walter et al. (2006) found a 14.7% increase in lake area around Chersky between 1974 and 2000 due to thermokarst processes, and therefore flooding from the draining of lakes or river changes might be of relevance for treeline transitions even at a broad scale.Crawford et al. (2003) concluded that in northern oceanic environments, ongoing climate warming is likely to result in paludification and consequently forest retreat, instead of advance.
Few studies have dealt with treeline dynamics in north-eastern Siberia (Harsch et al. 2009;compare Berner et al. 2013) and those available mainly focussed on growth trends or biomass accumulation patterns (e.g.Alexander et al. 2012;Berner et al. 2012Berner et al. , 2013)).Our analyses were conducted along a latitudinal gradient in the lower Kolyma river area, where Larix cajanderi Mayr is the only tree taxon.All studied sites were subjected to disturbances by fire or water at different intensities.We aim to investigate whether latitude (as a proxy for climate) or disturbances shape larch-stand patterns in this northern treeline ecotone.In this context we (i) applied uni-and bivariate point pattern techniques, to identify spatial patterns within and between the life stages of larches, (ii) identified changes in the age and size structure of larches and tested whether variances in stand structure values are best explained by latitude or by disturbance type and intensity and (iii) discuss how climate change might influence the disturbance regimes themselves.

Study area and field data collection
The study area is located in the lower Kolyma River area (68°23´N-69°08´N and 161°01´E-161°28´E, Fig. 1) and extends over an area of 1430 km².The climate is continental, with long, cold winters (mean January temperature -32.6 °C), and short, warm summers (mean July temperature 12.7 °C).Annual mean temperature is -10.6 °C and annual precipitation is low at 205 mm, half of it falling between May and September (Chersky meteorological station 1961-2015, station number 251230, http://www.ncdc.noaa.gov).As the meteorological record is relatively incomplete, only years with less than 10 missing days per month were used to calculate mean values.Mean annual temperatures are exceptionally warm since the new millennium, reaching -11.4 °C prior to 2000, and -9.4 °C from 2000-2015.The whole region is underlain by continuous permafrost, up to 300 m thick (Brown 1960).
Data on L. cajanderi stand structures were collected during three weeks in summer 2012.Sites were chosen to represent different disturbance regimes according to satellite images, correspondence with local inhabitants and site investigations.Altogether, seven sites (plots) of varying size were investigated in four areas.We recorded at least 30 trees (size-class definitions: seedling < 40 cm, sapling 40-200 cm, tree ≥ 200 cm) from a minimum area of 250 m².At all sites, height, basal diameter (right above the root collar), vitality (alive, dead), and, for individuals larger than 130 cm, diameter at breast height (DBH) were recorded for all individuals.Tree crown cover was estimated by measuring the longest crown axis, the length of its perpendicular axis and calculating the resultant area.Exact positions within each site of trees and saplings were mapped in the field.Seedlings were recorded in at least half of each site, but the positions were only recorded at a 4 m² subplot level.Upright dead trees and tree stumps were recorded including their exact positions and dead trees lying on the ground were recorded as dead wood ground cover.We collected tree discs at ground level or complete specimens of 126 individuals representatively across all sites and size classes (see definition above).Understorey vegetation was recorded and active layer depth, i.e. the depth of unfrozen ground above the permafrost table, measured at a 4 m² subplot level.

Site description
The Low Intensity Disturbance (LID) site (68°23´N, 161°27´E) was positioned on a dry hilltop, with a mean elevation of 30 m a.s.l. and active layer depth ranging from 16 to 80 cm.The soil was mainly covered by lichen (90%) and understorey vegetation consisted of e.g.mosses, Betula nana L., Ledum palustre L. and Poaceae.Larch trees and saplings were recorded in a 14 m × 20 m plot.Half of the trees on the plot were dead, and discs of two dead trees were dated to have died in the late 1950s.Although no fire scars were detected, the small crowns of the trees (see Supplementary file: Figure S2, available at https://doi.org/10.14214/sf.1666)hint to a low intensity crown fire at this site.
The Logging/Intense Disturbance 1 (Log/ID1, 68°31´N, 161°11´E) and Logging/ Intense Disturbance 2 (Log/ID2, 68°31´N, 161°12´E) sites were situated between two small channels of the Kolyma River, at an elevation of 22 m a.s.l.At Log/ID1, the active layer depth ranged from 30 to 73 cm; at Log/ID2 from 41 to 89 cm.The soil at both sites was densely covered by needles (>95%).Large Salix shrubs were growing in patches within the forest, dominating larchfree patches at Log/ID 2. The understorey vegetation included mosses and Pyrola spp.Larch trees and saplings were recorded from a 16 m × 20 m plot (Log/ID1) and a 16 m x 16 m plot (Log/ID2), but no age data were retrieved from the second site.Both sites lie in a region that was completely logged in the 1950s (local knowledge).Fire scars on cut edges of the tree stumps indicate a postlogging fire.No signs of fire were found on living trees or dead wood other than stumps.
The River Bank (RB) site (68°23´N, 161°28´E) was situated on a cut bank ~4 m above river level and in close proximity to a drained lake.Field and satellite image inspections led to the conclusion that water discharge into the river might have led to conditions which were dry enough to sustain tree growth, while the area surrounding the lake is too wet for larch.No signs of fire were discovered.The mean elevation is 17 m a.s.l. and the active layer depth ranged from 27 to 49 cm.The soil was almost completely covered by needles (>98%) and dead wood.Large Alnus shrubs were growing within the forest.Understorey vegetation included Pyrola spp., Poaceae and Ribes triste Pall.Larch trees and saplings were recorded from a 20 m × 20 m plot, but dendrochronological samples were only retrieved from adult trees.
The Meander1 (M1, 69°03´N, 161°13´E) and Meander2 (M2, 69°03´N, 161°12´E) sites are part of a wet, old river loop with polygonal structures.We found a few stumps, indicating moderate anthropogenic influence.At M1, the active layer depth ranged from 15 to 60 cm; mean elevation is 19 m a.s.l.At M2, active layer depth and mean elevation were not measured.The soil at both sites was mainly covered by mosses (>70%), and the understorey vegetation consisted of e.g.Cyperaceae, Alnus spp., Ledum palustre L. and Empetrum spp.Larch trees and saplings at M1 and M2 were recorded from 20 m × 20 m and 20 m × 24 m plots, respectively.Fire signs were detected on stumps only.
The Polygon (P) site (69°07´N, 161°01´E) was situated in a polygonal field with larches growing on polygon rims.As the tree density was sparse, no plot was demarcated, but we walked from tree to tree and searched for seedlings nearby, covering an area of 0.8 km².Twenty randomly placed 4 m² tree-free subplots were placed around the trees, to search extensively for seedlings.The active layer depth ranged from 13 to 57 cm and the soil was mainly covered by mosses (>76%).The understorey vegetation included Empetrum spp., Poaceae, Cyperaceae, Betula nana L. and Vaccinium vitis-idaea L. Increment cores were taken at breast height, but no basal rings were analysed.As no correction for missing piths has been made, only minimum ages can be inferred.

Dendrochronological approach
Discs were polished and ring width measured with a mechanical measuring table (Lintab, F. Rinn SA, Heidelberg, Germany) and WinDendro software.Samples were first cross-dated visually to detect missing rings and then checked with the help of Cofecha (Holmes 1983).Thin-sections of small individuals were dyed with astra blue and safranin, and rings were counted under a Zeiss light microscope and with WinDendro software.The ages of all trees were estimated with the help of age models, based on site-specific regressions (Suppl.file: Figure S1).

Statistical analyses
All analyses were conducted on larch only.Unless noted otherwise, only living individuals were used for statistical analyses.
Spatial point patterns, like the positions of trees within a forest patch, can be investigated using second-order statistics (Diggle 2003;Wiegand and Moloney 2013), for example Ripley's K-function (Ripley 1977) or the pair-correlation function (Stoyan and Stoyan 1994).We used univariate and bivariate point-pattern techniques to characterise spatial patterns within and between different life stages.The univariate Besag's L-function L 11 (r) (Besag 1977), a transformation of Ripley's K-function K11(r) (Ripley 1977), calculates the cumulative number of points within a distance r from a given point, divided by the overall intensity of the pattern and corrected for edge effects with isotropic edge corrections (Ripley 1977).Additionally, we used the pair-correlation function g 11 (r) (Stoyan and Stoyan 1994).This calculates the probability of finding trees at a given distance to another tree, using all inter-point distances, and is corrected for edge effects.
To determine spatial patterns of trees, saplings and dead trees/stumps and their relationship at different successional stages, we applied the bivariate L-function [L 12 (r)] and g-function [g 12 (r)] with isotropic edge correction and independence approach (Lotwick and Silverman 1982).At the LID site with 50% dead trees we used the random labelling approach (Besag and Diggle 1977), which does not shuffle tree positions, but tree labels (e.g."dead" and "alive") to gain information on the influence of the growing-position on a tree being dead or alive.
All spatial analyses were carried out up to ¼ of the shortest plot side, up to a maximum of 5 m.In the univariate analyses we tested against the null hypothesis of complete spatial randomness (CSR), which would result in L 11 (r) = r and g 11 (r) = 1, respectively.Values > r (or >1) indicate more inter-point distances at r than expected under CSR and thus clustering, values < r (or <1) indicate less inter-point distances at r than expected under CSR and thus segregation.For the bivariate cases, values > r (or >1) indicate positive spatial associations (attraction) of the two size classes, while values < r (or <1) indicate negative spatial associations (repulsion).
To create Monte Carlo simulation envelopes, we conducted 199 simulations (nsim) and used the fifth-highest and fifth-lowest simulated value for the confidence envelope (nrank), resulting in p = 0.05 (alpha = 2 × nrank/(1 + nsim); (Baddeley et al. 2015).The same number of points as found in the data was used for the simulated patterns.The spatial patterns were further assessed with a global goodness-of-fit test for maximum absolute deviation (Ripley 1977); mad.test in spatstatpackage (Baddeley and Turner 2005).At the RB site we found only 13 saplings 11 of which were dead, thus the sapling pattern at this site was excluded from spatial analyses.At site LID the positions of 1 sapling and 2 trees were not noted, nor, at site LOG/ID1, the positions of 6 trees and these individuals are not considered in the spatial analyses.
The biomass of living trees and saplings was calculated according to Alexander et al. ( 2012) with the function: total_biomass = 39.46 × basal_diameter 2.26 .
Linear mixed-effect modelling was used to investigate the influence of disturbing factor and disturbance intensity or transect position on stand, describing values such as tree size or density with site as random factor.For this purpose, the two areas where we have two replicated sites (Log/ID and M) were pooled into one site each to allow the analysis of individual tree data such as height and also stand data such as total biomass.As disturbance type and transect are highly correlated, we made two separate analyses and chose the better model with the help of Akaike's Information Criterion (AIC).The model with the smallest AIC and the lowest number of degrees of freedom was chosen as the best fitting model.
All statistical analyses were conducted in R (R Core Team 2015).For better comparability between different plot sizes, stand values were normalised to 1 ha.Growth rates of trees at the different sites were compared by age~diameter relationships.For spatial analyses we used the spatstat package (Baddeley et al. 2015); linear mixed-effect modelling was conducted with lme4 (Bates et al. 2015).

General stand characteristics and age structure
The highest mean tree height was found at the RB site, while the smallest trees were found in M1&2 and P (Table 1).In contrast, the mean DBH was large at LID, RB and P and small at M1. Large DBH values and small tree heights at P are probably a result of dead apical meristems found in many trees.
At all sites except P, seedlings occurred at high frequencies (Table 1).Highest densities of saplings and seedlings were found at the LID site where small individuals comprised 88% of all larch individuals.At the RB site, no recruits between 4 cm and ~150 cm height were found.Highest tree densities were found at the two Log/ID sites.In general, tree-to-recruit ratios were less than 1 on the older sites (LID, M1&2) and greater than 1 on younger sites (Log/ID1&2, RB), where the highest amount of total biomass was also found (Table 1).The lowest larch density and biomass amount were found at P.
While most stand variables can best be described by including disturbing factor and intensity as fixed factors, tree height is best explained by disturbing factor alone, with the smallest trees growing on water-affected sites.DBH and the number of stumps per ha are best explained by transect position (Table 2), with highest values at the southern and the northern sites.
Growth rates generally decrease with increasing stand age (Figs 2, 3) as trees at the younger sites (Log/ID1&2, RB) reached high diameter values much earlier than those at the older LID site.Comparing size and age distributions of the different sites suggests coherence for site history, with the four older sites (LID, M1&2, P) being similar, as well as the three younger sites (Fig. 3).Age, height and DBH distributions at LID have reversed J-shaped curves, as a result of much larger numbers of smaller and younger larch trees.In contrast, at M1&2, DBH and age distributions are skewed (Fig. 3), indicating increasing recruitment success until ~1980 with a sharp drop afterwards.Height distributions likewise progressively resemble a reverse-J shaped curve with low but increasing numbers of small individuals.Age and DBH distributions at the P site indicate sporadic recruitment before ~1900 but better reproductive success and recruitment afterwards.
Bimodal age distributions at the two Log/ID sites give rise to even-aged stands, where most trees established 30-50 years ago, but with a second peak of very young individuals (<10 years).While DBH distributions are fairly linear, height curves are bimodal with a pronounced peak of larches <50 cm tall and a broad, shallow peak of 5-10 m tall trees.Similar distribution patterns are seen in the RB data.

Univariate spatial patterns
According to the univariate pair-correlation function (g 11 (r)) and the L-function (L 11 (r)) (Table 3; Suppl.file: Figure S3), saplings at all sites show clumped patterns.Goodness-of-fit tests further indicate significant departures from random patterns.Inferred spatial patterns of trees vary between sites (Table 3; Suppl.file: Figure S4).Univariate g 11 (r) and L 11 (r) indicate clustering of trees at the young RB and Log/ID1&2 sites.No pattern is seen at the LID site and trees are randomly distributed at sites M1&2 at most scales analysed.Trees at the P site are significantly clumped.Dead trees and tree stumps show no deviations from CSR according to the goodness-of-fit (Table 3; Suppl.file: Figure S5).

Bivariate spatial patterns
While the g 12 (r) and L 12 (r) functions both indicate negative spatial associations between trees and saplings at the LID site, the overall pattern is only significant for L 12 (r) (Table 4; Suppl.file: Figure S6).Positive spatial associations between saplings and trees prevail at Log/ID1&2.While tree-sapling associations at M1 are random, positive associations are prevalent at M2.At site P, a positive association between saplings and trees is demonstrated at all scales.In contrast, overall patterns between saplings and dead trees/stumps (Table 4; Suppl.file: Figure S7), show no pattern according to the goodness-of-fit, despite L 12 (r) at Log/ID2 showing positive associations from one metre outwards.Positive associations are also indicated by g 12 (r) and L 12 (r) at LID. Associations of dead trees/stumps and trees are variable (Table 4; Suppl.file: Figure S8) and -despite a positive goodness-of-fit -random at site LID.At the Log/ID1&2 sites, they are clumped at/after ~1.5 and 2 m, respectively according to g 12 (r) and L 12 (r) (Suppl.file: Figure S8).

Spatial patterns vary between life stages and disturbance regimes
With our analyses, we identified uni-and bivariate spatial patterns of three larch life stages (saplings, trees and dead trees/stumps) under different types of disturbance regime.We were able to show that univariate patterns of trees differ between disturbance intensity, while sapling recruitment is clustered at small spatial scales across all sites and disturbance regimes.A similar result was noted in the Pyrenees, where recruitment patterns were similar between disturbed and undisturbed regions (Batllori et al. 2010).Furthermore, it has been shown that regeneration at alpine treelines often occurs in clusters (e.g.Mast and Veblen 1999;Zhang et al. 2009), while spatial patterns can be further influenced by disturbance type (e.g.fire intensity) and time since last disturbance (Kuuluvainen and Rouvinen 2000).All northern sites were located in wetlands, where larch growth is constrained by wet soil conditions.At P, the most extreme site, all size classes grow clumped at all scales.Clustering at this site can most probably be attributed to terrain, with larch recruitment limited to dry polygon rims (Zibulski et al. 2013;Frost and Epstein 2014).At such sites, seedlings are also in competition with a dense herbaceous layer (Bell et al. 2000).However, P is a special case, as underlying site conditions and the few larches at the larger spatial scales make it difficult to infer any second-order spatial patterns, and thus the results might be biased by random noise.At site M, the slightly clumped structure of saplings is most probably related to environmental heterogeneities in the wet environment.It is known that in unfavourable environments, competition may switch to facilitation, leading to tree islands -a pattern found at climatic and hydrologic treelines (Resler and Stine 2009 and references therein).Furthermore, once a tree has established, it usually improves its local environment, further enhancing tree recruitment (Bekker 2005).This might especially apply to saplings at sites M2 and P, which grow aggregated around trees.
The spatial analysis of the intensively disturbed sites Log/ID and RB show clumped patterns of trees.While saplings were virtually lacking at RB, they were positively associated with trees at Log/ID sites.In previous studies it was found that dead wood may serve as a seed bed (Johnson and Yeakley 2013).While we found a tendency of tree-clustering around stumps 10-20% of the seedlings growing at these two sites further germinated in dead wood remains (data not shown here).Bonnet et al.(2005), for example, showed that success of regeneration after fire was not only dependent on seed availability, which decreased with distance from unburned patches, but also on environmental conditions, including scorched litter on burned mineral soil and a high cryptogam cover, while exposed burned mineral soil and a dense understorey decreased recruitment success.While we found no fire marks at the RB site, large Alnus shrubs grew in treeless spots.Alnus is known to hamper larch recruitment (Goldammer and Furyaev 1996) and might be the underlying reason for tree clumping here.
At the comparably low disturbed site LID, recent recruits (saplings) grow clustered from 0-3 m and show a significant negative association to living trees.Although L. cajanderi is lightdemanding, seedlings can persist under closed forest canopies (Osawa et al. 2010) and can even grow to full tree height when a gap with sufficient light and growing space is available (Ban et al. 1998).However, the studied LID site did not have a dense tree cover, even when combining living and dead trees, so the light-gap theory does not satisfactorily explain our pattern.On the other hand, this result is in agreement with findings of (Kuuluvainen and Rouvinen 2000), who attributed this separation effect to interference zones around large trees.The LID site was positioned on a rather dry hilltop with dense understorey vegetation, so competition for water might have been the driving factor behind the development of this spatial pattern.

Stand structure is influenced by disturbance regime and intensity
Regional temperature is the most important factor in determining the natural treeline (MacDonald et al. 2008) and tree density decreases along the treeline ecotone (Šrůtek et al. 2002;Devi et al. 2008;Batllori et al. 2009).By analysing larch stand structure along a latitudinal transect, which could serve as proxy for temperature, with sites impacted by various disturbances, we found that these fire and water-related disturbances and their different intensities affect local stand characteristics to a greater degree than transect position.Furthermore, even at the northernmost site P, larches had ages over 160 years (at breast height), indicating that the position of the treeline is not currently shifting northwards.
All studied sites were impacted by disturbances and most stand variables best described by disturbance type and intensity.From local knowledge it is known that intensive logging in the 1950s led to large tree-free areas at Log/ID sites.Scorch marks on remaining logged stumps indicate that the sites must have been burnt after this logging event, but no other burned wood was found, thus we cannot make assumptions on the pre-fire stand structure or fire frequency.Berner et al. (2012), however, showed evidence for fires in the 1970s in the larger area.At the RB site, new land became available for forest colonisation ~60 years ago.Although the main processes behind this site development remain unknown, it seems likely that high water discharge from the cut bank into the river led to drier sites on the river bank, while areas around the lake remained too wet for Silva Fennica vol.51 no. 3 article id 1666 • Wieczorek et al. • Disturbance-effects on treeline larch-stands in… larch growth.The strong disturbances and influences at Log/ID and RB, which resulted in new land becoming available for colonisation, led to the highest tree densities of all studied sites.The capability of reaching high stand densities within a few decades after a disturbance is in agreement with observations by Alexander et al. (2012) and might be attributed to a lack of competition in the early stand-development phase (Goldammer and Furyaev 1996).Large establishment pulses are further known to follow stand-replacing fires under good conditions (Johnstone et al. 2004).At the Log/ID and RB sites, age curves show a sudden and high amount of recruitment, with a sharp drop afterwards.Diameter curves at these sites, however, rather resemble those of relatively unevenaged stands (Shorohova et al. 2009;Drössler et al. 2016), caused by highly variable growth rates of the single trees.The few recorded saplings at these sites might be attributed to the pronounced crown cover and resultant shading (Osawa et al. 2010).
Lowest densities were found at the water-disturbed sites M and P, which might be attributed water-logging and thus less favourable habitats for recruitment (Ohta et al. 2014) or to a lack of available seeds, as these sites had considerably fewer larch trees in the surrounding area than the more southerly stands (Hansen et al. 2013 and own observation).This finding is in concordance with (Kirdyanov et al. 2013), who found low growth rates in a wet riparian zone, in contrast to high tree growth rates on a river bank and terrace, within a short 100 m long transect.But while recruitment has been increasing in recent years at the polygonal site P, sites M1&2 experienced a drop in the last 20 years.Not only recruitment, but also growth was partly restricted, resulting in larches with an age of 38 years, but only 11 cm tall.Consequently, neither height nor diameter distributions can give any insight into the recruitment history of the different plots.The finding of decreasing recruitment success at sites M1&2 might be in concordance with satellite data, which indicate decreasing tree density in the region around Chersky, while shrub cover increases (Frost and Epstein 2014).Hansen et al. (2013) also estimated high losses of tree cover in our study area, as can be seen in their global map of forest cover change.From our observations in the field, however, it seems that forest cover from these satellite data is over-estimated for our study area, possibly being confounded by, for example, tall shrubs.Thus a loss of forest cover as determined from satellite images could be depicting places with intensive wildfires, as observed on the river bank opposite Chersky (own observation), but does not necessarily mean a loss of larch in general.
With respect to the intensity of disturbances, the LID site was probably the least disturbed.Here a fire event probably occurred in the late 1950s (which we infer from small tree crowns and two dated dead trees) that killed all recruits and half the adult trees.While 50% of the trees were dead, age and size distributions match well, resembling that of a natural, uneven-aged forest stand (Shorohova et al. 2009).Growth rates were much less variable than at the other sites, indicating common growing conditions for all individuals.

Climate change might influence disturbance regimes
Disturbances are the main driving factor in the landscape of the lower Kolyma in north-eastern Siberia, while the direct influence of climate seems to be of minor importance.However, climate change can impact the prevailing disturbance regimes.Wildfires in boreal forests, for example, have been found to increase over recent decades (Gillett et al. 2004) and were predicted to increase further in the future in parts of the circumboreal forest (Flannigan et al. 2009); this will have various consequences for forests at the treeline.Changing fire return intervals, for example, can lead to decreased seed availability (Brown and Johnstone 2012).While some studies found that fires led to decreasing recruitment success (Sirois and Payette 1991) or to a change in competition with increasing invasion of broad-leaved vegetation (Johnstone et al. 2010), others found that tree density after a fire led to high accumulations of biomass (Pare and Bergeron 1995).Ground fires following a stand-replacing fire, after which only the tallest trees survive, might even enhance the growth of larch trees (Schulze et al. 2012).A dense forest will then lead to a lower albedo and thus even higher temperatures (Otto et al. 2011).Although increasing tree density as a result of climate warming can lead to more intense fires, this seems not to apply to Russian wildfires, which might be because of the high incidence of ground fires (Wooster and Zhang 2004).As we found highest tree densities on those sites where new land was available for colonisation, it seems possible that increasing fire frequencies might be buffered by the strong recovery potential of these sites.If, however, only ground fires occur, tree density might strongly decrease, as can be inferred from low densities and slow regeneration processes at the LID site.
Higher temperatures might also lead to diverse reactions in periglacial wetlands: either to an increase in wet area due to permafrost thawing (Baltzer et al. 2014) or to a reduction of wetlands (Erwin 2008).A reduction in wetland can be manifest either as drained areas becoming too dry to sustain larch growth (Kharuk et al. 2015) or as an increase in favourable habitats for tree growth, which were previously too wet (e.g.site RB).Additionally, fires may lead to further permafrost degradation (Myers-Smith et al. 2008).We already saw a recruitment decline at the wet sites M1 and M2.As a result, the outcome of future stand patterns in the lower Kolyma area is highly speculative, due to the impact of climate change on the prevailing disturbance regimes.We have not seen a recent northward advance, and if recruitment and growth patterns stay as unfavourable as in recent years at site M, tree dieback and even treeline retreat are possible.

Conclusions
Our analyses of stand structure and spatial patterns highlight the complexity of the processes and factors influencing the formation of the treeline.We observed that tree regeneration is clustered at all sites, but not necessarily around living or dead trees.While high-intensity fires and drainage led to extremely high population densities and a subsequent decrease in recruitment intensity -probably due to competition -fires of lower intensity led to increasing but comparably lower recruitment success.Reproductive success at the northernmost sites in wetland areas might be influenced by climate change, changes in water disturbance regimes or a combination of both.These interactions of climate and the various disturbances make it difficult to predict future stand characteristics and the treeline position.

Fig. 1 .
Fig. 1.Locations of the study sites.A) Water disturbed sites in the northern study area, B) Sites influenced by logging and fire, C) Old, low disturbed site LID and water disturbed site in close proximity, D) Regional setting.Site codes: LID = Low Intensity Disturbance, RB = River Bank, Log/ID = Logging/Intense Disturbance 1&2, M = Meander 1&2, P = Polygon.

Fig. 2 .
Fig. 2. Age-basal diameter and height-age relationships of larches at differently disturbed sites.Note that age values for the P site (black triangles) are taken from breast height.Site codes: LID = Low Intensity Disturbance, RB = River Bank, Log/ID = Logging/Intense Disturbance 1&2, M = Meander 1&2, P = Polygon.

Fig. 3 .
Fig. 3. Stand values of living larches at the seven sites.Height values are given in 50 cm, basal diameter (DBB) in 1 cm and age values in 10 year bins.Grey bars for height and DBB indicate trees with height/DBB > site-specific median of trees.Grey bars for age data indicate measured individuals, while values for white bars are gained by regression analysis (see Suppl.file: Figure S1).Note that age values for the P site are taken from breast height.Note the different y-axes between sites.Site codes: LID = Low Intensity Disturbance, RB = River Bank, Log/ID = Logging/Intense Disturbance 1&2, M = Meander 1&2, P = Polygon.

Table 1 .
Larch stand parameters of the different sites.Stand values except for age, height, and diameter at breast height (DBH) have been scaled up to one hectare.Disturbing factor, intensity and transect position indicate the fixed effects used for mixed effect modelling (see Table 2).Site codes: LID = Low Intensity Disturbance, RB = River Bank, Log/ID = Logging/Intense Disturbance 1&2, M = Meander 1&2, P = Polygon.

Table 2 .
Linear mixed models for the different variables.Stand values except for age, height, and diameter at breast height (DBH) have been scaled up to one hectare.Disturbing factor, intensity and transect position indicate the fixed effects.All analyses were conducted on larches only.Site codes: LID =