Modelling the Abundance and Temporal Variation in the Production of Bilberry ( Vaccinium myrtillus L . ) in Finnish Mineral Soil Forests

Empirical models for the abundance and berry yield of V. myrtillus were constructed using generalized linear mixed model (GLMM) techniques. The percentage coverage of bilberry was predicted as a function of site and stand characteristics using the permanent sample plots of the National Forest Inventory (NFI) in 1995. The number of bilberries was predicted as a function of percentage coverage and stand characteristics using a sub-sample of the NFI plots in North Karelia. The between-year variation in the bilberry yield was analysed using the permanent experimental plots (MASI) established in different areas of Finland and measured in 2001–2007. The highest coverage of bilberry was found on mesic heath sites; on sub-xeric and herb-rich heath sites the values were 62% of that for mesic sites. The decreasing effect of deciduous trees (compared to spruce) was significant only on herb-rich heath sites. The coverage increased along with stand development up to certain stand ages and basal areas. The bilberry yields were higher in pine-dominated stands than in spruce-dominated ones. In spruce stands, the coverage of bilberry and stand basal area significantly affected the number of berries, whereas in pine stands only the coverage was a significant predictor. In the MASI data, the bilberry yield of pine stands was two times higher than that of spruce stands; however, the between-year variation in bilberry yield was higher in the spruce than in the pine stands. The estimated models were used to predict the bilberry coverage and yield along with stand development. On mesic heath sites in southern Finland (1200 dd.), the predicted annual yield of bilberry was about 25 kg ha–1 (95% confidence interval 9–73 kg ha–1) in a mature pine stand and about 10 kg ha–1 (3–35 kg ha–1) in a mature spruce stand. The models can be included in stand simulators, where they would facilitate the prediction of bilberry abundance and yields for silvicultural and forest planning purposes.


Introduction
Vaccinium myrtillus L. (bilberry, blueberry) is one of the most frequent and abundant vascular plant species in northern Europe, and is an important species in the understorey vegetation of conifer ecosystems (e.g.Ritchie 1956, Flower-Ellis 1971, Kardell 1980, Kuusipalo 1983, Hulten and Fries 1986, Storch 1993, Tolvanen 1995, Atlegrim and Sjöberg 1996, Mäkipää 1999, Kurki et al. 2000, Salemaa 2000, Coudun and Gégout 2007, Nielsen et al. 2007).In Finland, V. myrtillus is typical and abundant especially in spruce-and pine-dominated heath forests of medium fertility (Laakso et al. 1990, Salemaa 2000).Bilberry also occurs and produces yields in many marginal types of forest, and on pristine and drained peatland sites (Ruuhijärvi 1960, Sarasto 1961, Korpela and Reinikainen 1996, Salemaa 2000, Turtiainen et al. 2007, Salo 2008a).It has been noted that the coverage of V. myrtillus has decreased, on average, in Finnish mineral soil forests due to intensive forestry with clear-cuttings and soil preparation; the proportion of young forests has also increased, which has led to a more even distribution of the different age classes and had a negative effect on the abundance of V. myrtillus.On drained peatlands, the coverage of bilberry is increasing along with the post-drainage succession phases (i.e. recently drained, transforming phases, transformed phase) (Salemaa 2000).
As a food source and cover, bilberry is the major determinant of the selection of habitat by capercaillie (Tetrao urogallus) (e.g.Storch 1993).In addition, bilberry hosts several herbivorous invertebrate species that are an important food source for chicks (Lakka and Kouki 2009).The reduction in the average coverage of bilberry since the 1950s (Salemaa 2000) coincides with the beginning of intensive forestry in Finland as well as the decline of capercallie populations (Helle et al. 2003).Capercaillie as an old-growth forest specialist indicates also the well-being of many other forest dwelling species.Therefore, models for the abundance of bilberry in Finnish forests could be used in forest planning to enable more efficient game habitat management.
Many studies have been carried out on the fruit production of bilberry due to the large geographical area and wide range of habitats occupied by the species, as well as to the large number of factors (e.g.frost, precipitation, temperature, pollination success, drought, site type, tree stand structure) that influence the yield (e.g.Eriksson et al. 1979, Kardell 1980, Raatikainen and Raatikainen 1983, Raatikainen et al. 1984, Jäppinen et al. 1986, Kuchko 1988, Raatikainen and Vänninen 1988, Laakso et al. 1990, Belonogova 1993, Salo 1995and references therein, Wallenius 1999, Turtiainen et al. 2007).
Bilberry is one of the economically most important wild berry species in Finland, Sweden and Norway (e.g.Kardell 1980, Salo 1995, Saastamoinen et al. 2000, Kangas 2001).Bilberry is widely collected for both household consumption and sale.The picking of wild berries is of great economic significance to many regions in Finland due to the right of public access to all forest land (e.g.Kangas 2001).In addition, the income earned from the sale of e.g.wild berries, herbs and mushrooms is tax-free in Finland (Salo 1995(Salo , 2008b)).It has been estimated that the total bilberry crop in Finnish forests in a good year may be about 200 million kg, but in a poor year considerably lower (Turtiainen et al. 2005, Salo 2008b).In recent years, Finns have picked a total of 40-50 million kg lingonberries (Vaccinium vitis-idaea) and bilberries; i.e. about 10% of the total crop (Salo 2008b).The total value of the harvested wild berry crop was estimated at 115 million euros in 1997, and 95 million euros in 1998 (Saastamoinen et al. 2000).In Finland, the value of wild berries picked for domestic use, commercial sale, and direct and market sale in 2005 was 77.2 million euros (Salo 2008b).
The public and private forest landowners, both industrial and non-industrial, 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).In addition to timber production, forest uses include berry and mushroom picking, other forms of gathering (decorative lichen, herbs and other plants), Christmas tree harvesting, game management and hunting, reindeer husbandry, recreation and outdoor activities, landscape management and conservation (Salo 1995).Aspects such as the vitality of forests, scenic beauty and biodiversity are also regarded as important values of the forests in Finland (Kangas and Niemeläinen 1996).
Production functions suitable for planning calculations are needed when integrating non-wood forest goods and services in numerical forest planning.Some production functions have already been formulated for bilberry in Finland on the basis of expert or visual field estimates (Ihalainen and Pukkala 2001, Ihalainen et al. 2002, 2005).The empirical berry yield models developed so far have been based only on regional yield data from North Karelia (Ihalainen et al. 2003).Ihalainen et al. (2003) used stand age and site type as the predictors of the annual bilberry yield collected in 1981-1984. Furthermore, Turtiainen et al. (2005, 2007) calibrated the expert models of Ihalainen et al. (2005) with a set of measured yield data from various sources, and used the models to estimate preliminary national and regional yields of bilberry and lingonberry.However, the existing models do not consider the abundance (coverage) of bilberry as a predictor of bilberry yields.In addition, the annual variation in bilberry yields has not been included in yield predictions.
The aim of this study was to develop a model set, which enables predicting both the abundance of bilberry and bilberry yields including confidence intervals in mineral soil forests in Finland.The empirical models were designed to be included in stand simulators for silvicultural and forest planning purposes.Thus, only those site and stand characteristics were used as predictors that are usually known in forest planning systems.The performance of the estimated models was illustrated by simulating the development of the bilberry coverage in accordance with the development of Scots pine and Norway spruce stands in southern and northern Finland.The datasets on the bilberry coverage and the annual variation in bilberry yields covered the whole of Finland.Because the relationship between the coverage of bilberry and the number of berries was modelled using a dataset collected from North Karelia, bilberry yields were simulated only for southern Finland.Suggestions are given to improve the modelling datasets in order to obtain empirical yield models for bilberry in the whole of Finland.

