Boreal coniferous forest density leads to significant variations in soil physical and geochemical properties

At the northernmost extent of the managed forest in Quebec, Canada, the boreal forest is currently undergoing an ecological transition between two forest ecosystems. Open lichen woodlands (LW) are spreading southward at the expense of more productive closed-canopy black spruce– moss forests (MF). The objective of this study was to investigate whether soil properties could distinguish MF from LW in the transition zone where both ecosystem types coexist. This study brings out clear evidence that differences in vegetation cover can lead to significant variations in soil physical and geochemical properties. Here, we showed that soil carbon, exchangeable cations, and iron and aluminium crystallinity vary between boreal closed-canopy forests and open lichen woodlands, likely attributed to variations in soil microclimatic conditions. All the soils studied were typical podzolic soil profiles evolved from glacial till deposits that shared a similar texture of the C layer. However, soil humus and the B layer varied in thickness and chemistry between the two forest ecosystems at the pedon scale. Multivariate analyses of variance were used to evaluate how soil properties could help distinguish the two types at the site scale. MF humus (FH horizons horizons composing the O layer) showed significantly higher concentrations of organic carbon and nitrogen and of the main exchangeable base cations (Ca, Mg) than LW soils. The B horizon of LW sites held higher concentrations of total Al and Fe oxides and particularly greater concentrations of inorganic amorphous Fe oxides than MF mineral soils, while showing a thinner B layer. Overall, our results show that MF store three times more organic carbon in their soils (B+FH horizons, roots apart) than LW. We suggest that variations in soil properties between MF and LW are linked to a cascade of events involving the impacts of natural disturbances such as wildfires on forest regeneration that determines the vegetation structure (stand density) and composition (ground cover type) and their subsequent consequences on soil environmental parameters (moisture, radiation rate, redox conditions, etc.). Our data underline significant differences in soil biogeochemistry under different forest ecosystems and reveal the importance of interactions in the soil–vegetation–climate system for the determination of soil composition.

