Forestry, Faculty of
DSpace
UBCV
Tysor, Carolyn Susannah
Master of Science - MSc
University of British Columbia
The timing and duration of pollen shed and cone receptivity in lodgepole pine (Pinus contorta) influence its fecundity, levels of gene flow, and local adaptation. Accurate predictions of spring reproductive (“flowering”) phenology are essential for understanding local adaptation, gene flow, population structure, and climate change impacts. Using 16 years of data from lodgepole pine sourced from across British Columbia and grown at seven seed orchard sites in an approximate common garden design, I developed a multilevel Bayesian model to assess the forcing requirements for flowering onset and cessation, and the provenance effects on these requirements for both male and female strobili. The model reveals a countergradient pattern in flowering: populations from warm locations require more forcing by warm temperatures than those from colder locations. This countergradient effect should result in increased overlap among populations at the landscape level but decreased overlap within individual seed orchards. Pollen shed typically begins shortly before and ends slightly after receptivity. By applying the model to historical and predicted daily temperature records, I estimated flowering phenology for diverse provenances at multiple sites from 1945 to 2100, translating forcing predictions into specific days of the year. Historical flowering patterns are variable but show significant correlations across distant sites. Under moderate climate change scenarios, flowering is predicted to advance by two to three weeks without having appreciable impacts on synchrony.
Countergradient variation inlodgepole pine flowering phenologybyCarolyn Susannah TysorB.Sc., Northern Arizona University, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THEREQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Forestry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)January 2025© Carolyn Susannah Tysor, 2025The following individuals certify that they have read, and recommend to the Faculty of Graduate andPostdoctoral Studies for acceptance, the thesis entitled:Countergradient variation in lodgepole pine flowering phenologysubmitted by Carolyn Susannah Tysor in partial fulfillment of the requirements forthe degree of Master of Sciencein ForestryExamining Committee:Sally Aitken, Professor, Department of Forest and Conservation Sciences, UBCSupervisorElizabeth Wolkovich, Associate Professor, Department of Forest and Conservation Sciences, UBCSupervisory Committee MemberMichael Whitlock, Professor, Department of Zoology, UBCSupervisory Committee MemberSean Michaletz, Assistant Professor, Department of Botany, UBCAdditional ExamineriiAbstractThe timing and duration of pollen shed and cone receptivity in lodgepole pine (Pinus contorta) influence itsfecundity, levels of gene flow, and local adaptation. Accurate predictions of spring reproductive (“flowering”)phenology are essential for understanding local adaptation, gene flow, population structure, and climatechange impacts. Using 16 years of data from lodgepole pine sourced from across British Columbia andgrown at seven seed orchard sites in an approximate common garden design, I developed a multilevel Bayesianmodel to assess the forcing requirements for flowering onset and cessation, and the provenance effects on theserequirements for both male and female strobili. The model reveals a countergradient pattern in flowering:populations from warm locations require more forcing by warm temperatures than those from colder locations.This countergradient effect should result in increased overlap among populations at the landscape level butdecreased overlap within individual seed orchards. Pollen shed typically begins shortly before and endsslightly after receptivity. By applying the model to historical and predicted daily temperature records, Iestimated flowering phenology for diverse provenances at multiple sites from 1945 to 2100, translating forcingpredictions into specific days of the year. Historical flowering patterns are variable but show significantcorrelations across distant sites. Under moderate climate change scenarios, flowering is predicted to advanceby two to three weeks without having appreciable impacts on synchrony.iiiLay SummaryThe timing of flowering influences adaptation and where species occur because it affects plants’ ability toreproduce and which other plants they reproduce with. Lodgepole pine is an important species in BritishColumbia’s forests and forestry industry. After compiling a large data set, I modelled how much warmthlodgepole pine needs to “flower” (shed and receive pollen) in the spring. I found northern populationsrequire less warmth to flower than southern populations, so lodgepole populations tend to reproduce withothers in similar climates, which supports local adaptation. However, individual trees vary more thanpopulations, balancing local adaptation with genetic diversity within populations. With climate change,lodgepole are expected to flower several weeks earlier by 2100 than they do today without strongly alteringwhich populations can reproduce with each other. The predictive model developed here will enable moreresearch into how flowering timing affects lodgepole ecology and evolution.ivPrefaceThis project takes advantage of existing lodgepole pine phenology data collected over a decade and a halfby government and industry workers in seed orchards in British Columbia. The data was provided to meby Rita Wagner (Prince George Tree Improvement Station), Chris Walsh (Kalamalka seed orchards), andJoe Webber (Glynn Road Research Station). Rob Johnstone assisted in the transcription of handwrittendatasheets and I compiled, cleaned, and harmonized the phenology data into a usable data set.I was responsible for developing the research questions, designing and performing the modeling,analyzing the results, and writing the thesis under the supervision of Dr. Sally Aitken. My entire supervisorycommittee provided input into modelling and analysis decisions and feedback on interpretation of results aswell as detailed edits and comments on the thesis itself.Sarah Bisbing and Conor Phelan provided the shapefile for the Pinus contorta ssp. latifolia rangemap in Figure 2.1.The generative AI tool ChatGPT was used to:1) create custom legends for several figures and adjust figure aesthetics,2) help debug analysis code, and3) create a reverse outline of an earlier draft of this thesis which I used (without generative AI) to editthe structure of the document.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Study Species: Lodgepole Pine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Countergradient variation in lodgepole pine flowering phenology . . . . . . . . . . . . . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Study species . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.3 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.5 Analysis and prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.2.6 Data and code availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.3.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.3.2 Average forcing requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.3.3 Provenance MAT effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.3.4 Factor impacts and model error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.3.5 General growing degree day predictions . . . . . . . . . . . . . . . . . . . . . . . . . . 302.3.6 Day of year predictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.4.1 Comparison of dates and GDD with other studies . . . . . . . . . . . . . . . . . . . . 362.4.2 Relationship between pollen shed and receptivity . . . . . . . . . . . . . . . . . . . . . 372.4.3 Individual variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.4.4 Countergradient variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.4.5 Climate change impacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 422.4.6 Additional strengths and weaknesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Appendix A - Supplementary information for Chapter Two . . . . . . . . . . . . . . . . . . . . . . 57Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60viList of Tables2.1 Proportion of observations for each event interval censored or left or right end censored. . . . 132.2 Flowering state, symbols used to record them in the raw data, and phenophase definitions. . . 132.3 Table of intercepts 𝛼 and slopes 𝛽 for each site. . . . . . . . . . . . . . . . . . . . . . . . . . . 142.4 Comparison of model estimates to observations. Table entries are the proportion of meanmodel retrodictions contained in the observed range for each flowering event for forcing accu-mulation and day of year predictions, so 1 is perfect prediction and 0 is no overlap between theobserved range and modeled flowering GDD or day of year. In parentheses are the proportionof one 𝜎 intervals around mean model retrodictions that overlap the observed range. . . . . . 262.5 Mean forcing requirements and 95% HDI for beginning and end of pollen shed and conereceptivity for trees sourced from a provenance with a MAT of 0°C and ignoring the effectsof site, year, genotype, and tree. (Estimates correspond to thermal time model intercept 𝜇.)Based on 6000 draws from the posterior. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.6 Accumulated forcing requirements increase for warmer provenances. Slope (𝛽) estimates and95% highest density posterior intervals are from thermal time models for seed cone receptivityand pollen shed phenology in lodgepole pine. Forcing calculated in growing degree days withthreshold 5°C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28A2.1 Site metadata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60A2.2 Geographic distance is a significant predictor of similarity in flowering patterns across seedorchard sites. After stratifying by sex and flowering event, flowering day of year at each sitewas ranked from earliest to latest between 1945 and 2012 and pairwise correlations betweensites calculated as Kendall’s 𝜏 . R2 values and slopes for a linear model of Kendall’s 𝜏 as afunction of distance between sites illustrating the statistically significant negative relationshipbetween distance and rank correlations between sites. Slope in units of Kendall’s tau per km. 60viiList of Figures2.1 Site and provenance map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Observation dates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Site temperatures and forcing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4 Seed orchard design (factor replication). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.5 A conceptual diagram of how model retrodictions were compared to observations. . . . . . . . 212.6 Comparison of observations to model predictions. . . . . . . . . . . . . . . . . . . . . . . . . . 252.7 Model versus independent data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.8 Thermal time model offsets and variation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.9 Growing degree day predictions across provenance range. . . . . . . . . . . . . . . . . . . . . . 312.10 Day of year predictions 1945-2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.11 Provenance differences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.12 Geographic distance is a significant predictor of similarity in flowering patterns across seedorchard sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.13 Effects of countergradient variation on flowering overlap in seed orchards. . . . . . . . . . . . 352.14 Flowering phenology predictions under climate change. . . . . . . . . . . . . . . . . . . . . . . 36A2.1 Mean monthly temperatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57A2.2 Receptivity day of year predictions 1945-2011, all sites. . . . . . . . . . . . . . . . . . . . . . . 58A2.3 Pollen shed day of year predictions 1945-2012, all sites. . . . . . . . . . . . . . . . . . . . . . . 59viiiAcknowledgementsDepending on how you measure things, this thesis took 13 years or 3. I cannot thank enough the peoplewho understand crip spacetime who helped me build a path I could walk through an inhospitable institution,most especially Sandy and Caroline.My supervisor was patient when I needed it and her love for trees and forests kept me focused onthe big picture no matter how frustrated or confused I got. My entire supervisory committee - past andpresent - encouraged me, improved the research, and kept me engaged with their questions and insights. Ihad many wonderful lab mates and colleagues, whose interest and shared joy of discovery kept me motivatedand helped me pick through tricky problems. I am especially grateful for the friendship, understanding, andBayesian camaraderie of Joane Elleouet, Andrew MacDonald, and Dave Harris.Richard MacElreath, Michael Betancourt, and many contributors to the Stan forums taught methe statistics and programming I needed to build the phenology model at the center of this work.Seed orchard workers laboured for 16 springs shaking branches and inspecting seed cones, and theircareful records became an extraordinary phenology dataset that was a privilege to work with. I hope that I’veleft it in such a state that many others can benefit from their work. Jack Woods, Rita Wagner, Chris Walsh,and Joe Webber helped me obtain and interpret the seed orchard data this project relies on. I benefitedgreatly from their long years of work and hard-earned experience in BC forestry. Rob Johnstone transcribedtens of thousands of phenology observations from field sheets and scans before neural nets were up to suchcomplex OCR tasks. Sarah Bisbing and Conor Phelan provided the shapefile for the Pinus contorta ssp.latifolia range map.I would also like to thank GySgt T. Kerr for her example of integrity, grit, and good humour indifficult circumstances.Most of all, I am grateful for my partner and caregiver, Graydon. Without his financial, bodily,intellectual, technical, and emotional support, this thesis would not exist.ixDedicationTo my kith and kin who did not live to see this work finished, but walked part of this path with me.J.J.T., D.B., C.M., B.T., H.C., K.N., A.J., D.C., S.V., J.R.x1 IntroductionPlants must synchronize their life cycles with seasonal patterns of climate and resource availability to survive.Proper timing of life cycle events, such as flowering and budburst, allows plants to maximize resourceuse during favorable periods and avoid damage in harsh conditions. Reproductive phenology, particularlyflowering time, influences both gene flow and demography, shaping adaptation and species distributions(Chefaoui et al., 2019; Iler et al., 2021). In temperate regions, spring phenological events like budburstand flowering are primarily driven by temperature, and warming due to climate change is leading to theiradvancement (Piao et al., 2019).Temperate and boreal forest trees often exhibit strong clinal variation in traits such as growth andcold hardiness across their ranges (Savolainen et al., 2007). Clinal variation can be co- or counter- gradient(Conover & Schultz, 1995). Co-gradient variation occurs when the effect of selection on genotype aligns withthe environmental influence on phenotype, driving phenotypic differences along the cline. Counter-gradientvariation occurs when the effect of selection on genotype opposes the environmental effect on phenotype,driving phenotypic similarity along the cline. Since flowering phenology is primarily driven by temperature,populations along temperature gradients develop assortative mating, where individuals from similar thermalenvironments are more likely to mate. This restricted gene flow, due to assortative mating, leads to co-gradient variation even in the absence of selection (Soularue & Kremer, 2012). Selection can intensify theco-gradient variation or oppose the influence of gene flow to create countergradient variation (Soularue &Kremer, 2014). This clinal variation, in turn, influences how populations respond to environmental changeby either amplifying or mitigating its effects on phenology. As climate change causes phenological events,like flowering, to advance, it may disrupt existing patterns of gene flow and alter adaptive dynamics in foresttree species.As climate change shifts suitable habitats to higher elevations and latitudes, plant populations musteither migrate or adapt to persist (Aitken et al., 2008). Populations at the trailing edges of species’ rangesare expected to face increased mortality and regeneration challenges (Parmesan, 2006), while those in therange interior may struggle to keep pace with changing conditions (Davis & Shaw, 2001). Leading-edgeexpansion may be constrained by seed dispersal, maladaptation, and pollen limitation (Nathan et al., 2011;Zhu et al., 2012). While wind pollination enables long-distance gene flow, its success depends on the timingof flowering aligning with pollen dispersal, as gene flow only occurs when pollen reaches receptive flowers.Flowering phenology, by influencing both fecundity (Inouye, 2008) and gene flow (Franks & Weis,2009; Rivest et al., 2021), plays a key role in shaping spatial genetic structure, adaptation, and range limits1(Chuine & Beaubien, 2001; Kitamoto et al., 2006; Soularue & Kremer, 2014; Franks, 2015). Accurate predic-tions of flowering phenology are necessary for development of successful assisted gene flow and colonizationstrategies in conservation, as phenological mismatch can hinder population persistence and prevent gene flowbetween native and introduced populations (Wadgymar & Weis, 2017).Phenological stages in plants are difficult to observe directly because many developmental processesare hidden. Researchers rely on visible events, like budburst or flowering, to infer unseen stages. Calendardates for events are often converted into thermal time, which better reflects the effects of temperaturevariability on development (Trudgill et al., 2005). Temperature is the primary driver of spring phenology intemperate trees, with development proportional to the accumulation of heat over time (Sarvas, 1972). Oncea critical threshold of thermal time is reached, events such as flowering occur. However, accurately defininghow forcing accumulates is challenging, as developmental mechanisms remain unclear and vary among species(Chuine, 2000). A common approach uses a simple linear model where the daily contribution to thermaltime is the difference between the average daily temperature and a base threshold (Chuine & Régnière, 2017).These thermal time models are more reliable than calendar-based models for predicting phenology acrossyears and locations.1.1 Study Species: Lodgepole PinePinus cortorta (lodgepole pine) is an economically and ecologically important tree species facing multiplethreats from climate change (Hamann & Wang, 2006; Schneider et al., 2010; Sambaraju et al., 2012). Itsbroad geographic range, spanning 33∘ latitude and 31∘ longitude, covers diverse climates and soils withwidespread and significant local adaptation in many traits; while populations are often not growing in theiroptimum climate, survival and growth decline dramatically if trees are transferred too far from their homeclimates (Rehfeldt et al., 1999; Wang et al., 2010). Lodgepole pine are also adapted to geographic patterns offire throughout their range. Cone serotiny, an adaptive trait which allows for germination after fire events, ismuch more common in Rocky Mountain populations (Morgenstern, 1996) where fire is frequent. At smallerscales within the Rocky Mountain region, cone serotiny is predicted well by the type of disturbance thatinitiated the stand; serotinous cones are predominant in fire-initiated stands and open cones in other types(Muir & Lotan, 1985), with sharp boundaries between the two (Lotan, 1976).Pollen movement is a key vector of gene flow in lodgepole pine, as the species is predominantlyoutcrossing, with extensive pollen dispersal and relatively limited seed dispersal (Ennos, 1994). However,flowering phenology constrains population connectivity. Predicting pollen shed and receptivity would allowcharacterization of that constraint. Fortunately, high-quality reproductive and vegetative phenology data2have been collected from seed orchards across British Columbia, where lodgepole pine is cultivated forreforestation. This long-term data set is invaluable for modeling phenology across a wide geographic range.This study focuses exclusively on Pinus contorta ssp. latifolia, the most widespread subspecies. Thesmaller ranges and limited economic significance of the other subspecies (contorta, murrayana, and bolanderi)reduce the availability of phenology data. From this point forward, “lodgepole pine” refers specifically to P.contorta ssp. latifolia unless otherwise noted.1.2 Research questionsGiven the importance of gene flow in maintaining genetic diversity and shaping local adaptation, understand-ing flowering phenology in lodgepole pine is crucial not only for predicting future range shifts and adaptivepotential under climate change but also for explaining current patterns of spatial genetic structure and localadaptation. In this study, I used flowering phenology and temperature data to develop a multilevel Bayesianthermal time model for predicting pollen shed and cone receptivity in lodgepole pine. The data set spans16 years and includes 149 provenances grown in seven seed orchards across British Columbia from 1997 to2011. This study addresses the following questions:1) How much forcing is required for flowering events in lodgepole pine?2) Does provenance change forcing requirements?3) How will climate change affect lodgepole flowering phenology?32 Countergradient variation in lodgepole pine flowering phenology2.1 IntroductionForest tree species composing the vast temperate and boreal forests of the northern hemisphere often haveextensive ranges over large environmental gradients. Wind pollination is common in these trees, which canlead to gene flow over very long distances (Kremer et al., 2012; e.g. Campbell et al., 1999) if floweringtimes overlap. The timing of reproductive events influences both fecundity and gene flow, so understandingflowering phenology advances our understanding of range limits, spatial genetic structure, and adaptationin these species and has implications for conservation strategies like assisted gene flow (Chuine & Beaubien,2001; Inouye, 2008; Soularue & Kremer, 2012; Wadgymar et al., 2015; Gauzere et al., 2020).In addition to the evolutionary implications of flowering phenology, pollen shed from anemophiloustrees has a poorly studied but likely important role in forest nutrient cycling; forests produce enormousamounts of pollen that is high in plant macronutrients and decomposes rapidly, producing a reliable nutrientpulse in spring when many species are leafing out and in active growth (Lee & Booth, 2003). Floweringphenology is also relevant to forestry professionals making decisions related to seed production and treebreeding in orchards (Slavov et al., 2005), as well as reforestation on the landscape.While wind pollination results in extensive gene flow in trees (summarized in Savolainen et al.(2007), Kremer et al. (2012)), common garden experiments reveal extensive phenotypic and genetic clinesalong environmental gradients (reviewed in Alberto et al. (2013)). When the effect of selection on genotypealigns with the environmental influence on phenotype, co-gradient variation results in phenotypic differencesalong the cline. Countergradient variation, on the other had, can mask underlying genetic divergence byproducing similar phenotypes along the cline. While clinal variation can develop solely as a result of femaleflowering phenology filtering pollen along an environmental gradient (Soularue & Kremer, 2012), it mayalso result from selection (Soularue & Kremer, 2014) and lead to local adaptation. In northern hemisphereconifers, local adaptation is common (Leites & Benito Garzón, 2023) and often reflects strong trade-offsbetween avoidance of cold damage and competitive height growth (Howe et al. (2003)). Clinal variation inflowering is of particular interest because of its influence on gene flow and its correlation with and potentialimpact on variation in other traits.As the climate changes, the overall timing of spring phenological events will advance, especiallyfor plants active in rapid seasonal transitions and short growing seasons (Pau et al., 2011), like many highelevation and latitude trees. Significant advances in leaf out have already been observed (Linkosalo et al.,2009; Piao et al., 2019). Flowering phenology is largely driven by temperature, with flowering occurring4after a certain heat sum (accumulation of forcing temperatures) and is likely to be strongly affected bya warming climate. When populations occur along a temperature gradient, this temperature dependencedrives assortative mating that can influence evolution and lead to the development of co-gradient variationin other traits (Soularue & Kremer, 2012). Flowering phenology may also exhibit clinal variation itself whichcan amplify or counteract the environmental effect on flowering timing (Gauzere et al., 2020).Being able to predict flowering phenology improves evolutionary models and simulations and theirpredictions under climate change (Kremer et al., 2014; Duputié et al., 2015), but we lack knowledge ofvariation in reproductive heat sum requirements in many species. While there is extensive research on budburst from field observations and experiments (Alberto et al., 2013), flowering phenology research in treesis more limited because of time to and height at reproductive maturity. Research is especially scarce inconifers perhaps due to their less showy reproductive structures. However, bud burst timing is correlatedwith flowering phenology in many tree species and also primarily driven by temperature so that, on thelandscape, southern and low elevation populations break bud and flower earlier than northern and highelevation populations. While spring phenology in forest trees is often highly plastic, co- or counter- gradientvariation are also sometimes revealed in common gardens. Quercus petraea andQ. robur both show cogradientvariation, with budburst from southern populations occurring earlier than northern populations even in acommon garden, while environmental and genetic variation counter each other in Fagus sylvatica, leading toearlier bud burst of northern than southern populations in common gardens (Vitasse et al., 2009). Knowingwhether there are genetic effects or solely environmental effects is necessary for understanding the impactsof climate change on flowering phenology and their implications for evolution.Lodgepole pine (Pinus contorta) is an economically and ecologically important tree species facingmultiple threats from climate change (Hamann & Wang, 2006; Schneider et al., 2010; Sambaraju et al., 2012).It has a very large geographical distribution (across about 33∘ latitude and 31∘ longitude) encompassing awide range of climates and soils (Fig. 2.1) with widespread and significant local adaptation in many traits.For example, populations from both northern interior British Columbia and northern Idaho can survive inareas with mean annual temperatures between -4 and 6 ∘C, but the northern British Columbia populationsurvives and grows best where mean annual temperatures are ~ 1 ∘C and the Idaho population grows bestat ~ 4 ∘C (Rehfeldt et al., 1999). Local adaptation in lodgepole pine can be observed even at relativelysmall spatial scales when topographic variability is high: in a reciprocal transplant experiment, growthdeclines were observed when moving high elevation populations just 100m in elevation (Rehfeldt, 1983).Weak countergradient variation in growth initiation, which is likely to be correlated with flowering, has been5found in common garden and growth chamber experiments in lodgepole pine (Liepe et al., 2016; MacLachlanet al., 2017).While many traits of lodgepole pine show evidence of local adaptation, gene flow and within popula-tion genetic variation are high (Bisbing et al., 2021). Pollen is an important vector for gene flow in lodgepolepine because outcrossing is common, pollen dispersal is extensive, and seed dispersal is relatively limited(Ennos, 1994). There is evidence for spatially varying levels of gene flow in the species as populations fromareas with higher regional climate heterogeneity have higher genetic variance (Yeaman & Jarvis, 2006). Abetter understanding of flowering phenology could improve our understanding of past and future gene flowin lodgepole pine and the balance between selection and gene flow.Work on lodgepole pine flowering phenology is scarce and typically limited to presenting tablesof flowering dates (e.g. Owens et al. (2005)) or broad statements like “in May or June” (Owens, 2006).Nilsson (1981), however, reported on flowering abundance and phenology for fourteen provenances grown attwo sites observed over four years in Sweden, and O’Reilly & Owens (1988) observed one year of floweringphenology for seven provenances at a single site in British Columbia. While Owens (2006) reported thatlodgepole pine pollen shed and receptivity occur when degree days reach about 500, Nilsson (1981) foundthe flowering period ranged from about 161 to 247 GDD. Though Nilsson (1981) identified that provenancelatitude and altitude were associated with variation in flowering phenology and O’Reilly & Owens (1988)noted that southern provenances flowered after northern provenances at the same site, variation was notquantified in a way that was broadly applicable.To model pollen-mediated gene flow or how reproductive phenology could influence range limits,more precise estimates of lodgepole pine flowering phenology are needed.To address this need, I used flowering phenology and temperature data to fit a multilevel Bayesianthermal time model predicting pollen shed and cone receptivity in Pinus contorta ssp. latifolia. A thermaltime model estimates the amount of forcing required for a phenological event, which can then be translatedinto day of year. Phenological models that predict day of year directly are less able to make predictions dueto the variability of temperature on a yearly basis.I assembled and organized a 16-year flowering phenology data set from 149 provenances grown in 7lodgepole pine seed orchards in British Columbia. These seed orchards, which produce seed for reforestation,collected data on pollen shed and cone receptivity phenology from 1997 to 2012. Seed orchards mimiccommon gardens, with several to many provenances represented at each site and duplicated across sites. As6the same genotype is usually represented by multiple trees within a site and is often planted at multiple sites,it was possible to separate out genotype and provenance effects.The seed orchard data also presented certain challenges. Data was collected in many formats andphenological-event coding systems; I regularized such coding to enable meaningful comparison. Data wasalso collected in calendar time (day-of-year); I transformed it to thermal time (growing degree days). Finally,because trees were not observed daily, flowering dates are not known exactly; I explicitly accounted for thiscensoring in my analysis.Specifically, I answer: 1) How much forcing is required for flowering events in lodgepole pine? 2)Does provenance change forcing requirements? And 3) How will climate change affect lodgepole floweringphenology?2.2 MethodsMy aim was to estimate the thermal time required for flowering in lodgepole pine and to test for clinalvariation in those requirements. To determine the timing of pollen shed and cone receptivity, I modelledthe forcing requirements in growing degree days for flowering phenology in lodgepole pine as a functionof provenance mean annual temperature (MAT) while also accounting for heterogeneity due to site, year,genotype and individual tree. Using this thermal time model, I then predicted the timing and length of theflowering period for lodgepole pine at five sites across British Columbia historically and through 2100 undertwo climate change scenarios.2.2.1 Study speciesLodgepole pine is a wind-pollinated monoecious conifer with an extensive range and well documented localadaptation, especially in growth and cold hardiness-related traits (Rehfeldt et al., 2001; Wang et al., 2010;MacLachlan et al., 2017). While lodgepole pine has four subspecies, this work concerns only the mostwidespread subspecies, Pinus contorta ssp. latifolia. Pollen and seed cone buds differentiate in late summerand early fall, then go dormant (Owens, 2006). As temperatures warm the following year, buds resumedevelopment and strobili “flower”; receptive female strobili exude pollination drops between open bractsthat capture pollen shed from mature male strobili. Bracts on the female strobili swell closed after severaldays and pollen shed finishes in a similarly short period of time.2.2.2 Data2.2.2.1 Phenology Data72.2.2.1.1 Overview This project takes advantage of existing lodgepole pine phenology datacollected over a decade and a half by government and industry workers in seed orchards in British Columbia.Trees selected from across the province are grown in multiple, geographically distributed seed orchards aspart of regional tree breeding and seed production programs (Fig. 2.1). To plan for future seed productionand orchard management, seed orchard managers have monitored flowering phenology and seed productionat many of the orchards. Flowering phenology data was collected at the Prince George Tree ImprovementStation in British Columbia beginning in 1997. Collection at many other BC tree orchard sites began in2006 under the Forest Genetics Council’s Operational Tree Improvement Program 0722 (Webber, 2007) andcontinued intermittently through 2012. I contacted seed orchard managers and other forestry profession-als across British Columbia in 2012 and received data for a total of 743 individual trees representing 257genotypes grown at seven seed orchard sites.8Figure 2.1: British Columbia portion of the range of Pinus contorta ssp. latifolia in light green (after Bisbinget al., 2020). Parent tree origins marked with crosses. White circles represent orchard sites. Non-orchardsites Border and Trench shown with orange triangles. Four orchard sites are clustered near Vernon, BC.2.2.2.1.2 Provenance Scions and seeds were collected from trees in 143 natural stands acrossBritish Columbia in the 1970s and 80s and scions were grafted onto rootstock to create a genotype bank(Ukrainetz, 2011). The latitude, longitude, and elevation where each genotype was collected was availablefrom the British Columbia Ministry of Forests Seed Planning and Registry Application (Ministry of Forests(2014)). Breeding and selection programs in different seed orchards targeted better growth for different B.C.9regions (McAuley, 1996). Scions of selected genotypes from the bank were then grafted onto rootstock topopulate the seed orchards, which removes the possibility of confounding from abiotic and maternal effectsthat can occur when using seed-grown plants. Individual copies of these genotypes in the seed orchards arereferred to as trees. Genotype performance in progeny tests and in the orchards led to the removal of somegenotypes from the orchards and addition of others. Only genotypes representing first generation selectionsfrom natural stands were included in the analysis; offspring from crosses of these selections were excluded.Based on MacLachlan et al. (2017), selection for growth had little effect on vegetative bud break phenology,which should be correlated with reproductive phenology, so the seed orchard genotypes should provide agood indication of the phenology of natural populations.There is partial replication of most genotypes across years and sites as well as within sites (Fig. 2.4).While 51.4% of genotypes are represented by a single tree across all sites, 9.3% of genotypes are representedby trees at two different sites and 46.3% of genotypes are represented by 2 to 10 trees at the same site,primarily at the Prince George Tree Improvement Station (PGTIS) and Vernon. Additionally, genotypesand trees are often represented in multiple years. While incomplete, this crossing of factors allows for theestimation of genotype and tree effects in the model.10Figure 2.2: Dates phenological state observed in each orchard at each site for each year of observation.Multiple lines for a given site and year represent the different orchard divisions within each site.112.2.2.1.3 Sites and years Phenology data are from seven seed orchard sites. Model predictionswere extended to two additional sites, Trench and Border, that were further north and colder than the seedorchard locations (Fig. 2.1). Seed orchard sites span about 5∘ of latitude, although four of the orchardsare clustered near Vernon, B.C., and are at elevations from 464 to 632 m with mean annual temperaturesranging from 3.9 at PGTIS to 8∘C at Kalamalka.Each of the seven seed orchard sites is divided into multiple orchards. Genotypes are grouped intoorchards based on geographic origin, growth performance, and genetic diversity targets for seed planningzones, which are large areas that represent biogeoclimatic and administrative regions. PGTIS provided acontinuous 16-year record of its three orchards’ phenology and reliably observed multiple trees per genotype,but data from other sites are sparser.Data were provided by seed orchard workers in the British Columbia Ministry of Forests.2.2.2.1.4 Phenology scoring protocol Protocol C in Woods et al. (1996) was used as thebasis for collecting pollen shed and cone receptivity data as part of operational seed orchard management andgenetic diversity monitoring, though operational constraints led to some modifications. Workers monitoredseed orchards for the beginning of pollen shed and cone receptivity. Approximately 15 genotypes, usuallyrepresented by two trees each, were selected for detailed observations. When flowering was seen to bestarting, workers went into the orchard every few days to make observations on the selected trees. Theperiod of observation and number of days trees were observed were inconsistent and determined by workersat each site (Fig. 2.2).Stages of pollen and seed cone development described by Owens & Molder (1984) and updatedin Owens (2006) were used as a general guide for determining the phenological state of pollen and seedcones. Pollen cones are considered to be flowering (shedding) when tapping causes pollen to be released andseed cones are considered to be flowering (receptive) when there are gaps between the bract-scale complexes.Visible pollination drops are produced by female cones during flowering, though they recede midday if pollenis not present, so their absence is not a reliable indicator of flowering state (Owens, 2006) and was not usedto score phenology.Flowering states were recorded for both pollen shed and cone receptivity at the level of each tree.Protocol C recommends marking the dates when 20% of the cones on a tree have begun flowering and when80% of the cones on a tree have finished flowering. Operationally, there was some subjectivity and tree-levelstates for each cone type should be interpreted as “starting flowering” and “finished flowering.”12Table 2.1: Proportion of observations for each event interval censored or left or right end censored.Event_Label Sex interval left rightfirst_flowering FEMALE 0.63 0.37 NAfirst_flowering MALE 0.35 0.65 NAlast_flowering FEMALE 0.47 NA 0.53last_flowering MALE 0.48 NA 0.52Table 2.2: Flowering state, symbols used to record them in the raw data, and phenophase definitions.State Symbols Male Cones Female Cones0 0 none produced none produced1 1, 2.5, - not yet shedding not yet receptive2 3, 3.5, 4, 4.5, 5, pollenshed20, receptive20 shedding receptive3 -, receptive80, pollenshed80 finished shedding no longer receptiveAll observations are end or interval censored because trees were not observed each day of the surveyperiod (interval censoring) and not all survey periods were long enough to capture the full flowering period(end censoring) (Tab. 2.1). Rates of censoring for the end of receptivity and pollen shed were similar, butthe start of pollen shed was much more likely to be end censored than the start of receptivity.2.2.2.1.5 Data harmonization and cleaning There were some differences in how data wasrecorded at the Prince George Tree Improvement Station versus the other sites. For most of the Prince Georgedata, trees were marked as flowering or not flowering on each day of observation, so some observations wereredundant. At other sites, only the first day observed flowering and the first day observed after flowering wererecorded. At sites using a minimal recording schema, I inferred states for many trees for some observationdates at the orchard assuming that all trees in the same orchard were observed on each observation day. Forexample, if tree A and B are in the same orchard and A is observed to be receptive on Monday and tree Bon Tuesday, then I assumed that tree B was observed as not receptive on Monday. I cleaned and harmonizedthe data for analysis using R scripts so that, for pollen shed and receptivity, state 1 = not yet flowering,state 2 = flowering, and state 3 = finished flowering (Tab. 2.2).Flowering state in the field was recorded using different symbol sets and state resolutions. I sim-plified and harmonized these recording schemes by assigning each symbol to one of the three states above.The few trees that did not produce cones were assigned state 0 and were not included in the analysis.In total, there were 15662 usable observations of phenological state.2.2.2.2 Forcing units13Table 2.3: Table of intercepts 𝛼 and slopes 𝛽 for each site.Site 𝛼 𝛽Border -0.10 1.05Trench 0.16 0.98PGTIS 0.22 1.00KettleRiver 1.28 1.03Sorrento 1.92 1.02Tolko 0.14 0.99PRT 0.34 0.99Vernon 0.05 0.99Kalamalka 0.10 1.002.2.2.2.1 Site weather data Daily weather data at seed orchard sites was extracted fromPNWNAmet, a daily gridded meteorological data set at 1/16∘ (~ 6km) over northwest North America(Werner et al., 2019). The closest point in the PNWNAmet grid was used for each site. Mean dailytemperature was calculated as the average of the minimum and maximum daily temperatures, which addserror to the daily temperature data.Sites were located in areas with large topographic heterogeneity, so small distances between PN-WNAmet grid points and their corresponding sites led to significant elevation differences, which resulted inmisrepresentation of actual site temperatures by PNWNAmet data. More accurate temperature data wasavailable at monthly resolution from ClimateBC, which produces scale-free annual and monthly climate data(Wang et al., 2016), but daily temperatures were required for forcing unit calculations. So, daily grid pointtemperatures from the PNWNAmet data set were adjusted using site-specific monthly weather data. First,mean monthly temperatures for 1945-2012 at all sites were generated using ClimateBC v.7.30. Next, thismonthly site-specific data (𝑐𝑚𝑜𝑛𝑡ℎ,𝑠𝑖𝑡𝑒) was regressed on monthly PNWNAmet data (𝑝𝑚𝑜𝑛𝑡ℎ,𝑠𝑖𝑡𝑒) to generatesite-specific slopes (𝛽 and intercepts (𝛼)) (Eqn. (2.1)). The regression was performed with lm in R’s statspackage.𝑐𝑚𝑜𝑛𝑡ℎ,𝑠𝑖𝑡𝑒 = 𝛼𝑠𝑖𝑡𝑒 + 𝛽𝑠𝑖𝑡𝑒𝑝𝑚𝑜𝑛𝑡ℎ,𝑠𝑖𝑡𝑒 (2.1)Finally, daily mean temperature from PNWNAmet data (𝑝𝑑𝑎𝑡𝑒,𝑠𝑖𝑡𝑒) was corrected using the linearrelationship between monthly ClimateBC and monthly PNWNAmet data (Tab. 2.3, Eqn. (2.2))).𝑇𝑑𝑎𝑡𝑒,𝑠𝑖𝑡𝑒 = 𝛼𝑠𝑖𝑡𝑒 + 𝛽𝑠𝑖𝑡𝑒𝑝𝑑𝑎𝑡𝑒,𝑠𝑖𝑡𝑒 (2.2)14Corrected data, 𝑇𝑑𝑎𝑡𝑒,𝑠𝑖𝑡𝑒 was used for all forcing unit calculations.To verify that this was a reasonable correction, I compared ClimateBC mean monthly temps at thegridpoints to the PNWNAmet data (all slopes between 0.98 and 1.02), intercepts < 1.28, and 𝑅2 > 0.999).2.2.2.2.2 Forcing unit calculation Flowering is a developmental process whose rate is deter-mined by temperature. Forcing units describe the relative effect of temperature on development; warmertemperatures speed development and cooler temperatures slow development. Observable phenological eventsoccur only after a certain amount of forcing has accumulated. The mean daily temperature 𝑇 was mappedto degree day forcing units with a function that sets the threshold for forcing accumulation each day 𝑑 at5∘C.𝑅(𝑇 ) =⎧{⎨{⎩0 if 𝑇 < 5𝑇 − 5 if 𝑇 ≥ 5∘C(2.3)The 5∘C threshold is commonly used for calculating degree days. Experiments measuring thedevelopmental rate of several conifer and broad leaf tree species show that 5∘C results in a good linearapproximation to the true sigmoid response of development to temperature (Sarvas, 1972).Growing degree days 𝑓 are the accumulation of forcing units 𝑅(𝑇 ) over days 𝑑0 through day 𝑑.I chose ordinal day 1 (January 1) as the start time 𝑑0, but mean daily temperatures at our sites do nottypically rise above the 5∘C threshold before late March (Fig. 2.3).𝑓𝑑 =𝑑∑𝑛=𝑑0𝑅(𝑇𝑛) (2.4)Sites with very different mean annual temperatures, like Border and Trench, can have quite similarrates and amounts of forcing accumulation if most of the daily temperature differences are during the coldestmonths when temperatures are below the accumulation threshold of 5∘C.2.2.2.2.3 Typical years Some analyses are for a specific year, but others use the mean dailytemperatures for a “typical” year. The mean temperature for each day of the typical year is the averagetemperature that day across the years 1945 to 2012, the full temporal range of the PNWNAmet data set.15Figure 2.3: Mean daily temperature (A), daily forcing, (B) and forcing accumulation (C) for the first 180days of the year at both seed orchard (solid lines) and comparator (dashed lines) sites averaged across theyears 1945 to 2012.162.2.2.3 Provenance climates Mean annual temperature (MAT) for the 1961-1990 normal period forthe geographical source of each genotype were obtained from ClimateBC v7.30 (Wang et al., 2016) for useas the predictor variable in the thermal time model. Site MATs were also extracted for use when makingpredictions from the thermal time model (Fig. 2.4). MAT was chosen to represent provenance climatesbecause clinal variation typically occurs along temperature gradients (Conover et al., 2009) and MAT ishighly correlated with alternate temperature variables.Figure 2.4: Mean annual temperatures for sites (Panel A) and provenances grown at each site (Panel B). Inpanel A, circles denote seed orchard sites and triangles comparator sites. In panel B, point fill and shapeexpress levels of replication for each genotype grown at a site.2.2.3 ModelI fitted separate multilevel linear Bayesian models of the accumulated forcing 𝑓 for each of the four floweringevents of interest: start of receptivity (female begin), end of receptivity (female end), start of pollen shed(male begin), and end of pollen shed (male end). In the model, accumulated forcing for an event is predictedby the mean annual temperature (MAT) at the source location of a genotype (parent) while variation dueto site, year, genotype, and tree factors is captured with varying intercepts.17𝜙𝑖 is the sum of the intercept 𝜇, the effect of provenance as represented by MAT and quantifiedby the coefficient 𝛽, and offsets 𝛿 for each level of factors site, year, genotype, and tree. Offsets are GDDdifferences from the mean 𝜇, and an offset of zero means there is no effect of a factor level. Data were notsuitable for a nested design because trees were nested within site and genotype, but genotypes were notnested within site. Replication of genotypes within and across sites and years allowed for a more powerfulcrossed design (Fig. 2.4).𝜙𝑖 = 𝜇 + 𝛽MATi + 𝛿𝑆𝑖𝑡𝑒[𝑖] + 𝛿𝑌 𝑒𝑎𝑟[𝑖] + 𝛿𝐺𝑒𝑛𝑜𝑡𝑦𝑝𝑒[𝑖] + 𝛿𝑇𝑟𝑒𝑒[𝑖] (2.5)The likelihood function varies depending on the type of censoring the observation was subject to:left, right, or interval.If a tree is already flowering when it is first observed, the observation for start of flowering is leftcensored; I have an upper bound 𝑈 on the observation day’s accumulated forcing but no lower bound anduse the cumulative normal distribution for the likelihood (Eqn. (2.6)).Pr[𝑓𝑖 < 𝑈] = ∫𝑈∞Normal(𝑦 ∣ 𝜙𝑖, 𝜎)𝑑𝑦 (2.6)If a tree is flowering on the last day that it is observed, the observation is right censored; I havea lower bound 𝐿 on the day’s accumulated forcing but no upper bound and use the complement of thecumulative normal distribution for the likelihood (Eqn. (2.7)).Pr[𝑓𝑖 > 𝐿] = ∫∞𝐿Normal(𝑦 ∣ 𝜙𝑖, 𝜎)𝑑𝑦 (2.7)Interval censored observations are bounded by the flowering observation itself and the last dayobserved, 𝐿, before flowering (for start events) or the first day observed 𝑈 past flowering (for end events);so, the flowering observation is the upper bound 𝑈 for start events and the lower bound 𝐿 for end events.Pr[𝐿 > 𝑓𝑖 > 𝑈] = ∫𝑈𝐿Normal(𝑦 ∣ 𝜙𝑖, 𝜎)𝑑𝑦 (2.8)The residual standard deviation, 𝜎 , represents the variability in the observed data that is notcaptured by the rest of the model.182.2.3.1 Priors Priors were chosen to be regularizing and weakly informative. Priors limit parameterestimates to forcing accumulations possible at our sites during late spring, when lodgepole pine are knownto flower in British Columbia (Owens (2006)). Possible forcing accumulation through May and June wascalculated using the corrected daily PNWNAmet weather data from 1945-2012 across all of our sites. Forcingaccumulation is always positive and the overall distribution across reasonable ranges is strongly right skewed,so I used a gamma distribution with shape (𝛼) and rate (𝛽𝛾) parameters derived from the mean and standarddeviation of the May-June historical forcing accumulation.𝜇 ∼ Gamma(𝛼 = 2.7, 𝛽𝛾 = 0.012) (2.9)The prior on the provenance climate effect, 𝛽, is centered on zero as I do not know if there is aneffect of provenance climate, let alone its direction. A standard deviation of 25 GDD includes the mostextreme possible slope values given potential forcing differences between sites and different provenance meanannual temperatures.𝛽 ∼ Normal(0, 25) (2.10)Observation error is less than the largest accumulation of forcing possible in a two-week period inJanuary - June.𝜎 ∼ |Normal(0, 15)| (2.11)Factors included in the model were site, year, genotype, and tree. Factor offsets 𝛿𝑓𝑎𝑐𝑡𝑜𝑟𝑖 and theirstandard deviations 𝜎𝑓𝑎𝑐𝑡𝑜𝑟 were roughly constrained by both estimates of possible forcing accumulation andthe approximately known length of the flowering period from prior work (Owens et al., 2005). As site andyear have relatively few levels with which to estimate factor standard deviations, half normal priors withtheir relatively thin tails were used to help constrain offset parameter standard deviations.𝛿𝑓𝑎𝑐𝑡𝑜𝑟𝑖 ∼ Normal(0, 𝜎𝑓𝑎𝑐𝑡𝑜𝑟) (2.12)𝜎𝑓𝑎𝑐𝑡𝑜𝑟 ∼ |Normal(0, 9)| (2.13)19Priors on the factor offsets 𝛿𝑓𝑎𝑐𝑡𝑜𝑟𝑖 are adaptive and effects are estimated from the data with partialpooling.2.2.3.2 Model assumptions While trees need chilling temperatures to break endodormancy, I assumethat chilling requirements are always met and that chilling and forcing periods do not overlap. All trees havelikely received enough chilling to minimize forcing requirements. In an experiment where lodgepole pine froman Alberta population were grown in Ontario, chilling sufficient for normal growth initiation was reached inNovember and additional chilling in early spring had minimal impact on forcing requirements (Man et al.,2017). Forcing is a function of temperature only and does not incorporate photoperiod or chilling. Forcingincreases monotonically; thermal time cannot go backwards. Forcing in the model begins accumulating onJanuary 1, after chilling requirements are likely to have been met and before forcing (in reality or at our5∘C threshold) is occurring. Because of how forcing units are calculated (Eqn. (2.3)), this model does notaccount for very slow rates of development at temperatures under 5∘C or damage and abnormal developmentthat can occur at temperature extremes.2.2.3.3 Fit Models were fit in the probabilistic programming language Stan (rstan vers. 2.32.6 (StanDevelopment Team, 2020) via brms vers. 2.22.0 (Bürkner, 2017).) Stan uses the No-U-Turn Sampler, anefficient Markov Chain Monte Carlo method that extends the Hamiltonian Monte Carlo algorithm (Carpenteret al., 2017), to sample the joint posterior. Models were run with 6 independent chains for 4,500 iterations,discarding 2,250 warm-up iterations per chain for a total of 13,500 posterior samples for each parameter.Model convergence and performance were considered good; ?̂? was < 1.01 and bulk and tail effective samplesize was more than 100x the number of chains for all parameter values in all models (Vehtari et al., 2021).There were no divergences, and visual inspection of energy plots and rank plots showed acceptable samplingbehaviour.2.2.4 Validation2.2.4.1 Retrodictions To assess model performance, I tested whether the modelled start, end, andlength of flowering period were consistent with observations. Typically, model predictions are compareddirectly to point observations. However because flowering events are never observed exactly due to censoring,the model’s estimates were compared to observed ranges of forcing units and days instead (Fig. 2.5). Forexample, if a tree was observed to be not flowering on Monday and flowering on Thursday, then any meanretrodiction of Tuesday to Thursday is considered to overlap the observation interval. For left or rightcensored observations, the observation ranges were only bounded on one side. The retrodictions are not20just point estimates, but distributions, so a one standard deviation interval centered on the mean of eachretrodiction was also compared to the observation interval. Model performance was then gauged based onthe proportion of point estimates that fell within the observation range or whose one standard deviationrange around the estimate overlapped the observation range.Figure 2.5: A conceptual diagram of how model retrodictions were compared to observations. In this example,an observation of start of flowering is interval censored; sporadic observations determine the earliest date atree could have started flowering and the first day it was observed flowering. The model’s retrodiction of thisobservation could have a mean that lies within the observation interval (”mean in interval”), the mean couldfall outside of the observation interval but a one standard deviation interval around the mean could overlapthe observation interval (”one sd overlap”), or the observation interval may not overlap at all with the entireone standard deviation interval around the model mean (”no overlap”). Only the estimates coloured greenare counted as correct.Thermal time requirements for the exact combination of sites, years, genotypes, and trees observedwere retrodicted from the full posterior.2.2.4.2 Comparison of independent data with model predictions Model predictions were alsocompared to phenology data from three years of a field trial in central Sweden (Nilsson, 1981) and one yearof a provenance trial in central British Columbia (at the PGTIS site) (O’Reilly & Owens, 1988). Floweringdates and GDD were published in both studies only as provenance averages, not individual records. Sevenprovenances from the Swedish trial and seven from the British Columbia trial had enough data to include.These 14 total provenances were from the Yukon and British Columbia and included one Pinus contorta ssp.contorta provenance. The Swedish trial observed different trees for receptivity and pollen shed, while theBritish Columbia trial tracked male and female strobili on the same trees. The Swedish trial monitored 12-1521trees per provenance and sex, observing 10 strobili of each sex, while the British Columbia trial monitoredeight trees per provenance and five strobili of each sex. Observations were made daily in Sweden and everytwo to four days in British Columbia.Provenance MATs range from -4.5∘C (3.8∘C colder than the coldest provenance used to build themodel) to 5.4∘C. While the Swedish trial reported GDDs, the central British Columbia site reported dates,which were translated to forcing units (see Section 2.2.5.1). Only start data was included for the BritishColumbia trial as it reported flowering waning rather than termination dates. I combined these independentphenology data with estimates of provenance MATs (based on 1961-1990 normals) obtained with ClimateBCv7.30 then generated posterior predictions for each provenance MAT while ignoring random effects. Modelpredictions were plotted against the calculated forcing accumulations to test model performance with anindependent data set.2.2.4.3 Site effects Site offsets were compared to the mean annual temperature at each site to verifythat site effects were not likely to be climate effects.2.2.5 Analysis and prediction2.2.5.0.1 Forcing predictions To generate predictions for forcing requirements relevant to ex-isting seed orchards, I constructed a new data set across the full observed range of provenance MATs withhypothetical new levels for year, tree, genotype, and site. Posterior predictions were made from the fittedmodel while generating new levels from a normal distribution centered at zero with standard deviation sam-pled from the estimates for 𝜎𝑓𝑎𝑐𝑡𝑜𝑟𝑠, i.e. using the estimated priors for the random effects. Forcing predictionswere made across the entire range of provenance MATs observed. The full posterior was downsampled to6000 draws for these predictions.2.2.5.1 Translating forcing unit estimates into day of year A forcing unit estimate for a floweringevent can be translated into a specific day of year for any location and year and vice versa. A daily timeseries of accumulated forcing units was constructed for each desired year and location using Equation (2.4).GDD estimates for flowering events were then matched to accumulated forcing units in daily time series toidentify the corresponding day of year. Translations were made for each seed orchard site plus the morenortherly Border and Trench sites for a yearly time series between 1945 and 2012 and for 30-year historicaland future climate normal periods.222.2.5.2 Site synchronization Day of year estimates from 1945-2012 had patterns of early and lateyears that appeared to echo across sites. I simulated expected day of year of flowering events for provenancesgrown at their home sites with no random effects for all sites (including Border and Trench). I quantifiedpatterns of flowering across sites using rank correlation. Specifically, at each site I ranked the day of year onwhich flowering occurred across years, from the earliest to the latest, focusing on the coldest and warmestprovenances observed within the orchards. Then I calculated the Kendall rank correlation coefficient to assesspairwise correlations between each site for each flowering event. The choice of Kendall’s 𝜏 was motivatedby its superior ability to handle tied rankings relative to other methods like Spearman’s 𝜌, as tied rankingscould occur as often as 50% of the time in some comparisons.To understand the spatial extent over which these flowering patterns occur, I regressed rank corre-lations on distance between sites for each flowering event and provenance. Models were implemented withthe lm function in R. Great circle distances between points were calculated with the geosphere package inR.2.2.5.3 Countergradient variation impacts To understand the effects of the estimated countergradi-ent variation (𝛽), I compared the expected day of year for flowering events with and without the provenanceeffect. I translated the thermal time model’s accumulated forcing expectations to day of year for treessourced from all sites and grown either at 1) their home site or 2) at PGTIS in a typical year. To calculatethe expected day of year for flowering without the MAT effect, I translated the intercept 𝜇 estimated by themodel to day of year in a typical year at the site. The expectations exclude random effects and observationerror.2.2.5.4 Climate change I predicted flowering events at five locations along a temperature and latitudinalgradient in British Columbia for 30-year climate normal periods from 1951-2100. Kalamalka, Kettle River,and PGTIS are seed orchard sites; Trench and Border are additional sites added to make predictions acrossa wider range of climates, better reflecting both the provenances of trees grown at the seed orchards andcovering significantly more of the lodgepole pine range. Predictions are for these locations, not for thespecific seed orchard trees grown at some of these sites. Monthly climate normal data for historical andfuture periods was extracted from ClimateBC v7.30 (Wang et al., 2016). Future periods are based on an8-ensemble set of global climate model projections suggested by Mahony et al. (2022) when using a toollike ClimateBC (i.e. change-factor methods) to downscale climate data to high spatial and low temporalresolution. Shared Socioeconomic Pathways (SSPs) 2 and 5 tier one scenarios with radiative forcing levels in2100 of 4.5 and 8.5 W/m2, respectively, were chosen to show a range of possible shifts in phenology (O’Neill23et al., 2016). To obtain daily data for a future normal period, daily data for the typical year at each sitewas adjusted so that the mean monthly temperature of the typical year matched the ClimateBC projectionfor each climate normal period. The typical years were created by averaging the mean daily temperature ofeach day from 1945-2012 (the full set of daily data from the PNWNAmet data set).2.2.6 Data and code availabilityFlowering and climate data as well as model and analysis scripts are available on github and archived onZenodo (Tysor, 2025a,b,c). Model output is also archived on Zenodo with a guide for making new predictionsfrom the model (Tysor, 2025d).2.3 ResultsA thermal time model based on 16 years of flowering phenology data at seven seed orchard sites across BritishColumbia estimates the forcing accumulation required for lodgepole pine pollen shed and cone receptivity,which allows for flowering predictions at any location when daily temperature data and provenance MAT areavailable. Forcing requirements show a pattern of countergradient variation, with populations from warmerregions requiring more forcing to flower than those from colder regions. Phenology simulations under climatechange from the thermal time model show lodgepole pine flowering at least two weeks to more than a monthearlier by 2100 than they were in the mid-1900s.2.3.1 Validation2.3.1.1 Retrodiction Comparison between event observations, which are censored, and the median forc-ing estimate for those events shows that the thermal time model is responsive to censoring information (Fig.2.6). The model correctly estimates lower forcing requirements than measured for left censored observa-tions and higher forcing requirements for right censored observations. The model generally estimates higherforcing requirements for end events than observed, even for interval censored observations.24Figure 2.6: Comparison of observations to model predictions. The accumulated forcing at time of observationfor flowering events versus median estimates from posterior predictive distributions of the thermal time modelroughly follows a 1:1 line. The model attempts to correct for censoring in the data, so predictions for treesfirst observed after the start of flowering are lower than the accumulated forcing at the time of observation,and predictions for trees last observed before the end of flowering are higher.Because events were censored, model estimates for observations were compared to the observedaccumulated forcing intervals instead of the exact forcing on the day of observation (Conceptual diagramin Fig. 2.5). The mean thermal time estimate of flowering events was within the observed event rangefor forcing at least 60% of the time, and one standard deviation intervals around the mean overlapped theobserved range at least 96% of the time (Table 2.4).Length of flowering period was correctly retrodicted more often than the start or end of floweringevents because some poorly estimated events still have the correct length. The retrodicted mean floweringperiod length fell within the observed interval at least 90% of the time and a one standard deviation intervalaround the the retrodicted length overlapped the observed length interval at least 99% of the time.25Table 2.4: Comparison of model estimates to observations. Table entries are the proportion of mean modelretrodictions contained in the observed range for each flowering event for forcing accumulation and day ofyear predictions, so 1 is perfect prediction and 0 is no overlap between the observed range and modeledflowering GDD or day of year. In parentheses are the proportion of one 𝜎 intervals around mean modelretrodictions that overlap the observed range.Begin EndForcing Day of Year Forcing Day of YearReceptivity 0.6(0.97) 0.6(0.96) 0.71(0.97) 0.69(0.95)Pollen shed 0.76(0.99) 0.76(0.99) 0.7(0.96) 0.68(0.95)Receptivity (female) retrodictions were less likely than pollen shed to fall within the observed eventrange. As accumulated forcing was calculated at daily resolution, the proportion of day of year estimatesthat match observations is nearly identical to that for forcing.2.3.1.2 Independent data When confronted with independent phenology data, the model performswell (Fig. 2.7). Predictions for the P. contorta ssp. contorta provenance are not worse than for those forssp. latifolia. As random effects were excluded, model predictions are the same across years for a givenprovenance, but observations vary year-to-year.26Figure 2.7: Comparison between independent data and model predictions. Data is from provenance trialsin central British Columbia (O’Reilly & Owens, 1988) and Sweden (Nilsson, 1981) using provenances fromacross western Canada. One provenance was from ssp. contorta instead of ssp. latifolia. Measurements areyearly provenance averages; no uncertainty information was available. Model predictions show 50% and 95%highest density intervals (HDIs) around the estimate. Model predictions use 6000 draws from the posteriorand do not include random effects. A 1:1 line is shown for reference.2.3.2 Average forcing requirementsThe base heat sum requirement for the start of flowering for a population sourced from a climate with aMAT of 0 ∘C while ignoring effects of site, year, genotype, and tree is about 179 GDD for pollen shed and192 for receptivity (Tab. 2.5). Pollen shed occurs over a wider range of forcing than receptivity, beginning ata lower forcing accumulation and ending at a higher forcing accumulation than receptivity; the completionof pollen shed requires around 68 GDD above the requirements for pollen shed onset while receptivity endsafter approximately an additional 51 GDD. Estimates of forcing accumulations required to begin floweringare more certain than for ending flowering.27Table 2.5: Mean forcing requirements and 95% HDI for beginning and end of pollen shed and cone receptivityfor trees sourced from a provenance with a MAT of 0°C and ignoring the effects of site, year, genotype, andtree. (Estimates correspond to thermal time model intercept 𝜇.) Based on 6000 draws from the posterior.Receptivity Pollen shedBegin 192 (178-207) 179 (164-197)End 243 (224-260) 247 (226-267)Table 2.6: Accumulated forcing requirements increase for warmer provenances. Slope (𝛽) estimates and 95%highest density posterior intervals are from thermal time models for seed cone receptivity and pollen shedphenology in lodgepole pine. Forcing calculated in growing degree days with threshold 5°C.Receptivity Pollen shedBegin 5.21 (3.3-7.19) 5.85 (4.13-7.57)End 5.22 (3.06-7.28) 5.91 (3.96-8.02)2.3.3 Provenance MAT effectAs the MAT of the provenance grows warmer, accumulated forcing requirements also increase for all floweringevents (Tab. 2.6). The provenance MAT effect is higher for pollen shed than receptivity and higher for theend of flowering than the start. While the overall effect is small at between about 5 and 6 ∘C GDD of forcingrequirements per degree change in MAT, this translates to a difference from several days to more than aweek in expected flowering time given typical daily forcing at considered sites (Fig. 2.3 B).2.3.4 Factor impacts and model errorI estimated random effects for site, year, genotype, and tree (Fig. 2.8). Standard deviations among factorlevels are fairly similar for year, site, and genotype (ranging from about 12 to 21 GDD) and translate intoseveral days of difference during the flowering period (Fig. 2.4). Standard deviation parameter results arereflected in the median offset estimates (Fig. 2.8 B); greater spread in offset estimates corresponds to higherstandard deviation estimates for a factor. Site and year variation may reflect both year-to-year variation inclimate influences not captured by temperature, inaccurate mean daily temperature estimates, and differencesbetween sampling crews. There is no relationship between site MATs and their offset estimates (𝛿𝑠𝑖𝑡𝑒𝑖); siteoffsets are capturing random variation among sites, perhaps due to unmeasured weather effects or differencesin observers. Uncertainty in standard deviation for site (𝜎𝑠𝑖𝑡𝑒) and year (𝜎𝑦𝑒𝑎𝑟) is higher than for other factorsbecause there are relatively few sites (7) and years (16) (Fig. 2.8 A).Standard deviation of individual site effects were driven by the range of provenances grown at thatsite. The range of provenances grown at a site explains between 47% and 55% of the variation in standarddeviation in the site effect. The number of years a site was observed, however, was unrelated to the standard28deviation in site effects despite PGTIS having 16 years of observation compared to a maximum of 4 years atany other site. The larger the provenance MAT range at a site, the lower the uncertainty around a site effect.This suggests that at sites that sampled a smaller range of provenances, the model may have struggled topartition variation between the site and genotype effects.Genotype effects can be as large as 42.2 GDD (95% HDI 29.5-56.6) (for the start of receptivity ingenotype 1770), while tree-level effects are low at < 13 GDD, suggesting strong genetic control over flowering.Among-genotype variation is slightly higher for variation in receptivity than pollen shed, while the overalllow within-genotype variation (𝜎𝑇𝑟𝑒𝑒) is higher for variation in pollen shed than receptivity. This suggeststhat plasticity may be slightly higher for pollen shed than receptivity.Uncertainty in factors associated with the end of flowering is higher than for the beginning. Thisis likely because the end of flowering occurs at higher daily temperatures than the start of flowering, whichreduces the precision of GDD estimates; e.g. the forcing requirement can be estimated more precisely ifforcing accumulates at 10 GDD per day versus at 20 GDD.The residual standard deviation, 𝜎, is on the same scale as year, site and genotype standarddeviations. Error in the likelihood model could have several sources. Forcing estimates could be inaccurateeither because of error in daily temperature estimates or due to inaccuracies in the forcing model itself (e.g.the choice of threshold model or a different or even varying threshold). Other variables not included, suchas photoperiod or chilling, could affect phenology predictions as well.29Figure 2.8: A) Standard deviation parameter estimates from the thermal time model. Sigma (𝜎) estimatesuncertainty for the observation model, while other panels show standard deviations of the offsets for year,site, genotype, and tree factors (𝜎𝑓𝑎𝑐𝑡𝑜𝑟). B) Medians of offset estimates (𝛿𝑓𝑎𝑐𝑡𝑜𝑟).2.3.5 General growing degree day predictionsFigure 2.9 presents the thermal time model’s predicted forcing requirements for lodgepole pine flowering,which demonstrate a broad range of possible forcing accumulations, from approximately 105 to 370 GDD.Pollen shed typically spans a slightly broader GDD range than receptivity. Notably, on average, colderprovenances tend to require less forcing to flower compared to their warmer counterparts. Despite this, thevariation in forcing requirements is considerable, allowing for the potential overlap in flowering times amongpopulations from the coldest (MAT -0.7 ∘C) and warmest (MAT 6.8 ∘C) provenances when grown together.Usefully, these growing degree day estimates can be converted to specific days of the year for any locationwhere daily average temperatures are known, providing a practical tool for predicting flowering events.30Figure 2.9: Growing degree day requirements for pollen shed and cone receptivity in lodgepole pine acrossat different mean annual temperatures (MAT). Predictions are model estimates with offsets and uncertaintyfor an average site, year, genotype, and tree. Predictions are done at 100 MATs spanning the full range ofprovenances in the phenology data set and are based on 6000 draws from the posterior. Lines show medianestimates and ribbons show 95% HDI.2.3.6 Day of year predictionsGrowing Degree Day (GDD) predictions (Fig. 2.9) were converted into day of year estimates, illustratingsignificant interannual variability at all sites (Figs. A2.2 and A2.3). Flowering initiation at the warmestseed orchard site, Kalamalka, typically occurs in May, whereas at the coldest seed orchard site, PGTIS, itextends into June (Fig. 2.10). Within each site, the flowering period can fluctuate by over a month acrossdifferent years. While there is no discernible trend in flowering times across the 1945 to 2012 time period,notably earlier flowering periods were recorded during the early-1990s and the first decade of the 2000s.Warming trends at sites are small relative to year-to-year temperature variation and largest in the coldermonths of the year when forcing is not accumulating (Fig. A2.1), which is likely why flowering phenologydoes not appear to have advanced in recent decades. The majority of observed flowering events fell withinthe 95% highest density interval (HDI) of the model predictions (Figs. A2.2 and A2.3). Genotypes fromwarmer regions exhibited later flowering at a given site compared to those from cooler origins, highlightingthe influence of genetic differences on phenological responses to climatic conditions.31Figure 2.10: Predictions of flowering day of year 1945-2012 for trees from the coldest and warmest provenanceMAT growing at a cold and a warm seed orchard. Black horizontal lines show observed event ranges for alltrees grown at that site. Predictions use site-specific model estimates with offsets and uncertainty for anaverage year, genotype, and tree based on 6000 draws from full posterior. Lines show median estimates forall provenances and ribbons show 95% HDI.32Average differences between the coldest and warmest provenance at each site range from 4.2 to 6.2days (Fig. 2.11). The pattern of average differences within sites reflects the size of the provenance MATeffect (Tab. 2.6): the largest average differences are for male events and the smallest are for female events.Differences among the provenance differences are 0.8 to 1.4 days, but are swamped by the standard deviation.These differences between provenances are similar across sites and there is no relationship between them andsite MAT.Figure 2.11: Provenance flowering differences can be large. Mean predicted difference between flowering dayof year of the coldest and warmest provenance at each seed orchard site in a typical year. Error bars showstandard deviation.There is marked synchronicity in flowering patterns across the orchards: years that witnessed early(or late) flowering at one site generally corresponded to early (or late) flowering at others (Figs. A2.3and A2.2). Crucially, the similarity in flowering patterns was more pronounced among sites that weregeographically closer together (Fig. 2.12). A regression analyses, with results detailed in Table A2.2, showsstrong correlations between rank correlation and distance between sites. Distance explains a significantproportion of the variance in rank correlations (with R2 values between 0.93 and 0.95) and there is aconsistent negative relationship: rank correlations diminish by between 5.7% (± 0.23) to 5.9% (± 0.25) forevery 100 km increase in distance. All parameters in these models achieved statistical significance at a levelof p < 10-5.33Figure 2.12: Geographic distance is a significant predictor of similarity in flowering patterns across seedorchard sites. The y-axis gives the rank correlation (Kendall’s 𝜏) between pairs of sites, with distancebetween pairs on the x-axis. Solid lines represent the linear regression fits (and shaded areas denote the 95%confidence intervals) illustrating the statistically significant negative relationship between distance and rankcorrelations between sites.2.3.6.1 Effects of countergradient variation in seed orchards Our analysis revealed a counter-gradient variation in flowering phenology (Tab. 2.6), where populations across the landscape flower moresynchronously than expected in the absence of a provenance effect. In individual seed orchards, where manyprovenances are grown together, countergradient variation causes divergence in flowering times and an over-all longer flowering period for the orchard than a single provenance stand. When nine provenances weresimulated growing together at a single site (PGTIS), the countergradient effect of provenance MAT causesdifferent provenances to flower at different times (Fig. 2.13). If there were no effect of provenance, then allprovenances would flower at the same time in that orchard.34Figure 2.13: Day of year predictions for nine provenances (seed orchard sites plus Border and Trench) grownat PGTIS. Grey points show estimates without provenance effect (𝜇 from the thermal time model). Greenpoints show estimates with provenance effect. Point estimates (means) only; no uncertainty. Provenanceswere based on the seed orchard sites and two comparison sites and simulations ignored factor offsets.2.3.6.2 Climate change Simulations of historical and future flowering at nine sites spanning 9.4 ∘C ofMAT predicted significant advancements in expectations for both pollen shed and receptivity as temperaturesincrease throughout the century (Fig. 2.14). Advancements are reported as medians (95% quantile interval).In simulations for the 1951-1980 climate normal period, the earliest flowering event was the start of pollenshed at the warmest site on May 20 (May 18 - 22) and the latest was the end of pollen shed at the coldTrench site on June 26 (June 23 - 28). While 1951-1980 to 1981-2010 normal periods saw flowering advanceby 1 to 6 days on average, simulations showed even earlier flowering for future climate normal periods. UnderSSP2-4.5, expected flowering advanced by an average of 16 (16-17) days (at Border) to 23 (22-23) days (atKettleRiver). Under SSP5-8.5, a more pessimistic emissions scenario, expected flowering advanced by 24(24-25) days (at Border) to 32 (31-32) days (at KettleRiver).The length of the flowering period, overlap between male and female flowering events, and floweringat different sites do not appear to change much with climate warming this century.35Figure 2.14: Flowering events at five sites predicted by the thermal time model for 30-year climate-normalperiods 1951-2100 with two Shared Socioeconomic Pathway scenarios for 2012-2100. Flowering advances byseveral weeks to a month by the end of the century. 6000 draws, median model expectation with no factoroffsets, 50 and 95% HDI with 1961-1991 MAT normals.2.4 DiscussionI fit a thermal time model based on 16 years of flowering phenology data at seven seed orchard sites acrossBritish Columbia to estimate the forcing accumulation required for lodgepole pine pollen shed and cone re-ceptivity, which allows for flowering predictions at any location when daily temperature data and provenanceMAT are available. Forcing requirements show a pattern of countergradient variation, with populations fromwarmer regions requiring more forcing to flower than those from colder regions. Under a moderate climatechange scenario, lodgepole pine are predicted to flower at least two weeks earlier by the end of the centurywithout having large effects on flowering overlap across the landscape.2.4.1 Comparison of dates and GDD with other studiesWhile there is limited research on lodgepole pine flowering phenology, O’Reilly & Owens (1988) observedthe entire flowering period for seven provenances at three BC seed orchard sites in one year and Nilsson36(1981) did the same for 14 provenances grown at two sites over four years in Sweden. Results from a subsetof those provenances were reported in such a way that they could be compared to model results. This datais a challenging test of the model because it includes fourteen provenances not included in my data set,one of which was a different subspecies and five of which were from colder provenances than the model wasbuilt with. Additionally, the trial site in central Sweden was colder than any of the seed orchard sites inthis study. The trees were also younger and just entering maturity. Despite these differences, the modelpredicted flowering phenology well (Fig. 2.7). The end of flowering, especially for pollen shed, is generallyharder to determine and more variable than the start (Nilsson, 1981), which is reflected in the model’s higheruncertainty for parameters associated with this phase (Tabs. 2.5 & 2.6, Fig. 2.8) and poorer match betweenobserved and predicted flowering in the independent data.To my knowledge, the only other direct report of heat sum requirements for lodgepole pine floweringis from Owens (2006). That estimate of 500 GDD is much higher than those reported in Nilsson (1981) orassociated with flowering dates in other studies (O’Reilly & Owens, 1988; Owens et al., 2005). Unfortunately,the source of the phenology and climate data underlying the 500 GDD claim in Owens (2006) is unclearand could not be investigated. It is likely that a different weather data source led to a discrepancy inthe calculation of forcing units, a frequent source of difficulty in estimating and comparing GDD models(Chamberlain & Wolkovich, 2023).2.4.2 Relationship between pollen shed and receptivityThe models presented here predict that pollen shed starts slightly before and ends slightly after receptivity,which is congruent with previous findings at the stand level (Nilsson, 1981; Owens et al., 2005). Pollen shedand receptivity phenology are tightly correlated within individual trees (Yazdani & Fries, 1993), though theorder of flowering is not consistent year-to-year and flowering onset and termination within a tree can varyup to five days (Nilsson, 1981), with sunnier parts of the tree generally flowering earliest (Sarvas, 1962).Concerns expressed by Owens et al. (2005) about the potential for extreme protandry in seed orchards inwarm years are not supported by this work.2.4.3 Individual variationGenotypes have consistent effects across years; e.g. genotypes that flower early one year flower early the next.Consistent rankings among genotypes in leaf unfolding have been observed in other forest tree species suchas Quercus robur, Quercus petraea and Fagus sylvatica (Delpierre et al., 2017) and early and late floweringlodgepole pine genotypes are recognized in seed orchards (Yazdani & Fries, 1993; Owens et al., 2005). While37most genotype effects are small, they can be large relative to the provenance MAT effect. The largestestimated genotype effect was 42.2 GDD (95% HDI 29.5-56.6). This effect is equivalent to the expectedforcing accumulation gap between provenances with MATs that differ by 8.1 ∘C – slightly larger than theentire range of provenances examined in this study. Within-population variation has a larger effect onflowering phenology than countergradient variation. This is in keeping with previous work which found thatgenetic variation associated with growth initiation phenology is higher among individuals than populations(MacLachlan et al., 2021). Individual variation in flowering phenology can have direct and indirect effects onmean fitness, and large versus small amounts of variation can change the population response to non-linearselection (Zettlemoyer & DeMarche, 2022).Although flowering timing varies greatly from year to year, this plasticity reflects robust geneticcontrol of development in response to forcing temperatures, with variation predominantly occurring amongindividuals and, to a lesser extent, among populations. Differences in flowering forcing requirements forindividual trees of the same genotype year-to-year were typically much smaller and less variable than genotypeeffects, especially for female receptivity. This is in keeping with high heritability for spring phenologicaltraits identified in other conifers (Howe et al., 2003) and low population differentiation for spring vegetativephenological traits in lodgepole pine (Liepe et al., 2016). Larger within-genotype variation in pollen shedthan receptivity combined with larger among-genotype variation for receptivity than pollen shed suggestsgreater plasticity for male flowering than female. This may reflect the fact that pollen shed can be delayedby rainy weather; in Pinus sylvestris, pollen sacs cannot open if relative humidity is higher than 50% (Sarvas,1962); the thermal time model presented here does not account for this effect.2.4.4 Countergradient variationCountergradient variation requires reciprocal transplant or provenance experiments to quantify (Conover etal., 2009; Wadgymar et al., 2022). While seed orchards are not complete reciprocal transplant experimentsthat grow all provenances in all environments, by growing multiple provenances in orchards spanning a widerange of MATs, they resemble incomplete provenance trials. Combining this design with a model accountingfor the structure of the data, we were able to identify and quantify countergradient variation in flowering: forevery degree increase in provenance MAT, an additional 5.21 (3.3-7.19) to 5.91 (3.96-8.02) GDD are requiredfor flowering. This counteracts the primary environmental effect which leads to earlier flowering in the souththan the north (Owens et al., 1981).Two other studies found countergradient variation in lodgepole pine flowering (Nilsson, 1981;O’Reilly & Owens, 1988). While O’Reilly & Owens (1988) noted that receptivity varied markedly more38than pollen shed, Nilsson (1981) and this work show similar levels of countergradient variation in both. Ifound pollen shed showed a slightly stronger countergradient effect than receptivity, but there were notstrong differences. O’Reilly & Owens (1988) measured pollen catch rather than directly observing strobili,which may account for the difference. Interestingly, the closely related species Pinus sylvestris also showscountergradient variation in pollen shed (Luomajoki, 1993; Whittet et al., 2017), but receptivity was notexamined in that work.Unlike most other types of trees, pines do not have separate vegetative or reproductive bud breakevents, but instead have both types of bud together in compound long shoot buds that start elongating inspring. Due to this feature, spring vegetative and reproductive phenology are expected to be highly correlatedin lodgepole, and indeed growth initiation phenology in lodgepole has shown countergradient variation inboth growth chamber and field common garden experiments (Liepe et al., 2016; MacLachlan et al., 2017).In an experiment with populations sourced from more limited geographic areas along elevational rather thanlatitudinal gradients, however, countergradient variation in growth initiation was not observed (Green, 2007),and Liepe et al. (2016) found that overall earlier growth initiation in colder locations was lessened in wetter,higher elevation populations, perhaps due to the influence of heavy snowpack.Propagation and selective breeding for growth within seed orchards may have influenced our findingof countergradient variation in flowering phenology. Parent trees were selected from natural stands based onsize and stem form. While selection for superior growth traits in lodgepole pine did not lead to significantdifferences in growth initiation between natural and selected juvenile trees (MacLachlan et al., 2017), it ispossible that this selection could have amplified the flowering cline observed in this study.2.4.4.1 Implications and evolution of countergradient variation In lodgepole pine, countergradi-ent variation in flowering phenology increases overlap among populations in different climates. While thiseffect is weak, it combines with high within-population variation to increase synchronization among popu-lations. For instance, in Pinus sylvestris, countergradient variation resulted in populations 137 km apartwith a two-week difference in peak pollen shed still maintaining overlap in some years (Whittet et al., 2017).Despite both countergradient variation and high within-population variation in flowering acting to increaseoverlap, pollen mediated gene flow is strongly limited by differences in flowering timing across the landscape,especially at larger spatial scales.Countergradient variation may provide colder populations an advantage in a warming climate byfacilitating synchrony with warmer populations, which is particularly important during range expansionsdriven by rising temperatures. The lodgepole pine distribution is still expanding to its full climatic niche post-39glaciation (Johnstone & Chapin, 2003). At the extreme leading edge of the distribution, earlier receptivitymay enhance seed production by synchronizing with more central populations, mitigating local pollen scarcitythat can occur at the expanding edge (Elleouet & Aitken, 2019). Pollen scarcity is exacerbated not onlyby low population density but also by the age structure of trees; since northern provenances of lodgepolepine produce seed cones at earlier ages than they do pollen (Nilsson, 1981), fecundity in young stands isdependent on pollen from older stands. Conversely, at the trailing edge where local conditions are unsuitablefor seedling establishment, later pollen production could enable pollen-mediated gene flow from the warm anddeclining edge to the cooler, more favourable areas in the range center. Although existing local adaptationmight become maladaptive across lodgepole’s range under climate change (Davis & Shaw, 2001), increasedflowering overlap between cooler and warmer populations, along with substantial individual variation withinpopulations, may provide a buffer. Whether gene flow aids or hinders adaptation will depend on whetherpollen dispersal follows or opposes temperature gradients.While co-gradient variation can evolve without selection, countergradient variation requires it. Sim-ulations show that when flowering is driven by temperature, assortative mating along thermal gradients leadsto co-gradient variation, even in the absence of selection (Soularue & Kremer, 2014). This type of variationenhances population differences and can be intensified by divergent selection that aligns with the environmen-tal gradient. In contrast, countergradient variation mitigates differences along thermal gradients, reducingexpected phenological divergence. Countergradient variation in flowering indicates that plasticity alone maynot optimize phenotypes, with selection acting on flowering timing (or a linked trait) in opposition to theenvironmental effect.Selection may act directly on flowering timing or on linked traits such as growth initiation. Inlodgepole pine, female gametophyte development begins after pollination, with fertilization delayed until thefollowing year. However, the gametophyte must reach the free nuclear stage before dormancy (Owens, 2006).Short growing seasons in northern regions may constrain flowering phenology, as northern range limits areoften set by the ability to mature seeds within these short periods (Morin et al., 2007). Pollen and seedcone development in lodgepole pine is linked to growth initiation, and selection on growth may indirectlyinfluence flowering timing. In a model of budburst timing in trees with countergradient variation in thattrait, selection for earlier budburst was much stronger at higher elevations, where the optimal window forbudburst was narrower (Gauzere et al., 2020). Similar selective pressures on growth initiation could occurin lodgepole pine, which also exhibits weak countergradient variation in growth timing (Liepe et al., 2016).While flowering is correlated with growth initiation, it is unclear whether selection acts directlyon flowering, growth initiation, or another linked trait. Growth initiation is linked not just to flowering,40but also to other traits important to fitness like height growth and cold hardiness. In southern populations,selection for taller trees delays growth initiation and cessation, while in northern populations, it favoursearlier timing (MacLachlan et al., 2017). This may reflect historical pressures favouring competitive growthin the south versus survival and reproduction in harsher northern climates; earlier growth initiation innorthern populations is also linked to improved cold injury protection (MacLachlan et al., 2021). Due tothe intricate connections between growth initiation, flowering, and other traits, it is unclear what drivescountergradient variation in flowering.However, it is likely that strong selection for reproduction in northern climates drives the phenologi-cal differences observed in lodgepole pine. Flowering and reproduction are resource-intensive processes; tallertrees generally flower more than shorter ones, even in shaded conditions, and strobili, particularly female, arepredominantly produced in the sunniest parts of the crown (Sarvas, 1962; Nilsson, 1981). In northern envi-ronments, the harsh climate often limits the resources and time available for flowering. Southern provenancesthat manage to survive in these regions frequently fail to flower, and all provenances exhibit reduced flower-ing in the north, even when trees reach comparable heights (Nilsson, 1981). A study of 14 western Canadianlodgepole pine provenances grown at two sites in Sweden (Nilsson, 1981) revealed substantial variation inreproductive traits between northern and southern provenances, suggesting strong selection for reproductivesuccess in northern populations. Not only did both age at maturity and flowering abundance show counter-gradient variation but southern provenances were more likely to produce pollen cones alone, whereas northernprovenances tended to produce seed cones alone. (While mature trees typically produce both male and fe-male strobili, maturing trees often produce only one type.) The severe environmental constraints on bothgrowth and reproduction in northern climates likely drive countergradient variation in flowering phenology,reflecting adaptive trade-offs between growth and reproductive success in these challenging environments.2.4.4.2 Implications in seed orchards While countergradient variation enhances flowering synchro-nization in natural populations at their original locations, it has the opposite effect in seed orchards, leadingto increased divergence among populations (Fig. 2.13). Seed orchards are designed to ensure that selectedtrees contribute to genetic gains in growth and maintain genetic diversity above a target threshold. However,when flowering times vary among populations, it disrupts random mating, reduces the genetic diversity ofseed lots, and may even cause pollen limitation in some trees. Genotypic variation can either mitigate orexacerbate these issues. Because provenances have relatively small differences in phenology, differences inmating contributions within the orchards are more likely to be driven by individual variation in phenologyand provenance and individual variation in flowering abundance (Nilsson, 1981). To address synchronizationissues and disproportionate flowering abundance, seed orchards often rely on labour-intensive techniques like41supplemental mass pollination. Divergence in flowering times also extends the overall flowering period inseed orchards compared to natural stands, increasing the likelihood of pollen contamination from outsidethe orchard.2.4.5 Climate change impactsUnder current policies we are on track for warming of 3.1°C (range: 1.9–3.8) by 2100 (Environment Pro-gramme, 2024), congruent with expectations for the moderate climate change scenario SSP 2-4.5. Based onthe thermal time model presented here, flowering is expected to advance by 2-3 weeks under the moderateemissions scenario and 3-4 weeks under a more extreme emissions scenario. This change in flowering timecould lead to more frequent frost damage to strobili and reduce fecundity. Interestingly, even under extremewarming, the model does not predict large changes in the duration of flowering within a site or relative timingbetween sites, suggesting that gene flow may not be affected much by even large advances in flowering.At a broad geographic scale, climate is correlated across the landscape within a year. Flowering innearby populations may not be synchronized, especially in topographically variable environments, but thecloser populations are to one another, the more similar their patterns of early or late flowering will be. Evenat distances of 600 km, early and late years are still correlated half the time (Fig. 2.12). This correlation andany synchronization that depends on it could break down if warming occurs unevenly across the landscape.The role of pollen deposition in forest nutrient cycling is not well studied, but the research that doesexist suggests that changes to the timing of the nutrient pulse from pollen could have ecosystem impacts. Inmixed P. densiflora and P. rigida stands in South Korea, the nutrients provided by pollen were a significantfraction of those released by litterfall and occurred during a period of active growth in the forest (Lee & Booth,2003). Pollen that falls in peat bogs provides essential nutrients in a very nutrient-limited environment withpine pollen making an outsized contribution (Shumilovskikh et al., 2015). Mycorrhizae and plant growth inforests and nearby bogs could be affected by future changes in pollen deposition depending on the importanceof the timing of the nutrient pulse it produces.In some species, warming can reduce chilling accumulation and result in increased forcing require-ment to trigger developments, which could delay spring phenological events (Duputié et al., 2015). I assumehere that lodgepole flowering phenology in interior British Columbia has minimal variation due to differ-ences in chilling accumulation as all trees have likely received enough chilling accumulations to minimizeforcing requirements. In lodgepole pine from one Alberta population grown in Ontario, chilling sufficientfor normal growth initiation was reached in November and continued to accumulate into April; howeverafter February, additional chilling had only a very small impact on forcing requirements (Man et al., 2017).42A decline in chilling accumulation under climate change could lead to increased forcing requirements andslower advancement than projected here. On the other hand, chilling could accumulate faster under climatechange if warming occurs in colder months (as is occurring at the sites included here (Fig. A2.1)), whichcould advance phenology (Guy, 2014). The predictions for lodgepole pine presented here assume a linearresponse to warming, but interactions between chilling and forcing are complex (Flynn & Wolkovich, 2018)and can extend beyond the current growing season (Fu et al., 2014).2.4.6 Additional strengths and weaknessesThis thermal time model offers a straightforward and effective method for predicting flowering events inlodgepole pine. Its simplicity makes it easy to apply, yet it overlooks potential contributors to floweringphenology like correlations between phenological stages (Ettinger et al., 2018), the role of chilling (Flynn &Wolkovich, 2018), and photoperiod sensitivity (Zohner et al., 2016). The use of a constant forcing thresholdof 5 °C is common for forest tree species and is supported by empirical evidence (Chuine et al., 2003), but thisthreshold may vary by population (Chuine et al., 2001) and likely underestimates forcing at temperaturesbelow 5 °C while overestimating forcing at high temperatures (Sarvas, 1972). It also fails to account for theimpact of extreme events, such as frost or heat damage.A key contribution of this work is its ability to predict flowering for any lodgepole pine populationwith a known provenance. However, predictions for populations at the southern and far northern extremesof the species’ range should be approached with caution. The provenances used in this model cover only8 of the 27 degrees of latitude within the lodgepole pine (ssp. latifolia) range and cover MATs of just -0.7to 6.8 ∘C, and the seed orchard sites represent an even smaller fraction. Expanding the model to includemore provenances and sites would increase confidence in its predictions across the full geographic range ofthe species.433 ConclusionThis study uses a novel thermal time model to estimate forcing requirements for both receptivity andpollen shed in Pinus contorta ssp. latifolia (lodgepole pine), quantifies countergradient variation and within-population variation in flowering, and predicts significant advancements in flowering under climate change.Estimates are based on a model fit using my assembly of the most comprehensive data set of lodgepole pineflowering available to date, spanning 149 provenances grown at 7 sites and observed over 16 years. Unlikemany phenology models that rely on day of year, which limits their applicability, the model developed herecan be used across different locations and years.This model offers a powerful tool for predicting flowering timing across space and time, deepeninginsights into gene flow, adaptation, and climate change impacts on lodgepole pine. Its adaptability will allowflowering predictions across large regions to be assessed for overlap, with potential integration into windmodels to evaluate pollen-mediated population connectivity. Given flowering’s critical role in both fecundityand gene flow, these predictions will be valuable for investigating eco-evolutionary dynamics, particularly theinfluence of phenology on lodgepole pine range limits (Chuine & Beaubien, 2001). Additionally, phenologymodel outputs could be incorporated into larger process-based models to explore fitness consequences ofplasticity and adaptation under past and future climates (Duputié et al., 2015; Gauzere et al., 2020).The large within-population variation in flowering and modest countergradient variation associatedwith provenance mean annual temperature (MAT) that I found have important implications for gene flowand adaptation. Temperature-controlled flowering does restrict gene flow to populations from more similarclimates, but both the within-population and countergradient variation act to increase gene flow amongpopulations from more different climates. Phenology, however, is not the only factor controlling gene flow;studies that consider wind dispersal of pollen are necessary to understand the impacts of pollen mediatedgene flow in lodgepole pine and how it will affect adaptation under climate change (Kremer et al., 2012). Thedirection of pollen dispersal influences evolution along the gradient (Soularue & Kremer, 2014), but mod-elling pollen dispersal over large distances is a difficult problem. It could be simplified by using pollinationphenology to rule out gene flow between populations with unmatched flowering phenology.Advances in flowering phenology of several weeks are anticipated for lodgepole pine under climatechange. While the limited analyses here suggest advancement will not have a large impact on gene flow,further investigation is warranted using more sites or even considering the full range. Other climate changeimpacts could be incorporated as well. Importantly, gene flow could be impacted by extensive logging, rangeshifts likely under climate change (MacKenzie & Mahony, 2021), or other disruptions due to climate change44like increased fire and bark beetle outbreaks (Sambaraju et al., 2012; Parisien et al., 2023). Gene flow maysimply be lost from many populations because they are gone, temporarily or permanently.Data collected in seed orchards is invaluable for research, as demonstrated in this study. However,its use is limited by inconsistent collection methods, varied recording formats, and its fragmentation acrossmany sites. The flowering phenology dataset used here, sourced from both government and private entities,required substantial effort to obtain, standardize, and clean before analysis. Consolidating, standardizing,and publishing similar data from seed orchards would maximize the return on investment, making it moreaccessible for future research and broader applications.Although the model currently uses MAT to represent provenance, MAT may not capture the truedrivers of population differences — a hypothesis that could be tested by comparing phenology at sites withcontrasting MAT but similar spring warmups. MAT’s performance as a proxy variable could be tested byground truthing phenology at site pairs like Border and Trench, which have very different MATs due towinter differences but very similar spring warmups. When provenances with very different MATs have verysimilar spring warmups, their forcing requirements will be different according to the model here. This meansdivergence in flowering time and reduced gene flow potential between populations that are relatively similarin climate. Climate variables other than MAT may better represent provenance variation and should beinvestigated.The thermal time model presented here offers a straightforward yet powerful approach for predictinglodgepole pine flowering and provides a solid foundation for extension across the species range. To expandits use to more southerly regions of the range or to enhance climate change projections, the model couldbe tested and calibrated with data from additional southern provenances. This model was built usingobservational data, but experiments may provide better forcing estimates, clarify whether and how to includephotoperiod and chilling effects, and quantify forcing thresholds (Chuine et al., 2001; Hänninen et al., 2019).Additionally, incorporating both start and end events into a multi-stage model could yield more precisepredictions, especially for end events.This study’s flowering phenology model and data set for lodgepole pine offer a foundation forfuture research into climate-driven phenological shifts, gene flow, and adaptation. Building on this workcould reveal new insights into how phenology influences ecological and evolutionary dynamics under climatechange across lodgepole pine’s range.45BibliographyAitken SN, Yeaman S, Holliday JA, Wang T, Curtis-McLane S. 2008. Adaptation, migration orextirpation: Climate change outcomes for tree populations. Evolutionary Applications 1: 95–111.Alberto FJ, Aitken SN, Alía R, González-Martínez SC, Hänninen H, Kremer A, Lefèvre F,Lenormand T, Yeaman S, Whetten R, et al. 2013. Potential for evolutionary responses to climatechange - evidence from tree populations. Global Change Biology 19: 1645–1661.Bisbing SM, Urza AK, Buma BJ, Cooper DJ, Matocq M, Angert AL. 2021. Can long-lived specieskeep pace with climate change? Evidence of local persistence potential in a widespread conifer. Diversityand Distributions 27: 296–312.Bürkner P-C. 2017. Brms: An R Package for Bayesian Generalized Linear Mixed Models using Stan.Journal of Statistical Software 80.Campbell ID, McDonald K, Flannigan MD, Kringayark J. 1999. Long-distance transport of polleninto the Arctic. Nature 399: 29–30.Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, GuoJ, Li P, Riddell A. 2017. Stan: A Probabilistic Programming Language. Journal of Statistical Software76: 1–32.Chamberlain CJ, Wolkovich EM. 2023. Variation across space, species and methods in models of springphenology. Climate Change Ecology 5: 100071.Chefaoui RM, Serebryakova A, Engelen AH, Viard F, Serrão EA. 2019. Integrating reproductivephenology in ecological niche models changed the predicted future ranges of a marine invader. Diversity andDistributions 25: 688–700.Chuine I. 2000. A unified model for budburst of trees. Journal of Theoretical Biology 207: 337–347.46Chuine I, Aitken SN, Ying CC. 2001. Temperature thresholds of shoot elongation in provenances ofPinus Contorta. Canadian Journal of Forest Research 31: 1444–1455.Chuine I, Beaubien EG. 2001. Phenology is a major determinant of tree species range. Ecology Letters4: 500–510.Chuine I, Kramer K, Hanninen H. 2003. Plant development models. In: Schwartz M, ed. Phenology:An Integrative Environmental Science. Springer, 217–235.Chuine I, Régnière J. 2017. Process-Based Models of Phenology for Plants and Animals. Annual Reviewof Ecology, Evolution, and Systematics 48: 159–182.Conover DO, Duffy TA, Hice LA. 2009. The Covariance between Genetic and Environmental Influencesacross Ecological Gradients. Annals of the New York Academy of Sciences 1168: 100–129.Conover DO, Schultz ET. 1995. Phenotypic similarity and the evolutionary significance of countergra-dient variation. Trends in Ecology & Evolution 10: 248–252.Davis MB, Shaw RG. 2001. Range Shifts and Adaptive Responses to Quaternary Climate Change.Science 292: 673–680.Delpierre N, Guillemot J, Dufrêne E, Cecchini S, Nicolas M. 2017. Tree phenological ranksrepeat from year to year and correlate with growth in temperate deciduous forests. Agricultural and ForestMeteorology 234–235: 1–10.Duputié A, Rutschmann A, Ronce O, Chuine I. 2015. Phenological plasticity will not help all speciesadapt to climate change. Global Change Biology 21: 3062–3073.Elleouet JS, Aitken SN. 2019. Long-distance pollen dispersal during recent colonization favors a rapidbut partial recovery of genetic diversity in Picea sitchensis. New Phytologist 222: 1088–1100.Ennos RA. 1994. Estimating the relative rates of pollen and seed migration among plant populations.Heredity 72: 250–259.47Environment Programme UN. 2024. Emissions Gap Report 2024: No more hot air ... Please! With amassive gap between rhetoric and reality, countries draft new climate commitments. Nairobi.Ettinger AK, Gee S, Wolkovich EM. 2018. Phenological sequences: How early-season events definethose that follow. American Journal of Botany 105: 1771–1780.Flynn DFB, Wolkovich EM. 2018. Temperature and photoperiod drive spring phenology across allspecies in a temperate forest community. New Phytologist 219: 1353–1362.Franks SJ. 2015. The unique and multifaceted importance of the timing of flowering. American Journalof Botany 102: 1401–1402.Franks SJ, Weis AE. 2009. Climate change alters reproductive isolation and potential gene flow in anannual plant. Evolutionary Applications 2: 481–488.Fu YSH, Campioli M, Vitasse Y, De Boeck HJ, Van den Berge J, AbdElgawad H, Asard H,Piao S, Deckmyn G, Janssens IA. 2014. Variation in leaf flushing date influences autumnal senescenceand next year’s flushing date in two temperate tree species. Proceedings of the National Academy of Sciences111: 7355–7360.Gauzere J, Teuf B, Davi H, Chevin L-M, Caignard T, Leys B, Delzon S, Ronce O, Chuine I.2020. Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptivevalue of plasticity. A case study on tree phenology. Evolution Letters 4: 109–123.Green DS. 2007. Controls of growth phenology vary in seedlings of three, co-occurring ecologically distinctnorthern conifers. Tree Physiology 27: 1197–1205.Guy RD. 2014. The early bud gets too warm. New Phytologist 202: 7–9.Hamann A, Wang T. 2006. Potential effects of climate change on ecosystem and tree species distributionin British Columbia. Ecology 87: 2773–2786.48Hänninen H, Kramer K, Tanino K, Zhang R, Wu J, Fu YH. 2019. Experiments Are Necessary inProcess-Based Tree Phenology Modelling. Trends in Plant Science 24: 199–209.Howe GT, Aitken SN, Neale DB, Jermstad KD, Wheeler NC, Chen TH. 2003. From genotypeto phenotype: Unraveling the complexities of cold adaptation in forest trees. Canadian Journal of Botany81: 1247–1266.Iler AM, CaraDonna PJ, Forrest JRK, Post E. 2021. Demographic Consequences of PhenologicalShifts in Response to Climate Change. Annual Review of Ecology, Evolution, and Systematics 52: 221–245.Inouye DW. 2008. Effects of Climate Change on Phenology, Frost Damage, and Floral Abundance ofMontane. Ecology 89: 353–362.Johnstone JF, Chapin FS. 2003. Non-equilibrium succession dynamics indicate continued northernmigration of lodgepole pine. Global Change Biology 9: 1401–1409.Kitamoto N, Ueno S, Takenaka A, Tsumura Y, Washitani I, Ohsawa R. 2006. Effect of floweringphenology on pollen flow distance and the consequences for spatial genetic structure within a population ofPrimula sieboldii (Primulaceae). American Journal of Botany 93: 226–233.Kremer A, Potts BM, Delzon S. 2014. Genetic divergence in forest trees: Understanding the conse-quences of climate change. Functional Ecology 28: 22–36.Kremer A, Ronce O, Robledo-Arnuncio JJ, Guillaume F, Bohrer G, Nathan R, Bridle JR,Gomulkiewicz R, Klein EK, Ritland K, et al. 2012. Long-distance gene flow and adaptation of foresttrees to rapid climate change: Long-distance gene flow and adaptation. Ecology Letters 15: 378–392.Lee EJ, Booth T. 2003. Macronutrient input from pollen in two regenerating pine stands in southeastKorea. Ecological Research 18: 423–430.Leites L, Benito Garzón M. 2023. Forest tree species adaptation to climate across biomes: Buildingon the legacy of ecological genetics to anticipate responses to climate change. Global Change Biology 29:4711–4730.49Liepe KJ, Hamann A, Smets P, Fitzpatrick CR, Aitken SN. 2016. Adaptation of lodgepole pineand interior spruce to climate: Implications for reforestation in a warming world. Evolutionary Applications9: 409–419.Linkosalo T, Häkkinen R, Terhivuo J, Tuomenvirta H, Hari P. 2009. The time series of floweringand leaf bud burst of boreal trees (1846–2005) support the direct temperature observations of climaticwarming. Agricultural and Forest Meteorology 149: 453–461.Lotan JE. 1976. Cone serotiny-fire relationships in lodgepole pine. In: Proc. Tall Timbers Fire EcologyConference. Tallahassee, FL: Tall, 267–278.Luomajoki A. 1993. Climatic adaptation of Scots pine (Pinus sylvestris L.) in Finland based on maleflowering phenology. Acta Forestalia Fennica 0.MacKenzie WH, Mahony CR. 2021. An ecological approach to climate change-informed tree speciesselection for reforestation. Forest Ecology and Management 481: 118705.MacLachlan IR, McDonald TK, Lind BM, Rieseberg LH, Yeaman S, Aitken SN. 2021. Genome-wide shifts in climate-related variation underpin responses to selective breeding in a widespread conifer.Proceedings of the National Academy of Sciences 118.MacLachlan IR, Wang T, Hamann A, Smets P, Aitken SN. 2017. Selective breeding of lodgepolepine increases growth and maintains climatic adaptation. Forest Ecology and Management 391: 404–416.Mahony CR, Wang T, Hamann A, Cannon AJ. 2022. A global climate model ensemble for downscaledmonthly climate normals over North America. International Journal of Climatology 42: 5871–5891.Man R, Lu P, Dang Q-L. 2017. Insufficient Chilling Effects Vary among Boreal Tree Species and ChillingDuration. Frontiers in Plant Science 8: 1354.McAuley LH. 1996. Seed Orchards of British Columbia. Ministry of Forests, Silviculture Practices Branch:Province of British Columbia.50Ministry of Forests. 2014. Seed Planning and Registry Application (SPAR) - Province of BritishColumbia.Morgenstern EK. 1996. Genetic Variation in Some Major Species. In: Geographic Variation in ForestTrees: Genetic Basis and Application of Knowledge in Silviculture. UBC Press, 90–133.Morin X, Augspurger C, Chuine I. 2007. Process-based modeling of species’ distributions: What limitstemperate tree species’ range boundaries? Ecology 88: 2280–2291.Muir PS, Lotan JE. 1985. Disturbance History and Serotiny of Pinus-Contorta in Western Montana.Ecology 66: 1658–1668.Nathan R, Horvitz N, He Y, Kuparinen A, Schurr FM, Katul GG. 2011. Spread of North Americanwind-dispersed trees in future environments. Ecology Letters 14: 211–219.Nilsson J-E. 1981. Flowering in Pinus contorta: A study of flowering abundance and flowering phenologyof provenances from latitude 480N to 630N in a Swedish field trial.O’Neill BC, Tebaldi C, van Vuuren DP, Eyring V, Friedlingstein P, Hurtt G, Knutti R, KrieglerE, Lamarque J-F, Lowe J, et al. 2016. The Scenario Model Intercomparison Project (ScenarioMIP) forCMIP6. Geoscientific Model Development 9: 3461–3482.O’Reilly C, Owens JN. 1988. Reproductive growth and development in seven provenances of lodgepolepine. Canadian Journal of Forest Research 18: 43–53.Owens JN. 2006. The reproductive biology of lodgepole pine. Victoria, B.C.: Forest Genetics Council.Owens JN, Bennett J, L’Hirondelle S. 2005. Pollination and cone morphology affect cone and seedproduction in lodgepole pine seed orchards. Canadian Journal of Forest Research 35: 383–400.Owens JN, Molder M. 1984. The reproductive cycle of lodgepole pine. Victoria, B.C.: InformationServices Branch, Ministry of Forests.51Owens JN, Simpson SJ, Molder M. 1981. Sexual reproduction of Pinus contorta . I. Pollen development,the pollination mechanism, and early ovule development. Canadian Journal of Botany 59: 1828–1843.Parisien M-A, Barber QE, Bourbonnais ML, Daniels LD, Flannigan MD, Gray RW, HoffmanKM, Jain P, Stephens SL, Taylor SW, et al. 2023. Abrupt, climate-induced increase in wildfires inBritish Columbia since the mid-2000s. Communications Earth & Environment 4: 1–11.Parmesan C. 2006. Ecological and Evolutionary Responses to Recent Climate Change. Annual Review ofEcology, Evolution, and Systematics 37: 637–669.Pau S, Wolkovich EM, Cook BI, Davies TJ, Kraft NJB, Bolmgren K, Betancourt JL, ClelandEE. 2011. Predicting phenology by integrating ecology, evolution and climate science. Global ChangeBiology 17: 3633–3643.Piao S, Liu Q, Chen A, Janssens IA, Fu Y, Dai J, Liu L, Lian X, Shen M, Zhu X. 2019.Plant phenology and global climate change: Current progresses and challenges. Global Change Biology 25:1922–1940.Rehfeldt GE. 1983. Adaptation of Pinus contorta populations to heterogeneous environments in northernIdaho. Canadian Journal of Forest Research 13: 405–411.Rehfeldt GE, Wykoff WR, Ying CC. 2001. Physiologic plasticity, evolution, and impacts of a changingclimate on Pinus contorta. Climatic Change 50: 355–376.Rehfeldt GE, Ying CC, Spittlehouse DL, Hamilton DA. 1999. Genetic responses to climate in Pinuscontorta niche breadth, climate change, and reforestation. Ecological Monographs 69: 33.Rivest S, Lajoie G, Watts DA, Vellend M. 2021. Earlier spring reduces potential for gene flow viareduced flowering synchrony across an elevational gradient. American Journal of Botany 108: 538–545.Sambaraju KR, Carroll AL, Zhu J, Stahl K, Moore RD, Aukema BH. 2012. Climate changecould alter the distribution of mountain pine beetle outbreaks in western Canada. Ecography 35: 211–223.52Sarvas R. 1962. Investigations on the flowering and seed crop of Pinus silvestris. MetsäntutkimuslaitoksenJulkaisuja 53: 1–198.Sarvas R. 1972. Investigations on the annual cycle of development of forest trees. Active period. Metsän-tutkimuslaitoksen Julkaisuja 76: 1–110.Savolainen O, Pyhäjärvi T, Knürr T. 2007. Gene Flow and Local Adaptation in Trees. Annual Reviewof Ecology, Evolution, and Systematics 38: 595–619.Schneider RR, Latham MC, Stelfox B, Farr D, Boutin S. 2010. Effects of a Severe MountainPine Beetle Epidemic in Western Alberta, Canada under Two Forest Management Scenarios. InternationalJournal of Forestry Research 2010: 1–7.Shumilovskikh LS, Schlütz F, Achterberg I, Kvitkina A, Bauerochse A, Leuschner HH. 2015.Pollen as nutrient source in Holocene ombrotrophic bogs. Review of Palaeobotany and Palynology 221:171–178.Slavov GT, Howe GT, Adams WT. 2005. Pollen contamination and mating patterns in a Douglas-firseed orchard as measured by simple sequence repeat markers. Canadian Journal of Forest Research 35:1592–1603.Soularue JP, Kremer A. 2012. Assortative mating and gene flow generate clinal phenological variationin trees. BMC Evolutionary Biology 12: 79.Soularue JP, Kremer A. 2014. Evolutionary responses of tree phenology to the combined effects ofassortative mating, gene flow and divergent selection. Heredity 113: 485–494.Stan Development Team T. 2020. RStan: The R interface to Stan.Trudgill DL, Honek A, Li D, Straalen NMV. 2005. Thermal time – concepts and utility. Annals ofApplied Biology 146: 1–14.53Tysor CS. 2025a. Scisus/picolaDataFlowering: Interior lodgepole pine (Pinus contorta spp. Contorta) flow-ering phenology in British Columbia seed orchards 1997-2012 (1.0.0). Zenodo. DOI:10.5281/zenodo.14597484.Tysor CS. 2025b. Scisus/picolaDataClimate: Temperature data for seed orchard sites and lodgepole pineprovenances in British Columbia(1.0.0). Zenodo. DOI:10.5281/zenodo.14597516.Tysor CS. 2025d. Lodgepole pine flowering phenology brms model objects (1.0.0). Zenodo.DOI:10.5281/zenodo.14597582.Tysor CS. 2025c. Scisus/picolaModel: Thermal time model of interior lodgepole pine (Pinus contorta ssp.Latifolia) pollen shed and cone receptivity (1.0.0). Zenodo. DOI:10.5281/zenodo.14597656.Ukrainetz N. 2011. Program History: Lodgepole Pine Program Review.Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P-C. 2021. Rank-Normalization, Folding,and Localization: An Improved R^ for Assessing Convergence of MCMC (with Discussion). BayesianAnalysis 16: 667–718.Vitasse Y, Delzon S, Bresson CC, Michalet R, Kremer A. 2009. Altitudinal differentiation in growthand phenology among populations of temperate-zone tree species growing in a common garden. CanadianJournal of Forest Research 39: 1259–1269.Wadgymar SM, Cumming MN, Weis AE. 2015. The success of assisted colonization and assisted geneflow depends on phenology. Global Change Biology 21: 3786–3799.Wadgymar SM, DeMarche ML, Josephs EB, Sheth SN, Anderson JT. 2022. Local Adaptation:Causal Agents of Selection and Adaptive Trait Divergence. Annual Review of Ecology, Evolution, andSystematics 53: 87–111.Wadgymar SM, Weis AE. 2017. Phenological mismatch and the effectiveness of assisted gene flow.Conservation Biology 31: 547–558.54Wang T, Hamann A, Spittlehouse D, Carroll C. 2016. Locally Downscaled and Spatially Customiz-able Climate Data for Historical and Future Periods for North America (I Álvarez, Ed.). PLOS ONE 11:e0156720.Wang T, O’Neill GA, Aitken SN. 2010. Integrating environmental and genetic effects to predictresponses of tree populations to climate. Ecological Applications 20: 153–163.Webber J. 2007. Collection of Crop Statistics for Interior Lodgepole Pine Orchards. ProSeed Consulting.Werner AT, Schnorbus MA, Shrestha RR, Cannon AJ, Zwiers, Dayon G, Anslow F. 2019. Along-term, temporally consistent, gridded daily meteorological dataset for northwestern North America.Whittet R, Cavers S, Cottrell J, Rosique-Esplugas C, Ennos R. 2017. Substantial variation in thetiming of pollen production reduces reproductive synchrony between distant populations of Pinus sylvestrisL. In Scotland. Ecology and Evolution 7: 5754–5765.Woods J, Stoehr M, Webber J. 1996. Protocols for rating seed orchard seedlots in British Columbia.Research Report 06. Victoria, BC: Province of British Columbia, Ministry of Forests Research Program.Yazdani R, Fries A. 1993. Flower abundance, phenology, and pollination pattern in a Pinus contortaseed orchard. In: Pinus contorta: From untamed forest to domesticated crop. Umeå, Sweden: Institutionenför skoglig genetik och växtfysiologi. Swedish University of Agricultural Sciences. Department of ForestGenetics and Plant Physiology., 375–387.Yeaman S, Jarvis A. 2006. Regional heterogeneity and gene flow maintain variance in a quantitativetrait within populations of lodgepole pine. Proceedings of the Royal Society B: Biological Sciences 273:1587–1593.Zettlemoyer MA, DeMarche ML. 2022. Dissecting impacts of phenological shifts for performance acrossbiological scales. Trends in Ecology & Evolution 37: 145–157.Zhu K, Woodall CW, Clark JS. 2012. Failure to migrate: Lack of tree range expansion in response toclimate change. Global Change Biology 18: 1042–1052.55Zohner CM, Benito BM, Svenning J-C, Renner SS. 2016. Day length unlikely to constrain climate-driven shifts in leaf-out times of northern woody plants. Nature Climate Change 6: 1120–1123.56AppendicesAppendix A - Supplementary information for Chapter TwoFiguresFigure A2.1: Monthly mean temperatures for all sites for overlapping 30 year climate normal periods 1921-2020. Monthly data from ClimateBC v7.3057Figure A2.2: Day of year receptivity predictions 1945-2012 for trees from the coldest and warmest provenanceMAT growing at all sites. Black vertical lines show observed event ranges for all trees grown at that site.Predictions use site-specific model estimates with offsets and uncertainty for an average year, genotype, andtree based on 6000 draws from full posterior. Lines show median estimates for all provenances and ribbonsshow 95% HDI.58Figure A2.3: Day of year pollen shed predictions 1945-2012 for trees from the coldest and warmest provenanceMAT growing at all sites. Black vertical lines show observed event ranges for all trees grown at that site.Predictions use site-specific model estimates with offsets and uncertainty for an average year, genotype, andtree based on 6000 draws from full posterior. Lines show median estimates for all provenances and ribbonsshow 95% HDI.59TablesTable A2.1: Site metadata.Site Latitude Longitude Elevation MAT TypeBorder 59.53200 -127.7910 671 -1.4 comparisonTrench 56.15500 -124.1080 659 1.5 comparisonPGTIS 53.76839 -122.7158 598 3.9 seed orchardKettleRiver 49.22266 -118.9340 632 6.9 seed orchardSorrento 50.86859 -119.4517 506 7.0 seed orchardTolko 50.39683 -119.2099 522 7.3 seed orchardPRT 50.38709 -119.2918 464 7.6 seed orchardVernon 50.22995 -119.3201 525 7.8 seed orchardKalamalka 50.23570 -119.2766 486 8.0 seed orchardTable A2.2: Geographic distance is a significant predictor of similarity in flowering patterns across seedorchard sites. After stratifying by sex and flowering event, flowering day of year at each site was ranked fromearliest to latest between 1945 and 2012 and pairwise correlations between sites calculated as Kendall’s 𝜏 .R2 values and slopes for a linear model of Kendall’s 𝜏 as a function of distance between sites illustrating thestatistically significant negative relationship between distance and rank correlations between sites. Slope inunits of Kendall’s tau per km.Sex event r.squared estimate std.errorFEMALE begin 0.93 -0.00058 2.6e-05FEMALE end 0.95 -0.00058 2.3e-05MALE begin 0.94 -0.00059 2.5e-05MALE end 0.95 -0.00057 2.3e-0560
Thesis/Dissertation
eng
Forestry
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Aitken, Sally N.
Countergradient variation in lodgepole pine flowering phenology
Text