Sampling and Data on the Percentage Coverage of Bilberry
The study material consists of the measurements of site and stand characteristics and the coverage of plant species carried out on the permanent sample plots (300 m 2 ) of the Finnish National Forest Inventory (NFI) in 1995 (PSP3000 data).
The inventory in southern Finland was based on clusters of four permanent sample plots, arranged in a north-south direction, located systematically at intervals of 400 m.In northern Finland the corresponding layout was three plots at 600 m distance.The distance between the clusters in both the north-south and east-west direction was 16 km, and in northern Finland 24 and 32 km, respectively (Pysyvien koealojen… 1995, see also e.g.Kühlmann-Berenzon and Hjorth 2007).
The sample plots located on mineral soil sites (i.e.herb-rich forests and heath forests) on forest land (Tomppo et al. 1997) were used in the study.On mineral soil sites, the organic layer on top of mineral soil is not peat.Herb-rich forests (groves) have a fertile, brown soil with a mull layer, whereas heath forests have podzols with a grey leached horizon (e.g.Hotanen et al. 2008).The fertility of the sites was determined using the following site-type classification: site quality class I = herb-rich forests (groves) (1.7% of stands), II = herb-rich heath forests (Oxalis-Myrtillus type, OMT group) (19.8%),III = mesic heath forests (Myrtillus type, MT group) (46.5%),IV = subxeric heath forests (Vaccinium type, VT group) (27.8%) and V = xeric heath forests (Calluna type, CT group) (4.1%) (Frey 1973, Hotanen et al. 2006).Site quality type VI (barren heath forests, Cladonia type ClT group) was omitted because it contained only a few sample plots, and it is too infertile and dry for bilberry.Poorly productive land and waste land, i.e. exposed bedrock, cliff or sandy forestry land (class VII), as well as timber line forests including fell tops (VIII), were also excluded.
A total of 1888 sample plots were included in the study material.Most of these plots were completely located within one same stand, but 166 sample plots were divided into two stands and one sample plot into three stands due to the  1).The data had a hierarchical structure because the sample plots were located in 880 sample plot clusters that were located in 356 municipalities and 14 forestry centre regions.The altitude and mean effective temperature sum (threshold +5 °C) of the sample plot were used to describe the south-north variation in the growing conditions at the national level (Ojansuu and Henttonen 1983).The stand characteristics were measured and recorded according to the field instructions of the NFI (Pysyvien koealojen… 1995, see also e.g.Hotanen et al. 2006).On each sample plot, the projection coverage of individual species (including V. myrtillus) was, in most cases, estimated on four 2-m 2 sample quadrates (Reinikainen and Nousiainen 1995, Kühlmann et al. 2001), and the average values were used as the sample plotwise estimates of abundance.If the sample plot was bisected by a stand compartment boundary (total of 167 cases), the stand characteristics and species coverages were estimated separately for each stand (also Kühlmann-Berenzon and Hjorth 2007) and used in the modelling.As a result, the mean number of quadrates per stand in the whole material was 3.5.About 70% of the stands were naturally regen-erated, the remainder being planted or direct seeded.A small proportion of the stands (2.1%) were located on former agricultural land.Of all stands, 59% were dominated (according to the stand volume of the individual tree species) by Scots pine (Pinus sylvestris), 31% by Norway spruce (Picea abies) and 8% by deciduous trees (mainly birches).The rest of the stands were recently clear-cut, and the dominant tree species was not determined.During the last 10 years, 79% of the stands had not been treated by cuttings.Pre-commercial thinning, thinning, regeneration cutting and seed-tree removal had been carried out in 5, 8, 6 and 2% of the stands, respectively.