Abstract. At the northernmost extent of the managed forest in Quebec, Canada, the boreal forest is currently undergoing an ecological transition between two forest ecosystems. Open lichen woodlands (LW) are spreading southward at the expense of more productive closed-canopy black sprucemoss forests (MF). The objective of this study was to investigate whether soil properties could distinguish MF from LW in the transition zone where both ecosystem types coexist. This study brings out clear evidence that differences in vegetation cover can lead to significant variations in soil physical and geochemical properties.
Here, we showed that soil carbon, exchangeable cations, and iron and aluminium crystallinity vary between boreal closed-canopy forests and open lichen woodlands, likely attributed to variations in soil microclimatic conditions. All the soils studied were typical podzolic soil profiles evolved from glacial till deposits that shared a similar texture of the C layer. However, soil humus and the B layer varied in thickness and chemistry between the two forest ecosystems at the pedon scale. Multivariate analyses of variance were used to evaluate how soil properties could help distinguish the two types at the site scale. MF humus (FH horizons horizons composing the O layer) showed significantly higher concentrations of organic carbon and nitrogen and of the main exchangeable base cations (Ca, Mg) than LW soils. The B horizon of LW sites held higher concentrations of total Al and Fe oxides and particularly greater concentrations of inorganic amorphous Fe oxides than MF mineral soils, while showing a thinner B layer. Overall, our results show that MF store three times more organic carbon in their soils (B+FH horizons, roots apart) than LW. We suggest that variations in soil properties between MF and LW are linked to a cascade of events involving the impacts of natural disturbances such as wildfires on forest regeneration that determines the vegetation structure (stand density) and composition (ground cover type) and their subsequent consequences on soil environmental parameters (moisture, radiation rate, redox conditions, etc.). Our data underline significant differences in soil biogeochemistry under different forest ecosystems and reveal the importance of interactions in the soil-vegetation-climate system for the determination of soil composition.
tem (Richter and Yaalon, 2012;Van der Putten et al., 2013). Land-use changes or ecosystem shifts can have a wide range of impacts on soil properties such as nutrient availability, organic matter content, soil structure, erosion or soil water repellency (Li and Richter, 2012;Van der Putten et al., 2013;Willis et al., 1997). Concurrently, soil properties are of great importance for vegetation establishment and maintenance in space and time (Kardol et al., 2006).
In the central portion of Quebec's boreal zone, Canada, a southward expansion of open black spruce-lichen woodlands (hereafter LW) is currently being observed at the expense of more closed black spruce-moss forests (hereafter MF; Bernier et al., 2011;Girard et al., 2008Girard et al., , 2009Rapanoela et al., 2016). The current ecosystem shift could be due to a change in the regional fire regime (Ali et al., 2012;Rapanoela et al., 2016) that likely occurred several thousands of years ago (Richard, 1979;Asselin and Payette, 2005) and constitutes a hot stake raised by forest ecologists and forest managers in so far as open black spruce-lichen woodlands are less productive and, consequently, sequester less carbon than closed moss forests (Rapanoela et al., 2016;Van Bogaert et al., 2015). MF are characterised by dense stands mainly composed of black spruce (Picea mariana) with a ground layer dominated by feather mosses (Pleurozium schreberi and others) and sphagnum (Sphagnum spp.). In LW, black spruce also stands as the dominant and quasi-exclusive tree species, yet tree density cover is a lot scarcer and the predominant ground cover vegetation is composed of lichens, mostly Cladonia spp.
The main factors influencing soil formation, namely climate, organisms, topography, parent material and time, also determine soil properties (Jenny, 1994;Lundström et al., 2000;Mourier et al., 2010). Soil formation is polygenetic and sensitive to ecosystem and vegetation changes (Richter and Yaalon, 2012). Transformation can be perceptible within a few years or decades, as it has already been observed following land-use changes (e.g. Li and Richter, 2012;Richter and Yaalon, 2012). Some vegetation communities have commonly been observed in association with specific soil types (Lundström et al., 2000, Mourier et al., 2010, Willis et al., 1997. Podzols are typically found in boreal spruce forests under well-drained conditions (Sanborn et al., 2011;Ugolini et al., 1981). Yet it is still unclear (i) what part pre-existing soil properties play in the establishment and maintenance of the different ecosystems and (ii) how the persistence of a specific vegetation cover contributes to changing or maintaining specific soil properties. It is likely that both mechanisms coexist by relying on many feedbacks.
We have several reasons to hypothesise that ecosystems with different ground vegetation types, humus thicknesses and stand densities such as MF and LW could induce different soil geochemistry and soil horizon development. Differences in canopy openness and groundcover type may lead to variations in soil formation processes. MF soils should develop thicker accumulation (B) soil horizons as more or-ganic matter inputs are provided by the denser vegetation, leading to higher nutrient availability (Bonan and Shugart, 1989;Haughian and Burton, 2015). Conversely, the lowdensity canopy in LW should provide a lower nutrient supply, along with higher light and heat (Haughian and Burton, 2015;Sulyma and Coxson, 2001), leading to different soil microclimatic conditions. Variations in hydrological processes are to be expected beneath both ground covers as moss layers, which have a high water holding capacity, create a saturated environment that is less favourable to decomposition, whereas lichen layers may maintain soil moisture at lower levels (Bonan and Shugart, 1989). Water fluxes, insulation properties and soil acidity may have a great influence on chemical reactions during pedogenesis, particularly redox conditions, chemical element associations, and exchangeable cation circulation and mobilisation (Schwertmann, 1985;Brimhall and Dietrich, 1987). Lichens are also reported to be strong physical and chemical weathering agents of rock surfaces (Chen et al., 2000;Porada et al., 2014), while moss layer thickness may lead to a diffuse and weaker weathering. Differences in pedoclimate (Duchaufour, 1990) hydrological processes such as water fluxes, snowmelt rate, and drainage conditions (Buurman and Jongmans, 2005;Schaetzl and Isard, 1996), organic matter dynamics (Buurman and Jongmans, 2005) or weathering rate (Lundström et al., 2000) could deeply impact the podzolisation process.
Our main research objective was therefore to determine how soil biogeochemical properties and soil horizon development differ between two forest ecosystems, which are defined by different stand densities and ground vegetation compositions, but which both developed from glacial deposit in the transition region where both open-and closed-canopy forests occur as adjacent patches at the landscape scale. This is the first attempt, to our knowledge, in trying to identify a biogeochemical signature of boreal forest with open versus closed-canopy conditions -two states that are common to boreal forests of North America and Eurasia. In particular, this pedological investigation constitutes a first step to understand the implications of changes in boreal forest canopy openness that have happened and that will happen in the future.
We hypothesised that soils would be richer under dense forests with a moss cover that has higher C and N contents and higher base cation concentrations than lichen cover. Because iron complexation with organic compounds is a dominant process in podzolic B mineral horizons and is sensitive to properties such as soil pH, organic matter content or soil erosion (Holmgren, 1967;Li and Richter, 2012;McBride, 1987), we also hypothesised that the podzolisation process and hence iron reactive chemical species would be different depending on the local vegetation density. We assumed that the conditions found in soils covered by a lichen mat with low-density canopy would be more prone to iron and aluminium oxide accumulations than MF soils. Schaetzl et al. (2015) argued that water fluxes have a great influence on Biogeosciences, 14, 3445-3459, 2017 www.biogeosciences.net/14/3445/2017/ the intensity of podzolisation especially because they control the mobility of soluble organic complexes onto the soil profiles. Snow, snowmelt and deep percolation may thus vary between MF and LW because of tree density and lichen or moss cover. We therefore used an exploratory approach focusing on soil horizon thickness, carbon and nitrogen contents, base cation concentrations (Ca, Mg, K, Na), and species of iron (Fe) and aluminium (Al; organically bound and oxides) as elements of interest and investigated their relations with vegetation properties. Our study was first conducted at the pedon scale for a local approach (lichen vs. moss cover) on a 1 m 2 basis, and analyses were secondly run at the site scale (∼ 10 ha, LW vs. MF) to assess whether local observations could be generalised to the ecosystem.

Study sites
The study area was located in the Canadian northern region of the Manicouagan crater in Quebec in the north portion of the moss forest ecological domain (Fig. 1). Our study area covered approximately 900 km 2 . Climatic data from the nearest reported station in Wabush Lake (52 • 55 N, 66 • 52 W; 551 m elevation, located 150-200 km northwest of our study sites) show a mean annual temperature of −3.1 ± 3.3 • C and mean annual precipitation of 840 mm, with 51 % falling as snow, for the 1981-2010 period. The number of degree days above 5 • C for the period 1981-2010 was 817 per year (Environment Canada, 2013). Black spruce (Picea mariana (Mill.) BSP) was the dominant tree species in the area, but balsam fir (Abies balsamea (L.) Mill) was also scarcely found in MF. Groundcover vegetation was mostly composed of feather mosses (Pleurozium schreberi and others) with occasional patches of Sphagnum spp. in MF, while it was mostly composed of Cladonia spp. in LW. While vegetation in the area was essentially composed of black spruce-lichen woodlands north of the 52nd parallel and of black spruce-moss forests in the south, patches of both ecosystems could be found in a patchy distribution throughout the area. The study area was particularly relevant for our comparative study in so far as it presented a patchy distribution of plots covered by moss and lichen developed from glacial deposit, within regional sites of LW and MF ecosystems. Types of deposit were still slightly different from site to site (e.g. undifferentiated till, dead-ice moraine). We selected six independent sites (experimental units or EUs): three sites dominated by an MF ecosystem and three other sites dominated by an LW ecosystem (Table 1). Sites were preselected based on satellite images and maps of the vegetation ecological domains and of the surface deposit types to focus on ecological region 6P (côteaux du lac Caopacho) defined by the Ministère des Forêts, de la Faune et des Parcs du Québec. Site selec-tion was then validated through a field prospection. Each site was centred on a small headwater lake. Three transects were delineated per site; they extended from the lake shore out of the riparian zone up to 30 m. Sample plots were placed along each transect at distances of 10, 20 and 30 m from the end of the riparian area surrounding the lake, identified by the presence of Chamaedaphne calyculata. At each plot, soil samples were collected as described below. Stand characteristics were evaluated by measuring the basal area using a wedge prism relascope (factor 2) and by listing the measured tree species. Groundcover vegetation composition and abundance (in %) was estimated on a 1 m 2 area per plot prior to soil sampling; visual estimates were summed to 100 %. Our design for a given forest type was therefore composed of three EUs and nine sampling plots per EU, for a total of 27 soil sampling plots per ecosystem type (MF vs. LW, 54 altogether). This experimental design was conceived for further investigations linking soil and lake sediment composition at the watershed scale. Indeed, it corresponds to the first step of a palaeoecological investigation aiming to retrace the opening of the landscape over time using geochemistry analysis from lacustrine deposits. Before proceeding further, it was fundamental to test and demonstrate that nowadays the two types of ecosystem (LW and MF) display significant differences in soil properties at the watershed scale. Although the ecosystem type was determined at the lake watershed scale, plots within a site could display non-typical cover vegetation. For instance, local lichen patches could be found in some plots belonging to the MF ecosystem and vice versa. For practical purposes, lower case letters were therefore used when referring to vegetation type at the plot scale (lw and mf), while upper case letters were used when referring to the lake watershed/ecosystem (LW, MF). The most distant sites were 100 km apart. A general description of each site is given in Table 1. For statistical analyses including lichen and moss covers, binary values were attributed to plots depending on the dominant vegetation cover type.

Soil sampling and treatment
Soils in the study area were typical podzols whose development in northern Quebec dated back to the colonisation of bare grounds by the boreal forest after glaciation 9000 to 13 000 years ago (Lundström et al., 2000). Podzolisation is favoured by conifer stands and is formed by precipitation of litter organic matter that engages the mobilisation of soil Al and Fe (Buurman and Jongmans, 2005). The different methods used for sampling individual soil horizons are schematised on Fig. S1 in the Supplement. A 20 × 20 cm 2 template was used to sample the organic layer (FH horizon) after discarding the live green moss or lichen portion. Mineral soil fractions were sampled in each horizon. Soil samples of the mineral B and C horizons were collected in every plot (total N = 54 for both forest types) in June and September 2015 and in one out of the three plots per transect (10 m from the Table 1. General information and description of the six sites of study. LW is lichen woodland, and MF is moss forest. Values of basal area, ground cover composition and soil horizon thickness are given as means ± standard deviations to illustrate intra-site variability. Means where calculated from nine plots per site, regardless of the local ground cover type, but rather considering the regional ecosystem-type of the site they belonged to. NA means not determined during field sampling because of a fire that burned the forest in 2007. 18.9 ± 6.9 39 ± 13 57 ± 10 100 ± 0 0 ± 0 80.1 ± 11 10.8 ± 7.5 9 ± 4.5 68 • 0'45.545 W Figure 1. Overview of the study area and sites. (a) Geographical localisation of the study area and distribution of the study sites (experimental units). The study area is located in Quebec, Canada, at a latitude of 52 • N and a longitude of 67-68 • W. The right panel presents the sampling design as undertaken in each EU around the watershed lake with three transects and three sampled distances (plots) by transect using a basemap from google maps openlayers. (b) Open black spruce-lichen woodland. (c) Closed black spruce-moss forest. lakeshore) for the FH horizon (N = 18). In addition, the top 15 cm of the B horizon was collected volumetrically using a 5 cm diameter metal corer. Samples of the C horizon were retrieved with a soil auger as soon as the B-to-C horizon limit was reached. All soil samples were air dried, sieved at 2 mm and weighed. This fraction was used for conventional analysis (pH, texture). Samples were ground at 500 µm for geochemical analyses. C and N stocks as well as those of Al and Fe species were estimated and reported for the FH layer and the B top layer.

Geochemical analyses 2.3.1 Soil primary characteristics
Total C and N contents (%) were measured on all soil samples by combustion using an induction furnace (Leco ® Tru-Mac CNS Analyzer) following sieving at 2 mm, drying and grinding at 0.5 mm. The combustion is performed at 1350 • C under an oxygen gas atmosphere which turns C and N forms to CO 2 , N 2 and NO x . Gaz concentrations are then determined by thermal conductivity and infrared detection. pH values and effervescence to HCl indicated the absence of Biogeosciences, 14, 3445-3459, 2017 www.biogeosciences.net/14/3445/2017/ carbonate. pH was determined by potentiometric titrations in deionised water. Soil texture was determined by soil fractionation and grain size sedimentation following instructions from Carter (1993).

Exchangeable ions and extractable phosphorous
Soil samples were treated with an extracting Mehlich 3 solution (CH 3 COOH 0.2M, NH 4 NO 3 0.25M, NH 4 F 0.015M, HNO 3 0.013M, EDTA 0.001M) at a 1 : 10 ratio (Mehlich, 1984). Concentrations of exchangeable ions and extractable phosphorous (P) were then analysed by inductively coupled plasma atomic emission spectroscopy (ICP-AES). The Melich 3 solution mainly extracts exchangeable and soluble cations and phosphorus under aluminium, calcium and iron phosphate forms. Cationic exchange capacity (CEC) was calculated as the sum of exchangeable cation concentrations (K, Ca, Mg, Mn, Al, Fe, Na). Base saturation was calculated as the sum of main base cations (Ca + Mg + K + Na) divided by the CEC.

Reactive Fe and Al chemical species
Podzols are characterised by a B horizon predominantly made up of amorphous material constituted of organic matter bounded at different degrees to Al and Fe oxides and hydroxides (Canadian Soil Survey Committee, 1978;Lundström et al., 2000). Chemically bound species of Fe and Al were extracted from mineral soil samples by means of three chemical methods that rely on different reduction reagents: sodium pyrophosphate (Na 4 P 2 O 7 ), ammonium oxalate (C 2 H 8 N 2 O 4 ) and dithionite-citrate-bicarbonate (Na 2 S 2 O 4 , Na 3 C 6 H 5 O 7 , NaHCO 3 ; Pagé and Kimpe, 1989). We followed the extraction protocols of Mehra and Jackson (1960) and McKeague (1978). Concentrations of extracted ions were then analysed by ICP-AES. Species extracted by pyrophosphate, oxalate and dithionite-citrate reagents will be respectively referred to as "pyro", "oxa" and "dit" hereafter. Pyrophosphate is the weakest extractor, known to specifically isolate organically bound iron (Fe pyro ). Oxalate removes Fe pyro as well as the inorganic amorphous iron (Fe oxa ). Dithionitecitrate removes Fe pyro , Fe oxa and crystalline iron (Fe dit ), i.e. most Fe species (Mehra and Jackson, 1960;Blume and Schwertmann, 1969;McKeague et al., 1971). Al species extracted by all three methods are reported to be of similar nature as Fe species, although extractions may be less specific and result in some overlaps (McKeague et al., 1971). In particular, quantities of Al extracted by oxalate (Al oxa ) may be higher than quantities extracted by dithionite-citrate (Al dit ) in some cases such as in acid soils or podzols (Johnson and Todd, 1983;Pagé and Kimpe, 1989). McKeague et al. (1971) showed that these two extractants are less useful for distinguishing species of Al in soils than for Fe species. Relying on iron extraction specificity, relative quantities of crystalline iron (e.g. goethite, hematite) can be obtained by subtracting the quantities of iron extracted with oxalate from those extracted with dithionite-citrate (Fe CRI = Fe dit −Fe oxa ). Amorphous inorganic iron, also designed as short-range order (SRO) mineral phases, is calculated by Fe SRO = Fe oxa − Fe pyro (Johnson and Todd, 1983;Pagé and Kimpe, 1989). Finally, we computed the commonly used iron crystallinity index (CI = Fe dit : Fe oxa ) in order to asses soil development rate (Arduino et al., 1984;Blume and Schwertmann, 1969). Amorphous iron species gradually aggregate into crystalline forms during the pedogenesis process, which makes them good indicators of soil age (Arduino et al., 1984;Johnson and Todd, 1983).

Statistical analyses
To test the influence of spatial scale on associations of soil properties and belowground geochemistry between forest types, we performed statistical analyses at two different spatial scales: one at the soil pedon scale (plot scale) and another at the lacustrine watershed scale (site scale). At the plot scale, we tested for differences in each soil property between lw and mf using a mixed analysis of variance (ANOVA) with the variables "site" as random intercept, "transect" nested into "site" as random intercept, and the dominance of moss vs. lichen as the fixed effect. All mixed models were fitted using the nlme R package (Pinheiro et al., 2016) in the R environment (R Core Team, 2013). We then assessed covariances among aboveground attributes of vegetation (stand basal area and percent cover of moss and lichen) and chemical variables using partial least squares canonical analysis (PLSCA) and cross-correlation tests. PLSCA identifies latent variables ("orthogonal canonical component") that maximise correlations between two sets of variables given symmetric roles. PLSCA also makes it possible to calculate the proportion of variance explained by each orthogonal canonical component of one set into the other and to determine the "intragroup community index" (Tenenhaus, 1998) that represents the weight of one variable in a set in explaining variations of variables in the other set. To this effect, we grouped the vegetation variables estimated at the plot scale (e.g. basal area, moss and lichen cover) into one set of variables and the soil geochemical variables into another set. PLSCA was performed using the plsdepot R package (Sanchez, 2012). Then, to determine whether association patterns between vegetation attributes and geochemical soil properties at the soil profile level could be scaled up at the site scale, we performed a permutational multivariate analysis of variance (PERMANOVA) to test whether geochemical soil properties differed among ecosystem types (LW vs. MF). Following the parsimony principle, we restricted the number of soil chemical variables in our analyses to the five elements (Al oxa , Fe SRO , Al SRO , Ca and Fe) that had the highest intragroup community indexes relative to vegetation components in the PLSCA. The justification for this choice is that variables that have a low intragroup community index are loosely linked to www.biogeosciences.net/14/3445/2017/ Biogeosciences, 14, 3445-3459, 2017 variables of the other set (Tenenhaus, 1998). PERMANOVA performs a permutational MANOVA procedure (Anderson, 2001;McArdle and Anderson, 2001) on similarity matrices (here obtained using Euclidean distances) taking the hierarchical structure of the experimental design into account in order to test whether the locations of centroids differ among the groups of interest in multivariate space. The three tested factors were (1) ecosystem type (LW vs. MF) as the fixed effect, (2) site as the random effect (N = 6) and (3) transect as the random effect nested within site (N = 18, or 3 per site). Most importantly, we avoided any pseudo-replication issue by testing the effect of ecosystem type (LW vs. MF) with the random variable "site" as the error term (ddl of error term = 4; see Table 5). Moreover, to ensure valid results, we normalised our data prior to analysis and used the Monte Carlo method to maximise the number of permutation combinations when estimating the p value (Anderson, 2001;McArdle and Anderson, 2001). PERMANOVA was performed using Primer 6 software (Anderson, 2001;McArdle and Anderson, 2001). We also tested whether multivariate dispersions around each group's centroid were homogeneous between the LW and MF ecosystem types using the vegan R package (Oksanen et al., 2016).
3 Results and discussion

Soil profile analysis
Soil profiles from the two vegetation cover types showed visually prima facie pedological differences at the plot scale. All horizons were thicker under the mf cover than under lw cover (Fig. 2a, Table 2). Colours, although not defined with a soil handbook, also differed among soil profiles: while we observed a reddish to light yellow B horizon in lw plots, they were darker and browner in mf plots (Fig. 2b and c). Colour hues of mineral horizons have been reported to reflect various species and concentrations of Fe(III) oxides (Arduino et al., 1984;Schwertmann, 1985), known as indicators of soil age and soil weathering rate. The more important the pedogenic process of Fe reduction and subsequent removal are, the less colourful the soils are, displaying mostly the grey colours of the silicate matrix (Schwertmann, 1985). It is thus likely that pedogenic processes differ between lw and mf plots. Because all sites developed from a similar geomorphological base (glacial deposit), we thought that the differences observed could portend different developmental or formation processes depending on the vegetation cover. Regarding texture prospection, we found no significant differences (p value > 0.1) in the mean percentages of sand, clay and silt between mf and lw samples of the C horizon nor in the mean percentage of sand and clay of the B horizon (Table 1).

Soil chemical properties
Regarding the FH humus layer, C and N concentrations were 30 to 50 % higher in mf plots than in lw plots (Fig. 3, Table 2). Concentrations of the two dominant exchangeable base cations (Ca, Mg) were respectively 2.5 and 1.5 times higher in mf humus than in lw humus (Fig. 3). This translated into direct consequences on humus cation exchange capacity, which was ∼ 50 % higher in mf than in lw. However, no difference in humus pH was observed between mf and lw ( Table 2). With regard to the B horizon, contrary to the FH horizon, we detected no difference in C or N concentrations between mf and lw plots (Fig. 4, Table 2). As previously reported in acid soils of North American boreal coniferous forests (Bonan, 1990), the N concentration in mineral soil was very low (0.5-0.8 %). The B horizon appeared to have higher Ca and Mg concentrations in mf than in lw, while other ion concentrations were similar between plot types (Fig. 4). Regarding (c) Picture taken in the field of a soil profile representative of mf. Soil profiles were composed of a thin eluvial pale grey horizon (Ae), a mineral B horizon whose colour varied between sites from light, reddish to dark brown and faded to paler colours with depth, and a grey mineral C horizon.  the pH, mf B horizons were ∼ 5 % more acidic than lw B horizons (Table 2).
Concerning the C horizon, as expected, C percentage was lower in this deep horizon than in the FH and B horizons (Table S1 in the Supplement). However, the percentage of C was ∼ 2 times higher in the C horizons of mf plots than in that of lw plots (Table S1). This organic enrichment can only have a biological origin coming from the upper layers as no C is provided by the mineral parent material in these acidic soils evolved from a granitic bedrock. Similarly, extractable phosphorous (P) was 2.5 times higher in the lw C horizon than in that of mf. The accumulation of products of mineral weathering as well as the migration of organic P compounds could explain this difference.

Fe and Al reactive species in mineral horizons
Our results for the different Fe and Al species in the B horizon show that ranges of organically bound metal concentrations (Fe pyro , Al pyro ) were similar in lw and mf plots (Fig. 5). However, lw B horizon exhibited 1.5 to 2.5 times higher concentrations of Fe and Al extracted by oxalate and dithionitecitrate than mf B horizon. In compliance with other studies performed in acid forest soils (Johnson and Todd, 1983), we observed higher concentrations of Al oxa than Al dit (Fig. 5). The main difference between lw and mf B horizons lay in the proportion of inorganic amorphous Fe and Al (Fe SRO , Al SRO ; Fig. 5, Table 2) while crystalline iron concentrations (Fe CRI ) were similar (Table 2).
In the C horizon, no major variations in Al and Fe species concentrations could be observed between lw and mf plots (Table S1). This deeper horizon was also less concentrated Concentration in g kg −1 Figure 5. Distribution of Fe and Al species concentrations in the B horizon of lw and mf plots. Boxplots represent the distribution around the median values. "Pyro", "oxa" and "dit" subscripts mention the type of extractions used for isolating Fe and Al species (respectively pyrophosphate, oxalate and dithionitecitrate). "SRO" stand for short-range order species (Al SRO = Al oxa -Al pyro , Fe SRO = Fe oxa -Fe pyro ). Asterisks indicate a significant difference using mixed models (see Materials and methods), with * p value < 0.05, * * p value < 0.01 and * * * p value < 0.001. than the B horizon in all reactive Fe and Al species, which is consistent with the findings that, in podzolic soils, metallic elements are more concentrated in the upper centimetres of the B horizon (Lundström et al., 2000). Because no differences in Al and Fe species concentrations were found in the C horizon as opposed to our observations in the B horizon, our result suggest that the B horizon's chemical structure and metal oxides composition may have been influenced by different pedogenetic development under mf and lw cover rather than by the mineralogical origin of their parent materials. In particular, the absence of any noticeable difference in the iron crystallinity index of B horizon between all studied sites indicated that they had identical ages (Table S2). Structural and composition variations are thus derived from other drivers than soil origin and instead depend on factors influencing horizon formation processes.
Our results diverge from the observations made by Ugolini et al. (1981), who found no morphological variations nor geochemical differences in the soil profiles of lichen tundra and spruce forest in a boreal zone of Alaska. In particular, they found no difference in reactive Fe species concentrations (Fe pyro , Fe dit ) between the two ecosystems. They concluded that time and climate were stronger drivers of soil formation in their study area than the vegetation type. We suggest that because of the proximity of our study sites and the patchy plot distribution within sites, climate could be considered homogeneous in the present study. Our results are rather in line with those of Li and Richter (2012), who showed that land-use changes between old hardwood forests, cultivated agricultural fields and old-field pine forests can induce transformation and redistribution of soil iron oxides over relatively short timeframes (Fe oxa , Fe dit ). They inferred that some differences could be due to the different erosion levels and biological activities impacting soil iron oxides and organometallic compound transformation (Li and Richter, 2012).

Phosphorus distribution
We observed low concentrations of extractable phosphorous in the B horizon, in contrast to the FH and C horizons (Figs. 3 and 4, Table S1). These results could be explained by P sorption properties of Fe SRO and Al SRO species. As a matter of fact, in the B horizon, most P is bound within Fe-P and Al-P complexes (Grand and Lavkulich, 2015;Li and Richter, 2012). Grand and Lavkulich (2015) showed that P sorption to short-range order Al and Fe mineral phases decreased the availability of labile P. In the C horizon, in addition to its enrichment in organic C, the smaller amounts of Fe and Al oxides may be the reason why more labile P is available (Table S1). Reactive Fe and Al species are known to play an important geochemical role in acidic forest soils due to their sorption properties that influence carbon and nutrients bioavailability through coupling reactions, which makes them good predictors of nutrient availability (Grand Biogeosciences, 14, 3445-3459, 2017 www.biogeosciences.net/14/3445/2017/ Ratios were calculated by dividing for each soil type concentration measured in the B horizon by that measured in the C horizon. "(e)" indicates exchangeable or extractable elements. "Pyro", "oxa" and "dit" subscripts mention the type of extractions used for isolating Fe and Al species (respectively pyrophosphate, oxalate and dithionite-citrate). "SRO" stand for short-range order species (Al SRO = Al oxa -Al pyro , Fe SRO = Fe oxa -Fe pyro ). Asterisks indicate statistical significance using mixed models (see Materials and methods), with * p value < 0.05, * * p value < 0.01 and * * * p value < 0.001. and Lavkulich, 2015; Li and Richter, 2012). Here, P sorption to short-range order Al and Fe mineral phases in the B horizon may reduce nutrient availability for plants, resulting in limited P supply. P limitation is commonly observed in acidic soils of boreal forests (Giesler et al., 2002).

Relations between B and C horizons
Soil properties were analysed by considering the local variability of the parent material and, therefore, by studying the ratio of B to C horizons for various properties. We found that B : C ratios were different between mf and lw soils for many chemical properties (e.g. Ca, Al oxa , Al dit , Al SRO , Fe dit , Fe SRO ; cf. Fig. 6), which suggests that dissimilar biogeochemical processes and vertical transfers occur locally in the soil of the two vegetation types. These differences in B : C ratios between mf and lw soils also confirm that C horizon composition is unlikely to drive most of the variation observed between lw and mf B horizons and that it is the influence of vegetation that impacts most of the soil biogeochemistry. Furthermore, the low concentration of chemical elements in the C horizon compared with B and FH horizons also invalidates the hypothesis of a deeper mineralogical influence explaining the main differences in geochemical composition between lw and mf plots. Finally, variations in soil conditions such as temperature, pH and soil hydrology could play a role in differentiating horizon composition. The thickness of the organic layer and its higher water retention capacity in mf forests could greatly affect soil processes, and this is also reflected in the thickness of the B horizon.

Covariance between vegetation and soil geochemical variables
Results of the multivariate PLS canonical analysis conducted using the five most significant compounds (according to Figs. 4 and 5, Table 2) support the hypothesis of a different biological influence of vegetation on soil chemical composition and structure that discriminates between the two forest types. Indeed, our results revealed that the five soil geochemical variables with the greatest intragroup community indexes were Al oxa , Fe SRO , Al SRO , Ca and Fe (Fig. 7, Tables 3 and 4), meaning that they were highly linked to variables in the vegetation set. Vegetation variables, on the other hand, all had an important weight in explaining the variability of soil physicogeochemical variables (Table 4). This result was consistent with stand basal area (representative of forest production) and cover type vegetation being tightly correlated (Fig. 7). Fe and Al complex species were positively correlated with each other, positively correlated with lichen cover, and negatively correlated with moss cover and basal area (Fig. 7). This result was consistent with our aforementioned quantitative observations at the plot scale. Exchangeable Fe, Ca and Mg were positively correlated with each other as well as with the vegetation characteristics of dense moss-covered stands (Fig. 7). Exchangeable Fe, Ca and Mg behaviours were very much alike, displaying negative correlations with Fe and Al oxides   Axes correspond to principal orthogonal canonical components. Positive, null or negative correlations between variables are indicated by acute, right or obtuse angles, respectively, between the corresponding vectors. "(e)" indicates exchangeable or extractable elements. Regarding the chemical elements, only variables that showed the greatest differences in mean values between lw and mf plots in the B horizon were included in the PLSCA.
(Figs. 5, S2 and S3). The different behaviours of exchangeable Fe and bound Fe could be explained by their different mobility properties and abilities, in particular since fluxes could vary under different soil environmental conditions and soil thicknesses between lw and mf plots. The conversion reactions of iron oxides depend to a large extent on pedoenvironmental factors (pH, water activity, temperature, etc.;Schwtermann, 1988). These factors vary with depth and depend on the groundcover. Furthermore, organic matter seems to have an influence on iron oxides by inhibiting their crystallinity (Borggaard et al., 1990): in mf plots, the thicker and denser organic matter layer could explain the lower concentrations of Fe and Al oxides in B horizons. Fe and Al oxides species could also differ between mf and lw plots because soil temperature and moisture are also responsible for different goethite : hematite ratios (Schwtermann, 1988). Hematic soils develop in warmer conditions and are characterised by reddish brown colours, while goethitic soils develop under colder environment and turn yellowish-brown. This is consistent with our observations of clearer red to yellow soils under lw cover where little organic matter accumulates as  opposed to mf soils overlaid by a thick dark brown organic layer which could lead to warmer temperatures. We also found positive correlations between C and N concentrations, moss coverage, and Mg and Ca concentrations in the FH organic horizon (Fig. S3). The higher base cation bioavailability in mf plots could be explained by greater inputs of organic matter to the soil surface and by a higher coniferous basal area cover. Decomposing organic matter and litter are known to be important sources of base cation supply such as Ca and Mg (Finzi et al., 1998;Grand and Lavkulich, 2015).

Differences in soil geochemistry at the site scale
We scaled up the effect of ground cover type (lw vs. mf) on Al oxa , Fe SRO , Al SRO , Ca and Fe at the site scale (lake watershed scale) to test explicitly whether the same chemical elements differed between ecosystem type (LW vs. MF) using PERMANOVA. Our results showed significant differences between MF and LW ecosystem types (P(MC) = 0.0005, pseudo-F = 12.939). The type of ecosystem alone explained 47.6 % of the total variance in the studied set of geochemical variables (Table 5). In addition, our test of multivariate  Table S3), which reinforces the conclusion that the difference observed between MF and LW originates from differences in mean values of Al oxa , Fe SRO , Al SRO , Ca and Fe rather than differences in variance (Fig. 8). This finally highlights interactions between ecosystem structure and geochemical composition of the soil at the site scale.

Total element stocks in the soil
Because LW and MF displayed variations in soil composition and thickness, we hypothesised that they should hold different total amounts of chemical species. Total net stocks of C and N were scaled up from the plot to the site scale (layer thickness × organic matter concentration) in FH and B horizons (top 15 cm). On average, the FH horizon in MF held 3.5 times higher amounts of C and 4 times more N than that in LW (Table 6). Similarly, in the B horizon, both C and N stocks were 2 times higher in MF than in LW. These results suggest that, in addition to C sequestration in a greater aerial biomass, closed moss forests also hold more C in their soil. However, regarding Fe and Al species, total net stocks were higher in LW than in MF, despite their thinner B horizons (Table S4).

Biological influence
Overall, because we did not find any difference in the geochemistry and texture of the C horizon between ecosystem types and because our sites developed from surficial deposits (undifferentiated till, dead-ice moraine) of similar origin (glacial), our results, which are in line with those of other studies, suggest a biological influence of vegetation on soil profile development and soil chemistry (Finzi et al., 1998;Haughian and Burton, 2015;Wood et al., 1984). Haughian and Burton (2015) showed that more variability in soil composition could be explained by vegetation functional groups (e.g. mosses vs. lichens) than by abiotic characteristics (soil texture, topography). They found that variations in nutrient availability were the result of differences in vegetation types (lichens, feather mosses and vascular plants) rather than the opposite (i.e. nutrient availability as a cause of vegetation patterns). In their study of podzol biogeochemical vertical stratification in hardwood forests of New Hampshire, Wood et al. (1984) concluded that the B horizon is subject to strong geochemical control that is under biological influence rather than of mineralogical origin. Finzi et al. (1998) observed different distributions of Ca, Mg, Fe and Al exchangeable cations in 0-7.5 cm mineral soils under various tree species. They found an association between tree species and soilspecific chemical properties at the tree scale and suggested that vegetation influenced soil acidity and cation cycling in the forests studied (Finzi et al., 1998). Here, we suggest  Figure 9. Schematic illustration of feedback processes between stand biomass (basal area and/or density) and soil biogeochemistry as a consequence of climatic conditions based on the present study interpretations.
that the differences we observed in soil geochemical structure and composition between MF and LW had a biological origin that resulted in repercussions of the local vegetation on soil environmental conditions, which in turn influenced soil biogeochemical processes. We suggest the following potential causes for the difference observed in soil between LW and MF ecosystems: (i) the abundance of tree cover could influence soil formation through its direct consequence of canopy openness on soil micro-environmental conditions and drainage; (ii) snowfall rate and snowmelt duration in turn may influence reductive conditions and cause differential dissolution of iron oxides that accumulate in the B horizon (Giesler et al., 2002); (iii) thickness of organic soil layers could create differential insulating properties (Lawrence and Slater, 2007); (iv) because water flow also affects the distribution of chemical elements such as labile P, Fe and Al (Giesler et al., 2002), differences in drainage conditions, slope, moisture and hydrochemical processes under moss and lichen covers (Brown et al., 2010;Haughian and Burton, 2015;Price et al., 1997) could also explain the observed variations in lw and mf soil chemical properties; and, finally, v) lichen surface weathering capacity (Chen et al., 2000;Porada et al., 2014) may be responsible for the higher concentrations of Fe and Al species in LW upper mineral horizons, while a dense moss cover may be less aggressive in mineral weathering; however, because this environment is more productive, it could generate more organic acid and favour a deeper profile development.

