Disturbances can control fine-scale pedodiversity in old-growth forests : is the soil evolution theory disturbed as well ?

Biota–soil interactions in natural ecosystems are the subject of considerable research. Our hypothesis is that individual trees play a significant role through biomechanical and biochemical disturbances affecting soil formation in temperate forests, resulting in a complex spatial pattern of disturbance regimes and a close relationship between disturbance histories and soil units. In Žofínský Prales (Czech Republic) – the fourth oldest, continuously protected reserve in Europe and the first site of global research network SIGEO (Smithsonian Institution Global Earth Observatories) in continental Europe – we compared extensive dendrochronological, soil and pit– mound microtopography data both temporally and spatially from an area of anthropogenically unaffected 42 ha collected from 2008–2012. These data sets differ in terms of information complexity and length of memory: tree cores contain complex information about the disturbance history of the past 350 years, footprints of disturbances from the uprooting of a specific tree can persist 1700 years, and soils represent an extensive composite phenotype that has been developing for at least the entire postglacial period (10 500 years). On average, 6.18–13.41 % of the canopy on individual soil units was disturbed per decade. Even though the “backbone” of key events in the development of the forest ecosystem remained the same (e.g. the 1870s, 1880s and 1980s), the internal structure of disturbance history often differed among soil units; the most exceptional were Gleysols and Histosols, where important feedback from soil to trees was expected. However, the characteristics of treethrow dynamics as well as the frequencies of stronger releases in core series also significantly differed along a gradient of terrestrial soil weathering and leaching (Haplic Cambisols – Dystric Cambisols – Entic Podzols – Albic Podzols). These results suggest the existence of several disturbance regimes within the forest, controlling fine-scale pedodiversity.


Soil evolution theory
Recent studies on the diversity and development of soils in the landscape have revealed the limits of traditional pedogenetic theories created from the middle of the 19th century, which emphasized the role of geology (Morton, 1843) or climate (Glinka, 1914) in pedogenesis.The traditional concept of "zonal soils" (Hartemink and Bockheim, 2013) or the recent concept of "endemic soils" (Guo et al., 2003;Bockheim, 2005) came from early views on landscape evolution (Johnson et al., 1990) and Clementsian concepts of vegetation succession (Cline, 1961).Terms such as "soil maturity" or "climax soil" have been used in scientific studies (e.g.Leirós et al., 2000).Although the idea of multifactorial soil evolution originated in the 1930s (Kellogg, 1936) -including the role of biotic factors -the role of individual organisms has been underestimated at coarser spatial scales.

P. Šamonil et al.: Disturbances can control fine-scale pedodiversity
In traditional theory, steady-state thickness is an important component of "developed zonal soils" (Lisetskii, 1999).
The assumption of a maximum weathering rate under a soil of mean (critical) thickness, with a decrease in both shallower and thicker soils (Gilbert, 1877;Davis, 1892), can even be considered the pillar of the theory.If the theory of steady-state thickness holds true, this implies that within an area of constant or minimally variable climate and geology, thickness should be closely related to topography, and soil taxonomical units should be formed in homogeneous large patches closely related to major topographical forms.Even now, "Gilbert's humped soil production function" has supporters (see Small et al., 1999;Anderson 2002) as well as critics, and its validity is much discussed (Gabet and Mudd, 2010;Phillips, 2010).Carson and Kirkby (1972) concluded that "Gilbert's steady-state soils" are actually metastable and the model does not operate at the landscape scale.Phillips (1993) andD'Odorico (2000) described significant variability in soils as a result of the unstable magnification of small perturbations, and Phillips et al. (2005) revealed considerable variability in the thickness of regolith in Arkansas, which was not connected to topography but was rather the result of local lithological variation and the effects of individual trees.Based on an erosion-sedimentation model, Gabet and Mudd (2010) concluded that treethrow pit-mound topography creates an uneven distribution of soils on slopes.It is not known to what extent the biomechanical and biochemical impacts of individual trees can influence not only the soil thickness but also the trajectory of soil formation and local pedodiversity (Huggett, 1998;Ibañez and Bockheim, 2013).
An evaluation of complex biota-soil interactions is difficult due to the need for extensive terrain data on coarser spatial scales and differences between timescales in vegetation and landform development.Whereas vegetation may vary continuously and rapidly, landforms sometimes vary episodically and relatively slowly (Phillips, 1995).In natural forest ecosystems, the dominant tree species can be considered to be "ecosystem engineers" (sensu Wilby, 2002), and determine the environment for other taxa through positive or negative niche constructions (Odling-Smee et al., 2003).Nonetheless, even the development of the populations of ecosystem engineers is determined -through feedback -by changes in the environment.The biotic component of ecosystems does not merely adjust to the environment but is an active participant in ecosystem creation.Modern pedogenetic theory needs to incorporate these effects.
As far as complex biota-soil interactions are concerned, so far only partial processes have been described.In the direction trees → soil environment, a surprisingly strong and spatially non-random impact of selected biomechanical and biochemical effects of trees on soil morphology and chemistry has been observed (e.g.Phillips and Marion, 2006;Šamonil et al., 2008a, b, 2010a, b); biogenous creep could even have at least the same importance as inorganic creep in the soil and regolith flux on slopes (see Caine, 1986;Pawlik et al., 2013).shows the border between the (probably) historically anthropogenically affected (32 ha) and unaffected (42 ha) zones; the grey solid lines represent contour lines; the black solid lines represent streams; grey areas show places significantly affected by water (semihydromorphic and hydromorphic soils).The study area was overlain with a rectangular network of 353 sample plots (grid 44.25 × 44.25 m); five regularly located soil profiles were studied in terms of soil properties within each plot (crosses in the detail); water-affected sites were evaluated separately from a rectangular network (for details see Šamonil et al., 2011).One rectangle of 16.04 ha served for the spatial analysis of releases in core series; eight rectangles of 1.50 ha served for the spatial analysis of pit-mounds.
In the other direction, several studies have revealed the significant and species-specific role of treethrow dynamics in the reproduction of forest tree species (Schaetzl et al., 1989a, b;Šebková et al., 2012).Research has often used mathematical models for tree-soil interaction evaluations (Yoo et al., 2007;Gabet and Mudd, 2010;Constantine et al., 2012), but as far as we know studies connecting the disturbance history with current soil diversity are still lacking for temperate forests.The absence of detailed observations on tree-soil interactions at the level of stands hinders a deeper understanding of the dynamics of natural forests.Our paper deals with this topic.
Cambisols, Dystric Cambisols, Entic Podzols and Albic Podzols -soils that, according to the traditional theory of zonal soils, should follow distinct zonality.Soil units differ significantly in terms of soil morphology as well as soil chemical properties (Šamonil et al., 2011, unpublished).Šamonil et al. (2011) also documented a surprisingly short range of spatial autocorrelation of the genetic depth of soils as well as a decrease of ranges of soil horizon thickness with increasing depth, which they connected to treethrow dynamics.The exceptional diversity and variability of soils at the site cannot be explained by traditional pedogenetic factors because the geological bedrock and climate are almost homogeneous, and geomorphology explains only a small part of the variability in data.These results therefore led to a hypothesis in which individual trees play a significant role in the development of soils, creating a complex pattern of disturbance history.