Bilberry Yield Data
In June 2006, the number of unripe berries was recorded in a total of 32 stands in North Karelia.Stands (sample plots) that were, on the basis of the measured site and stand characteristics (NFI10_NK data), potential growing sites for bilberry were selected from the permanent sample plots of the 10th NFI (Valtakunnan metsien… 2005).If the sample plot was divided into two or more stands, only one stand per a sample plot was selected.In each stand, the number of unripe berries and the coverage of bilberry were recorded on four 1-m 2 quadrates (Table 1).The mean number of bilberries and the mean coverage of bilberry on the quadrates were used in the analyses.The site and stand characteristics were obtained from the database of the 10th NFI (Valtakunnan met-sien… 2005).Only pine-and spruce-dominated stands were selected for the study material (18 pine stands and 14 spruce stands).Of 32 stands, 6 stands were growing on herb-rich (II), 20 stands on mesic (III) and 6 stands on sub-xeric heath (IV) sites.About half of the stands had been artificially regenerated.
The annual variation in the yield of bilberry was studied during 2001-2007 using the so-called MASI data collected from permanent sample plots in different areas of Finland (Salo 1999(Salo , 2005)).The Finnish name MASI literally means berry (MArja) and mushroom (SIeni) information system, and it consists of a number of different databases.During the years 2001-2002 the MASI network consisted of 250 stands, and after 2002 of 220 stands (Salo 2005).Bilberry yield has been recorded in stands that were considered to be potential growing sites for bilberry.In each stand, the number of bilberry ripe berries has been counted annually on five 1-m 2 quadrates.The mean number of bilberries on the quadrates was used in the analyses.Stands on mesic heath sites (III) were included in this study.In addition, there had to be at least two observations per stand, i.e. the annual number of berries on the stand had to be recorded at least twice during 2001-2007.This resulted in a total of 51 stands with a total of 221 observations on the annual number of berries (Table 1) (on average, 4.3 annual yields per stand during the seven-year period).Half of the stands were pine-dominated and half sprucedominated stands.The development classes of the stands were as follows: young thinning stands 14%, advanced thinning stands 63%, and mature stands 23%.

