The study site is located in Valsaín (Segovia, Central Spain, 40° 51′ 35′′ N, 4° 3′ 52′′ W) in a north-west facing slope (15%) at an elevation of 1350 m a.s.l. It is found at the transition zone between P. sylvestris- and Q. pyrenaica-dominated woodlands and, thus, at P. sylvestris local dry elevational limit. The stand is composed of patches of mature P. sylvestris trees (sampled trees 100 ± 7 years old), younger P. sylvestris trees (sampled trees 43 ± 4 years old) and Q. pyrenaica trees (sampled trees 42 ± 1 years old). According to records obtained from the Spanish Meteorological Agency (AEMET) nearby stations, mean annual precipitation was 728 mm and mean annual temperature 9.1°C during the period 1950–2014.

To investigate the impact of drought on the physiology of the two studied species, three rainfall exclusion plots of 65–150 m² were established in April 2012 based on the structure of the stand: one for Q. pyrenaica and two for P. sylvestris, one for each size class to assess the species response to the treatment regardless of tree size (Table 1). Rainfall was excluded from reaching the soil by fully covering it with corrugated PVC boards, which remained year-round throughout the experiment. Additionally, a 0.5 m ditch was dug and covered with plastic on the upper side of each plot to minimize downward subsurface water flow. Soil humidity was monitored in all rainfall exclusion plots and the control area with sensors set at 25 and 50 cm depth.

Micro-cores were sampled from six trees per species and treatment, three per age class in the case of P. sylvestris, every 7–15 days between April and November during 2012–2014, with a total of 24 trees being monitored per year. In order to minimize inter-annual bias due to previous year punching injuries, we generally avoided using the same trees on consecutive years. Most of sampled trees subjected to rainfall exclusion were selected from the inner part of the plots to reduce the border effect. Samples were extracted with a Trephor (Rossi et al., 2006) and stored in a 70% ethanol solution at 4°C. P. sylvestris micro-cores sampled in 2012 were cut in sections of 10–18-μm in thickness with a sliding micro-tome. Since this technique was unsuccessful for Q. pyrenaica, the rest of the samples were dehydrated with increasing concentrations of ethanol and xylene, embedded in paraffin and cut with a rotary microtome in sections 8–12 μm thick. All sections were stained with Safranine O and Astra Blue, dehydrated with increasing concentrations of ethanol and xylene and mounted with either Canada balsam or a mixture of distyrene, tricresyl phosphate and xylene (DPX).

Stained sections were observed under visible and polarized light to distinguish the different xylem cell developmental stages: cambium, cell enlargement, cell-wall lignification and mature cells. For P. sylvestris, the number of tracheids in each phase was counted along three radial segments per date and then averaged. For Q. pyrenaica, photographs of the sections were taken with 10x magnification with an Olympus Colorview camera attached to a microscope. The number of cambium cells was counted and the width of earlywood and latewood increments at each cell development phase were measured along three segments with ImageJ (Schneider et al., 2012). Due to the complex anatomy of Q. pyrenaica’s xylem, the end of the wall-thickening phase could not be assessed on this species. Growth rates (r, in cells d^-1 or μm d^-1) were calculated as the ratio between the final number of tracheids or earlywood and latewood widths and the length of the cell enlargement phase, which represents a good approximation of the period of cell production.

Crown phenology was recorded on the same dates of wood sampling. In 2012, leaf phenology was recorded by species at the stand level, whereas in 2013 and 2014 it was recorded individually on the trees used to study cambial phenology. The dates of bud-burst, end of leaf and needle elongation and, in Q. pyrenaica, the beginning of autumn leaf discolouration, were identified.

