We studied the timing of budburst of valley oak (Quercus lobata Née) at Hastings Reservation, central coastal California. Similar to other taxa, budburst was advanced by warmer temperatures. Over the 30-year study period, however, there were no significant trends in either air temperature or the timing of budburst, except during the 2014–2016 drought, during which the earliest budburst dates were advanced. Several individual tree characteristics correlated with budburst timing, including access to ground water, soil available phosphorus, and elevation, the effects of which were in turn correlated with winter microclimatic conditions of individual trees. Budburst timing was significantly related to both subsequent acorn production and radial growth; trees leafing out on or near the population mean for the year experienced greater radial growth and produced larger acorn crops than trees leafing out earlier or later than the mean. Differences in acorn production were due to both differences in phenology among trees and plasticity in the phenology of individual trees across years, while differences in radial growth were primarily due to plasticity in individual tree phenology. Valley oak phenology exhibits considerable variability; the extent to which this plasticity will help this keystone California species adapt to future climate change remains to be seen.
The connection of vegetation phenology to meteorological conditions has been recognized since at least the first half of the 18th century, when the French naturalist Reaumur (1735) proposed the method of cumulative thermal summation as a means of predicting the onset of budburst in temperate trees. Since then, the study of phenology has become an industry unto itself due to its role in reconstructing and predicting primary productivity in this era of global climate change (Chmielewski and Rötzer 2001; Badeck et al. 2004; Menzel et al. 2006, Piao et al. 2019). Nonetheless, many questions remain. In particular, although a close relationship between environmental conditions—primarily temperature—and budburst is well established in many taxa, other factors influencing phenology are poorly studied, despite evidence demonstrating considerable variability in phenological landscapes over small spatial scales both within and between species (Denéchère et al. 2021). Studies also indicate that genetics (Faticov et al. 2020; Papper and Ackerly 2021; Wright et al. 2021), elevation and soil type (Cole and Sheldon 2017), and even anthropogenic light pollution (ffrench-Constant et al. 2016) can affect tree phenology. Analysis of the factors influencing phenological differences both within populations and across years has the potential to improve our understanding of a wide range of ecological factors affecting plant populations, including seed production, disease susceptibility, and insect damage, as well as the likely effects of future climate change (Koenig et al. 2015; Pearse et al. 2015a; Bogdziewicz et al. 2017; Faticov et al. 2020).
We studied budburst phenology in a population of valley oak (Quercus lobata Née) over a 30-year period, 1991–2020. Prior work on this population has focused on the relationships among budburst timing, phenological synchrony, and acorn production (Koenig et al. 2012, 2015), and the relationship between budburst timing and herbivore damage (Pearse et al. 2015a). Here we look specifically at within-population and among-year variability in budburst. We address three issues: (1) variability in annual budburst across years and among individuals within the population; (2) factors correlating with budburst phenology; and (3) the relationship between budburst and subsequent acorn production and radial growth, including the relative importance of within- vs. between-tree effects.
Methods
Species and Study Site
Quercus lobata is a “white” oak in the subgenus Quercus (Pearse and Hipp 2009). It is winter deciduous, matures acorns in one year, and, like other members of this genus, wind-pollinated and an obligate outcrosser. Its range, entirely within California, encompasses foothill regions ringing the Central Valley below ∼1,800 m in elevation (Griffin and Critchfield 1972). Throughout this range, stands typically consist of relatively few large trees (White 1966; Bolsinger 1988). Our study site, Hastings Natural History Reservation (36°12′N, 121°33′W), is in the outer coast range of Monterey County in central California. This locality consists of typical coastal oak savanna dominated by valley, blue (Q. douglasii Hook & Arn.), and coast live (Q. agrifolia Née) oaks interspersed within a mixed deciduous forest and grassland matrix (Griffin 1974). The climate is Mediterranean with warm, dry summers and cool, wet winters. Based on 30 years of records taken between 1990 and 2019 at the Reserve headquarters, mean annual rainfall is 534 mm and mean annual temperature is 14.0°C.
Individual Tree Variables
Budburst for 25 individual Q. lobata was measured each year from 1991–1996, inclusive. The survey was restarted and expanded in 2003 to encompass 86 trees (including the original 25), chosen as part of a study of acorn production. Thus, the data analyzed here covered the 30 years (24 of which were surveyed) of 1991–1996 and 2003–2020. Diameter at breast height (dbh) and elevation were recorded for each tree.
Each spring we surveyed trees weekly for budburst starting between mid-January and 1 March, depending on the year. Budburst was used as a proxy for flowering and pollen production, which follow budburst by an average of 12.3 days (Koenig et al. 2012). Analyses were based on the day-of-year when a tree was determined to have undergone budburst, defined as the first date on which at least 5% of the tree had leafed out and turned green. Occasionally a tree leafed out prior to the first survey of the year. If 25-75% of the tree had leafed out prior to the first survey, we estimated budburst as having occurred 7 days prior to the first survey date (2.7% of tree-years). If >75% of the tree had leafed out as of the first survey, budburst was estimated to have been 14 days prior to the first survey date (0.9% of tree-years).
For each of the 86 trees we measured water availability, soil available nitrogen (N), soil available phosphorus (P), and winter microclimatic temperature. Water availability was measured by predawn xylem water potential in September 1991 using a pressure bomb (Knops and Koenig 1994). Although values vary from year to year depending on rainfall, differences in xylem water potential values are concordant among trees (Knops and Koenig 2000). Thus, although only measured once, these values provide a good index of overall differences in water stress among trees. Soil available N and P were estimated using four ion-exchange resin bags placed 5–10 cm in the soil underneath each tree from October 1992 to April 1993. Resin bags were subsequently analyzed for NO3, NH4, and PO4 with values expressed and mg/l effluent (Knops and Koenig 1997).
Acorn production was estimated each year in September using visual surveys. Two observers scanned different parts of the canopy of each tree and counted as many acorns as they could in 15 s (Koenig et al. 1994, 1996). Values were added together and ln-transformed to reduce non-normality. Radial growth was measured using dendrometers (Cattelino et al. 1986) established in 1994 and subsequently measured annually at the time of the autumn acorn survey. Annual dendrometer growth values were scaled within individual trees to a mean of 0 and standard deviation of 1.
Environmental Variables
Daily temperature (maximum and minimum) and rainfall were obtained from the Hastings Reserve weather station located within 3.5 km of all trees. As rainfall in the Mediterranean climate of the study area falls almost entirely between early November and late March, we summed rainfall for each year from 1 November to the mean budburst date of the survey trees; we refer to this variable as winter rainfall. Mean daily temperature in our study site is rarely below 0°C; consequently, using a cumulative degree day model was not feasible. Instead, we used an iterative linear predictor (Bush and Mosteller 1955) to quantify temperature, which mitigates problems associated with measuring values over a fixed period (Gienapp et al. 2005). The linear predictor weighted values according to the formula:
where λ(t) is the weighted mean temperature on day t and α the weighting factor. We started calculations on 1 December of the prior year and ended at the mean (for the population analyses) or individual tree budburst date. Calculations started on 1 December because valley oaks often retain acorns and leaves into November, but only rarely longer (Koenig et al. 2014).
λ weights each day's temperature by αn, where n was the number of days prior to budburst. Preliminary analyses comparing daily maximum, daily minimum, and daily mean ([maximum + minimum] / 2) temperature and values of α between 0.01 and 1.0 (α = 1.0 is tantamount to averaging mean daily values) indicated that the best predictor of budburst was using mean daily temperature and α = 0.98, which predicted budburst significantly better than averaging mean daily temperatures (ANOVA comparing α = 0.98 and α = 1; χ2 = 24.3; P < 0.001). Thus, each day's mean daily temperature prior to budburst was weighted 98% as heavily as the subsequent day. As an example, the λ-weighting factor 30 days prior to budburst was 0.9830 = 0.545. Statistical analyses used the weighted mean temperature values, but (unweighted) mean temperature values were used in figures to simplify interpretation of observed relationships.
Winter microclimate was measured using small, automated temperature recorders (iButtons; Maxim Integrated Products, Sunnyvale, CA) located on the north side of each tree ∼1.5 m above ground (Koenig et al. 2015). iButtons were programmed to record temperatures at 4-hr intervals each day starting at 04:00 hours. For the analyses performed here, we calculated the mean of the daily maximum and minimum temperatures for each tree between 1 December and 31 March for the years 2004–2019; the 08:00 temperature reading was not included due to logistic issues. We then averaged mean values across all years to yield an overall index of the winter microclimate of each tree. Preliminary analyses indicated that correlations with mean minimum values were much stronger than with mean maximum values, and thus we used only mean minimum values (referred to as winter microclimate) to analyze relationships between trees experiencing different microclimates and the other tree characteristics measured.
Data Analysis
Analyses were conducted in R 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria) using Cox proportional hazards models (procedures coxph and coxme in package coxme) and Pearson correlation coefficients. For the Cox proportional hazards models, we list the exponential of the regression coefficients for the fixed effects (exp(β)). exp(β) is the ratio of the hazards between two individuals whose values differ by one unit of the fixed factor when all other covariates are held constant. Thus, positive values indicate that the variable positively affects the response variable, while negative values indicate the reverse. In order to render effect sizes comparable within models, all variables were scaled to a mean of 0 and standard deviation of 1 prior to analysis.
Cox proportional hazards models in which the response variable was the mean, earliest, or latest budburst date of the year included weighted mean temperature, winter rainfall, and time (year) as fixed factors. These models tested the effects of temperature and rainfall on budburst of the population and whether there was a temporal trend in budburst timing over the length of the study. For individual trees, we calculated Pearson correlation coefficients between variables. We did this because some variables were measured only once for each tree (xylem water potential and soil nutrients) or were characteristics of the trees that stayed the same (elevation) or were concordant (dbh) during the study, and because correlations between microclimate and the tree variables helped interpret their relationships with budburst date. For these individual tree analyses, we used mean budburst date of trees calculated by averaging the day-of-year budburst date for each individual across the years when all trees were surveyed (2003–2020).
To test the effects of budburst date on radial growth and the subsequent acorn crop, we calculated budburst for each tree relative to the annual mean of all trees ([budburst of tree i in year x] – [mean budburst date of all trees in year x]) and used this value (both linear and quadratic) as fixed factors in mixed-effects models (procedure lmer in package lmerTest) in which the acorn crop or radial growth of the tree in year x was the response variable and tree ID was a random effect.
Finally, we performed analyses distinguishing whether the observed effects of phenology on growth and reproduction were due to differences among trees (between-tree effects) or plasticity in the phenology of individual trees (within-tree effects). For each tree, we calculated the mean budburst date across all years and the difference between a tree's observed budburst date in a particular year and its mean budburst date (centered budburst date). Including these two variables in a linear regression model as fixed effects distinguishes within-tree effects (centered budburst date) from between-tree effects (mean budburst date) (van de Pol and Wright 2009).
Results
Annual Budburst Timing
A summary of the earliest, latest, and mean budburst dates illustrates the variability observed among years, with mean budburst ranging from Day 48 (18 February 2015) to Day 93 (3 April 2012) (Fig. 1a). The earliest budburst date detected during the study was Day 14 (14 January 2015). Mean (unweighted) temperature from 1 December to the mean budburst date and winter rainfall were also highly variable (Fig. 1b). Cox proportional hazards models testing the effects of weighted mean temperature, rainfall, and year on budburst dates indicated that the mean, earliest, and latest budburst dates all were significantly advanced when temperatures were warmer (Table 1; Fig. 2a). Budburst date was generally delayed in wetter years, but the effect of rainfall was not significant (Table 1; Fig. 2b). Over the course of the study, there was no significant trend in the mean or the latest budburst date as indicated by the ‘year’ term. There was, however, a significant positive trend in the earliest budburst date (Table 1).
Fig1.
(a) Summary of annual mean budburst dates for valley oaks at Hastings Reservation, 1991–1996 and 2003–2020. Plotted are the mean (marked by X), standard deviation (thick lines), and earliest and latest (dashed lines) of budburst dates, 1991–1996 and 2003–2020. (b) Mean temperature (1 December of the prior year to mean budburst date; plotted values unweighted) and seasonal rainfall (1 October of the prior year to 30 March).