Statistical Modelling
Models were prepared for bilberry coverage (Model 1), bilberry yield (Model 2) and the annual variation in bilberry yield (Model 3) using the PSP3000, NFI10_NK and MASI datasets, respectively.The mean percentage coverage of bilberry and the mean number of berries in the stand were used as the response variables and modelled using generalized linear mixed models (GLMMs) with a binomial response (logit-link function) and a Poisson response (log-link function), respectively (e.g.McCullagh and Nelder 1989).The response variables were rounded to the nearest integer due to discrete (binomial and Poisson) distributions.The multilevel hierarchy of the data (stand, sample plot, sample plot cluster, municipality, forestry centre region), and subsequently correlated observations, was taken into account by including random effects at different levels in the variance component models, and by allowing the intercept to vary randomly across the levels (e.g.Searle et al. 1992, Snijders and Bosker 1999, Goldstein 2003).As suggested by Browne et al. (2005), overdispersion in the GLMMs was taken into account by adding a random term as being at the bottom level ("pseudo" level).The GLMMs were estimated with the pseudo quasi-likelihood (PQL) method using the glmmPQL function of the R software (R Development Core Team 2007).
The percentage coverage of bilberry was treated as a proportion, and the general multi-level binomial model was assumed as follows (cf.Browne et al. 2005): where y is the observed coverage of bilberry in the stand (i.e. the mean percentage coverage of 2-m 2 quadrates in the stand); binomial(n, p) denotes the binomial distribution with parameters n (binomial sample size, in our case n ijklm are all equal to 100) and p (proportion of successes i.e. the expected coverage of bilberry); ln(p/(1p)) is a logit-link function; and f(•) is a linear function with arguments X ijklm (i.e.fixed predictors) and β (i.e.fixed parameters).Subscripts i, j, k, l and m refer to forestry centre region, municipality, sample plot cluster, sample plot and stand, respectively.u i , u ij , u ijk , u ijkl and u ijklm are random, normally distributed between-forestry centre region, between-municipality, between-sample plot cluster, between-sample plot and betweenstand effects with a mean of 0 and constant variances.Random terms at different hierarchical levels were assumed to be uncorrelated.
The general multi-level Poisson model was as follows: ) ) ( ) where, in the NFI10_NK data y is the mean number of unripe berries on four 1-m 2 quadrates in the stand, and in the MASI data the mean number of ripe berries on five 1-m 2 quadrates in the stand; the conditional distribution of y, given the expected value π, is the Poisson distribution; ln(π) is a log-link function; and f(•) is a linear function with arguments X ijklt (i.e.fixed predictors) and β (i.e.fixed parameters).Subscripts i, j, k, l and t refer to forestry centre region, municipality, sample plot cluster, sample plot (or stand) and year, respectively.u i , u ij , u ijk and u ijkl are random, normally distributed between-forestry centre region, between-municipality, betweensample plot cluster and between-stand with a mean of 0 and constant variances.Because the glmmPQL function did not allow the random cross-effects, the effects of the years 2001-2007 were included using dummy variables in fixed predictors.Random terms at different hierarchical levels were assumed to be uncorrelated.The temperature sum, altitude, stand age, stand basal area and bilberry coverage were continuous fixed predictors considered in the models.The squared transformation of the predictor was used to describe the possible non-linear relationship between the response and predictor variables.The categorical predictors (site type, dominant tree species, regeneration method and history of the stand) were included using dummy variables so that the most common situation was selected as the basis of the models (mesic heath forest, Norway spruce, natural regeneration and forest land, respectively).In Eq. 2 (MASI data), the reference year of the dummy variables for the fixed year effects was 2006, since the NFI10_NK dataset was collected in 2006 (see 2.4 Simulations).All the fixed predictors included in the models had to be logical and significant at the 0.05 level, and no systematic errors were observed in residuals.

Simulations
The performance of the estimated models was illustrated by predicting the percentage coverage of bilberry and the bilberry yield and its annual variation in both Scots pine and Norway spruce stands.The development of the coverage of bilberry was presented as a function of stand characteristics for several site types in both northern and southern Finland.However, because the data on the number of unripe berries mainly covered mesic heath (III) sites in North Karelia, the bilberry yield and its annual variation were presented only for mesic heath sites in southern Finland (1200 dd.).
Stand development was simulated using the Motti stand simulator (Hynynen et al. 2002, Salminen et al. 2005).Simulations were con-ducted for northern and southern Finland, with the temperature sum set to 900 dd. and 1200 dd., respectively, and the altitude to 100 m.The development of pine stands was simulated on mesic (III), sub-xeric (IV) and xeric (V) heath sites, and the development of spruce stands on herb-rich forest (I), on herb-rich (II) and mesic (III) heath sites, which are typical growing sites for pine and spruce.In the Motti tree-level growth simulator, the initial stands on different site types and at different geographical locations were generated with a biological age and stem number of 10 years and 2000 stems ha -1 , respectively.The stands were assumed to be naturally regenerated.In the simulations, stand management followed the silvicultural recommendations applied in privately owned forests in Finland (Hyvän metsän-hoidon…2006).Depending on the tree species and site type, the stands were thinned 2-3 times before regeneration felling at the age of 59-83 years in southern Finland and 91-137 years in northern Finland.
The percentage coverage of bilberry was predicted as a function of the site and stand characteristics (Eq. 1, Model 1, PSP3000 data), and the number of unripe bilberry berries was predicted as a function of the coverage of bilberry and site and stand characteristics (Eq.2, Model 2, NFI10_NK data).The number of ripe berries was calculated by assuming that 80% of the unique berries become ripe.The ratio between ripe and unripe berries has been calculated on the basis of the MASI data collected in 2000 (Ihalainen et al. 2003).The number of ripe berries was converted into the bilberry yield (kg ha -1 ) by multiplying the mean fresh weight of one ripe berry.Using a sub-sample of the MASI data from North Karelia, Eronen (2004) calculated that the mean fresh weight of one ripe berry in mesic heath forests is 0.35 g, and this value was used in our study.
Model 2 was based on the observations made in 2006, and thus did not take into account the annual variation in the bilberry yield.The year 2006 was also used as the reference year in the model for the annual variation in bilberry yield (Eq. 2, Model 3, MASI data).Therefore, the average annual yield of bilberry was predicted using the sample mean of the fixed year effects of Model 3 as a year effect in Model 2. The variance of the year effects was estimated from the sample variance of the fixed year effects by subtracting the additional variance caused by the estimation errors (see e.g.Miina 1993).The lower and upper bounds of the 95% confidence interval for the average annual yield of bilberry were calculated using the 0.025 and 0.975 quantiles of the normal distribution and the standard deviation of the year effects, respectively.

