Modelling the coverage and annual variation in bilberry yield in Finland

The coverage of bilberry (Vaccinium myrtillus L.) was modelled as a function of site and stand characteristics using the permanent sample plots of the National Forest Inventory (NFI) (Model 1). The sample sites consisted of mineral soil forests as well as fells and peatland sites. Annual variation in the bilberry yield (Model 2) was analysed based on measurements over 2001–2014 in the permanent sample plots (so-called MASI plots) in various areas of Finland. We derived annual bilberry yield indices from the year effects of Model 2 and investigated whether these indices could be used to estimate annual variation in bilberry crops in Finland. The highest bilberry coverage was found in mesic heath forests and fell forests. On peatlands the coverage was, on average, lower than on mineral soil sites; the peatland sites with most bilberry coverage were meso-oligotrophic and oligotrophic spruce mires and oligotrophic pine mires. Our bilberry yield indices showed similar variation to those derived from the mean annual berry yields reported and calculated earlier using the MASI plots; the correlation between the indices was 0.795. This approach to calculating annual berry yield indices is a promising way for estimating total annual bilberry yields over a given period of time. Models 1 and 2 can be used in conjunction with the Miina et al.’s (2009) bilberry yield model when bilberry coverage, average annual yield and annual variation in the yield are to be predicted in forest planning.


Introduction
Bilberry (Vaccinium myrtillus L.) is one of the most abundant wild berry species in Finland (Salemaa 2000).It has adapted to a wide range of different site and land types in coniferous ecosystems and, as a result, is widely distributed across different parts of Europe and northern Asia (Ritchie 1956;Coudun and Gégout 2007).In North America, there are several species closely related to V. myrtillus (Ritchie 1956).
In Finland, bilberry is typical and abundant especially in conifer heath forests of medium fertility (Laakso et al. 1990;Salemaa 2000).It is also present and producing crops in many marginal forest types (e.g. fell forests), and on pristine and drained peatland sites (Sarasto 1961;Salemaa 2000).It has been estimated that Finnish forests and peatlands could produce an average of about 184 million kg of bilberries annually (Turtiainen et al. 2005(Turtiainen et al. , 2007)), however wild berry yields vary greatly from year to year due to e.g.frost, pollination success and variation in precipitation and temperature (Raatikainen 1993;Wallenius 1999).Turtiainen et al. (2011) calibrated the total bilberry yield in an average crop year for different crop years using the inventory data on wild berries (referred to as MASI data; Salo 1999) which was collected from permanent berry sample plots in 1997-2008.The calibration was a two-stage process.First, mean annual national bilberry yields (kg ha -1 ) were calculated from the numbers of bilberries counted annually on the MASI sample plots and site-specific figures for the fresh weight of a bilberry.The arithmetical mean berry yield was calculated on the basis of mean annual berry yields.Second, the minimum and maximum values from the series of mean annual berry yields for the study period (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) were divided by the arithmetical mean berry yield and multiplied by 184 million kg (i.e. total berry yield during an average crop year).Turtiainen et al. (2011) thus calculated that total annual bilberry yields in Finland would vary between approximately 90 and 310 million kg.
Picking wild berries is an important forest custom in the northern hemisphere, particularly in the Nordic countries and former territories of the USSR (Lee et al. 2007;Turtiainen and Nuutinen 2012).It is also so that today both the public and private forest landowners place a high value on the multiple-use aspects of forests and objectives other than wood production have got increasing weight in forestry decision-making (Kangas 1998;Rämö et al. 2013;Filyushkina et al. 2016).Consequently, production functions are needed for non-wood forest products and services so that the various aspects of forest use -both timber and non-timber -can be integrated into forest planning calculations (Borges et al. 2014).In these functions, site type and stand characteristics are the most reasonable predictors as they are known in forest planning systems.
A number of models for the species Vaccinium have been developed in different countries.For example, Coudun and Gégout (2007) used climatic and edaphic factors to help predicting the abundance of bilberry in French's forests.In Canada, models for lowbush blueberries (Hall et al. 1982) and Vaccinium vitis-idaea L. (Krebs et al. 2009) were developed to predict the yield using climatic variables.In Sweden, various ecosystem services (including bilberry production) were modelled using both stand characteristics and climatic and edaphic factors as predictors (Gamfeldt et al. 2013).However, none of the models mentioned above are suitable for forest planning calculations because of the predictors they employ.
In Finland, several models for the yields of bilberry and cowberry (V.vitis-idaea L.) have been developed (Ihalainen and Pukkala 2001;Ihalainen et al. 2002Ihalainen et al. , 2003Ihalainen et al. , 2005;;Miina et al. 2009;Turtiainen et al. 2013).All these models have been prepared for forest planning purposes.So far only Miina et al. (2009) have taken into account the annual variation in bilberry yields using the MASI sample plot measurements from 2001-2007.The models of Miina et al. (2009) predicted bilberry coverage as a function of site fertility and stand characteristics, and used these results to predict berry yield and annual variation in berry yield as a function of bilberry coverage and stand characteristics.The coverage model can be applied to mineral soil sites throughout the country but the yield model is based on regional data from North Karelia, Finland.
In this study, the work of Miina et al. (2009) was continued so that the model for the coverage of bilberry considers not only mineral soil sites but also peatland sites and fells (Model 1).Further, the annual variation in bilberry yield during 2001-2014 (Model 2) was modelled using a considerably longer time series of data from the MASI sample plots than in the earlier study by Miina et al. (2009).Finally, we investigated whether the annual yield indices derived from Model 2 could be used to estimate total annual bilberry crops in Finland for the 2001-2014 period.

