Drivers for Annual Cork Growth under Two Understory Management Alternatives on a Podzolic Cork Oak Stand

: Understory management practices and stand density characteristics allow one to distinguish a cork oak traditional silvopastoral system (known as a montado ) from a cork oak forest system. Although understanding the manner in which different management practices affect cork growth is imperative, there are still only a few outputs from experimental research that contribute to this knowledge. The effect of potential drivers on annual cork growth was analyzed using a linear mixed model approach. Two dimensions of drivers were considered: intraspeciﬁc competition, assessed by tree level distance-dependent indices; and interspeciﬁc competition, assessed by variables characterizing understory management. The present dataset was collected from an experimental trial established on a cork oak stand in Podzolic soil on the Tagus river basin, covering two different cork growth cycles over the period from 2003 to 2015. The adjusted models considered two understory management alternatives: spontaneous shrubs maintenance and forage application. In both models, annual precipitation displayed a positive effect on annual cork growth, as expected. However, no signiﬁcant effect of intraspeciﬁc competition was found. Additionally, there was a positive effect on annual cork growth associated with the spontaneous shrubs growth and a negative effect associated with lupine presence; both effects linked to different cork ring ages’ thresholds. The study main contributions are the following: (i) the introduction of the interaction between cork growth cycle stage and understory management practices, only possible with cork sample collections from different cork rotation cycles; (ii) the ﬁnding that there was no signiﬁcant effect of intraspeciﬁc competition on cork growth.


Introduction
Cork oak (Quercus suber L.) is an important species in the Mediterranean forest, as it is provides raw material suitable for cork stoppers and other high value manufactured products. This species mainly occurs in a traditional silvopastoral system called a montado, characterized by low density and sparse trees. The species also occurs in a forest system with higher tree density, usually dominated by spontaneous vegetation. The Portuguese National Forest Inventory characterizes cork oak stands by an average density of around 66 trees per hectare [1], making no distinction between these two systems: montado and cork oak forest. Since there is no clear definition of a tree density threshold to differentiate between these two systems, they can only be distinguished by their management practices and production aims: the first is oriented for cork production, crop production and/or grazing; the second focuses mainly on cork production.
In a montado, the standard understory management implies mechanization practices, which may influence several ecosystem components, such as tree regeneration, soil biota, soil compaction, and/or root development (e.g., [2][3][4][5]). These practices have several aims: facilitate access to trees during the debarking operation; reduce phytovolume of a specific native or invasive shrub species; reduce fire hazard; allow permanent and improved pastures to increase forage availability for the livestock; contribute to soil fertility (e.g., [6][7][8][9]).
There are still few studies focusing on the impact of stand density and management practices on cork growth; namely, there is a clear lack of knowledge, when considering simultaneously the effect of intraspecific and interspecific competition. Caritat et al. [10] also looked for tree radial growth under distinct understory management alternatives, but did not separate wood and cork growth. Paulo et al. [11] explored the effect of intraspecific competition using distance-dependent tree competition indices, but focused on tree diameter and crown variables. Sánchez-González et al. [12] developed a total cork thickness model considering relative tree dimension variables and one single distance-independent competition index, but did not test the interaction with alternative management options. Faias et al. [13] compared the annual cork pattern growing under two distinct understory management alternatives along one cork debarking rotation cycle of 9 years, but did not consider tree intraspecific competition effect.
The dataset in this study is an extended version from the one used by Faias et al. [13], with additional cork samples extracted in the same trial from a different cork growth cycle, covering the period from 2006 to 2015. Tree coordinates, as well as understory evolution information, were also collected for the present study. This extended dataset allows us to explore not only the effect of both intra-and interspecific competition drivers on cork's annual growth, but also its interaction with the ring age due to the existence of cork samples from two different cork debarking rotation cycles. The study intends to answer the following questions that are relevant to support management decisions, such as deciding the best timing to remove the understory: (i) Is there an impact of intraspecific competition on cork growth?; (ii) What is the impact of understory (interspecific competition) management practices on cork growth?; (iii) What is the interaction between cork age and understory management?