To avoid bias in inter-annual comparisons due to the selection of different trees each study year, a micro-core including the tree rings belonging to the period 2011-2014 was taken from every monitored tree (Table 1) in November 2014 to assess the cumulative effect of the rainfall exclusion treatment on wood anatomy. Anatomical sections were prepared according to the abovementioned protocol and photographed with 10x magnification and a resolution of 2.2 pixels/μm. Cell-wall thickness, lumen area and earlywood and latewood widths were measured along three radial rows in P. sylvestris. For Q. pyrenaica, earlywood and latewood widths and the area of earlywood vessels with a lumen area > 7500 μm² were measured. Earlywood vessel area was additionally measured in all the section photographs taken for the analysis of cambial phenology. All measurements were carried out with ImageJ (Schneider et al., 2012).

Given the different functionality of earlywood (EW) and latewood (LW), the ratio between latewood width and ring width (LW/RW) was calculated for both species to assess whether rainfall exclusion influenced the proportion of each type of wood. Lumen areas were used to estimate vessel and tracheid radii, which were in turn used to calculate the hydraulic diameter (D_H). D_H, computed as Σd⁵/Σd⁴, where d is conduit diameter, is proportional to the hydraulic conductance given by the Hagen–Pouseuille law and, thus, a proxy of xylem hydraulic conductivity (Sperry et al., 1994).

Morning net photosynthetic rate (A; μmol m^-2 s^-1) and stomatal conductance to water vapor (g_s; mol m^-2 s^-1) and midday leaf water potential (Ψ; MPa) were recorded for P. sylvestris in late July or early August between 2012 and 2014. Six trees per treatment (three per age class) were measured. The inaccessibility of leaves of Q. pyrenaica trees at the rainfall exclusion plot precluded the possibility to carry out measurements in this species. Gas exchange measurements were performed on current-year needles exposed to direct sunlight with an LCpro+ portable photosynthesis system (ADC BioScientific Ltd., Hoddesdon, Herts, United Kingdom). Leaf water potential was measured with a Scholander pressure chamber (PMS Instrument Co. 7000, Corvallis, OR, United States). For each tree, a subsample of the measured needles was scanned to measure needle area. Subsequently, needles were dried and weighed to obtain their specific leaf area (SLA; cm² g^-1), calculated as the ratio between needle area and dry mass. SLA was used to recalculate net photosynthetic rate and stomatal conductance on a leaf area basis.

Additionally, collected needles were ground to analyze carbon and nitrogen elemental concentration and stable carbon isotope ratio (δ¹³C). Samples of 3–4 mg were analyzed using a PDZ Europa ANCA-GSL elemental analyzer interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, United Kingdom). The average δ¹³C analytical precision was 0.1‰. δ¹³C data were transformed to discrimination (Δ) following Farquhar et al. (1989):

Source air δ¹³C (δ¹³C_a) data were obtained from NOAA Earth System Research Laboratory (White et al., 2015).

Needle carbon and nitrogen concentrations on a mass basis were used to calculate C/N ratio and photosynthetic nitrogen use efficiency (PNUE; μmol g^-1 s^-1), which is the ratio between net photosynthetic rate and leaf nitrogen content.

For each sampling date, phenological observations were converted into presence/absence (binary) data. Logistic functions were used to estimate the day of the year (DOY) at which each cambial and foliar developmental phase began and ended. Preliminary analysis showed no differential effect of the treatment between P. sylvestris tree sizes. Hence, data were pooled by species to increase sample size and, thus, statistical robustness. Differences in phenology between treatments were assessed fitting generalized linear mixed models (GLMM) with a binomial distribution and a random intercept associated with tree. The assessed fixed factors were DOY, year (to account for inter-annual climatic variability), treatment and the interaction ‘treatment × year.’ The Akaike Information Criterion (AIC) was used to compare the different GLMMs fitted. Lower values of AIC indicate a better fit. A model was considered to have a significantly higher explanatory power than another when the difference in AIC (ΔAIC) between models was ≥ 2 (Burnham and Anderson, 2002). If ΔAIC < 2, the difference among models was considered non-significant and, therefore, the more parsimonious model (i.e., with fewer variables) was selected.