Disturbance regimes
Dendrochronological studies of disturbance history in natural forests often imply the uniformity of the local disturbance regime, thereby ignoring local variations in soil development such as specific tree-soil interactions.Typical publications include overviews of canopy disturbances in individual decades as far into the past as possible (e.g.D 'Amato and Orwig, 2008;Fraver et al., 2009).In European areas settled since prehistory, much recent research has dealt with the disturbance history of old-growth forests of a few to hundreds of hectares.Based on these results, the dynamics of temperate, boreal and Mediterranean forests has been discussed (e.g.Panayotov et al., 2011;Trotsiuk et al., 2012).However, the size of sites is almost always determined negatively (by the absence of human impact); therefore there is reason to doubt its ecological relevance in forest ecosystem dynamics.Some sites can be too small to study disturbance regimes; other sites can be too large and include several different regimes.Even though some studies observed significant spatial nonrandomness in the occurrence of disturbances (e.g.Panayotov et al., 2011), researchers usually ignore this issue.Spatial relationships are usually studied only in the search for synergy in the occurrence of individual disturbance events and in the evaluation of disturbance area and severity (e.g.Fraver et al., 2009).
Detailed analyses of spatial interdependence in disturbance impacts are also not available because of the limitations in data and methods.Disturbance history is still dominated by patch-model studies, in which the original disturbance record of individual trees is extrapolated for the entire study plot, which is of constant size and often internally heterogeneous.This method leads to general data that loses its variability.Only in exceptional cases has disturbance history been studied using spatially explicit, individual-based models at the level of individual trees (see Shimatani and Kubota, 2011;Šamonil et al., 2013a), which makes it possible to study biota-soil interactions in deeper detail.As a result, the possible co-existence of several disturbance regimes in natural forests that interact with soils remains an unresolved scientific issue.Researching this topic may lead to an understanding of key relationships between ecosystem components at the stand and landscape scales, including the creation of pedodiversity and its role in forest dynamics.Results can also help identify suitable sizes of sites and data sets for studies of disturbance history.A description of local patterns in disturbance history and their relationship to soils is one of the aims of the multidisciplinary study presented here.

Hypotheses and preconditions
In this study we examine the following hypotheses: (i) multiple disturbance regimes exist in a natural beech-dominated forest.At the same time, we hypothesize that (ii) specific disturbance records are connected to the local occurrence of soil taxonomic units (see below).To test these hypotheses, we use the extensive dendrochronological record of disturbance history (Šamonil et al., 2013a), the occurrence of treethrow disturbances (microtopographical features caused by a single uprooted tree, with or without the uprooted trunk, Šebková et al., 2012) and a taxonomic classification of soils (Šamonil et al., 2011).
We presume that two-way connections exist between the tree layer and the soil, through which they are mutually determined (sensu definition of geoecosphere by Huggett, 2001).This complicates the causal argumentation in the studied processes.Although we primarily discuss disturbance history as a factor of soil formation, disturbance phenomena and the preservation of their traces could be at the same time partially influenced by soil.The specific features of each data set are the length of its "memory", and the level of information complexity it includes (Fig. 3).Tree cores contain compact information about the past 350 years and the maximum length of a core series is 450 years; the duration of treethrow pit-mounds can be more than 1700 years (Sect.3.2); and soils at the site have been developing for at least the entire postglacial period (circa 10 500 years).
As far as information complexity is concerned, soil taxa represent an extended composite phenotype (Phillips, 2009) with often non-linear memory.They develop under the complex influence of many factors, of which one is disturbance.Dendrochronological records contain relatively solid information about the youngest phase of the disturbance history.They implicitly record a set of biomechanical and biochemical influences of trees on soils, but these influences cannot be separated.The longer memory of pit-mound microtopography also records rare disturbance events.The disadvantage of these otherwise key data (see Sect. 3.2) is that they represent only one of the biomechanical effects of trees.This is why we consider treethrow dynamics to be a qualitatively less complex record of disturbance history than dendrochronological data.
We test hypothesis (i) (see above) through a spatial analysis of the number of disturbance events on released trees and through a spatial analysis of treethrow pit-mound shapes; significant spatial non-randomness in the long-term disturbance record should suggest the existence of more disturbance regimes; H 0 supposes a random distribution.We test hypothesis (ii) through an analysis of various characteristics of the disturbance regime with respect to corresponding soil units; H 0 supposes that the characteristics of the disturbance regime are not dependent on the soil.

Study site
The research was carried out in the Novohradske mountains, in the Žofínský Prales National Nature Reserve (Fig. 1).This reserve is situated along an altitudinal gradient ranging between 735 and 830 m a.s.l.; gentle NW slopes predominate.Bedrock is almost homogeneous and consists of fineto medium-grained porphyritic and biotite granite (CGS, 2014).Annual average rainfall is 704 mm (CHMU, 2014).According to Tolasz et al. (2007), snow cover in the locality lasts on average 100-110 days per year.The average seasonal total depth of new snow is 150-200 cm, and the average seasonal maximum snow cover depth is 50-75 cm in the locality.The average seasonal maximum of water content of the snow cover is 100-150 mm.Annual average temperature is 4.3 • C. Plant communities can be classified in the Galio odorati-Fagetum (most frequent), Mercuriali perennis-Fagetum, Calamagrostio villosae-Fagetum, and Luzulo-Fagetum associations; spring-area plant communities can by classified in the alliance Caricion remotae or in the association Equiseto-Piceetum (Boublík et al., 2009).

Natural disturbance factors and stand structure patterns
The predominantly fine-scale mosaic of forest patches contrasts sharply with the few hectares of blowdowns created by a storm on 18 January 2007 (named Kyrill).The uneven temporal occurrences of treethrow pit-mounds (Šebková et al., 2012) as well as the multimodal record of disturbance history (Šamonil et al., 2013a) suggest that such strong events are unlikely to be unique in the history of Žofín.At the same time, Šamonil et al. (2013a) found an unclear relation between disturbance intensity and the corresponding spatial autocorrelation range of disturbance events, which suggests different disturbance agents in the past.In addition to strong winds and snow storms, Chadt (1908), Kruml (1960) and Brázdil et al. (2004) also found biotic disturbance factors (Ips typographus L.) in the region in the past.Strong wind events originate mainly in connection with convective storms that generally occur during the summer half-year (April to September) or in the form of gales along conspicuous horizontal pressure gradients that usually occur during the winter half-year (October to March, Dobrovolný and Brázdil, 2003); tornadoes as well as tropical cyclones are generally absent in Central Europe (but see Setvák et al., 2003).Tectonic, volcanic and flood disturbances are not relevant in the locality.
In the core zone, Fagus sylvatica reached a maximum age of 430, Abies alba 420 and Picea abies 439 years (ring counting after cross-dating by Šamonil et al., 2013a).Repeated tree censuses over a period of almost 40 years have shown that approximately two-thirds of all dead trees perished standing and subsequently fell with minimal soil physical disturbance.The remaining one-third was uprooted.Pit-mound microtopography can persist for more than 1700 years after a tree uprooting event (Šamonil et al., 2013b).Treethrow pit-mounds currently cover 11.65 % of the reserve (74.2 ha), with 7.7 % consisting of mound areas and 4.0 % of pit areas (Šamonil et al., 2011).The overall rotation period (how often an area equivalent to the entire study area is disturbed, sensu Pickett and White, 1985) of tree uprooting dynamics in Žofín was calculated to be approximately 1380 years.However, if we take into account the postdisturbance spread of mound materials by water soil erosion onto adjacent, originally undisturbed sites, the rotation period for disturbance is reduced to approximately 870 years (Šamonil et al., 2013b).This is an average number for the entire site, which does not take into account the local nonrandomness in the occurrence of disturbance events.

Tree census and dendrochronological survey
Individual data sets were completed based on a regular network of square plots with a lateral of 44.25 m established in 2008.The centres of squares were located geodetically with an accuracy of circa 0.05 m within the 74.2 ha area (Fig. 1).The basis for the whole-scale dendrochronological investigation was a tree census that began in 2008 and was completed in 2013.The locations and selected properties (tree species, diameter at breast height -DBH, health condition -dead standing or lying tree, broken tree, stump, etc.) of all trees with DBH ≥ 10 cm were recorded within the 74.2 ha area, and a stem position map including 23 881 standing or lying tree trunks was created (for details, see Král et al., 2010a).Precise stem positions, which were measured from the centres of the regular network with a precision of circa 0.5 m, were connected to soil maps and served as the basis for spatial analyses of releases in core series (see Sect. 3.4).
We applied a widely accepted approach to tree core sampling (e.g.Lorimer andFrelich, 1989, 2002).Tree cores were taken at a height of 1.3 m (recruitment age) on each of the 353 square plots, one from each cored tree.Six cores per plot (nine in cases of the occurrence of gaps) were extracted from non-suppressed trees closest to the plot centre.Apart from living trees, dead trees disturbed during the Kyrill storm were also cored.A historical stem-position map from 1997 enabled us to eliminate recent, dendrochronologically irrelevant disturbances from the sampling (Nowacki and Abrams, 1997); thus we cored only trees released before 1997.Šamonil et al. (2013a) examined the more recent disturbance history with the help of a tree census.Cores visually accepted in the field were evaluated in the laboratory.Eliminated samples were replaced by coring a neighbouring tree on the plot.Three rounds of field sampling and laboratory analyses were carried out from 2008-2011, with 3020 individual trees cored from the 74.2 ha area.
To fully evaluate the disturbance history, we cored an additional 212 juvenile trees with DBH ≤ 10 cm in gaps or suppressed state in order to calculate gap origin thresholds.We measured 636 crown areas to calculate the DBH-crown-area relationship, and we evaluated the social status of 2293 randomly selected trees (e.g.Lorimer and Frelich, 1989;Frelich, 2002;Splechtna et al., 2005).A detailed description of these data sets and methodological steps is given in Šamonil et al. (2013a).
The cores were dried and smoothed with fine sandpaper and processed using standard dendrochronological techniques (Schweingruber et al., 1990).We measured the widths of the tree-rings using the PAST 4 program (SCIEM, 2007), with 0.01 mm accuracy.Damaged cores, cores without subbark tree-rings and cores that were missing more than 30 mm of the pith were rejected.A total of 1986 cores were accepted for further disturbance history analyses.Cores were crossdated and we used a pith locator to estimate the number of rings to the pith (see details in Šamonil et al., 2013a).

Soil survey
In 2008, 1765 soil cores (for the evaluation of deeper soil horizons and properties) as well as shallow soil profiles (for the evaluation of upper horizons) were taken on 353 square plots within the 74.2 ha study area.We took five profiles per plot using a strict sampling design (see details in Šamonil et al., 2011).Each profile was classified based on the World Reference Base (Driessen et al., 2001, IUSS Working Group WRB 2006).Because of the sharp transition between terrestrial and (semi)-hydromorphic soils, we separately mapped soils strongly influenced by water (Fluvisols, Stagnosols, Gleysols, Histosols).The stem position map was used for precise mapping (see above).We then used these results to determine locally predominant soil taxonomic units (Fig. 2 Places with the occurrence of three or more soil taxonomic units were considered to be mixed-dry (terrestrial) or mixedwet ((semi)-hydromorphic) sites.

Survey of pit-mound microtopography
Pit-mounds were field mapped in 2008.All features whose height or depth reached 20 cm (for pit-mounds without a fallen trunk), or when the fallen trunk had a DBH ≥ 10 cm were included in the analyses.Moreover, additional characteristics of pit-mounds on 353 circular plots of 23 m in diameter, which were established around the centres of a regular network, were also recorded (Fig. 1, Šebková et al., 2012).
Overall, detailed characteristics of 1733 pit-mounds were collected in 2009 and 2010.
The following properties were recorded for each pitmound form within the circles: length and width (to calculate the area of the pit-mound, an approximation to an ellipse was used), the height of the mound and the depth of the pit relative to the adjacent undisturbed ground surface, an estimation of the proportion of pit and mound areas, the presence of a fallen trunk (see details in Šamonil et al., 2009).The spatial pattern of treethrow pit-mound positions was studied (see below), while the detailed characteristics of treethrows were used to create descriptive statistics for the individual taxonomic units.

Data processing
To eliminate the historical impact of humans in the study of tree-soil interactions, we focused on the core zone of the reserve.Because of indications of possible human influences, the 50.0 ha core zone established by Šamonil et al. (2013a) was reduced for this study to 42.0 ha (Fig. 1).As a result of their small spatial extent, Fluvisols, Stagnosols, Stagnoglej (a specific unit distinguished in Czech soil taxonomy by Němeček et al. (2011), but not in WRB) and mixed water-affected soils were excluded from the analyses.In the end, the study plot was reduced to 40.4 ha, where we distinguished locally dominating Haplic Cambisols, Dystric Cambisols, Entic Podzols, Albic Podzols, Gleysols + Histosols or mixed terrestrial soils.These units covered between 1.36 (Albic Podzols) and 13.79 ha (Haplic Cambisols, Table 1, Šamonil et al., 2011).Corresponding data on disturbances were associated with these units.The amount of data within all units was sufficient for statistical analysis (see below).

The dendrochronological record and its relationship to soils
Disturbance history was studied on an irregular network of points represented by precisely located cored trees.We separately evaluated (1) the initial growth of trees -whether it occurred in a gap or under the canopy; and (2) disturbance events during subsequent growth on the basis of detected releases in the radial tree growth.
The threshold ring width of gap origin was calculated to be 1.45 mm for F. sylvatica and other broadleaved species, and 1.58 mm for P. abies and other conifers (Šamonil et al., 2013a).The release threshold was detected using the boundary line approach (BL, Black and Abrams, 2003).According to these authors boundary line represents the maximum percent growth change that is physiologically possible at a given level of prior growth.We considered growth changes (GC) at the level ≥ 100 % of BL as major releases, GCs in the range 50-99 % of BL as moderate releases, and GCs in the range 20-49 % of BL as weak releases.Pulses that were < 20 % of BL or pulses with duration under 10 years were considered to be inconclusive.To prevent overestimating the proportion of canopy area disturbance by lateral releases, the calculation of disturbance history was limited only to trees that did not reach the tree canopy.The threshold was calculated by Šamonil et al. (2013a) according to Lorimer and Frelich (1989) and D'Amato and Orwig ( 2008) at DBH = 75 cm for broadleaved species and DBH = 60 cm for conifers.Thicker trees could be overtopped with ≤ 5 % probability, hence they were rejected from canopy disturbance area calculations.
We evaluated the decadal disturbance history.As a measure of disturbance intensity, we used the proportion of disturbed tree canopy (e.g.Fraver et al. 2009).For that purpose, we modelled the relation between the area of exposed crown and DBH for particular woody species (Šamonil et al. 2013a).According to this relation and the calculation of cumulative radial increment of DBH the canopy areas were estimated for every year a cored tree was in the records.Finally, we calculated the decadal relation between the disturbed area and the summary canopy area (i.e.sample depth) separately according to the type of release in core series (weak, moderate, and major; for details see Šamonil et al., 2013a).
On the basis of their precise position, cored trees were attributed to particular soil taxonomic units.For the data sets that emerged from this process, we calculated partial disturbance histories (Fig. 4) and basic descriptive statistics (Table 1).Chronologies were truncated when sample depth dropped below 15 trees (Fraver and White, 2005;D'Amato and Orwig, 2008).The shortest chronologies started in 1790 (Albic Podzols, Gleysols + Histosols); the longest in 1660 (Haplic Cambisols, Entic Podzols).
Individual time series of canopy area disturbance and the number of juvenile trees were compared between different soil types.Time series (beginning in 1790) were filtered using a central moving average with a window width of 10 years.To get a clearer trend, the moving average was applied twice.The resulting trend curves that correspond to individual soil units were compared by computing distances between soils based on the slopes of curves in corresponding years.More specifically, for each soil unit a series of differences (representing curve trend slopes) between values in successive years was calculated, and a Euclidean distance matrix based on these differences was then computed and used in cluster analysis (complete linkage method, Fig. 5).This variable is dependent on the number of cored trees within soil unit and relation between DBH and crown area (Table 1, Šamonil et al., 2013a).The canopy areas were estimated for every year a cored tree was in the record, crown areas gradually decreased to the past.Chronologies were truncated when the sample depth (refers to the second y-axis) dropped below 15 trees (for details, see Sect.3.4.1.).
The statistical significance of each distance (similarity) was assessed via Monte Carlo permutation tests with 10 000 permutations (Table 2).At the same time, we evaluated the spatial relationships between trees with differing numbers of releases.Up to the moment of reaching canopy closure, 0-9 releases ≥ 20 % of BL were detected in tree-ring series (Šamonil et al., 2013a).Marked differences in the data refer to the possible existence of various disturbance dynamics within the locality.On the one hand, pulsating gaps with many releases existed; on the other hand some places experienced closed or open tree canopy for prolonged periods.We chose only those trees that had already reached the canopy, i.e. broadleaved trees with DBH ≥ 75 cm and conifers with DBH ≥ 60 cm (Fig. 6, n = 682).Releases were monitored up to the point of reaching of the canopy.Spatial analyses were carried out on a rectangle of 16.04 ha (Fig. 1).
To describe spatial patterns of core series as well as uprootings (see below) we used various types of pair correlation function (Stoyan and Stoyan, 1994).Stoyan and Penttinen (2000) define the pair correlation function as follows: consider two infinitesimally small discs of areas dx and dy at distance r.Let p(r) denote the probability that each disc contains a point of the process.Then p(r) = λ 2 g(r)dxdy Eq. ( 1), where λ is density.In a different way, the pair correlation function g(r) is the probability of observing a pair of points separated by a distance r, divided by the corresponding probability for a Poisson process (Baddeley, 2008).
To describe the tree spatial patterns of core series, we used the "i-to-any" summary g i• (r), which is the corresponding analogue of the pair correlation function (Baddeley, 2008, see below).A randomization test was used to test the null hypothesis of random labelling: given the locations X, the marks are conditionally independent and identically distributed.In randomization tests, the simulated patterns X are generated from the data set by holding the point locations fixed and permutating the marks randomly.Under random labelling, g i• (r) = g(r) Eq. ( 2), we thus used g i• (r) − g(r) to make a test statistic for random labelling.We generated 199 simulations of our null model to obtain pointwise critical envelopes for this model.Edge effects were corrected using Ripley's isotropic correction (Ripley, 1988).

Pit-mound microtopography and its relation to soils
The position of each pit-mound (Fig. 6) was attributed to particular soil units.We also placed eight rectangles of 1.5 ha in terrestrial stands to study spatial relationships among treethrows.To describe the uprooting density variability at sites dominated by various soils, the pair correlation function g(r) was estimated for each plot at steps of 0.5 m for r values up to 25 m.A null model of complete spatial randomness (hereafter CSR) was used on the assumption that the first-order intensity λ is constant within particular plots.We used 199 Monte Carlo simulations of CSR to obtain pointwise critical envelopes for a test of our null hypothesis (Besag and Diggle, 1977).Edge effects were corrected using Ripley's isotropic correction (Ripley, 1988).All spatial analyses were conducted using the package "spatstat" (Baddeley and Turner, 2005) in the statistics software R (R Development Core Team, 2006).

Dendrochronological disturbance history (circa 350 years of memory)
Partial disturbance histories provided a complex perspective on the most recent disturbance history of the site.(Fig. 4).On average, 6.18-13.41% of the canopy on individual soil units was disturbed.Even though mean values of decadal canopy disturbance area decreased along a gradient of weathering and leaching processes (Haplic Cambisols → Dystric Cambisols → Entic Podzols → Albic Podzols), because of the rather wide confidence intervals the differences between soils were not statistically significant (Table 1).Mean values also did not significantly differ for either all releases ≥ 20 % BL or stronger releases ≥ 50 % BL.The analysis of temporal series (Table 2, Figs.4-5) showed some specific details of partial disturbance histories.Even though the "backbone" of key events in the development of the forest ecosystem remained the same (e.g. the 1870s, 1880s and 1980s), the internal structure of the disturbance history often differed among soil units.Those series including only stronger releases (≥ 50 % BL) demonstrated a higher degree of dissimilarity.From this perspective, all soils were mutually dissimilar; the most exceptional were Gleysols and Histosols .
Including weak releases (20-49 % BL) in the analysis of temporal series from 1790 led to an increase in the statistical similarity of partial soil disturbance histories.Significantly similar were mainly Haplic Cambisols, Dystric Cambisols and Entic Podzols (Table 2,.We also found similarities between Dystric Cambisols and Albic Podzols, and also between mixed terrestrial soils and Haplic Cambisols and Entic Podzols.Gleysols + Histosols again showed a unique disturbance history. As opposed to data on canopy disturbance area, juvenile individuals of gap origin recorded the presence of gaps rather than their origin.In particular decades, they reflected the appearance of new gaps as well as the long-term presence of open gaps.The presence of juvenile individuals of gap origin was again significantly similar between Haplic Cambisols, Dystric Cambisols, Entic Podzols and Albic Podzols, and we also observed significant similarity with p ≤ 0.05 towards Gleysols+Histosols (Table 2).By contrast, no similarities between soils were observed when the regeneration of all juvenile individuals (growing in gaps and under closed canopy) was analysed.This indicates the different age structure of the forest on various soils while preserving partial similarities of the disturbance record.
We analysed the number of releases ≥ 20 % BL in individual core series (Table 1).Although the fewest releases were experienced by trees on Gleysosls + Histosols and mixed terrestrial soils, the differences exceeded the level of statistical significance α = 0.05 by a narrow margin.On average, 2.09-2.80releases tree −1 were recorded.By contrast, the number of releases ≥ 50 % BL statistically significantly increased along a soil leaching and weathering gradient from Haplic Cambisols (0.09 release tree −1 , after Box-Cox transformation) to Albic Podzols (0.50 release tree −1 , Table 1).A significant gradual increase in releases ≥ 50 % BL was observed particularly in the turbulent period between 1870 and 1889 along a gradient of soil leaching intensity from Haplic Cambisols to Albic Podzols (Fig. 4).
The analysis of the spatial pattern of trees with varying numbers of releases ≥ 20 % BL in tree-ring series brought surprising results (Fig. 7).Individuals with ≥ 5 releases tree −1 had a marked tendency towards clustering at a distance of 5-30 m, while at ca. 20 m their occurrence was even statistically significantly clustered.On the other hand, trees with fewer releases at these distances were rather regularly distributed.

Tree uprooting dynamics (circa 1700 years of memory)
Treethrow dynamics showed clear connections to soil taxonomy.The average number of uprooted trees increased along a gradient of soil weathering and leaching, from almost 9 trees ha −1 on Haplic Cambisols to more than 21 trees on Entic and Albic Podzols (Table 1).In the case of older treethrow pit-mounds without the fallen trunk, this trend was observed in the direction Haplic Cambisols (32.7 pit-mounds ha −1 ) → Dystric Cambisols (53.6 pitmounds ha −1 ) → Entic Podzols (61.9 pit-mounds ha −1 ), while on Albic Podzols the number of pit-mounds decreased (48.4 pit-mounds ha −1 ).Gleysols + Histosols formed a unique unit in which huge numbers of uprooted trees were recorded (70.4 pit-mounds ha −1 ).In many cases these were freshly fallen trunks downed by the Kyrill storm on 18 January 2007.Because of the much shorter longevity of mostly organic treethrow features, the number of older treethrows without the trunk on Gleysols + Histosols was extremely low (10.9 pit-mounds ha −1 ).The highest number of treethrow pit-mounds among all soil units was observed in terrestrial stands with high soil diversity (86.1 pit-mounds ha −1 , Table 1).The average size of treethrow microtopographical shapes and the spatial proportion of mounds within microsites significantly differed according to soil units.Once again Gleysols + Histosols turned out to be exceptional, where the smallest average size of treethrow shapes (7.59 m 2 ) and the lowest average proportion of mounds (35.28 %) were observed.This was undoubtedly connected to the specific "plate-like" shape of root-balls in waterlogged places, where treethrow pit-mounds have different dynamics.In terrestrial stands the average size of pit-mounds significantly decreased along a gradient of soil leaching from 11.13 m 2 on Haplic Cambisols to 8.60 m 2 on Albic Podzols.The proportion of microsites -mound vs. pit -did not significantly differ in terrestrial areas (Table 1).
The analysis of the spatial pattern of treethrow pit-mounds on eight plots of 1.5 ha showed a general tendency towards clustering between ca. 3 and 17 m.This tendency markedly grew along a gradient of soil leaching from Haplic Cambisols to Albic and Entic Podzols (Fig. 8); on Podzols, the pattern of pit-mounds was already statistically clustered.

Discussion
Disturbance events in the natural temperate forest studied here occurred non-randomly in space and time, and there were significant differences between disturbance histories of various soil units.While periods of intensive disturbances were clearly visible in all partial dendrochronological records, weaker events were site specific.This indicates that the traditional sum of disturbance history (e.g.D 'Amato and Orwig, 2008;Fraver et al., 2009) can mask significant differences in the pattern of weaker events and encumber the ecological interpretation of results.At the same time, weak disturbances are known to be of key importance in the dy-namics of beech-dominated temperate forests (Splechtna et al., 2005;Šamonil et al., 2013a).
The observed clustered occurrence of trees with higher disturbance events in the dendrochronological record at a distance of 5-30 m may be connected to the existence of a selection phase of the forest cycle (Král et al., 2010b).This phase is characterized by significant vertical layering, and, according to some authors, can be stable for decades or even hundreds of years (see Rubin et al., 2006;Westphal et al., 2006;Král et al., 2010b).Dying individuals in the canopy layer gradually release the lower tree layer, which is indicated in radial growth as a (weak) release (see Šamonil et al., 2013a).In this case, the long-term dynamics of lower layers is determined by higher layers rather than by the relationships among individuals in the lower tree layers.Such relationships are crucial for individuals with a low number of releases in their tree-ring series.Such trees can be found more often in lighter edaphic gaps or regenerated after intensive disturbances that remove several canopy trees.Competition among such individuals could lead to their regular spatial distribution.Temporal changes in the spatial pattern of forest cycle phases are a promising subject for future research.
The dendrochronological record (particularly releases ≥ 50 % BL) was generally in accordance with the occurrence of treethrow pit-mounds.Lucas et al., 1993) could locally intensify or suppress podzolization processes.Weaker manifestation of podzolization process can be seen as the typical course of pedogenesis at terrestrial sites corresponding to other pedogenetic factors (see the assessment of soil taxonomic adjacency by Phillips, 2013b).Gradual changes of some disturbance regime properties (e.g.treethrow dynamics) along a gradient of soil weathering and leaching support this assumption.However, our results from the forest stand scale do not prove a causal relationship between disturbance histories and soil units.Our findings should be subsequently verified by a detailed geophysical measurement stratified according to the soil units (Neal, 2004, verification of biomechanical effects) and by pedoanthracological analyses (Robin et al., 2013, verification of biochemical effects) accompanied by radiocarbon dating (Walker, 2005).The co-formation of soil units by spatially non-random disturbances could lead to an increase of taxonomic pedodiversity and self-organization over time, which was empirically demonstrated by Phillips (2001) and Saldaña and Ibañez (2004) in chronosequences of fluvial terraces.In a rain forest, Lobo and Dalling (2013) revealed a significant relation between soil type and canopy disturbance regime.The importance of past disturbances in pedodiversity can be quantified by a calculation of soil complexity (Phillips, 2013a).On the theoretical level, our findings support rather an evolutionary view of pedogenesis (with a variety of pathways or multidirectional soil succession) than the developmental view (unidirectional pedogenic changes, Huggett, 1998).
The biomechanical effect of trees on soils has repeatedly been described on the fine scale of disturbed microsites (reviews by Šamonil et al., 2010a;Pawlik, 2013).Studies have revealed changes in the rate as well as the trajectory of pedogenesis (e.g.Veneman et al., 1984).Tree uprooting, the crucial biomechanic process in studied forest ecosystems, re-sults in an overturning of the soil profile as well as the emergence of the C substrate horizon (see Šamonil et al., 2011, 2013b).The coarse-grained material of the mound represents specific material (e.g.Schaetzl et al., 1989a), in which soil horizons are newly formed; similar material is uncovered in the pit as well.The pit-mound microtopography is gradually levelled over time (with a maximum pit-mound longevity up to circa 1700 years, Šamonil et al. 2013b), but the unique soil features persist for a longer period (e.g.Pawluk and Dudas, 1982).Our latest soil analyses showed an increasing proportion of the coarse soil fraction (0.1-2.0 mm) in the upper mineral A horizon (unaffected by pedogenesis in deeper horizons) along the gradient of soil weathering and leaching from Haplic Cambisols to Albic Podzols, which could be the footprint of past disturbances.
On the forest stand scale the overall clustered arrangement of treethrows was connected with the simultaneous uprooting of several trees within one gap and also with local historical contingency (Phillips, 2004), which influences new disturbance events.Most tree species at the site significantly prefer treethrow mounds to other microsites (Šebková et al., 2012).This is caused by more favourable microclimatic conditions (Beatty and Stone, 1986) and the lower thickness of organic horizons on mounds (Šamonil et al., 2008a, b).The affinity of regeneration for mounds (and gaps after treethrows) increases the probability of the occurrence of further pitmounds in the given microsite, which is further supported by unstable mounds (see Mayer et al., 1989).The significantly higher tendency of treethrows towards clustering on Albic Podzols is most probably an extreme example of this connection.A more frequent occurrence of the humic Hh-and residual Hr-form of the humification horizon (Klinka, 1997) and an increase in their thickness was observed on Albic Podzols (Šamonil et al., 2011).An increase in the thickness of the forest floor on Podzols can facilitate the affinity of tree species for bare treethrow mounds and for gaps after treethrows.On Podzols, the higher density of treethrows probably leads to a shorter rotation period than the usual 1380 years (Šamonil et al., 2013b).The significant tendency towards clustering on Podzols can eventually result in a further shortening of the rotation period.However, because of the significant spatial non-randomness of disturbance events, one cannot exclude the existence of places that have not been disturbed for the entire Holocene.Our latest, still unpublished findings arising from pedoanthracological, geophysical and radiocarbon analyses of currently undisturbed areas within Žofín support such an idea.
Not only the biomechanical but also the biochemical effects of trees likely influence pedodiversity in Žofín.The characteristics of the top organic horizons on Albic Podzols may be related to a higher proportion of needle-leaved trees in such stands in the past.Albic Podzols are often found near Gleysols and Histosols, which are ecologically favourable for conifers, especially Picea abies.A constant source of diaspores could influence nearby terrestrial stands in the long run, which in general are more suitable for Fagus sylvatica and other broadleaved trees.The hypothesis that conifers locally dominated on Albic Podzols in the ancient past and had a significant effect on pedogenesis (through the effect of decomposing needles) still needs to be tested (see above).A short-term tree-census in Žofín since 1975 apparently supports the hypothesis.At the same time, it supports the concept of trees as ecosystem engineers (Wilby, 2002;Corenblit et al., 2011).
However, tree-soil interactions operate reciprocally and the disturbance history of soil units may be partly determined by soils.This element of tree-soil interactions was the most obvious on Gleysols and Histosols; the occurrence of these soils does not depend on the effects of trees but on topographic conditions (platforms and concave forms).Disturbance histories of these hydromorphic soils were completely unique.Mutual connections between the environment and tree layer disturbances move the dynamics of forests from the level of one-way biogeomorphology to ecogeomorphology (Corenblit et al., 2011).So far this topic has been studied mostly in fluvial systems (Fisher et al., 2007).

Figure 1 .
Figure 1.The study area: Žofínský Prales National Nature Reserve in the Novohradské mountains, Czech Republic.The dashed lineshows the border between the (probably) historically anthropogenically affected (32 ha) and unaffected (42 ha) zones; the grey solid lines represent contour lines; the black solid lines represent streams; grey areas show places significantly affected by water (semihydromorphic and hydromorphic soils).The study area was overlain with a rectangular network of 353 sample plots (grid 44.25 × 44.25 m); five regularly located soil profiles were studied in terms of soil properties within each plot (crosses in the detail); water-affected sites were evaluated separately from a rectangular network (for details seeŠamonil et al., 2011).One rectangle of 16.04 ha served for the spatial analysis of releases in core series; eight rectangles of 1.50 ha served for the spatial analysis of pit-mounds.

Figure 3 .
Figure 3. Memory length and complexity of individual data sets.The dash-dotted line represents the zone of our main interest, where we hypothetically expect significant relations between data sets.

Figure 4 .
Figure 4. Disturbance history of the Žofín core zone based on dendrochronological data in relation to the predominant soil taxonomical units (Fig. 2).Columns represent the number of juvenile trees (lower part of the figure) and the canopy area disturbances (upper part of the figure; the height of white, grey and black bars indicate the canopy area disturbed separately for weak, moderate and major release events in radial growth).Sum of canopy is the relevant area of all exposed crowns in the dendrochronological record, i.e. sample depth in individual decades.This variable is dependent on the number of cored trees within soil unit and relation between DBH and crown area (Table1,Šamonil et  al., 2013a).The canopy areas were estimated for every year a cored tree was in the record, crown areas gradually decreased to the past.Chronologies were truncated when the sample depth (refers to the second y-axis) dropped below 15 trees (for details, see Sect.3.4.1.).

Figure 5 .
Figure 5. Statistical similarities of individual soil taxonomical units with respect to their disturbance histories.Cluster analyses were calculated separately for all releases ≥ 20 % BL (left) and only stronger releases ≥ 50 % BL (right).Heights refer to Euclidian distances between STUs; for details see Sect.3.4.1.

Figure 6 .
Figure 6."Footprints" of disturbance events in the Žofín core zone.The left panel represents the number of releases per canopy tree; only trees which have already reached the tree canopy are shown (n = 682).Locations of trees were used for evaluating the spatial pattern of releases (for details see Sect.3.4.1).The right panel shows the occurrence of treethrow pit-mounds with or already without the uprooted tree.Locations of pit-mounds served to calculate spatial relations between soil disturbances (see Sect. 3.4.2).

Figure 7 .
Figure7.Test statistic for random labelling of releases per canopy trees.We used the "i-to-any" summary g i• (r), which is the analogue of the pair correlation function.Under null hypothesis of random labelling (g r• (r) − g(r)) theo (r) is valid g i• (r) = g(r).We generated 199 simulations of our null model to obtain pointwise critical envelopes for this model.If the value of statistic (g r• (r)−g(r)) obs (r) is larger than the value of the (g r• (r) − g(r)) hi (r) then the releases show significant clustering.If the value of (g r• (r) − g(r)) obs (r) is smaller than value of the (g r• (r) − g(r)) lo (r) then the releases show a regular distribution.In the grey zone, we cannot reject the null hypothesis of random labelling.The variable r refers to distance (for details see Sect.3.4.1).

Figure 8 .
Figure 8. Spatial patterns of uprootings at sites dominated by various soils.Significant spatial distributions were determined through 199 Monte Carlo simulations of the null model of complete spatial randomness.If the value of g hi (r) is larger than the value of the g hi (r) then the uprootings show significant clustering.If the g obs (r) is smaller than value of the g lo (r) then the recruits show a regular distribution.In the grey zone, we cannot reject the null hypothesis of complete spatial randomness.The variable r refers to distance (for details see Sect.3.4.2).

Table 1 .
Descriptive statistics of selected characteristics of disturbances in relation to main soil taxonomical units in the Žofín core zone.Box-Cox transformation and Kruskal-Wallis test were applied (within a row, numbers with different letters significantly differ at α = 0.05).All confidence intervals are 95 %.

Table 2 .
Euclidean distances (under diagonal) between dendrochronological disturbance histories of forests on individual soil taxonomical units (STUs) and the statistical significance of their similarities (above diagonal, for details see Sect.3.4.1 and Fig.5).