Soil total phosphorus and nitrogen explain vegetation community composition in a northern forest ecosystem near a phosphate massif

The relationship of the community composition of forest vegetation and soil nutrients were studied near the Sokli phosphate ore deposit in northern Finland. Simultaneously, the effects of the dominant species and the age of trees, rock parent material and soil layer on these nutrients were examined. For this purpose, 16 study plots were established at different distances from the phosphate ore along four transects. Phosphate mining may take place in Sokli in the future, and the vegetation surveys and soil sampling conducted at the plots can be used as a baseline status for following the possible changes that the mining may cause in the surrounding ecosystem. The total phosphorus (P) and nitrogen (N) contents of the soil humus layer were positively related with species number and abundance of the understorey vegetation, and the correlation was slightly higher with P than N. This is interesting, as N usually has the most important growth-limiting role in boreal ecosystems. The spatial variation in the content of soil elements was high both between and within plots, emphasizing the heterogeneity of the soil. Dominant tree species and the soil layer were the most important environmental variables affecting soil nutrient content. High contents of P in the humus layer (maximum 2.60 g kg⁻¹) were measured from the birch-dominated plots. As the P contents of birch leaves and leaf litter were also rather high (2.58 and 1.28 g kg⁻¹, respectively), this may imply that the leaf litter of birch forms an important source of P for the soil. The possible mining effects, together with climate change, can have an influence on the release of nutrients to plants, which may lead to alterations in the vegetation community composition in the study region.

1 Introduction

Climate and availability of soil nutrients are important factors controlling the species composition of tree stand and understorey vegetation in boreal forests (Cajander, 1909, 1949; Kuusipalo, 1985; Økland and Eilertsen, 1996). High-latitude forest ecosystems are characteristically cold, have a short growing season and are nutrient poor. Organic matter decomposition and nutrient release are usually slow in cold climates (Hobbie et al., 2002). The edaphic conditions are reflected in the growth and chemical composition of plant species, as well as in species composition of vegetation (Vinton and Burke, 1995; Salemaa et al., 2008). In addition, tree cover affects the species composition and abundance in the understorey by shading (Verheven et al., 2012; Tonteri et al., 2016) and regulating nutrient input in throughfall precipitation (Salemaa et al., 2019) and litterfall (Ukonmaanaho et al., 2008).

Nitrogen (N) and phosphorus (P) are generally the main growth-limiting nutrients for plants (Koerselman and Meuleman, 1996). Boreal forests are mostly N limited (Tamm, 1991), and fertilization with N usually speeds up forest growth (Saarsalmi and Mälkönen, 2001). Nitrogen is bound in organic material, and only a little is directly available for plants as inorganic ammonium (NH⁴⁺) and nitrate (NO³⁻) (Marschner, 1995) or as organic forms like amino acids (Näsholm et al., 2008, and references within). Phosphorus deficiency occurs in temperate and tropical forest ecosystems, but P is rarely a limiting factor in boreal forests on mineral soil (Augusto et al., 2017). However, P can be growth limiting on boreal peatlands (Moilanen et al., 2010; Brække and Salih, 2002). The ratio of soil N to soil P is significant for forest growth on a global scale (Augusto et al., 2017). Hedwall et al. (2017) found that the species richness of vascular plants in a temperate forest doubled with combined NP fertilization in southern Sweden but not when either of the nutrients was added alone. This positive effect was strongest in grass species. In boreal N-limited forests, the number of vascular plant species (grasses and forbs) increased with increasing N concentration of the organic layer (Salemaa et al., 2008). Hofmeister et al. (2009) noticed that in a temperate forest the species richness of the herb layer was higher in P-rich than P-poor soils but only if strong N limitation occurred simultaneously in the P-rich soils. However, in many regions where humans have enhanced atmospheric N deposition, the number of plant species has decreased (Dirnböck et al., 2014). For instance, high soil N was related to decreased herb-layer species richness in deciduous forests in Sweden (Dupré et al., 2002).