Materials
Data were collected for the percentage coverage of bilberry on the permanent sample plots (300 m 2 ) of the Finnish National Forest Inventory (NFI) in 1995 (so-called PSP3000 data) (Finnish Forest Research Institute 1995).Miina et al. (2009) and Turtiainen et al. (2013) used the same data in their studies and have described the inventory method in detail.Sample plots representing seven different categories (a-g) were used in this study.The categories were initially developed for modelling the coverage of cowberry and were selected on the basis that cowberry occurs frequently and abundantly on all of them (Turtiainen et al. 2013).As bilberry and cowberry very often grow side by side, at the same sites, it seemed reasonable to use the same categories in this study.
Categories (a-e) pertained to forest land and were as follows: a) mineral soils, b) spruce mires -transforming phase, c) spruce mires -transformed phase, d) pine mires -transforming phase, and e) pine mires -transformed phase.In the case of poorly productive land and unproductive (waste) land, only fell forests (f) and summits (g) (i.e.site quality class VIII; Finnish Forest Research Institute 1995; Tomppo et al. 2011) were considered in the modelling (Supplementary file 1, available at http://dx.doi.org/10.14214/sf.1573).It is worth noting that when talking briefly about fells, both fell forests and summits are meant.
A total of 2515 sample plots from the PSP3000 data were included in the dataset for this study.For all sample plots bilberry coverage was estimated from four 2 m 2 quadrates; the averages were used as the sample plot-wise (stand-wise) estimates of abundance.Due to the heterogeneity of the sample plot, 11% of the plots were divided into two or three stands.As a result, 2801 stands were included in the analyses (Suppl.file 1).The stand characteristics and mean bilberry coverage were calculated for all stands and used in the modelling.A description of stand characteristics by categories is given by Turtiainen et al. (2013, p. 6).The data had a hierarchical structure because the sample plots were located in 983 clusters, which in turn were located in 367 municipalities and 14 forestry centre regions.
The annual variation in bilberry yield was analysed using MASI data for the 2001-2014 period (Salo 1999).The MASI inventory reports annual counts of ripe bilberries on five permanent 1 m 2 quadrates located in forest stands which had been found to be good growing sites for the species.The mean annual number of ripe bilberries and the mean bilberry coverage for these quadrates were used in modelling.The dataset included 306 observations of mean annual bilberry numbers in 50 stands (Supplementary file 2, available at http://dx.doi.org/10.14214/sf.1573).Stands were characterised in terms of site type, dominant tree species, stand age and basal area.The majority (60%) of stands in the dataset were dominated by pine (Pinus sylvestris L.), 32% were dominated by spruce (Picea abies (L.) Karst.) and the rest by deciduous tree species, mainly birches (Betula pendula Roth, B. pubescens Ehrh.).The stands were located in 29 municipalities and in 13 forestry centre regions.

Statistical modelling
Models were prepared for the mean percentage bilberry coverage (Model 1) and the mean number of bilberries in the stand (Model 2) using generalised linear mixed models (GLMMs) (Supplementary file 3, available at http://dx.doi.org/10.14214/sf.1573).Model 1 was expressed using a logit-link function with a binomial response, and Model 2 using a log-link function with a Poisson response.
In Model 1 all the continuous predictors (e.g.stand age and basal area) and categorical predictors (e.g.site quality classes in different categories a-e) had to be logical and statistically significant at the 0.05 level, with no systematic errors in residuals.The coverage was modelled and the categories were combined as described by Turtiainen et al. (2013) and hence categories (b) and (c) were combined and are later jointly described as 'spruce mires'; similarly categories (d) and (e) were combined and are jointly described as 'pine mires'.
In Model 2, bilberry coverage, altitude and mean effective temperature sum were used as potential fixed predictors.The year effects in 2001-2014 were considered using indicator variables as fixed predictors.
We adjusted for the multi-level hierarchy and unbalanced structure of the datasets by including random effects at different levels.We adjusted for over-dispersion in the response variables by adding an additional random term at the bottom level ('pseudo' level).The GLMMs were estimated using the glmmPQL function of the R software (R Development Core Team 2012).

Simulations
The performance of Model 1 was illustrated by predicting the bilberry coverage in various stands whose initial characteristics, development and management were simulated using the Motti stand simulator (Hynynen et al. 2005).The simulated stands were located in southern Finland, with the temperature sum and altitude set at 1200 dd. and 100 m, respectively.We used the same pine and spruce stands representing several site types on mineral soils as Miina et al. (2009), and the same pine stands on pine mires as Turtiainen et al. (2013).We also simulated spruce stands on spruce mires (sites II and III) using the Motti simulator.

Berry yield indices
Bilberry yield indices were calculated from the year effects in Model 2. The year effects were transformed into exponential form, multiplied by one hundred and scaled to have a mean of 100.The indices obtained in this study were compared with berry yield indices derived from the mean annual berry yields presented by Salo (2015).Salo (2015) calculated a series of mean annual bilberry yields for 1997-2013 using MASI data and the method described by Turtiainen et al. (2011).For the purpose of comparison we also scaled Salo's (2015) mean annual berry yields to have a mean of 100.

Model for bilberry coverage
In the case of mineral soils, Model 1 for bilberry coverage produced results very similar to those described by Miina et al. (2009); the highest coverage was found on mesic heath forests (site III, Table 1).On sub-xeric heath forests (site IV) and herb-rich heath forests (site II), the coverage was 63% and 57% of that on site III respectively.On rocky and sandy soils (VII) and xeric heath forests (V) the coverage was considerably lower, and on herb-rich forests (I), it was scarce.
On mineral soils bilberry coverage increased with stand age and stand basal area up to the age of 175 years and a density of 26 m 2 ha -1 , after which coverage gradually decreased.It is worth noting that stand basal area was a significant predictor not only on mineral soils but also on spruce  and pine mires.On mineral soils and spruce mires, the dominant tree species was also a significant predictor of bilberry abundance.The coverage was lower in spruce-dominated mineral soil forests than pine-dominated forests.The negative effect of deciduous trees on bilberry coverage applied on mineral soils (sites I and II) and spruce mires.Meso-oligotrophic (III) and oligotrophic spruce mires (IV) did not differ significantly from mesic heath forests with respect to bilberry abundance (Table 1).Coverage at the other sites (I-II) on spruce mires was, on average, one third of the coverage on site III, which is on mineral soils (the reference in Model 1).Of the pine mire sites, site IV had the highest bilberry coverage, although not as high as site III (on mineral soils) and sites III and IV (on spruce mires).At the sites on meso-oligotrophic and more fertile pine mires (I-III), the coverage was two thirds of that on site III (mineral soils), whilst on poor oligo-ombrotrophic pine mires (V), it remained considerably lower.
Model 1 indicated that coverage on poorly productive land, i.e. on fell forests at high altitude, was higher than on site III, on mineral soils (Table 1).In addition altitude had a positive effect on bilberry coverage in all categories (a-g), suggesting that coverage is, on average, higher in eastern and northern Finland than in western and southern Finland.However, both altitude and stand age and basal area are factors in the largescale variation in bilberry coverage in Model 1.
In the simulated stands the effect of site fertility on predicted coverage was clear and the coverage increased with stand age and basal area (Fig. 1).In most cases coverage was temporarily affected by thinnings, but always decreased significantly following regeneration felling.On average, coverage was higher on mineral soils than on pine and spruce mires.

Annual variation in bilberry yields
During the study period of 2001-2014, predicted year effect in Model 2 for the number of bilberries was highest in 2005 and 2012, and lowest in 2002 and 2010 (Fig. 2, Table 2).The sample mean and variance for the fixed year effects of Model 2 were 0.1374 and 0.0367 respectively.2) and mean annual bilberry yields (kg ha -1 ) presented by Salo (2015).The former index series covers years 2001-2014 and the latter one years 1997-2013.In both series, the mean value of the indices is 100.
The bilberry yield indices in this study show similar variation to those derived from Salo's (2015) study; the correlation between the indices in the two studies was 0.795.According to Salo (2015), the best and worst bilberry crops were obtained in 2012 and 2004 respectively.For the 2001-2013 period there is agreement between the index series on the above-and below-average years, with the exception of 2009.However, the range of the indices calculated in this study was smaller than that in Salo's (2015) study (69-150 and 50-180 respectively).

Discussion
This study represents a continuation of work by Miina et al. (2009) and Turtiainen et al. (2011).Our model of bilberry coverage (Model 1) considered a wider range of sites, including other sites where bilberry occurs in the field layer (categories b-g) as well as the mineral soil forests (category a) considered by Miina et al. (2009).Covering a wider range of sites was needed, because one third of forestry land in Finland is growing on peatlands (Natural Resources Institute Finland (Luke) 2015) which are not considered by the earlier model of Miina et al. (2009).Model 1 was prepared for forest planning purposes and it can be used together with the yield model of Miina et al. (2009: Table 3).a The number of observations at each level is given in parentheses.A random term at 'pseudo' level accounts for the overdispersion.
In the case of mineral soil forests (category a), our results on bilberry coverage were very similar to those of Miina et al. (2009).Miina et al. (2009) considered their results in the context of earlier empirical findings and concluded that they were consistent with these findings (see also Gamfeldt et al. 2013;Hedwall et al. 2013).It follows that our results are also consistent with the empirical findings.
Bilberry coverage on peatlands was, on average, lower than on mineral soil sites (Salemaa 2000).Of the peatland sites, sites III and IV on spruce mires and site IV on pine mires had the highest bilberry coverage.These findings are in line with earlier studies (Sarasto 1961;Salemaa 2000).According to Model 1, deciduous trees decreased bilberry coverage on both spruce mires and the most fertile mineral soil sites.To the best of our knowledge this is the first study to consider the effects of tree species on bilberry coverage on peatlands.
Bilberry was highly abundant on fell forests, and altitude also had a positive effect on bilberry coverage.In northern Finland bilberry grows on poorer sites than in southern Finland owing to the higher relative humidity and higher soil moisture in the north.Bilberry is widespread on the lower slopes of fells and mountains, but it does not occur at such high altitudes as Vaccinium vitis-idaea L. or Empetrum nigrum L. (Salemaa 2000).
Our modelling of bilberry coverage was based on an extensive dataset collected from the permanent sample plots of the Finnish NFI in 1995 (Suppl.file 1).The MASI dataset was considerably smaller (Suppl.file 2) and comprised mainly data from mature stands of medium site fertility which had been found to be good for bilberry.Therefore, it was unfeasible to fit a new yield model in this study, and, although it is based on regionally limited data, we recommend the bilberry yield model of Miina et al. (2009) for future use (Kilpeläinen et al. 2016).We used the MASI data from 2001-2014 to model annual variation in bilberry yields (Model 2).Our Model 2 is based on data from twice as many years as the model by Miina et al. (2009) and we therefore recommend that it is used in conjunction with the Miina et al.'s (2009) bilberry yield model when both average annual yield and annual variation in the yield are to be predicted.
In this study, we first derived annual bilberry yield indices from the year effects of Model 2 and then investigated whether these indices could be used to estimate annual variation in bilberry crops in Finland.The bilberry yield indices varied quite similarly with those of Salo (2015) (Fig. 2).There are several possible reasons for the differences between these two index series.Firstly, we used, on average, fewer MASI stands per year than Salo (2015) (22 vs. 69 stands per year).We used only those MASI stands for which data on stand characteristics and bilberry coverage were available in addition to data on bilberry numbers and site fertility.Secondly, by using a mixed modelling approach we were able to adjust for the unbalanced, cross-correlated structure of the MASI data.Data were not available for all sample plots for every year and observations from a given year were cross-correlated due to annual variation in conditions.Thirdly, Turtiainen et al. (2011) and Salo (2015) calculated mean annual bilberry yields (kg ha -1 ) using the method described by Turtiainen et al. (2011), which used average weights of ripe bilberries on different site fertilities.The weights were obtained from earlier studies by Kuchko (1988) and Ihalainen et al. (2003).However, one can question whether the weights used were accurate enough (cf.Eronen 2004).It should also be noted that Turtiainen et al. (2011) did not consider annual variation in bilberry weight.In this study we derived bilberry yield indices directly from the number of bilberries without converting berry counts into fresh weights.
In conclusion, the modelling approach introduced in this study is a promising way for estimating total annual bilberry yields over a given period of time.The approach can also be applied to cowberry.However, models of annual variation in berry yield should be based on more comprehensive data than were used in this study.Access to a comprehensive dataset is required to determine how to adjust berry yields for an average crop year to predict total yield in a given year.The results of this study suggest that a very simple calibration method will suffice: the berry yield index for the year in question is divided by 100 and multiplied by the total yield for an average crop year (cf.Turtiainen et al. 2011: Eq. 2); however, this must be verified in future studies.

Fig. 1 .
Fig. 1.Predicted coverage of bilberry in pine and spruce stands of different site fertilities (i.e.sites I-V; see the definitions in Suppl.file 1).The development of stands, representing mineral soils (A, B), pine mires (C) and spruce mires (D) in southern Finland, was simulated using the Motti simulator (arrows indicate thinnings).

Fig. 2 .
Fig. 2. Bilberry yield indices calculated on the basis of Model 2 (Table2) and mean annual bilberry yields (kg ha -1 ) presented bySalo (2015).The former index series covers years 2001-2014 and the latter one years 1997-2013.In both series, the mean value of the indices is 100.

Table 1 .
The multi-level binomial model (Model 1) estimated for the mean percentage coverage of bilberry on the 2 m 2 quadrates in the stands of the PSP3000 data.Sites I-V and VII-VIII refer to different site quality classes (see Suppl.file 1).Mineral soils, spruce mires and pine mires pertain to forest land, i.e. categories (a-e) of this study.The number of observations at each level is given in parentheses.A random term at 'pseudo' level accounts for the overdispersion.
a The parameter estimates of site variables "site III, spruce mires", "site IV, spruce mires" and "site VIII, waste land" were not statistically significant.b ArtificialRegen (artificial regeneration) is an indicator variable for the regeneration method (ref.natural regeneration).c FormerAgrLand (former agricultural land) is an indicator variable for stand history (ref.former forest).d An indicator variable for the dominant tree species (the reference is other tree species).e In this context, forest land refers to categories (a-e) of this study.f

Table 2 .
The multi-level Poisson model (Model 2) estimated for the mean number of bilberries on five 1 m 2 quadrates inMASI stands, measured in 2001MASI stands, measured in  -2014.   .