Model for the Percentage Coverage of Bilberry
The site type was a significant predictor in Model 1, estimated for the percentage coverage of bilberry using the PSP3000 data covering the whole of Finland (Table 2).The highest coverage was found on mesic heath (III) sites.Compared to III sites, the coverage of bilberry decreased in the following order: sub-xeric heath (IV) (62% of that on III), herb-rich heath (II) (62%), xeric heath (V) (22%) and herb-rich forest (I) (11%) (cf.Odds ratio, Table 2).The coverage of bilberry in direct seeded and planted stands was lower than that in naturally regenerated stands.If the stand had been afforested (i.e.former agricultural land), the coverage of bilberry was lower.
The stand characteristics significantly affected the coverage of bilberry.The coverage was the higher, the higher the age and basal area of the stand up to the age of 191 years and a density of 25 m 2 ha -1 , after which the coverage gradually decreased.Because dominant tree species was an important predictor of bilberry yields, it was used also as a predictor of the abundance of bilberry though it would not be significant at the 0.05 level.In pine-dominated stands, the coverage of bilberry was 1.13 times higher (p = 0.07) than that in stands dominated by spruce or deciduous trees (Table 2).The decreasing effect of deciduous trees (compared to spruce) was significant (p = 0.02) only on herb-rich heath sites (II).Interaction terms among the stand and site characteristics were nonsignificant in the study material (p > 0.1).
The temperature sum was not a significant predictor but, in contrast, altitude had a positive effect on the coverage of bilberry; i.e. the cover-age was, on the average, higher in east and north Finland than in south and west Finland.Due to the fact that the stands are older and sparser in the north than in the south, stand age and basal area also explain the large-scale geographical variation in the coverage of bilberry in Finland.Thus, the effects of temperature sum and geographical variation were partly reflected in the modelling through altitude and stand characteristics.
Most of the unexplained variation (46%) of the model for the coverage of bilberry was found at the stand level.The sample plot-, sample plot cluster-and municipality-level residual variations almost equally accounted for the residual variation (15-18%).The forest centre region-level residual variation accounted for only 2% of the total residual variation.

Models for Bilberry Yield
In Model 2, estimated using the NFI10_NK data measured in 2006 from North Karelia, the site type was not a significant predictor for the number of unripe bilberries (Table 3).In the spruce stands, the number of unripe berries was predicted as a function of the percentage coverage of bilberry and stand basal area, whereas in the pine stands the coverage of bilberry was the only significant predictor.The number of unripe berries increased when the coverage of bilberry increased; the maximum was achieved at a coverage of 43% in pine stands and 38% in spruce stands.In the PSP3000 and NFI10_NK data, the maximum coverage of bilberry was 71% and 54%, respectively (Table 1).The number of unripe berries was the higher, the higher basal area of the spruce stands a) The number of observations at the level is given in parentheses.A random term at "pseudo" level accounts for the overdispersion.
up to a density of 14 m 2 ha -1 , after which the berry yield strongly decreased with increasing stand basal area.The between-year variation in the number of ripe berries was studied using the MASI data measured in 2001-2007 in different parts of Finland (Model 3).Because the NFI10_NK data on the number of unripe berries were recorded in 2006 and the predictions of the count models had to be comparable, year 2006 was selected as the reference year (category) in the model for the number of ripe berries (Table 4).Also in the MASI data, the tree species had a significant effect on the yield of bilberry, and thus separate models were estimated for pine-and sprucedominated stands.According to the models, the number of ripe berries was low in 2002, 2004 and 2006, and high in 2001, 2003, 2005 and 2007.The sample means of the fixed year effects of Model 3 were 0.1422 for pine stands and 0.5450 for spruce stands.These values were used as the year effects when the average annual yield of bilberry was calculated using Model 3. The variances of the year effects of Model 3 were 0.2907 for pine stands and 0.4504 for spruce stands.According to the models estimated using the MASI data, the bilberry yield in pine stands has been two times higher than that in spruce stands.