Table 1.
Results of Cox Proportional Hazards Models of Annual Budburst Dates of Valley Oak at Hastings Reservation, 1991–1996 and 2003–2020. Log-likelihood ratios and regression coefficients for three models in which the mean, earliest, and latest budburst date of the population were the response variables. All models included the mean weighted temperature, rainfall, and year as fixed effects, scaled (mean = 0, standard deviation = 1) prior to analysis. A regression coefficient larger than 1 means that a variable had a positive effect on the hazard; that is, it increased the probability of budburst. * P < 0.05; *** P < 0.001; other P > 0.05.

Fig. 2.
(a) Mean budburst date of valley oak at Hastings Reservation (day-of-year) vs. mean temperature (1 December to mean budburst date), 1991–1996 and 2003–2020. (b) Mean budburst date of valley oak vs. winter rainfall (1 November to mean budburst date). In both cases, each dot represents a year (n = 24 years).

Budburst Timing of Individual Trees
Correlations of mean budburst date and winter microclimate with water availability, soil nutrients, elevation, and dbh are summarized in Table 2. Trees growing in sites with warmer winter microclimates experienced earlier budburst (Fig. 3a). Budburst was also earlier among trees that had lower water availability and higher soil P (Fig. 3b); there were no significant relationships with soil N. Trees growing in warmer winter microclimate sites had less water availability and higher soil P (Fig. 3c).
Table 2.
Pearson Correlation Coefficients Between Mean Budburst Date (Left) and Winter Microclimate (Right) and Individual Tree Variables of Valley oak at Hastings Reservation, 1991–1996 and 2003–2020. * P < 0.05; ** P < 0.01; *** P < 0.001; other P > 0.05, n = 86 trees.