Because data were not normally distributed and measures were repeated within individuals, differences between treatments in the number of cambial cells, total growth, growth rates and wood anatomical features were assessed using GLMMs with a Gamma distribution and an identity link function, except for Q. pyrenaica growth rates, for which a log link function was used. We fitted GLMMs with the fixed effects treatment, year and the interaction ‘treatment × year’ and a random intercept to account for correlation within trees. If the best model, given by its AIC, included the interaction ‘treatment × year,’ we tested the differences in specific years using Tukey contrasts. Differences in gas exchange and leaf chemical properties between treatments were assessed using Mann–Whitney–Wilcoxon test.

All analyses were carried out with R version 3.1.1 (R Core Team, 2014). GLMMs were fit with the packages “lme4” (Bates et al., 2014) and “glmmADMB” (Fournier et al., 2012). Multiple mean comparisons (Tukey contrasts) were performed with the package “multcomp” (Hothorn et al., 2008).

In Pinus sylvestris, the number of cambial cells during dormancy was 6 ± 2, whereas active cambium had 10 ± 3 cells (Supplementary Figure S1). The number of cambium cells, both active and in dormancy, was significantly lower in trees subjected to rainfall exclusion only in 2013. Xylem formation started between mid-April and early-May (DOYs 105-122), with cell wall-thickening starting on average 2–5 weeks later (Figure 2 and Supplementary Figure S1). The first cells completed maturation between the end of May and the end of June (DOYs 142-176), depending on the year (Figure 2). Xylem enlargement ended between mid-July and early-September (DOYs 203–234), depending on the year. Ring maturation was completed in October to early-November (DOYs 283–299). No significant differences between treatments were found in the beginning of the enlargement or wall-thickening phases (Table 2), whereas rainfall exclusion significantly influenced the end of both phases, with trees subjected to rainfall exclusion finishing growth one to 3 weeks earlier than control trees (Table 2 and Figure 2). The strength of the effect depended, however, on the year (Table 2).

Moreover, one of the old pine trees subjected to drought died between the 2013 and 2014 growing seasons. During 2013, cambium was mostly inactive and formed at most four tracheids, which in many cases did not lignify (Supplementary Figure S2). Cambial inactivity occurred before crown damage was evident. Since the growth of this tree before the rainfall exclusion was applied did not differ from that of other trees in the stand, its death was likely a direct consequence of the intensification of drought by the rainfall exclusion treatment. Furthermore, the presence of mistletoe (Viscum album L.) on this tree might have exacerbated the rainfall exclusion effect.

In Q. pyrenaica, the number of cambium cells during dormancy was 8 ± 1, whereas active cambium had 10 ± 2 cells, with no significant differences between treatments (Supplementary Figure S3). Xylem growth presented a bimodal pattern, corresponding to the peaks in earlywood vessel and latewood formation (Supplementary Figure S3). Earlywood growth onset occurred between late-March and early-April (DOYs 91–94) and vessel lignification started 1–8 weeks later, depending on the year (Figure 2). Latewood development started between June and early-July (DOYs 160–182) and its lignification began up to 2 weeks later (Figure 2). Latewood enlargement ceased at the end of July-mid-August (DOYs 211–231). No significant differences were observed between treatments in Q. pyrenaica cambial phenology (Table 3 and Figure 2).

Budburst occurred synchronously among individuals of each species within each study year; in P. sylvestris, between DOYs 132 and 147, depending on the year (Figure 2). In Q. pyrenaica, budburst occurred between DOYs 97 and 133, depending on the year, with a second leaf sprouting occurring between DOYs 173 and 189 (Figure 2). Oak leaf discoloration started between DOYs 233 and 279, depending on the year. Leaf phenology was unaffected by rainfall exclusion in both species (Tables 2, 3).

No significant differences were found between treatments in P. sylvestris growth rates for any of the studied years (Figure 3). For both treatments, growth rates were highest in 2012, the driest year. Control trees had significantly more tracheids per ring than trees subjected to rainfall exclusion only in 2012 (Figure 3). None of the other anatomical features assessed (LW/RW ratio, D_H and cell-wall-thickness) showed differences between treatments (Figure 3). However, the ratio LW/RW was lowest in 2012, the driest year, and highest in 2014, the wettest. Cell-wall thickness behaved similarly, whereas D_H decreased in 2014, although this drop was only significant for control trees.