Simulation Results
The performance of the estimated models was investigated by predicting the coverage and yield of bilberry in pine and spruce stands whose development and management were simulated using the Motti simulator.The predicted coverage of bilberry followed the development of the stand basal area and age (Fig. 1).The coverage increased with increasing stand age and basal area, and was temporally affected by thinnings.The site type had a clear effect on the coverage of bilberry, but there were no major differences between pine and spruce stands on the same site type (i.e.mesic site).The coverage was slightly higher in the stands in northern Finland (900 dd.) than in southern Finland (1200 dd.) due to the lower stand density in northern Finland.
The average annual yield of bilberry, as well as the 95% confidence interval for the between-year  a) The number of observations at the level is given in parentheses.A random term at "pseudo" level accounts for the overdispersion.
10 kg ha -1 in a spruce stand.The 95% confidence interval of these figures was 9-73 kg ha -1 in a mature pine stand and 3-35 kg ha -1 in a mature spruce stand.
The predicted bilberry yields were compared to the yields computed using the yield model and the correction factor 0.28 presented by Turtiainen et al. (2005: Eq. 2 and Table 4).The comparison was made in pine and spruce stands on mesic heath sites in southern Finland (1200 dd.) so that our models were replaced by the model of Turtiainen et al.In a pine stand, the patterns of the predicted bilberry yields were corresponded very closely to each other.However, compared to the model of Turtiainen et al. our models gave higher yields in a pine stand and lower yields in a spruce stand, especially at the end of the rotation.variation in the yield, were simulated in pine and spruce stands on mesic heath sites in southern Finland (1200 dd.).The annual bilberry yields increased with increasing stand age and reached the level of 20 kg ha -1 after the first thinning in a pine stand and 10 kg ha -1 in a spruce stand (Fig. 2).In a mature spruce stand, the bilberry yield decreased heavily when the stand basal area was above the optimal level, and increased temporarily when the stand was thinned at the ages of 43 and 58 years.In a pine stand, stand basal area affected the bilberry yield only slightly through its effect on the percentage coverage of bilberry (i.e.stand basal area had no direct effect on bilberry yield in pine stands).Five years before the end of rotation, the annual yield of bilberry was about 25 kg ha -1 in a pine stand and about a) The number of observations at the level is given in parentheses.A random term at "pseudo" level accounts for the overdispersion.

Abundance
The site type significantly affected the coverage of bilberry.The highest coverage was found on mesic heath sites, i.e. in the Myrtillus type (MT) group.The result is similar to that reported in other studies (e.g.Kalela 1970, Raatikainen et al. 1984, Jäppinen et al. 1986, Salemaa 2000, Tonteri et al. 2005).Bilberry is also abundant in sub-xeric and herb-rich heath stands; the coverage values were 62% of that for mesic heath stands.On sub-xeric, and especially on xeric sites, bilberry probably suffers slightly from a shortage of water, whereas in fertile herb-rich heath and especially herb-rich forests, competition from grasses and herbs limits its abundance (Kuusipalo 1983).
The bilberry abundance increased along with secondary succession; i.e. the coverage was the higher, the higher the age and basal area of the stand up to certain limits; after which the effect of stand age and basal area gradually decreased.Bilberry suffers from clear-cuttings and subsequent soil preparation in extensive areas of Fennoscandia (e.g.Ferm and Sepponen 1981, Tolvanen 1994, Atlegrim and Sjöberg 1996, Salemaa 2000, cf.small clear-cuts in very humid conditions; Nielsen et al. 2007).Consequently, bilberry is sparse in seedling and sapling stands as well as in dense and shaded young thinning stands (e.g.Kardell 1980, Raatikainen et al. 1984, Salemaa 2000, Tonteri et al. 2005).In direct seeded and planted stands, bilberry is not as abundant as in naturally regenerated stands, where the light conditions are more suitable for bilberry and soil preparation is not so intensive.
The tree species composition also had an effect on the coverage of bilberry.In the pine-dominated stands, the coverage of bilberry was higher than that in stands dominated by spruce or deciduous trees.The decreasing effect of deciduous trees compared to spruce was significant only on herb-rich heath sites.In such fertile, well-lighted deciduous forests, herb and grass vegetation is the predominant understorey vegetation (Tonteri et al. 2005).Most probably due to the abundant herb and grass vegetation (caused by the site history and relatively high nitrogen concentration in the soil) (Wall and Hytönen 2005), the coverage of bilberry was low in afforested stands, i.e. on former agricultural land.The proportions of pine, spruce and deciduous trees affect the light conditions in the field vegetation layer.It seems evident that a moderate supply of light is an essential factor affecting the abundance of bilberry (e.g.Raatikainen et al. 1984, Salemaa 2000, Nielsen et al. 2007).Hence it follows that, in pine-dominated forests in Finland, bilberry may be the predominant species and a good competitor, and have a strong effect on all the vegetation (Kuusipalo 1983).In spruce forests, more shadetolerant herbs and other plants (e.g.Laakso et al. 1990), and in well-lighted deciduous forests, herbs and grasses may gain benefit from the lower competitive ability of bilberry.Thus, even if the site quality is equal, the abundance relationships between bilberry and other components of the field layer vary over wide ranges.
The temperature sum was expected to be a significant predictor of bilberry coverage.However, the possible effect of temperature sum and subsequent geographical variation appeared through the altitude and stand characteristics.It is also known that bilberry thrives on poorer soils in northern Finland than in southern Finland (Kalela 1970, Tonteri et al. 2005).The coverage on less fertile sites should increase towards the north due to the more humid climate and increasing soil moisture (Kardell 1980, Salemaa 2000).