Fig. 3.
Scattergrams of variables measured for 86 individual valley oaks at Hastings Reservation. (a) Winter microclimate (2004–2019) vs. mean budburst date (2003–2020). (b) Mean budburst date vs. xylem water potential (xwp) and soil available phosphorus. (c) Winter microclimate vs. xylem water potential and soil available phosphorus. * P < 0.05; ** P < 0.01; *** P < 0.001.

Relationship Between Budburst and Fitness
Mixed-effects models revealed significant nonlinear effects of budburst on both a tree's subsequent acorn crop and radial growth (Table 3, top). In both cases, trees undergoing budburst relatively early and relatively late produced smaller acorn crops and experienced less radial growth than trees whose phenology matched the mean of the population for the year (Fig. 4).
Table 3.
Results of Models Testing the Effects of Relative Budburst Date (Linear and Squared terms) (Top) and Within- vs. Between-tree Effects of Phenology (Bottom) on Reproduction and Growth of Valley Oak at Hastings Reservation, 1991–1996 and 2003–2020. Acorn crop was ln-transformed; radial growth was standardized within trees. Tree ID included as a random effect in the models testing for the effects on reproduction and growth. Effect sizes multiplied by 100. ** P < 0.01; *** P < 0.001; other P > 0.05.

Fig. 4.
(a) Scattergram of the (ln-transformed) acorn crop vs. the prior spring's relative budburst date (budburst relative to the mean budburst date of valley oaks that year). (b) Scattergram of standardized radial growth vs. the prior spring's relative budburst date. Both lines plotted using predicted values from regressions including both linear (budburst date) and quadratic terms (budburst date squared). For statistical tests, see Table 3.