No significant differences were found between treatments in Q. pyrenaica growth rates, neither for EW nor for LW, although LW rates tended to be lower under rainfall exclusion (Figure 4). Both growth rates were lowest in 2013, although this reduction was only significant for EW growth rates between 2012 and 2013 in control trees. Ring width and the LW/RW ratio decreased in trees subjected to rainfall exclusion, with differences between treatments becoming significant in 2014 (Figure 4). D_H did not differ between treatments or years in this species.

No significant differences were found between treatments or years in needle C/N ratio (Figure 5). However, needle Δ significantly increased along the study period in both treatments and was significantly lower in trees under rainfall exclusion than in control trees in 2012 and 2013. No differences in PNUE, A or g_s were observed between treatments, although g_s was marginally significantly lower (P = 0.067) in trees subjected to rainfall exclusion than in control trees in 2013 (Figure 5). PNUE was highest in 2013. A did not suffer significant inter-annual variation, although it was marginally significantly lower (P = 0.065) in 2012, the driest year of the three studied, particularly under rainfall exclusion. g_s greatly varied among years, especially in control trees, being lowest in 2012 and highest in 2013. Midday Ψ was significantly more negative in treated than in control trees in 2013. Leaf water potentials barely differed among measurement dates under the rainfall exclusion treatment but were significantly less negative in control trees in 2013 (Figure 5).

Xylem formation was significantly shortened under the rainfall exclusion treatment in P. sylvestris (Figure 2). Conversely, the onset of xylogenesis was not affected by rainfall exclusion in either of the two species studied, indicating that other factors could be involved in the resumption of cambial activity, such as temperatures or photoperiod (Rossi et al., 2008; Pérez-de-Lis et al., 2016). Thus, as hypothesized, the length of the growing season was shorter under drought, particularly in P. sylvestris, which confirmed its higher sensitivity to drought than the studied sub-Mediterranean oak. Concurring with our results, Klein et al. (2005) and Eilmann et al. (2011) found that trees had shorter growing seasons under drought due to an earlier cessation of growth rather than a later onset in Pinus halepensis and P. sylvestris, respectively. This earlier cessation of xylogenesis is consistent with a shift in the allocation of resources under stress, favoring carbon storage, bud formation or fine root growth over stem secondary growth (Waring, 1987; Eilmann et al., 2011; Wiley et al., 2013). On the other hand, Belien et al. (2012) and D’Orangeville et al. (2013) did not observe an effect of rainfall exclusion on cambial phenology in Picea mariana and Abies balsamea located in boreal forests.

Conversely, no differences were found in leaf phenology as a result of rainfall exclusion for either of the study species (Figure 2). Likewise, Adams et al. (2015) found that Juniperus monosperma was unaffected by a drought experiment, although they observed a delay in Pinus edulis leaf phenology under the treatment during an already dry year. Ogaya and Peñuelas (2004) found a delay in Arbutus unedo crown phenophases as a result of a drought experiment, but not in Quercus ilex and Phillyrea latifolia. Concurring with our study, Klein et al. (2005) found no changes in P. halepensis needle phenology due to irrigation, while the treatment strongly affected stem growth period length. Therefore, the severity of drought necessary to trigger shifts in leaf phenology seems to be species-specific (Peñuelas et al., 2004), but stronger than that modifying cambial phenology.