Berry Yield
The site type was not a significant predictor of the number of unripe bilberries in the dataset collected from North Karelia in 2006 (NFI10_NK), although the number of unripe berries increased with increasing site fertility (Table 1).This may be due to the low number of observations (stands) in herb-rich and sub-xeric heath forests and the poor bilberry year.Note that the coverage of bilberry was also higher on fertile sites (Table 1), which may cover the effect of site quality on the berry yield in the small dataset.In general, forest stands of medium fertility (i.e.mesic sites) produce the highest yields (e.g.Ihalainen et al. 2003, 2005 andreferences therein).Bilberry production on sub-xeric sites is also good (e.g.Eriksson et al. 1979, Raatikainen et al. 1984, Jäppinen et al. 1986, Belonogova 1988, Kuchko 1988, Ihalainen and Pukkala 2001).This is often also the case with mature coniferous, herb-rich heath forests in moderate light conditions, and with xeric heath forests in northern Finland (Raatikainen et al. 1984, Turtiainen et al. 2005).
Tree species had a significant effect on the number of bilberry in both the NFI10_NK and the MASI data, and separate models were estimated for pine-and spruce-dominated stands.The yields were higher in pine-dominated than in spruce-dominated forests, which is in agreement with the findings of earlier studies (e.g.Eriksson et al. 1979, Ihalainen andPukkala 2001).In addition to the percentage cover of bilberry, the stand basal area also significantly affected the number of berries in spruce stands (NFI10_NK).In pine stands, on the other hand, only the percentage coverage was a significant predictor.The number of berries increased when the coverage of bilberry increased; the maximum was reached at a coverage of 38% in spruce stands and 43% in pine stands.The number of berries was at its maximum when the basal area of spruce stands was 14 m 2 ha -1 .Spruce-dominated stands which are not too dense produce relatively good yields (Raatikainen et al. 1984, Ihalainen et al. 2005).Spruce shades the ground vegetation much more than pine; and this strong shading has a negative effect on berry production (Raatikainen et al. 1984, Laakso et al. 1990, Ihalainen and Pukkala 2001).Clear-cutting areas, seedling and sapling stands, dense and shady young thinning stands, produce the lowest yields (e.g.Kardell 1980, Raatikainen et al. 1984, Ihalainen and Pukkala 2001).According to Raatikainen et al. (1984), the dominant tree species was a more important factor than the site type in explaining the bilberry yield; surprisingly, they found no significant correlation between the bilberry coverage and yield, which was probably due to poor bilberry years.
In this study, a mean fresh weight of 0.35 g for one berry (Eronen 2004) was used for all years.Ihalainen et al. (2003) applied a mean weight of 0.31 g for ripe berries; this value was obtained from the bilberry yield data collected in years 1981-1984.It is known that the weight may vary between the years, especially on sub-xeric and xeric sites, due e.g. to variation in the moisture conditions (Raatikainen et al. 1984).In addition, berries are assumed to be slightly heavier on moist and fertile sites than on dry and infertile ones (e.g.Eronen 2004).
The bilberry yields in the MASI data are higher than the average due to the sampling method.The bilberry populations were chosen such that the bilberry plants were pickable, i.e. in addition to a high coverage, bilberry had been observed to produce good crops in earlier years (Salo 1999).If a random or systematic sample of the Finnish forests was used, there would also be sample plots with few berry plants and, consequently, lower yields.The stands of the NFI10_NK data were assumed to be potential growing sites for bilberry on the basis of the NFI measurements (e.g.dominant tree species, stand age and site fertility).In the NFI10_NK data collected from North Karelia, where was a poor bilberry year in 2006, the mean number of unripe bilberries was 9.2 (SD 19.1) berries per m 2 (Table 1).Assuming that the proportion of unripe and ripe berries is 0.8 (Ihalainen et al. 2003) and the mean weight of one berry is 0.35 g (Eronen 2004), then the mean yield of the NFI10_NK data was 26 kg ha -1 (SD 53 kg ha -1 ).The corresponding mean yield from different parts of Finland in the MASI data, was much higher, i.e. 180 kg ha -1 (SD 154 kg ha -1 ).