In models including both between-tree (mean budburst date) and within-tree (centered budburst date) effects, both contributed significantly to the effects of phenology on the subsequent acorn crop (Table 3, bottom). Only the mean within-tree effect (centered budburst date) exhibited a significant relationship with radial growth.
Discussion
Budburst is a relatively easily measured indicator of climate effects (Cleland et al. 2007). In order to interpret the role that climate plays in affecting budburst, however, it is important that other factors influencing phenology are recognized. Although a variety of genetic, environmental, and anthropogenic factors have been found to affect phenology in various tree species, few studies have examined the effects of multiple factors on budburst. Phenology can also have important effects on the fitness of plants, affecting the timing of species interactions (Yang and Ruidolf 2010), resource acquisition (Nord and Lynch 2009), insect herbivory (Pearse et al. 2015a, 2015b), and reproduction (Galen and Stanton 1991; Koenig et al. 2012). The role of phenology as both an indicator of change and as a driver of fitness has made it the focus of a wide range of ecological studies.
We measured the timing of budburst over three decades on a population of 86 Quercus lobata for which we had data on temperature, water availability, soil nutrients, size, growth, and reproduction. As found in numerous prior studies, temperature had a strong effect on budburst, with warm temperatures advancing the budburst date of the population (Fig. 2a) and of individual trees (Fig. 3a). Using the weighted mean temperature variable, each experienced unit of λ(t) increased the probability of population budburst by over an order of magnitude (Table 1). Rainfall was not significant in a multivariate analysis including temperature, although the correlation between annual rainfall and mean budburst date was positive and nearly significant (Fig. 2b). This most likely reflects rainfall being inversely correlated with maximum temperatures: over the same 24 years for which phenology was surveyed, there was a significant inverse correlation between mean maximum temperature and rainfall during the winter (November through March) (Pearson correlation r = –0.51; P = 0.01; n = 24 years). Thus, wetter conditions correlate with cooler, more moderate temperatures at the study site and elsewhere in the Mediterranean regions of California, making the effects of temperature and rainfall on phenology difficult to distinguish (Gerst et al. 2017; Armstrong-Herniman and Greenwood 2021).
Over the time period of the study there was little overall trend in the timing of budburst, with only the earliest budburst date advancing significantly (Table 1). Examination of Fig. 1a suggests that this was most likely due to the extreme drought years of 2014–2016 (Griffin and Anchukaitis 2014; Luo et al. 2017), during which winter rainfall was low and temperatures high.
Overall, annual mean temperatures at the study site did not increase significantly during the years of the survey (Pearson correlations between year and annual mean maximum, annual mean minimum, and annual mean average temperatures = 0.20 [P = 0.37], 0.24 [P = 0.27], and 0.36 [P = 0.09], respectively). Models, however, suggest that climate change is likely to have significant effects on California oak distributions in the future, including valley oak, primarily due to warming temperatures rather than changes in precipitation (Kueppers et al. 2005; Loarie et al. 2008). Thus, although we detected no clear trend over the course of this study, valley oak phenology will potentially advance in the future. Although valley oaks are winter deciduous, the relatively moderate climate occasionally allows acorns of this species to remain on the trees throughout the winter (Koenig et al. 2014). In addition, the length of time that valley oaks are without leaves at our study site is as short as 3 months (mid-November/early December through mid-February/early March). We did not record timing of leaf senescence in this study, but both budburst and senescence will potentially be influenced by climate change (Firmat et al. 2017; Zani et al. 2020), affecting the length of the growing season and fitness via seed production (Journé et al. 2021) and susceptibility to insect herbivory the following spring (Karban 2007).
Other than temperature, we detected several variables correlating with budburst of individual trees, including water availability, soil P, and elevation (Table 2). In general, valley oaks avoid drought by tapping into ground water on alluvial terraces (Knops and Koenig 1994). Nonetheless, trees with relatively poor access to ground water, and thus under greater water stress, leafed out significantly earlier than trees with greater access to water (Fig. 3b). This is likely a consequence of trees with less access to ground water corresponding to warmer sites (Fig. 3c). Such trees will be more water-limited during the dry, Mediterranean summers that limit carbon gain (Hollinger 1992). Consequently, they will be more resource limited, and may leaf out earlier in order to use more surface soil water left over from the wet winter in order to compensate for lower carbon gain in the summer. Gene transcriptional responses to water stress in valley oaks affect over half its genome, indicating considerable opportunity for local adaptation that is likely to be the proximate driver of this relationship (Gugger et al. 2016).
Soil available phosphorus also correlated significantly with mean budburst date of individual trees. Trees growing in soils with more P, and thus enhanced availability of soil nutrients, leafed out earlier than trees with less access to P (Fig. 3b). This is also likely to be a consequence of microclimatic differences among sites. Soil microbial activity that causes litter decomposition is temperature dependent and warmer sites have faster rates of litter decomposition (Meentemeyer 1978). Although we do not have litter decomposition data, our results are consistent with previous work demonstrating that phosphorus is largely released from decomposing litter, whereas nitrogen is not (Knops et al. 2010). Thus, warmer sites with faster litter decomposition should have higher P, but not necessarily higher N. Mean budburst date of individual trees also correlated with elevation, with trees at higher elevation growing in sites with warmer winter microclimate and leafing out earlier. Tree size (dbh) did not correlate with budburst; previous work has shown considerable variability among species in terms of the effect of tree size on phenology (Marchand et al. 2020).
We found strong effects of budburst phenology on subsequent growth and reproduction of trees (Table 3). Confirming the convex-shaped trend found in earlier analyses (Koenig et al. 2012), trees leafing out relatively early and relatively late compared to the population mean had significantly reduced acorn crops the following autumn (Fig. 4a), a relationship consistent with pollen limitation playing an important role in acorn production (Pesendorfer et al. 2016). More surprising was that radial growth was also greater among trees leafing out at or near the population mean budburst date (Fig. 4b). Thus, trees leafing out early in the season apparently did not benefit by acquiring more resources allowing them to either grow more or produce more acorns, as might be predicted if such trees are able to photosynthesize for a greater length of time. Countering this line of reasoning, previous work found that trees in this population leafing out earlier experienced greater rates of folivore damage than trees setting leaves later (Pearse et al. 2015b). In addition, early budburst entails a greater risk of frost harming flowers and subsequent seed production (Fowells and Schubert 1956). Clearly there are multiple factors affecting the relationship between phenology and fitness. In any case, we found strong selection against trees leafing out early (or late) within a year as evidenced by their relatively poorer radial growth and smaller subsequent investment in reproduction. The explanation for the greater radial growth in trees leafing out around the mean population budburst date remains to be determined.
In an analysis quantifying the relative importance of between-tree effects and within-tree plasticity, both significantly contributed to differences in acorn production (Table 3). Trees that leafed out at or near the population mean produced larger acorn crops than trees leafing out relatively early or late (between-tree effect). Among years, however, individual trees produced larger acorn crops when they leafed out at or near the population mean for the year compared to when they did not (within-tree effect). Only the latter within-tree effect, indicative of individual plasticity, contributed significantly to the effect of budburst timing on radial growth. These results support analyses using common garden plantings that have found considerable environmental plasticity in phenotypic variation in both valley oak (Wright et al. 2021) and the largely sympatric blue oak (Papper and Ackley 2021). Budburst phenology is yet another character of valley oaks, along with acorn size (Koenig et al. 2009), sex allocation (Knops and Koenig 2012), and a range of other life-history traits (Barringer et al. 2013), that exhibits wide variability both across years and among individuals.
Valley oaks have been aptly described as “. . .the monarch of California oaks by virtue of [their] size, age, and beauty” (Pavlik et al. 1991). Unfortunately, this species has suffered inordinately over the past 150 years due to agricultural and residential development and stands to lose an estimated 54% of its modern potential range as a consequence of climate change by the end of the 21st century (Kueppers et al. 2005). Phenology is thought to play a key role, not only as an indicator of climate change, but as a major determinant of plant species ranges (Chuine and Beaubien 2001). It remains to be seen whether the considerable variability in valley oak phenology documented here will help this species adapt to the changing conditions of California's future.
Acknowledgments
We thank the reviewers for comments, the field assistants who helped gather data, and Dick Sage for talking me into studying phenology more seriously. Financial support came from the University of California's Integrated Hardwoods Range Management Program, the California Department of Fish and Game, and National Science Foundation grants DEB-0816691 and DEB-1256394 to WDK. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.