Soil, climate and vegetation dynamics
It has previously been suggested that the two vegetation types considered could be two ecological states and that LW could be an alternative stable state resulting from regional disturbance history (Jasinski and Payette, 2005). We propose here a diagram ( Fig. 9) that synthesises our interpretations of the possible feedback processes between climate, vegetation and soil biogeochemistry that can be considered to result from successive fires. Indeed, the progression of opencanopy forests could be due to a greater fire frequency resulting from the changing climate (Rapanoela et al., 2016). Differences in fire events may lead to direct and indirect consequences at the soil level arising from fire impacts on vegetation structure and soil properties (Certini, 2005). Fire effects on soil properties have been shown to range from negative short-term effects (removal of organic matter, erosion, loss of nutrients through volatilisation, alteration of microbial communities, etc.) to long-term consequences (enhanced productivity, impact on forest successions, etc.; Certini, 2005). While negative effects on soil properties seem to be short lived (detectable some years post-fire, at most a decade) and restricted to a few top centimetres of superficial layers (Certini, 2005), we suggest that indirect effects may be wider and could have longer-term consequences. Frequent fires may hinder the accumulation of a thick top layer of organic matter, thus leading to direct aftereffects on soil physical and microclimatic conditions. The low stem basal area (likely associated with lower tree density) inherited from frequent fires results in both higher radiative insolation and precipitation reaching the soil surface, thus setting an environment more suitable to the establishment of light-tolerant lichen rather than colonisation by moss species. Repeated fires may spearhead lichen dominance in LW by maintaining preferential environmental conditions for its colonisation and establishment (Girard et al., 2009). In return, lichen establishment could also maintain a specific soil composition that is low in nutrients because of its low primary productivity (Moore, 1980). Altogether, the disturbance regime in boreal forests could determine the ground cover vegetation type and impact soil development through both direct and indirect effects, by generating poorer soils and sustaining the establishment of less productive forests.