Model Evaluation
The models estimated for the abundance and yield of bilberry were based on the datasets collected from the forests, which partly have a different management history than tomorrow's forests will have in Finland.Today's mature forests have been regenerated naturally without site preparation (excluding controlled burning) and they most likely were not fully stocked in their young stages.During the last decades, artificial regeneration methods with site preparation have been used commonly and young stands have been recommended to be grown dense.Changing forest management may retard bilberry to recolonise regeneration areas and disperse as stand ages.Therefore, the estimated relationship between bilberry yields and stand characteristics may not be just the same in the future.In our models, the effects of changing forest management on bilberry coverage and yield were considered to some degree by stand age, stand basal area and regeneration method.Moreover, the high morphological plasticity and rapid vegetative recovery of the deciduous bilberry allow it to respond rapidly to changed environmental conditions (Tolvanen 1995, 1997, Tolvanen and Laine 1997), e.g. to typical and possibly suitable light conditions (for bilberry) in future mature forests.
The simulations showed that the performance of the estimated models was logical and the predicted bilberry yields corresponded well with the bilberry yields presented earlier (see Turtiainen et al. 2005).Especially in the spruce stands, the bilberry yields in dense stands were considerably lower than those in thinned stands (Fig. 2).In spruce stands, the effect of thinning on the abundance and yield of bilberry was reflected through stand basal area.Bilberry yields decreased strongly when the stand basal area exceeded the optimum stand basal area for the berry yield (i.e. 14 m 2 ha -1 ), but increased temporarily after thinnings at the stand age of 43 and 58 years (at the stand basal area of 26 m 2 ha -1 and 29 m 2 ha -1 ) (Fig. 2).According to the simulations, thinnings were also needed to maintain the percentage coverage of bilberry in both pine and spruce stands (Fig. 1).In southern Finland stand basal area was kept near to the optimal level for the coverage of bilberry (i.e.24 m 2 ha -1 ), but in northern Finland the post-thinning stand basal area was much less than 24 m 2 ha -1 , and thus the predicted coverage of bilberry decreased temporarily after thinnings.
Our models may underestimate the bilberry yields, especially in mature spruce stands (Fig. 2).The coverage of bilberry and the number of berries correlated significantly (0.78, p < 0.001) in the NFI10_NK data collected during a poor berry year (cf.Raatikainen et al. 1984).However, the relationship between the coverage and number of berries may not be fixed, but varies from year to year depending on the quality of the berry year.In the simulations, the parameter estimate of the coverage was fixed, and only the intercept of Model 2 was calibrated to correspond to an average berry year using the year effects of Model 3 (Table 4).Because our models predict the mean yield of stands having the same site and stand characteristics, the predictions are also lower than the yield of those stands where bilberries are usually collected (cf.Eriksson et al. 1979, Ihalainen et al. 2003).
The models were estimated for use with stand simulators where they would predict the abundance and yield of bilberry for silvicultural and forest planning purposes.Because the site and stand characteristics are known in forest planning and forest management greatly affects the stand characteristics, site and stand characteristics are the most reasonable predictors in the models (Ihalainen et al. 2003).The bilberry yields also depend on the state of the understorey vegetation and weather conditions, information about which is, however, not available in practise.Due to the large annual variation in bilberry yields, yield data should be collected from the same stands for many years (see Ihalainen et al. 2003).Thus, we first predicted the abundance and then the annual yield of bilberry as a function of site and stand characteristics, including bilberry coverage.
The abundance and the annual variation in the bilberry yield were predicted using extensive datasets (the sample plots of the NFI and the permanent sample plots in different areas of Finland, respectively), but the linkage between the abundance and berry yield was based on the regional dataset collected in North Karelia in a single year.The modelling data, and consequently the models, could be improved by determining the percentage cover of bilberry and stand characteristics (e.g.dominant tree species and stand basal area) for the stands of the MASI data.Using these additional measurements and the modelling approach of this study, more reliable empirical models could be obtained for the bilberry yield in the whole of Finland.So far, empirical models for the bilberry yield have been prepared only for North Karelia (Ihalainen et al. 2003).In order to promote the exploitation of wild berries, bilberry yields predicted by the estimated models could be utilized, for example, in berry-yield-oriented forest management planning, or by producing maps of potential stands for bilberry picking.So far, the Finnish Forest Research Institute (Metla) has boosted wild berry picking by developing national yield forecasts based on sample plots of the MASI dataset.

Fig. 1 .
Fig. 1.Simulated coverage of bilberry in Scots pine and Norway spruce stands of different site type in northern (900 dd.) and southern Finland (1200 dd.) using Model 1 (Table 2).The stand development was simulated using the Motti simulator, and thinnings and final cutting were simulated according to the silvicultural recommendations.Site type: I = herb-rich forest, II = herb-rich heath (Oxalis-Myrtillus group), III = mesic heath (Myrtillus group), IV = sub-xeric heath (Vaccinium group), and V = xeric heath forest (Calluna group).Arrows indicate thinnings.

Table 4 .
The multi-level Poisson models (Model 3, Eq. 2) estimated for the mean number of ripe berries on five 1-m 2 quadrates in Scots pine and Norway sprucedominated stands of the MASI data, measured in 2001

Fig. 2 .
Fig. 2. Simulated average annual yield of bilberry and its 95% confidence interval in a Scots pine and Norway spruce stand on a mesic heath site (Myrtillus group) in southern Finland (1200 dd.) using Models 1-3 (Tables 2, 3 and 4).Predictions were compared to the yields computed by the model of Turtiainen et al. (2005; Eq. 2 andTable 4).Arrows indicate thinnings.

Table 1 .
Main characteristics of bilberry, sites and stands in the study materials.

Table 2 .
The multi-level binomial model (Model 1, Eq. 1) estimated for the mean percentage coverage of bilberry on the 2-m 2 quadrates in the stands of the NFI data.Site type: I = herb-rich forest, II = herb-rich heath (Oxalis-Myrtillus group), III = mesic heath (Myrtillus group), IV = sub-xeric heath (Vaccinium group), and V = xeric heath forest (Calluna group).

Table 3 .
The multi-level Poisson models (Model 2, Eq. 2) estimated for the mean number of unripe berries on four 1-m 2 quadrates in Scots pine and Norway spruce-dominated stands in North Karelia, measured in 2006.