In this study, we analysed whether plant species composition and nutrient levels of tree leaves indicate soil total N and P at a northern boreal (Hämet-Ahti, 1981) research site in Sokli, Finland. At this site, the soil contains naturally large variations in P content. In Sokli, there is a large deposit of phosphate rock, a carbonatite complex mainly consisting of apatite (Ca₅(PO₄)₃F), which was discovered by the Mining and Steel Company of Rautaruukki Oy in 1967 (Vartiainen and Paarma, 1979). Plans to open a phosphate mine in Sokli have been on display for decades and will possibly be realized in the future. The vegetation at the carbonatite complex differs from that of the typical forests of the region (Talvitie, 1979; Pöyry Environment, 2009). Downy birch (Betula pubescens) is dominant and often the only tree species, whereas more typical forests in the region are dominated by Scots pine (Pinus sylvestris) or Norway spruce (Picea abies). Understorey vegetation at Sokli is slightly richer in herb and grass species compared with surrounding forests, where dwarf shrubs, bryophytes and lichen dominate the understorey. However, vegetation similar to that in Sokli can be found as patches elsewhere in the region.

The general aim of this study was to determine the undisturbed baseline status of the forest ecosystem in terms of soil, understorey vegetation and tree layers in the Sokli area in case there is a need to monitor the effects of phosphate mining. Phosphate mining can cause, for instance, aerial deposition of heavy metals and phosphate onto the surroundings of the mine (Reta et al., 2018), which can lead to changes in the abundance and species composition of the understorey. Vegetation, soil and foliage chemistry surveys provide data on the current state of the ecosystem (from the year 2015) that can be used as a reference level for the changes. Our specific aim was to identify which factors in the soil and tree layer explain the composition and abundance of plant species. In addition, we studied which environmental variables could explain soil nutrient contents, especially total P content.

We hypothesize that there are positive relationships between the following factors:

• a.

  N and P contents of the soil humus layer and the abundance and species composition of the understory vegetation,

• b.

  N and P contents in the topmost soil layers and the N and P contents of needle and leaf biomass,

• c.

  N and P contents in the topmost soil layers and the occurrence of birch trees in the research plots.

2 Material and methods

2.1 Site description

We established 16 study plots along four transects (A–D) around the planned Sokli mining district (67^∘48^′ N, 29^∘16^′ E) in Savukoski, eastern Lapland, in 2014 and 2015 (Fig. 1). The plots were located different distances from the phosphate ore in four transects, enabling evaluation of the possible effects of the mine in the future. No plots were located inside the mining district, as accessing and doing research at the mining district would have required a permit from the mining company. The carbonatite massif of Sokli belongs to the Devonian Kola Alkaline Province (Tuovinen et al., 2015). Nine of the plots were located in Natura 2000 conservation areas. Plots A4, A5, and A6 were in Värriö; A1 and A2 were in Yli-Nuortti; B1, B2, and B3 were in Törmäoja; and D5 was in the UK-puisto–Sompio–Kemihaara Natura area. In terms of topography, Törmäoja is a valley, reminiscent of the form of a kettle (kattilalaakso in Finnish) (Natura 2000: Standard Data Form FI1301512 and FI1301513). The central parts of the valley are treeless or the trees are at sapling stage because cold winds blowing through the valley kill the new buds in the spring. Thus, our plots at Törmäoja were on the edge of the less steep western part, where some mature trees grow. The majority of the plots had a mixed composition of at least two tree species, while in some plots the tree cover consisted of only one species. An additional factor affecting vegetation cover and species composition at our research site is reindeer herding. Since all the plots were located in areas where reindeer roam freely, we assume that the pressure caused by grazing and trampling is equal in all plots. The plots of this study, together with SMEAR 1 (Station for Measuring Ecosystem-Atmosphere Relations) at the Värriö Subarctic Research Station (67^∘46^′ N, 29^∘35^′ E) (Hari et al., 1994), serve as a gradient type network for monitoring the current status and the possible, mining-induced, changes of the environment in the future.