Conclusions
We identified clear relationships between soil and vegetation structures that are reflective of a whole integrative system relying on feedback interactions. Although the correlation patterns between the ecosystem's biological components and the soil variables seem complex, our results suggest that in comparison to closed-canopy forests, open forests with a lichen ground cover are associated with a soil impoverished in C and available nutrients that develops a thinner B horizon characterised by high concentrations in amorphous species of Fe and Al. Ecosystem productivity and carbon sequestration are affected twice in LW compared with MF: through a lower density of trees (lower basal area and fewer stems) and through nutrient-and organically limited soils. This emphasises the current economic and climatic stakes that forest opening represents in so far as it could have important consequences in terms of carbon sequestration capacity. Our study of soil compartments confirms that the current opening of black spruce forests is an ecological, economic and climatic stake with underlying long-term consequences, notably in the terrestrial carbon budget. A good understanding of the processes governing soil biogeochemistry and feedback interactions between soil and vegetation remains fundamental for forest management, es-pecially in areas of ecological transitions representing major challenges, such as the northern boundary between productive boreal forests and open lichen woodlands. Our study highlights that natural disturbances may influence landscape remodelling and ecosystem heterogeneity in much more ways than through their direct impacts on seedling regeneration and soil nutrient short-term depletion. Repeated disturbance events could have long-term consequences on soil formation and development. Disturbance history and inheritance could promote the establishment and maintenance of specific vegetation-soil systems (or ecological states). If so, soil science and biogeochemistry could become interesting proxies in disturbance ecology, notably for palaeoecological investigations aiming to reconstruct changes in vegetation. Further investigations should explore vegetation structure (basal area and tree density) and soil relations in other ecosystems and focus on microclimatic and drainage conditions as well as disturbance history as explanatory drivers of variability.