Study Area and Trial Description
An experimental trial to compare two understory management alternatives was implemented, as a complete randomized design with two blocks, in uneven-aged cork oak stands near Montargil village (39 • 3.24 N, 8 • 10.58 W). Each block contains two quadrangular plots, with 2 hectares each, including a 20-m border to ensure no impact of non-treated areas on the monitored trees ( Figure 1). Soils were classified as Podzols, according to the IUSS (International Union of Soil Sciences) Working Group [14], available at the Portuguese Agency for the Environment website (http://sniamb. apambiente.pt/webatlas/). Soil samples, with pH 5.5 to 6.1 and organic matter percentage of 1.5% to 1.8%, indicated similarity between plots. Monthly precipitation was gathered from the nearest available meteorological station at the SNIRH (Sistema Nacional de Informação de Recursos Hídricos; www.snirh.pt) network for the period included in the dataset (Table 1).
Between 2003 and 2012, two different understory management operations were performed in one plot randomly selected in each block ( Figure 1): (i) periodic removal of the spontaneous shrubs species, with incorporation of organic matter into the soil, followed by the seeding of Lupinus luteus L. (RUL-understory removal with lupine); (ii) spontaneous shrubs species maintenance, during a cork debarking rotation cycle of 9 years (NUR-no understory removal). The lupine presence was confirmed by visual assessment in the field until 2009, and after its first application in 2003 the need for seeding again was identified in 2007 and 2009 (Table 1). In 2007, lupine seeding had a low germination rate, due to unsuitable environmental conditions. Table 1. Yearly operations calendar carried out for each understory management alternative. with corresponding annual precipitation, computed from October of the year before the growth period to September of the growth period year, shrubs height (HShrubs) estimation with equation 3, and cork ring ages for each cork year sampling. RUL represents the understory removal followed by lupine seeding and NUR represents the maintenance of spontaneous understory vegetation.

Period
Precipitation

Understory Evaluation
Stand understory layer was composed of spontaneous shrubs species dominated by Cistus salvifolius L. (25%), alongside with sparsely distributed Ulex airensis M.D. Esp.Santo, Cubas, Lousa, C.Pardo & J.C.Costa (15%), Lavandula pedunculata (Mill.) Cav. (10%), and Calluna vulgaris (L.) Hull (5%). Until 2003, the understory was mechanically removed with a forestry mulcher, in the complete stands area, with an interval of 3 to 4 years. In order to assess understory evolution with time, since 2012, for each of this shrub species, additional measurements of shrub height (m) were performeda minimum of two samples per species-in the following years after clearing (TShrubs): 2, 4, 6 and more than 9 years in a similar neighboring stand. These measurements were used to develop a simple quadratic model, dependent of TShrubs, as in Navarro et al. [15], allowing the estimation of shrub height in intermediate years.

Tree and Cork Measurements
One important feature of the trial that allowed the present study to go forward in comparison to the one used by Faias et al. [13] was the presence of two different cork rotation periods within the stand understory management period: 2003 to 2012 and 2006 to 2015 (Figure 1). This detail allows the interaction between cork ring age with understory removal and periodic lupine application to be studied. The diameter at breast height without cork (du) of each debarked tree was measured at each debarking year. Mean stand density between 2003 and 2015 had a slight trees per hectare variation, from 130 to 120 on block 1 and 93 to 74 on block II. Individual tree coordinates were taken, considering the Cartesian coordinate system, in order to allow the computation of distance-dependent competition indices.
Cork samples, approximately 20 cm × 20 cm, were taken at breast height from the debarked trees in 2003, 2006, 2012, and 2015. The samples were boiled for 1 h in water at 100 °C and atmospheric  Hull (5%). Until 2003, the understory was mechanically removed with a forestry mulcher, in the complete stands area, with an interval of 3 to 4 years. In order to assess understory evolution with time, since 2012, for each of this shrub species, additional measurements of shrub height (m) were performed-a minimum of two samples per species-in the following years after clearing (TShrubs): 2, 4, 6 and more than 9 years in a similar neighboring stand. These measurements were used to develop a simple quadratic model, dependent of TShrubs, as in Navarro et al. [15], allowing the estimation of shrub height in intermediate years.

Tree and Cork Measurements
One important feature of the trial that allowed the present study to go forward in comparison to the one used by Faias et al. [13] was the presence of two different cork rotation periods within the stand understory management period: 2003 to 2012 and 2006 to 2015 ( Figure 1). This detail allows the interaction between cork ring age with understory removal and periodic lupine application to be studied. The diameter at breast height without cork (du) of each debarked tree was measured at each debarking year. Mean stand density between 2003 and 2015 had a slight trees per hectare variation, from 130 to 120 on block 1 and 93 to 74 on block II. Individual tree coordinates were taken, considering the Cartesian coordinate system, in order to allow the computation of distance-dependent competition indices.
Cork samples, approximately 20 cm × 20 cm, were taken at breast height from the debarked trees in 2003, 2006, 2012, and 2015. The samples were boiled for 1 h in water at 100 • C and atmospheric pressure, and were polished with a belt sander machine (FAI™ Modelo LCU, Paredes, Portugal). Cork rings were marked on a vertical axis in samples where rings were unmistakably identified. Annual cork growth or cork ring width (rw) was measured for the eight complete rings after digital scanning, using the image analysis software ImageJ [16,17]. These measurements include different cork ring ages within the same period (Table 1), having a range of mean values between 1.5-4.5 in 2006 (n = 98 samples), 2.1-3.9 in 2012 (n = 114 samples), and 1.7-4.7 in 2015 (n = 101 samples). Cork samples from 2003 were used only in the previous study from Faias et al. [13], as a baseline to characterize both stands before the trial establishment.

Annual Cork Growth Modeling
Due to the nested structure of the data-trees inside plots and plots inside blocks-the data analysis was carried out using a linear mixed modelling approach [18]. Comparing the large number of cork ring measurements to the number of repeated cork rings measured, with the same age at the same tree, in different cork samples (from 2006 and 2015), an absence of serial correlation was assumed, as in Paulo et al. [19].
The first step for the model development was the definition of a tree distance-depending intraspecific competition index (CI) more correlated with cork ring width. The indices tested were those that were demonstrated by Paulo et al. [11] to be significantly correlated with tree diameter ( Table 2). The competitors were selected from a 20-m circle around the subject tree. These indices were fitted, one at a time, using the complete dataset with the following equation: where rw is the cork ring width (mm) for tree i within plot j and block k; CI is a competition index; β 0 , β 1 are the fixed effects parameters; µ 0jk is the random effect; and ε ijk is the error term. Table 2. List of competition indices tested, respective parameter estimates (β 1 ), significance (p-value), and Akaike information criterion (AIC) values, where ij is the distance between subject tree i and its competitor tree j; du is the diameter at breast height under cork; Nc is the number of competitors. Hyperbolic cork ring age (tc) and annual precipitation (Pr) were explanatory variables included in the model, due to the known decreasing pattern of the annual cork growth with cork age and the positive correlation with annual precipitation (e.g., [20][21][22][23][24]). The precipitation was computed, as defined in Faias et al. [13], for the period between October of the year before the growth period and September of the growth period year (Table 1). Lastly, two final sets of models were separately fitted with the dataset for each understory management alternative, RUL or NUR, according to the following expression: rw ijt = (µ 0j + β 0 ) + β 1 tc ijt + β 2 Pr t + β 3 CI ij + β 4 Tr jt + β 5 tc ijt × Tr jt + ε ijt (2) where rw is the cork ring width (mm); tc is the cork ring age (years); Pr is the annual precipitation (mm); Tr represents the variables of each understory management alternative described in Table 1; β 0 , β 1 , β 2 , β 3 , β 4 , β 5 , are the fixed effects parameters; µ 0j is the random effect; i is the tree index; j is the plot index; t is the cork growth year index; and ε ijt is the error term. Considering the RUL management, two alternative variables were tested (Table 1): LSeeding, which takes value 1 in the following year after lupine seeding in Autumn; and PLupine, which identifies the visual presence of lupine, by taking value 1 in the two years following its seeding. For the NUR management, three variables were considered (Table 1): Shrubs, which identifies the presence of shrubs species in understory observed two years after clearing; TShrubs, which identifies the time without understory removal; and HShrubs, which gives the understory mean height (m) estimated with the previously developed model as a function of time.
All models were fitted using the procedure PROC MIXED from software SAS 9.4 (SAS, Cary, NC, USA) [25]. The RANDOM statement, from this procedure, was applied to specify the random effects associated with nested structure (trees inside plots). The random effect parameter was tested in all the fixed parameters and the selection criteria for its inclusion in the model was the lowest Akaike information criterion (AIC) value. Variance components used the covariance structure. Significance of the parameter estimates was evaluated considering α = 0.05.

Understory Model
The main shrub species within the trial showed an increasing pattern of total shrub height (HShrubs) in the following years after understory removal (Figure 2), which potentially reached a maximum total height at about 7 years, before starting to decline afterwards. The following equation was fitted for the mean total height (R 2 = 0.93), as a function of time without understory removal: Forests 2018, 9, x FOR PEER REVIEW 6 of 13 identifies the visual presence of lupine, by taking value 1 in the two years following its seeding. For the NUR management, three variables were considered (Table 1): Shrubs, which identifies the presence of shrubs species in understory observed two years after clearing; TShrubs, which identifies the time without understory removal; and HShrubs, which gives the understory mean height (m) estimated with the previously developed model as a function of time.
All models were fitted using the procedure PROC MIXED from software SAS 9.4 (SAS, Cary, NC, USA) [25]. The RANDOM statement, from this procedure, was applied to specify the random effects associated with nested structure (trees inside plots). The random effect parameter was tested in all the fixed parameters and the selection criteria for its inclusion in the model was the lowest Akaike information criterion (AIC) value. Variance components used the covariance structure. Significance of the parameter estimates was evaluated considering α = 0.05.

Understory Model
The main shrub species within the trial showed an increasing pattern of total shrub height (HShrubs) in the following years after understory removal (Figure 2), which potentially reached a maximum total height at about 7 years, before starting to decline afterwards. The following equation was fitted for the mean total height (R² = 0.93), as a function of time without understory removal:

Annual Cork Growth Models
Testing each of the distance-dependent competition indices as unique explanatory variable for the cork ring width model, the only index presenting a coefficient significantly different from zero ( Table 2) was the one expressed as the number of competitors, with tree diameter higher than the subject tree diameter, within a 20-m radius (CI2). This index estimated parameter showed a negative

Annual Cork Growth Models
Testing each of the distance-dependent competition indices as unique explanatory variable for the cork ring width model, the only index presenting a coefficient significantly different from zero (Table 2) was the one expressed as the number of competitors, with tree diameter higher than the subject tree diameter, within a 20-m radius (CI2). This index estimated parameter showed a negative effect, demonstrating that an increasing number of larger trees within a 20-m radius is associated with a reduction of the cork ring width. Later, the CI2 index was included while fitting the full models for cork ring width (Equation (2)). Table 3 shows the main results for the two alternative models (original and final) for the RUL model, showing the results of four alternative models. The LSeeding variable, as well as its interaction with cork age, was not statistically significant (RUL1A). However, the PLupine variable was statistically significant, presenting a negative coefficient, as well as, its interaction with the cork age with a positive coefficient (RUL2A). In both RUL models, the CI2 index was not significant and was removed from the final model (RUL1B, RUL2B).
Considering the NUR model, Table 4 shows three alternative models (original and final), where the coefficients for both variables, Shrubs and TShrubs (Table 1), and their interaction with cork age were statistically significant, with similar AIC values (NUR1A, NUR2A). The coefficient for HShrubs (NUR3A) was positive and the coefficient of its interaction with cork age was negative. For these models, the CI2 index was not significant and therefore was removed from the final models. Simultaneously, the coefficient for the cork ring age was not statistically significant, and the final model was fitted excluding this variable (NUR1B, NUR2B, NUR3B).
To understand either the dynamics of PLupine or HShrubs interactions with cork ring age on the cork ring width, the first partial derivative was computed for each understory management variable (tc > β 4 /β 5 ). The value obtained for the PLupine interaction suggests a negative effect on the ring width, up to 5.4 years of cork growth (β 4 /β 5 < 0); and then it is converted into a positive effect (β 4 /β 5 > 0). For the HShrubs interaction there was a positive effect on cork ring width up to 6.6 years of cork growth (β 4 /β 5 > 0), representing the point where the shrubs potentially reached their maximum height, and then this effect was no longer observed (β 4 /β 5 < 0).
The significance of the plot's random effects was tested in all the fixed parameters through the refitting of both models. Results showed that the random effect parameter convergence, significance and a lower AIC value was obtained when it was added to the intercept parameter of each model. Table 3. Parameter estimates for the two alternative cork ring width (rw) models (with and without the CI variable) regarding the RUL management variables listed in Table 1 (n = 1008), where: tc is the cork ring age (years); Pr is the annual precipitation (mm); CI2 is the number of competitors, within 20 m, with duj > dui; LSeeding takes value 1 in the year after lupine application (Autumn); PLupine takes value 1 in the two years following its seeding.  Table 4. Parameters estimates for the three alternative cork ring width (rw) models (with and without the CI variable) regarding the NUR management variables listed in Table 1 (n = 1006), where: tc is the cork ring age (years); Pr is the annual precipitation (mm); CI2 is the number of competitors, within 20 m, with duj > dui; Shrubs is the presence of shrubs in understory observed two years after clearing; TShrubs is the time without understory removal; HShrubs estimates understory mean height (m) as a function of time.

Intraspecific Competition Driver: Tree Level Distance-Dependent Indices
A previous study by Paulo et al. [11] revealed a significant effect of distance-dependent tree competition indices on tree crown characteristics, when looking for pure cork oak stands from 98 to 238 trees per hectare. For the present study, the number of competitors within a 20-m radius with tree diameter larger than the subject tree diameter presented a negative effect on annual cork growth, when fitted as single explanatory variable. However, when this tree competition index was tested jointly with the cork age, the annual precipitation, and the understory management related variables, it did not reveal any further significant effect on annual cork growth for this specific stand density range. Hence, there should be further research on the intraspecific competition effect on the tree and cork growth for different stand densities and structures.

Interspecific Competition Driver: Lupine Application Characteristics
The periodic application of a leguminous species (lupine) beneath cork oak trees as in the RUL plots from the present trial represent a common understory management practice of the montado system, but in this particular case without cattle pasture. Faias et al.'s [13] results indicated that periodic lupine application did not show any effect on annual cork growth, as well as, tree radial growth, considering one cork debarking rotation cycle. The same conclusion could be obtained from the present study if only considering the effect of the lupine application in the autumn on the cork ring of the following year. However, when considering the lupine visual presence on the field-the two years following its application-a negative effect on the annual cork growth was revealed until cork ring ages around 5 years. At this point, it changes to a positive effect on cork ring width. More research on the effect of other applied forage species on the annual cork growth is required, for instance, with natural and sown rainfed grasslands [8], which have been associated with ecosystem biodiversity increase, pasture productivity, and resilience [26].

Interspecific Competition Driver: Understory Growth Characteristics
The measured total shrub height on these specific spontaneous species in the understory showed a pattern over period of a cork debarking rotation cycle, in line with the reported pattern noted by Santana et al. [27], which had been described as having a sharp increase with time after clearing followed by a leveling off. It could also be seen that Cistus and Lavandula species potentially reach a maximum cover, and then start declining, as noted in Díaz Barradas et al. [28]. However, the Ulex species pattern was distinct from the increasing pattern showed in several studies, where it was considered a period longer than the one observed in this study (e.g., [29]). The present results revealed a positive effect of understory height on the cork ring width until 6 years of cork growth after debarking. From this cork ring age, the understory revealed a negative effect on annual cork growth. These results are in line with conservative understory management strategies that minimize soil intervention to reduce the negative impact on tree root development [30] and on tree regeneration patterns and viability [31]. It is also in line with the findings of Mendes et al. [4], which suggest that making a shrub cut every 5-6 years supports the normal natural regeneration of the understory shrubs, with a minor effect on local species richness of soil biota. In contrast, Caritat et al. [10] showed no significant effect of the presence of shrubs on radial growth on 10 cork oak trees growing with or without shrubs. Martin et al. [32] also indicated no significant effect when comparing understory and soil management practices on holm oak growth in sites with different soil and understory characteristics. These results should also not be seen as contradictory to the experimental study by Caldeira et al. [6], which suggests that oak trees are more vulnerable to hydraulic failure in a stand with shrub invasion, due to two main differences. First, the latter not only addresses extreme drought events, but also an extreme stand invasive shrubs condition, dominated by Cistus ladanifer L., which is not similar to any spontaneous shrub species in the present trial. Second, the NUR plots from the present study represent a common management alternative for the cork forest systems, whereas the other represents an extreme case of the absence of management.

Conclusions
Linear models developed in this study were used to contribute to the knowledge about the effect of intra-and interspecific drivers, such as distance-dependent competition indices or understory management operations, on annual cork growth. Their application should not be extrapolated for understory conditions and compositions distinct from the ones described-Podzolic soils with stand densities up to 132 trees per hectare. The analyses undertaken do not show any impact of intraspecific competition on cork growth, expressed through a set of competition indices. In fact, competition indices are not significant when added to a cork growth model that already includes cork age, precipitation of the previous growing year, and understory management.
Under these trial conditions, regarding interspecific competition, this study suggests that after debarking, the understory can be maintained during some years without affecting cork growth (6 years in this trial). Not taking into consideration fire risk, which would be higher for a cork forest system, this conclusion may lead to a decreasing number of clearings within a cork debarking rotation cycle that may be favorable to local biodiversity. On the other hand, lupine application immediately after debarking seems to have a negative impact on cork growth, but may be considered after some years (5 years in this trial).
Answering the last question research question originally proposed for this study, there is a clear interaction between the understory management and cork age, with a different impact of the first depending on whether its application takes place at the beginning of the cork debarking period or in another year within the period.