Figure 1 Geological map of the research area. Plots are marked with black stars. The easternmost plot is located at the SMEAR 1 Station. (source: Hakku Service, https://hakku.gtk.fi/en/locations/search, last access: 28 February 2019.)

Meteorological parameters from the years of data collection and for the climatological normal period of 1980–2010 are presented in Appendix A. The wind blows almost equally from the south-west and north-east during spring and summer, whereas in winter and autumn the prevailing wind direction is south-west (Ruuskanen et al., 2003). The growing season, when daily average temperature exceeds 5 ^∘C, lasts from June to September. Soils are haplic Podzols with sandy tills (FAO, 1988).

Table 1 Tree species composition of the research plots (dbh = diameter at breast height).

2.2 Plot set-up and vegetation characterization

The distance between two plots depended on the topography and existing roads, but generally it was about 2 km. A plot consisted of four clusters, each including three 1 m² subplots for observations and sampling (Fig. 2). The size of the whole plot was 30 m × 30 m. We recorded all tree species growing on the plots, measured their heights and diameter at breast height (dbh) (equivalent to a height of 1.3 m) (Table 1). Stem volumes were estimated using the equations of Laasasenaho (1982). We estimated tree age by measuring dbh and examining the existing approximated tree age from plot A6 at SMEAR 1, where the mature trees are about 70 years old. We considered trees with dbh 1–9.9 cm as young, 10–14.9 cm as middle-aged and > 15 cm as old. We visually assessed the cover (percent surface area) and counted the number of plant species in the understorey vegetation per plot in all 12 subplots in the summers of 2014 and 2015 (Appendix B). We used a 1 m² frame to delineate the subplot (Salemaa et al., 1999). All species in the bottom layer (bryophytes and lichens) and field layer (dwarf shrubs, tree seedlings, grasses, sedges and forbs of height < 50 cm) were included.

Figure 2 Set-up of each research plot with clusters and subplots within clusters. Trees were measured from the whole 30 m × 30 m area.

2.2.1 Sampling of soil

Altogether, 256 soil samples were collected from 16 plots using a soil corer (inner diameter 5 cm) in June 2015. The soil was sampled within a 1 m distance from the subplots (see Liski, 1995). The soil cores were separated by visual criteria into four soil horizons: the top layer, which is a mixture of litter and decomposing organic layer (F); the humus layer (O); the eluvial layer (A); and the illuvial layer (B) (see Köster et al., 2014). The rocky soil and shallow humus layer made it impossible to sample the mineral soil layers in some clusters. The soil samples from each horizon were combined into composite samples in each cluster in the field. The composite samples were air-dried, except for the organic F and O horizons, which were dried at 60 ^∘C for 48 h. Dried mineral soils were sieved with a 2 mm sieve, and the F and O horizons were milled before storing in a dry place for further analyses.

2.2.2 Sampling of needles and leaves

Five pines and five spruces per plot were chosen for needle sampling in September 2015, when the needle growth had ended. If less than five trees per species were present, all of them were chosen. Three branches (length approximately 50 cm) were taken from the upper third of the canopy using a branch saw. We took only second-order branches because cutting of first-order branches would have been too destructive to the trees (see Helmisaari, 1990). Needle age classes (C = current year; C + 1 = 1-year-old needles; C + 2 = 2-year-old needles) were separated from each branch and dried at 65 ^∘C for 48 h, milled and stored in a dry place for further analyses. The samples were combined so that there was one C, one C + 1 and one C + 2 composite needle sample per tree.

We sampled green birch leaves in July 2015 and sampled leaf litter in September 2015. Approximately 10 green leaves from 10 different trees were picked and combined (Rautio et al., 2010). Only mature, undamaged leaves were chosen. Birch litter was collected under the same tree canopies from which the green leaves had been taken and in approximately the same number as the green leaf samples. We aimed to take litter leaves shed in the current year, so that they were decomposed as little as possible. Green and litter leaves were dried at 65 ^∘C for 48 h and manually cleaned of extra material, such as soil particles and needles. The needles and the few soil particles attached on the litter leaves were removed with tweezers. The green leaves did not need cleaning. The litter leaves were also rather clean, as it had rained at the time of sampling. After cleaning, the leaves were milled and stored in a dry place for further analyses. Needles and leaves were sampled at a different time than the soil. Both needle (e.g. Helmisaari, 1990) and soil nutrient contents vary between the seasons. However, as all soil and all needle sampling was conducted at the same time of the season, the comparison between the plots was not hindered.

Figure 3 Soil total phosphorus (P) contents of the research plots (based on dominant tree species) in different soil horizons. Species 1 = pine, 2 = spruce and 3 = birch.

2.3 Laboratory analyses

Total element contents of potassium (K) and P were analysed from soil and foliar samples by inductively coupled plasma optical emission spectrometry. For this analysis, the samples were first wet combusted. A 1 g amount of mineral soil sample and 0.3 g of organic sample were combusted with 1 mL of H₂O₂ and 10 mL HNO₃ and heated in a microwave oven. The samples were then filtered with Whatman Grade 589/3 filter paper and stored in plastic bottles in a cooler until they were analysed.

Total carbon (C) and N were analysed directly from dried and milled foliar samples and from the F and O soil layers. Samples of 2–3 mg were measured and analysed with an element analyser, which uses a high-temperature combustion method with subsequent gas analysis of CN (VarioMax, Elementar Analysensysteme GmbH, Germany). Soil pH was measured from two O layer samples per plot, and their average value was used. A total of 20 mg of dried sample was mixed together with ultrapure water (50 mL). The suspension was covered and left standing for 24 h, and pH was measured with a glass electrode.

2.4 Statistical analyses

We used one-way analysis of variance (ANOVA) and Tukey's honest significant difference post hoc test for analysing the plot-wise differences in needle element contents. Plot-averages of needle elemental contents were calculated across all needle age classes and both conifer species. One-way ANOVA was also used in analysing differences between the needle age classes. We grouped the plots based on their dominant tree species into pine, birch and spruce plots and calculated the average soil nutrient contents in each horizon in these plots. We then compared the nutrient contents in each soil horizon with one-way ANOVA.

We tested the effects of environmental variables on soil total P and N contents and C:N with linear mixed-effect models. We used dominant tree species, estimated age class, rock parent material (Fig. 1), and soil horizon as fixed effects and plot as random effect. Soil total P needed to be log-transformed, while for N and C:N the visual inspection of residual plots (Fig. C1 in Appendix) did not reveal obvious deviations from homoscedasticity or normality. We obtained p values for the fixed effects by likelihood ratio tests, where the full model with all the fixed effects was tested against a model where each fixed effect was removed in turn. We used package lme4 (Bates et al., 2015) in R programme 3.4.3 (R Development Core Team, 2017) for building the models. Pseudo R² values for the models were calculated by using package r2glmm (Jaeger, 2017). The models took the following form:

$$\begin{array}{}\text{(1)}& {\mathrm{SC}}_{\mathrm{P},\mathrm{N},\mathrm{CN}}=B_{\mathrm{0}}+B_{\mathrm{dt}}+B_{\mathrm{ta}}+B_{\mathrm{g}}+B_{\mathrm{h}}+\in ,\end{array}$$

where ${\mathrm{SC}}_{\mathrm{P},\mathrm{N},\mathrm{CN}}$ is the soil nutrient content (total P or N or C:N ratio), B₀ denotes a fixed intercept parameter, B_dt denotes the fixed unknown parameters associated with the dominant tree species, B_ta denotes the fixed unknown parameters associated with the age of the dominant tree species, B_g denotes the fixed unknown parameters associated with the rock parent material and B_h denotes the fixed unknown parameters associated with the soil horizon. The random effect ∈ is assumed to take the following form:

$$\begin{array}{}\text{(2)}& \in ={\propto }_{\mathrm{p}}+u,\end{array}$$

where ∝_p denotes the random parameters related to the research plot and u is an unobservable error term. Random effect parameters and the random error term are assumed to follow normal distributions ${\propto }_{\mathrm{p}}\sim N(\mathrm{0},{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}})$ and $u\sim N\left(\mathrm{0},{\mathit{\sigma }}_u^{\mathrm{2}}\right)$.