Unlike our hypothesis, no cumulative effect was observed in P. sylvestris anatomical traits during the 3-year rainfall exclusion; they were, however, influenced by year-to-year climatic variability (Figure 3). Control trees had a significantly higher number of tracheids than treated ones only in 2012, the driest year. Given the lack of differences in growth rates between treatments, this is probably due to tree-intrinsic characteristics, such as the number of cambial cells, and treatment-derived effects, such as the shorter growing season found in trees subjected to drought. Moreover, the number of tracheids was significantly lower in 2012 and 2013, the driest years. Experiments performed in different conifer species have found an increase in cell production in irrigated compared to non-irrigated trees (Abe et al., 2003; Eilmann et al., 2009; de Luis et al., 2011; Balducci et al., 2013). On the other hand, Belien et al. (2012) did not find differences in the number of tracheids formed under drought and D’Orangeville et al. (2013) only found differences in tree ring size on the second year of the experiment. This variability may suggest that the response of cambial production to drought is highly dependent on its severity and duration. Drought episodes have been found to strongly influence secondary growth even after several years, either by cumulative hydraulic damage (Anderegg et al., 2013) or by affecting the recovery of the canopy and its indirect effect on carbon reserves due to reduced photosynthetic carbon uptake (Galiano et al., 2011). Thus, a stronger effect of drought on tracheid production may be expected after prolonged exposure to rainfall exclusion.

Similarly to Balducci et al. (2013), we did not observe differences in P. sylvestris wood anatomy between treatments (Figure 3). Nonetheless, hydraulic diameter tended to be higher and cell-wall thickness lower in the driest year, particularly in control trees, which resulted in a greater proportion of earlywood, indicating that anatomical traits respond to climatic variability. By contrast, Belien et al. (2012) and D’Orangeville et al. (2013) found that the tracheids of trees under drought had smaller lumen size and thicker cell walls than those of control trees. Eilmann et al. (2009, 2011), however, also found greater lumen diameters under drought in P. sylvestris. Therefore, this species may tend to produce fewer tracheids under dry conditions without reducing water transport efficiency due to the greater proportion of large-lumen cells. This enhancement in hydraulic diameter may not, however, necessarily increase its vulnerability to cavitation as drought-induced cavitation mostly depends on pit membrane properties, rather than conduit diameter (Hacke and Sperry, 2001; Martínez-Vilalta et al., 2009).

On the other hand, the more-drought tolerant Q. pyrenaica did not show differences in earlywood vessel hydraulic diameter between treatments or years (Figure 4), whilst other studies in sub-Mediterranean oaks have found both wider (Corcuera et al., 2006) and narrower (Eilmann et al., 2009) earlywood vessels under drought. However, consistent with these studies, we found that Q. pyrenaica trees subjected to rainfall exclusion suffered an increasing reduction in latewood width along the study period (Figure 4). This hypothesized cumulative effect of the rainfall exclusion treatment might have been more evident in this species than in P. sylvestris due to the heavy reliance of ring-porous species on stored carbohydrates to form their annual rings (Michelot et al., 2012). This could have intensified the rainfall exclusion influence as reserves may have not been fully replenished throughout the experiment. Although rainfall exclusion did not affect the size of earlywood vessels and, thus, the capacity of the main conductive system remained unaltered, the reduction in latewood under drought could ultimately challenge hydraulic performance in this species as latewood vessels provide water transport for several years once earlywood vessels become non-functional due to cavitation (Granier et al., 1994).

The reduction in the number of tracheids formed by P. sylvestris in drier years suggests that the allocation of resources to xylem was reduced under drought. The observed increase in hydraulic diameter under drought could have, however, partially offset this reduction in carbon allocation to xylem formation by likely enhancing tracheid water transport efficiency. By contrast, Q. pyrenaica maintained earlywood vessel size but reduced latewood width under rainfall exclusion. Therefore, this species invested less carbon in xylem formation but at a greater risk of hydraulic failure, particularly in extremely dry years when latewood functionality could be critical to maintain hydraulic sufficiency (Granier et al., 1994; Corcuera et al., 2006). Thus, although with different effects on their hydraulic system, our observations suggest that both species followed carbon-saving strategies in xylem formation under drought. This is consistent with studies that indicate that trees adjust to drought by first limiting secondary growth rather than carbon assimilation (Waring, 1987; McDowell et al., 2008; Oberhuber et al., 2011; Fatichi et al., 2014; Adams et al., 2015). This limitation was more evident in the old P. sylvestris tree that died during the experiment, which barely had any cambial activity prior to its death while the crown was still functional. Albeit an isolated event within our plots, water stress-related mortality is not unprecedented within the study area (Gea-Izquierdo et al., 2014). The death of this tree under the rainfall exclusion treatment could, therefore, portend high vulnerability to drought-induced mortality in P. sylvestris at its dry limit (Martínez-Vilalta and Piñol, 2002; Vilà-Cabrera et al., 2011).