We calculated plot-wise averages from the percentage covers of the plant species in the subplots. We ordinated this vegetation data by global non-metric multidimensional scaling (Minchin, 1987) using the vegan package (Oksanen et al., 2018) in R programme 3.4.3 (R Development Core Team, 2017). Ordination pattern of the plots based on the Bray–Curtis dissimilarity indices in floristic composition was analysed to find the main environmental gradients behind the vegetation variation. We analysed the data in three-dimensional space but present the results in one vs. two and one vs. three dimensions (the results in two vs. three dimensions did not give any new information). We then fit the plot-wise data of soil elements (O horizon), needle element contents, volume of birch (percent of total tree volume), species cover (percent of the surface area) and plot distance from the phosphate ore as linear vectors to the ordination pattern of the plots. The correlation between the environmental variables and the ordination was calculated by a linear vector procedure (envfit in vegan). The total P in the soil O horizon was also fitted as a smooth surface to the ordination pattern in order to analyse the form of the relationship (linear or non-linear). The fit was done by a generalized additive model (Gaussian distribution error).

Figure 4 Soil total N content (a, b) and soil C:N ratio (c, d) of the research plots (based on dominant tree species) in soil horizons F and O (species 1 = pine; species 2 = spruce; and species 3 = birch.

Figure 5 Correlation (Pearson) including soil elements (K, P, C:N, pH, C, N:P), number of species (with _n in the end of the name) and total percentage cover of plant species in different layers (moss and lichen, grasses, herbs and sedges, dwarf shrubs and trees), needle elements (N, P, K), plot distance from Sokli, green birch elements (N, P, K) and birch litter elements (N, P, K). Levels of significance: ^* = 0.05; ${}^{**}$ = 0.01; ${}^{***}$ = 0.001. Positive correlations are displayed in blue, and negative correlations are displayed in red. Colour intensity and the size of the circle are proportional to the correlation coefficients.

3 Results

3.1 Soil element contents

The outlying points in Fig. 3. and the high standard deviations of P showed rather high variation between and within plots (Table C1, Fig. C2 in Appendix). Birch-dominated plots had the highest P and N contents and lowest C:N ratio compared with coniferous plots in all soil layers where the elements were measured (Figs. 3 and 4), and these differences were mostly statistically significant (Table C2). We found no statistical evidence for differences in soil N:P ratio or total C content based on the dominant tree species of the plots. There was only one spruce-dominated plot and thus only four soil samples from each horizon from the spruce plot, which may have affected these results. In general, topsoil had the highest P content, but in many plots deeper layers also had high P content. Certain plots (A1, A3, B1, D2) had clearly different P and N contents and C:N ratios to most other plots. The N contents and C:N ratios were in most cases higher in the F horizon than the O horizon. The majority of the plots had higher N:P ratio in the F horizon than the O horizon.

Table 2 Average foliar element contents (g kg⁻¹) of the three major nutrients and C, and the relationships of C:N and N:P with standard deviations. Needle age classes: C = current year; C + 1 = 1-year-old needles; C + 2 = 2-year-old needles.

3.2 Needle and leaf element contents

Needle P contents were highest in the C needles and significantly different from other age classes in both pine and spruce (Table D1). Against our expectations, the needle P contents of both conifer species were rather similar across plots (Table D2). On the other hand, N and C contents, as well as the C:N ratio of the conifers, showed some between-plot variation (p<0.05), but no significant variation was found in the foliar N:P ratio in either species. Spruce had slightly higher needle P contents than pine in all age classes, whereas N contents were higher in pine than in spruce needles (Table 2). Birch had higher P contents of green leaves than the conifers. Leaf litter of birch also had quite high P contents, and in general litter leaves showed more variation in element contents than the green leaves. Nitrogen contents were lower in birch leaf litter than in green leaves, but the contents of C increased slightly from green leaves to litter. However, no differences between the plots were detected in either of the elements. According to the correlation matrix between the elements in the soil O horizon, tree leaves and needles, and number of species in the understory (Fig. 5) birch K (green leaves) correlated with soil K and pH, birch litter N correlated with soil N:P, and birch litter K correlated with soil N. Otherwise, no significant correlations between foliar element contents and soil element contents were found. The number of mosses and lichen correlated negatively with soil total C (p<0.01) and C:N (p<0.05), the number of sedges and grasses positively correlated with soil total P (p<0.001) and pH (p<0.01) and negatively with soil C:N (p<0.01) and N:P (p<0.001), and the number of dwarf shrubs and trees positively correlated with soil total K (p<0.05). No significant correlations were found between number of species and foliar element contents. The cover (percent of surface) of grass and sedge species correlated positively with soil P (p<0.001) and number of sedges and grasses (p<0.05), while the cover of dwarf shrubs and trees correlated positively with the P content of green birch leaves (p<0.05).

Figure 6 Ordination pattern of the research plots in dimensions 1 and 2 (a, b) and dimensions 1 and 3 (c, d). Panels (a) and (c) give plot ordinations, and (b) and (d) give weighted averages of the most abundant species (highest cover percentage of the surface). Less abundant species are marked with light-coloured crosses. The names of species are combinations of the first four letters of their genera and species names (e.g. Solivirg stands for Solidago virgaurea). The tree species mentioned in the figure are at seedling stage. In (a) plots A4, B3, and A2 and in (c) plots A1 and B1, B2 and B3, and A4 and A2 were located on top of each other.

3.3 Mixed-effect model results

We used linear mixed-effect models for determining which environmental factors can best explain soil total P and N contents and the C:N ratio. The dominant tree species and soil horizon explained 45 % of the total P of soil, and the soil horizon explained 20 % of the total N of soil and the C:N ratio of soil (Table C3). The other fixed effects tested had p values > 0.05 and were for that reason excluded from the models. The highest estimates of P were produced with birch as the dominant tree species and F as the soil layer, and the highest estimates of N were produced with F as the soil layer in the final models.

Figure 7 Ordination pattern with smooth surface fit and linear vector fit of soil phosphorus (P) in the O layer (a), linear vector fits of soil element contents in the O layer (b), linear vector fits of foliar data and plot distance from Sokli phosphate ore (c), and linear vector fits of number of species in different layers of understorey (d). Moss_lichen includes moss and lichen species; Grass_sedge includes forb, grass, and sedge species; and D.shrubs_trees includes dwarf shrubs and tree seedlings. Plots A4, B3 and A2 were located on top of each other, and only A4 is shown.

3.4 Ordination analysis of understorey vegetation

The closer the plots were to each other in the ordination space, the more similar their vegetation was (Fig. 6a, b). Plots positioned more on the left-hand side of each panel had a higher number of forbs and grasses growing on them than the plots positioned on the right-hand side of each panel (Fig. 6a, b). Species such as Calamagrostis epigejos, Carex spp., Rubus arcticus and Luzula pilosa had relatively high coverage on the left of each panel. Plots further to the right in each panel had more species that tolerate poor and dry growing conditions, such as Cladonia lichens. The tree species also changed from right to left, as the plots on the right were dominated by pine, whereas furthest on the left in plot D2, birch was the only tree species present. In general, the fertility trend in the vegetation followed the first dimension, while the moisture gradient followed the second dimension. Moisture-demanding species, such as Equisetum sylvaticum, are located in the upper part of Fig. 6b, and those tolerating drier conditions, such as Stereocaulon tomentosum, are located in the lower part of the ordination space in Fig. 6b. Another moisture gradient, expressing specific paludified conditions, seemed to follow the third dimension. Peatland species like Sphagnum angustifolium and Aulacomnium palustre are located in the upper part of Fig. 6d, and species preferring dry conditions, such as Cetraria islandica, are in the lower part of Fig. 6d. Considering all three dimensions of ordination space, the generalist species, such as Pleurozium schreberi and Vaccinium myrtillus are located in the middle of Fig. 6b, d.

The vector arrows fitted to the ordination space (Fig. 7c, d) depict the maximum correlations between environmental variables and plot ordination. The length of an arrow indicates the magnitude and direction of the polarity (plus–minus) of the correlation. The highest correlations occurred between the plot-wise average P content of the soil O horizon and the ordination pattern of the plots (Table C4). The isocline gradient of soil P in relation to the ordination pattern was almost linear (Fig. 7a). Vectors of soil pH, N and P content all increased towards the more fertile plots, but the vectors of soil C:N and N:P went in opposite directions (Fig. 7b) indicating poor soil conditions. The average total number of grass, forb, and sedge species and their coverage in the study plots also increased towards the more fertile plots (Fig. 7d).

4 Discussion

All the plant species growing in the study plots were common forest species in Finland (e.g. Reinikainen et al., 2000, Finnish Biodiversity Information Facility https://laji.fi/en, last access: 20 March 2019) (Appendix B). However, in some plots the structure and abundance of species in the understorey clearly differed from the surrounding, more typical northern boreal forests. We found evidence that the number of species in the group of grasses and sedges, as well as the cover (percent of the surface) of the same plant group, had a higher positive correlation with humus P content than N content (Fig. 5). However, both of these nutrients were important factors explaining the vegetation composition in the ordination configuration (Fig. 7), which supports our first hypothesis. We also found that the humus C:N ratio correlated negatively with the abundance and species composition in the understorey. Additionally, Salemaa et al. (2008) have observed that total N and the C:N ratio of the humus layer explained most large-scale vegetation variation across several forest sites in Finland. They also measured extractable soil P, which seemed to have more power to explain vegetation patterns in northern Finland than in southern Finland. Soil P availability was one of the key factors in plant community variation in alpine habitats in Troms, northern Norway (Arnesen et al., 2007), where a higher variety of lichen species and the frequency of occurrence of Salix herbacea and certain sedge and grass species were explained by higher availability of P in soil. We conclude that the possible aerial deposition of phosphate from the mine in Sokli could lead to changes in plant species composition and abundance if high amounts of P are deposited into the ecosystem surrounding the mining region.

Our second hypothesis stated that the N and P contents of the topmost soil layers correlate with the N and P contents of foliar biomass, but our results (Fig. 5, Table D2) did not support this hypothesis. The reason could be that we measured total contents of N and P in soil instead of the plant-available contents of these nutrients. The plant-available contents of these nutrients might have given different results. Perhaps plant species composition in ground vegetation is sensitive to even small additions of available N and P in the upper soil layers where the roots occur, whereas higher contents of these elements are required for there to be any effect on the needles. The P and N levels of our needle samples were similar to those previously measured in Finland (Helmisaari, 1990; Merilä and Derome, 2008; Moilanen et al., 2013). The higher P contents of C needles compared with older needles is common for conifers and occurs because the dry weight in recently matured needles increases faster than the transportation of P to the needles (Helmisaari, 1990). The N contents of both green birch leaves and leaf litter agreed with those reported by Ferm and Markkola (1985). The P contents of the green leaves were higher than measured in that study (approximately 2.0 g kg⁻¹). Although the foliar N and P contents were not reflected in the uppermost soil layers, our results support the third hypothesis, and the occurrence of birch correlates positively with the N and P content of the top layers of soil (Table C2). The plots dominated by birch had significantly higher total P content in all but the B layer compared with plots dominated by the conifers (Figs. 3 and 4). Birch leaves were a major source of litter in the plots where soil P was high. These findings are supported by the study of Lukina et al. (2019), which found that the extractable P content of organic soil layers was significantly higher in birch- and spruce-dominated forest sites than in sites dominated by pine in north-western Russia. Viro (1955) found that the leaf litter of birch had remarkably high P content compared with other Finnish tree species. The litter P contents in this study were near the approximately 1.50 g kg⁻¹that Ferm and Markkola (1985) measured from a 40-year-old forest but much less than those reported from younger forests. In a litter experiment in Abisko (northern Sweden), the addition of birch litter increased both the total P (Sorensen and Michelsen, 2011) and the available P (Rinnan et al., 2008) contents in the organic soil layer in those subarctic heaths, where Hylocomium splendens dominated the moss layer. These results imply that birch is an important factor in recycling and providing P to the soil in certain types of northern forest sites.

In general, the spatial variation in soil element contents between plots was high, emphasizing the heterogeneity of soil fertility level (Figs. C2 and C3). As our results showed, this heterogeneity can partially be explained by the dominant tree species of the research plot, which especially affects the topmost soil layers. According to the nutrient-uplifting hypothesis (Jobbágy and Jackson, 2004), trees and other vegetation can transport minerals such as P and K from the deep soil layers to the surface of soils. The P contents of soil samples (Table C1) in our study (1.80–2.60 g kg⁻¹ in the O horizon) fell mostly in the category we could expect based on the literature. The P content of the humus layer in southern Finnish forest soil has been observed to vary between 0.80 and 2.10 g kg⁻¹ (Mäkipää, 1999), whereas different studies in northern Finland have found the P contents of 0.39–3.00 g kg⁻¹ in the organic topsoil (Mikkola and Sepponen, 1986; Reimann et al., 1997). Most of our plots had the highest P content in the organic soil layers, implying that decaying plant parts were a major source of P added to the soil. Low Arctic soils tend to have organic P as the primary form of P (Weintraub, 2011). The content of organic P usually gets smaller in the deeper soil (Achat et al., 2009). Thus, if P content is high in deep soil layers, as it was in some of our plots, the source of P in these plots is most likely to be the underlying bedrock. The plot-wise average pH of our soil samples agreed with that measured by Köster et al. (2014), who conducted their study at the same site, albeit not in the same plots. The pH of the soil humus layer correlated positively with the number of grass, herb and sedge species, which is reasonable, since higher pH usually implies a more fertile site. The soil N contents from our plots agreed with the reported values from Finnish forest sites (Merilä and Derome, 2008; Salemaa et al., 2008), ranging between 9.8 and 12.8 g kg⁻¹. Salemaa et al. (2008) reported a soil C:N ratio of 40 from a northern Finnish forest site, which is higher than we measured.

Our study area does not represent typical northern boreal forest, as it was located near the phosphate massif, the effect of which needs to be considered. Talvitie (1979), who used remote sensing for a geobotanical survey of the Sokli massif, found that the density of occurrence of birch, juniper and grass species increased when carbonatite was the underlying rock material. The surveys related to Natura 2000 (Standard Data Form FI1301512 and FI1301513) stated that the Törmäoja and Yli-Nuortti areas, where plots B1–B3 and A1–A2 were, have a high occurrence of grass species and a sparse birch-dominated tree cover due to carbonatite in the soil. According to the geological map (Fig. 1), only small parts in the western ends of both the Törmäoja and Yli-Nuortti Natura areas are located on top of carbonatite rock. Similarly, the map shows that those of our plots where the vegetation community was reminiscent of Sokli have something other than carbonatite as the rock parent material. However, all of our plots have metamorphic (tonalitic migmatite and amphibolite) or igneous (mafic volcanic and ultramafic) rock as the parent material, and phosphate mineral apatite can occur in such rocks (Walker and Syers, 1976). It is likely that these types of rock materials leach more phosphate than other types of bedrock (Arnesen et al., 2007). Thus, the rocks outside of the carbonatite massif may also have locally high P content, which affects the P content of the soil. The mixed-effect model factor “geology” did not consider this, which could be the reason why it was not important in explaining soil P content.

The baseline status and the current vegetation composition of our research site was worth studying for several reasons. We conducted our study in a region which has for decades been under more or less heated discussion related to whether mining activities will begin or not. The site is very remote and the plan is to move the material from the mine to locations of further production by trucks (Pöyry Environment, 2009). In addition to the aerial deposition from the mine, this could increase the dust and pollution caused by transportation, the amount of which is currently minimal. The effects of mining on the surrounding ecosystem and its vegetation composition can be unpredictable when combined with the changes caused by climate change. High-latitude regions are considered more vulnerable to climate change than more southern regions (Hartmann et al., 2013). Soil microbial activity may change due to a warmer climate, and therefore N may become more available from organic sources (Rustad et al., 2001). This, together with high soil P, may induce growth and affect vegetation dynamics. Climate change has already caused variation in the vegetation at high latitudes, as deciduous shrub coverage has expanded in the Arctic region (Sturm et al., 2001; Park et al., 2016). Greater deciduous shrub cover causes increased leaf litter input, which in turn may bring more nutrients that are recyclable to the ecosystem.

5 Conclusions

We found that the total P content of the soil humus layer was an important factor explaining the community composition of forest understorey vegetation near the Sokli phosphate ore in Finnish Lapland. The plots with high soil total P in the humus layer had birch as the dominating tree species. As green birch leaves and leaf litter also had high contents of P, we suggest that the litter caused the high total P contents in the humus layer. As climate change and the possible mining activities may affect the nutrient and vegetation dynamics in the studied region, the research that we carried out has an important role in both clarifying the current situation and forming a baseline for evaluating the magnitude of changes in the future.

Appendix A

Table A1 Meteorological parameters from Värriö. Values for the climatological normal period are from Pirinen et al. (2012). Growing degree day sum was calculated as the average daily temperature (average of daily maximum and minimum temperatures) above the 5 ^∘C base temperature, accumulated on a daily basis over the year. Negative values are treated as 0 and ignored. SWE stands for snow water equivalent. ^a Data from SMEAR 1 Station. ^b Data from SMEAR 1 Station (only 2009–2015). All other data are collected from the Värriö Subarctic Research Station by the Finnish Meteorological Institute.

Appendix B

Table B1 Average coverage (percentage of surface area) of understorey plant species per plot.

Appendix C: Soil element contents within and across plots and statistical analyses

Figure C1 The fitted values vs. residuals: q–q plots and histograms of residuals from the mixed-effect models (the upper row shows values for phosphorus, and the lower row shows values for nitrogen).

Figure C2 Soil total P content within plots in different soil layers.

Figure C3 Soil total N content and C:N ratio within plots in the F and O layers.

Table C1 Average total contents of elements (g kg⁻¹) in soil layers and their standard deviations in parentheses. All plots are included.

Table C2 Statistical differences of soil elements in each soil layer of different plots, grouped by their dominant tree species. Levels of significance: ^* = 0.05; ${}^{**}$ = 0.01; ${}^{***}$ = 0.001.

Table C3 Results from the mixed-effect models, testing the effects of environmental variables on soil total P and N content and C:N ratio. The tested variables were dominant tree species of the research plot, estimated tree age, rock parent material (geology) and soil layer. Random effect was related to plot number. Pseudo-R² was calculated based on Nakagawa and Schielzeth (2013), Johnson (2014), and Jaeger et al. (2016).

Table C4 Linear correlations of element contents of soil, needles, and leaves; number of species in different vegetation layers; and plot distance from Sokli with the non-metric multidimensional scaling ordination pattern. The group of “grasses and sedges” includes forb, grass, and sedge species and “d. shrubs and trees” includes dwarf shrubs and tree seedlings. Levels of significance: ^* = 0.1; ${}^{**}$ = 0.05; ${}^{***}$ = 0.01. The first seven rows are soil values.

Appendix D: Needle and leaf nutrient contents per plot

Table D1 Statistically significant differences between needle age group by species (C = youngest needles; C + 1 = 1-year-old needles; C + 2 = 2-year-old needles) and p<0.05 one-way analysis of variance, with Tukey's honest significant difference post hoc test.

Table D2 Statistically significant differences of needle nutrient contents between plots, p<0.05.