PNUE was unaltered by the rainfall exclusion treatment; it suffered, however, high year-to-year variability, suggesting a climatic influence on this variable (Figure 5). PNUE inter-annual changes followed those in g_s and, to a lesser extent, A. These variables were highest in 2013, despite 2014 was the wettest year, which may be due to REW was higher on the dates measurements were taken in 2013 than in the other two study years (Figure 1). The pattern followed by these variables is consistent with studies that show that decreases in g_s reduce PNUE, while increase intrinsic water use efficiency (iWUE; Limousin et al., 2015).

Needles of trees under the rainfall exclusion treatment had lower Δ and, thus, higher iWUE, than control trees (Figure 5). This effect of rainfall exclusion on the more integrative δ¹³C-derived iWUE was, however, less evident in instantaneous gas exchange measurements, which only showed marginally significant differences between treatments in A and g_s. Due to the non-linearity of the relationship between A and g_s, stomatal closure causes a proportionally greater decrease in transpiration than in photosynthesis (Limousin et al., 2015). This may explain the higher response to the rainfall exclusion treatment and inter-annual variability observed in g_s and Δ, both lowest in the driest year, than in A. Nonetheless, photosynthetic rates are tightly coupled to stomatal control of water use and, therefore, long-term stomatal limitations could significantly alter the tree’s carbon reserves (Galiano et al., 2011).

Midday leaf water potentials followed changes in g_s in control trees whereas they remained more negative in trees subjected to rainfall exclusion, being the difference between treatments significant in 2013 (Figure 5). Nevertheless, midday Ψ observed at our site were above those reported to cause significant xylem cavitation (Martínez-Vilalta and Piñol, 2002; Poyatos et al., 2008) and the absolute minimum values reported for the species (-2.5 MPa; Martínez-Vilalta et al., 2009). Hence, thanks to the tight control exerted by stomata to avoid excessive water loss and, consequently, prevent severe damage to the hydraulic system (Irvine et al., 1998; Poyatos et al., 2008), P. sylvestris may afford increasing current year tracheid size under dry conditions, as observed in this and other studies (Eilmann et al., 2009, 2011), without increasing the risk of drought-induced cavitation.

Rainfall exclusion and irrigation experiments have been central for improving our understanding of how tree physiology responds to drought (Beier et al., 2012). Our rainfall exclusion treatment was successful in decreasing soil water content (Figure 1). However, most of the studied traits showed non-significant differences between treatments, despite they responded to inter-annual rainfall variability. This lack of significant differences might be due to trees being able to cope with imposed drought by absorbing water from deeper soil layers (Belien et al., 2012), from roots reaching beyond the exclusion area or from intercepted rainfall (Breshears et al., 2008), minimizing the effect of the treatment. Additionally, the field and laboratory intensive requirements characteristic of cambial phenology, wood anatomy and gas exchange studies resulted in small sample sizes, which probably obscured statistically significant effects of the rainfall exclusion treatment (Metz et al., 2016). Despite these shortcomings, our study provides valuable information on how these traits respond to an intensification of natural droughts. Although the focus of this study was the response of cambial phenology and wood anatomy to drought in the two study species, with the additional insight provided by leaf phenology in both species and gas exchange measurements in P. sylvestris, complementing these measurements with other variables, such as leaf chlorophyll, partial dehydration or leaf area index, would provide a more integral understanding of tree response to drought. Nonetheless, the results reported showed that secondary growth responded more strongly to the treatment than the canopy variables assessed, which would be consistent with the shift in resource allocation under drought stress observed in other studies (Waring, 1987; Oberhuber et al., 2011; Fatichi et al., 2014).