INVESTIGATION INTO CLIMATIC EFFECTS ON THE GROWTH AND GENETIC S TRUCTURE OF SKY ISLAND PONDEROSA PINE By Paula E. Marquardt A DISSERTATION Submitted to Michigan State University i n partial fulfillment of the requirements for the degree of Plant Biology Doctor of Philosophy 2018 ABSTRACT INVESTIGATION INTO CLIMATIC EFFECTS ON THE GROWTH AND GENETIC STRUCTURE OF SKY ISLAND PONDEROSA PINE By Paula E. Marquardt In the desert southwest, significant variations in moisture and temperature occur along steep altitudinal gradients. Ponderosa pine forests in the Santa Catalina Mountains , southeast Arizona, USA, consist of two partially sympatric species , the variable needle Pinus ponderosa var. brachyptera (3 - and mixed - needle) that pr efers cool and moist conditions , and the 5 - needle P. arizonica that prefers hot and dry conditions . The objective of this research program is to determine how disturbance that is variations in climate contributes to shifts in population structure . Da ta w as collected along two south - facing slopes (of similar elevation and aspect) for average needle number per fascicle, tree - ring widths, and gene frequencies, to examine limiting factors to growth, climatic shifts in grow th , and population genetic structure. Our data show that c o - occurring ponderosa pine needle - types (3 - , mixed - , and 5 - needle) have filled different ecological niches ; their growth is limited by seasonal water availability, which also controls the length of the growing season . Three distinct genetic groups are present in the Santa Catalina Mountains. Genetic variability is reduced for the 5 - needle type suggesting a possible bottleneck or founding e vent. Results from this study about moisture limitations will help land managers to determine range limits for seed planting zones. iii This dissertation is dedicated to my husband, Jerry Rayala, for all your support. iv ACKNOWLEDGMENTS I thank my major professor, Frank Telewski, and other committee members Kim Scribner, Jim Smith and Nate Swenson for mentoring and guidance on the project , and for their helpful comments o n the final version s of the manuscript s . The USDA Forest Service, and t he Paul Taylor Fellowship, Department of Plant Biology and Beal Botanical Garden, both of Michigan State University, supported the project. Departmenta l staff from Plant Biology (Michigan State University) and Northern Research Station (United States Forest Service ) aided with purchasing and logistics. The Coronado National Forest provided access to samples and the soils data. I thank the following people for help on the project: J. Rayala, J. Grimes, G. Friedlander, J. Kilgore, and B. Epperson assisted with field data collection. A. Foss , R. Williging, J. Lund, M. Meisenheimer, H. Jenson, H. Stricker, and S. Lietz , provided lab , data processing, o r GIS support. J. Stanovick provided statistical support. S. Chinn supplied expertise in the use and interpretation of dendrochronology software. J. and L. Griffith, and A. and T. Harlan (both deceased) provided logistic support for the field research. R. Cronn contributed the DNA sequences for the microsatellite marker development project. D. McKenney provided the spatial climate data sets for the Santa Catalina Mountains. B. Sturtevant provided co nstructive feedback on the project . D. Donner, J. Smith, and K. Scribner reviewed papers. L. Rayala edited paper s . Reprinted by permi ssion from Springer , Trees - Structure and Function ; Variable Climate Response Differentiates Growth of the Sky Island Ponderosa Pines ; Marquardt , P. E., Miranda, B. R., Jennings, S., Thurston , G., and Telewski F. W.; © This is a U.S. v government work and its text is not subject to copyright protection in the United States; however, its text may be subject to foreign copyright protection 2018 ; a dvance online publication, 7, Oct., 2018, doi.org/10.1007/s00468 - 018 - 1778 - 9 , TREES . The final publication is available at link.springer.com. vi TABLE OF CONTENTS LIST OF TABLES ix LIST OF FIGURES xi CHAPTER 1 MONITORING FOREST RE SPONSES TO CHANGING CLIMATE 1 AB S TRACT 1 INTRODUCTION 2 RATIONALE AND SIGNIF ICANCE 6 METHODS 9 CONCLUSIONS AND FUTURE WORK 1 0 APPENDICES 1 4 APPENDIX A: Manuscripts 1 5 APPEN DIX B: Tables and Figures 1 7 REFERENCES 2 0 CHAPTER 2 VARIABLE CLIMATE RES PONSE DIFFERENTIATES THE GROWTH OF SKY ISLAND PONDEROSA PINES 2 5 ABSTRACT 2 7 INTRODUCTION 2 8 METHODS 3 1 Study area and Sampling 3 1 Climate data 3 2 Trend and C limate - growth analysis 3 3 RESULTS 3 7 A nnual and Monthly climate 3 7 Seasonal climate 3 7 Series length, Chronology statistics, and Growth patterns 3 8 Correlations 3 9 Simple correlation and Response function 3 9 PCP correlation and TAVG partial correlation 40 PDSI correlation and TAVG partial correlation 4 1 DISCUSSION 4 1 S eries length, Chronology statistics, and Growth patterns 4 1 Response Function Analysis 4 4 Correlations to Climate 4 6 Sensitivity to PDSI 4 7 vii Respiration 4 8 C ONCLUSIONS 4 9 APPENDICES 5 2 APPENDIX C : Validation of Climatic Data 5 3 APPENDIX D: Supplementary F igures 5 6 APPEN DIX E: Tables and Figures 69 REFERENCES 8 2 CHAPTER 3 SHIFTS IN CLIMATE - GROWTH RELATIONSHIPS OF THE SKY ISLAND PINES 9 1 ABSTRACT 9 2 INTRODUCTION 9 3 METHODS 95 Study area and Sampling 95 Climate data 96 Tree growth chronologies 96 Climate - growth relationships 97 RESULTS 99 Seasonal climate - growth relationships 99 Temporal stability of climate - growth relationships 99 DISCUSSION 10 0 Seasonal climate - growth relationships 10 0 Temporal stability of climate - growth relationships 10 3 C ONCLUSIONS 1 04 APPENDICES 1 06 APPENDIX F : Supplementary Tables and Figures 1 07 APPEN DIX G: Tables and Figures 1 12 REFERENCES 1 1 9 CHAPTER 4 GENETIC DIFFERENTIATION AMON G SKY ISLAND POPULATIONS OF PONDEROSA PINE 1 2 4 ABSTRACT 1 2 5 INTRODUCTION 1 2 6 METHODS 1 2 8 Study s ite 1 2 8 Sampling and Marker analysis 1 2 9 Diversity and Population structure 1 3 2 RESULTS 1 3 3 viii I dentification of groups 1 3 3 Genetic variation among haplotypes and populations 1 3 4 Population structure 1 3 5 DISCUSSIO N 1 3 6 Population structure and Differentiation 1 3 6 C ONCLUSIONS 1 3 8 APPENDICES 1 4 1 APPENDIX H : Supplementary Tables and Figures 1 4 2 APPE N DIX I : Tables and Figures 1 52 APP ENDIX J: Taxonomy of ponderosa p ine 1 6 4 REFERENCES 1 6 8 ix LIST OF TABLES Table 1 .1 Geologic time scale modified from White et al. (2007; Table 1) 1 8 Table 2. 1 Sampling design used for two south - facing slopes, Santa Catalina R ange (sites MTL and BIG) 7 0 Table 2 .2 Summaries for climate data obtained from gridded data sets (McKenney et al. 2011) for two study sites: MTL and BIG 7 1 Table 2. 3 General statistics for the common interval (1925 - 2009) of the de - trended chronologies for P. arizonica and P. ponderosa var. brachyptera from transition zones at two sites (MTL and BIG) 7 2 Table 2. 4 Correlation coefficients of seasonal PD SI with tree - ring c hronologies 7 3 Table 2. 5 Effect sizes ( r ) of PDSI as a measure of drought sensitivity ( 7 5 Table 3. S1 Ponderosa pine tree - ring data obtain ed from the ITRDB, except where noted. Three needle types were analyzed [5 - , 3, and mixed (M)] at two sites: BIG and MTL. N = the number of trees sampled. (pre - 1950) than Marquardt et al . ( xxxxx chronologies will be deposited with ITRDB within 1 yr.) 1 08 Table 3.1 Average needle counts for individual trees sampled for tree - r i ng correlation analyses at two sites: MTL and BIG 1 13 Table 3.2 Percent difference (% DIFF) in growth between one dry winter (1961) and one normal winter (1962) for 3 - , mixed - , and 5 - needle types 1 14 Table 4.S1 Glossary of genetic terms 1 43 x Table 4.S2 F luorescently labeled unique sequences used in 3 - primer PCR 1 45 Table 4.S3 Summary information for the 42 chloroplast haplotypes analyzed 1 46 Table 4.S4 Hapl otype frequencies (42 haplotypes; 1 locus; 427 samples; N, number o f individuals per 21 population s ) 1 48 Table 4.1 Summary information for the 21 pine populations analyzed in this study 1 53 Table 4.2 Sites, needle type, populatio ns, sample sizes, and haplotype polymorphism characterized for 21 populations of ponderosa pine at 5 sites 1 55 Table 4.3 Analysis of molecular variance of plastid SSRs for 9 populations of ponderosa pine at contact zones in southwest Arizona, USA 1 56 Table 4.4 Pairwise PhiPT values between sympatric and distant populations of ponderosa pine 1 57 Table 4.5 Pairwise PhiPT values between populat ions of ponderosa pine from contact zones in the Santa Catalina Mountains 1 58 Table 4.A1 Authorities for three taxa of the genus Pinus growing on Mountain Islands in southwestern USA: Sky Island 7 ; P. arizonica and P. brachyptera were obtained from The Plant List (2013) 1 66 xi LIST OF FIGURES Figure 1 .1 Establishing the links between t ree - ring traits, genotypes, and the climatic constraints limitin g to tree growth. (a) Possible consequences of environmental heterogeneity on fitnes s, population structure, and distribution of taxa. (b) Benefits of linking basic ecology to applied research through the study of population structure 19 Fig ure 2. S 1 Figures of soil AWHC (%) used to validate PDSI data analyzed in study. 5 7 Figure 2. S 2 Correlations of modeled and locally collected climate data are collinear and validate the accuracy of the gridded climate dataset (i.e. same month correlations with high positive r - value, e.g. May - May). The McKenney et al. (2011) modelled variables are o n the y - axis and corresponding local NOAA variables are on the x - axis . (A) Total PCP: modeled PAL (y) to local Pal (x) 59 Figure 2. S 3 Plots of average annual P DSI values estimated for two st u dy P lots MTL (top) and BIG (bottom), for the reference period 1925 - 2009. Yearly values are depicted along the x - axis and average PDSI index along the y - axis. Negative values denote dry conditions 6 2 Figure 2. S 4 Climate graphs of mean monthly climate data analyzed for BIG (Mt. Bigelow) research site ( - 110.71495, 32.41378; 2535 m elevation ) for the reference period of 1925 - 2009 (McKenney et al., 2011). 6 3 Figure 2. S 5 Correlation analysis for P. arizonica (white) , P. ponderosa var. brachyptera (grey) growing at two transition zones within the Santa Catalina Mountains: MTL (hatched bars) and BIG (non - hatched bars). 6 5 Figure 2. S 6 Scatterplots of the mean tree - ring index (ARSTAN or Standard chronology) against the average primary PCP variable xii (summing across months) for the most highly correlated periods in the Seascorr analyses (Meko et al. 2011). These periods are defined as the three season - lengths with the highest correlat ion . (A) Plot of Mt. Lemmon P.arizonica tree - ring index (y) by PCP (x); all p - values are highly significant 67 Figure 2. 1 Locations of two study sites sampled for tree - ring analysis: MTL (Mt. Lemmon) and BIG (Mt. Bigelow). Black triangles mark the study sites northeast of Tucson, AZ, USA 76 Figure 2. 2 ARSTAN or Standard chronologies of P. arizonica (A,C) and P. ponderosa var. brachyptera (B,D) when Expressed Population Signal (EPS) is greater than 0.85, and sample depth is > 14, for the reference period of 1925 2009. Two sites were evaluated: MTL (A,B) and BIG (C,D). The corresponding sample depth (number of tree c ores) is indicated by shading 77 Figure 2. 3 Climate - growth relationships determined by principal components multiple regression for two species of ponderosa pine [ P. arizonica = white , P. ponderosa var. brachyptera = grey] growing at two transition zones within the Santa Catalina Mountains: MT L (hatched bars) and BIG (non - hatched bars). 78 Figure 2. 4 Climate - growth associations obtained by relating total p recipitation (PCP) as primary climate variable (top in each panel) and mean a verage temperature (TAVG) as secondary climate variable (bottom in each panel) to the ARSTAN ring - width chronologies of P. arizonica (A,B) and P. ponderosa var. brachyptera (C,D) at two sites MTL (A,C) and BIG (B,D) [Note: Standard chronology for (A)]. PCP was summed and TAVG was averaged for a period of three months whose ending months are shown on the x - axes (previous August to current September). Significant correlations and partial correlations (p < 0.05) are shown by dark bars 8 0 Figure 2. 5 Climate - growth associations obtained by relating average PDSI a s primary climate variable (top in each panel) and mean average temperature (TAVG) as secondary climate variable (bottom in each panel) to the ARSTAN ring - width chronologies of P. arizonica (A,B) and P. ponderosa var. brachyptera (C,D) at two sites MTL (A,C) and BIG (B,D) [Note: Standard C hronology for (A)]. PDSI was summed and TAVG was averaged for a period xiii of three month whose ending months are shown on the x - axes (previous August to current September). Significant correlations and partial correlations (p < 0.05) are shown by dark bars 8 1 Figure 3.S1 Average yearly standard chronology for 3 - , mixed - , and 5 - needle type in the early BIG tree - ring study (Dodge 1963; Fritts 1963). Each chronology represents the mean of 12 cores from 6 trees for the reference period of 1881 - 1960. A = 3 - needle pine, B = mixed - needle pine, and C = 5 - needle pine 1 09 Figure 3.S2 Average yearly standard chronologies for 3 - , mixed - , and 5 - needle types in the recent BIG tree - ring study . Each chronology represents an EPS > 0.82 and sample depth > 12, for the reference period of 1950 2007. A = 3 - needle type, B = mixed - ne edle type, and C = 5 - needle type 1 10 Figure 3.S3 Average yearly standard chro nologies for 3 - , mixed - , and 5 - needle types in the recent MTL tree - ring study. Each chronology represents an EPS > 0.84 and sample depth > 15, for the reference period of 1950 2007. A = 3 - needle type, B = mixed - needle type, and C = 5 - needle type 1 11 Figure 3.1 Locations of the two populations sampled for growth analysis: Mt. Lemmon (MTL) and Mt. Bigelow (BIG). Black symbols locate study plots northeast of Tucson, AZ where 3 - , mixed - , and 5 - needle pines are sympatric 1 15 Figure 3.2 Total PCP - growth relationshi ps were determined by principal components multiple regression for three needle types [3 - needle = white , mixed - needle = blue, 5 - needle = grey] growing at two transition zones within the Santa Catalina Mountains. W ithin a plotted color, standard index composite chronologies analyzed were BIG - E (double hatched bars) for the early period of 1895 1960, and BIG - R (hatche d bars) and MTL - R (non - hatched bars) for the recent period 1950 2007 . Climate variables analyzed were seasonal PCP: PrSummer (Jul - Sep), Winter (Nov - Mar), Spring (Apr - Jun), Summer (Jul - Sep)] 1 16 Figure 3.3 TAVG - growth relationships were determined by moving window correlations. From top to bottom are correlations for xiv 3 - , mixed - and 5 - needle types. A = BIG - E data, B = BIG - R data, and C = MTL - R data. Color scale = correlation coefficients; significant correlations (p < 0.05) are denoted by *. Intervals during which a variable was not significantly correlated to tree growth are shaded pale red or pale blue ( c. + 0.2). 1 17 Figure 4. 1 Black circles mark the locations of five study sites comprising 21 populations that were sampled for genetic analysis. 159 Figure 4. 2 Climate diagram for the Santa Catalina Mount ains, for the 84 year reference period spanning 1925 - 2009 for (a) average annual temperature (°C; red line) and (b) total annual precipitation (m m; blue line), respectively. The average temperature increased by 1.3°C over the course of the study (from 10.6°C to 11.9°C). Site specific climate data sets (MTL, BIG) were generated using the methods of McKenney et al. (2011) and averaged to gether to construct this diagram. The dotted lines are the linear trend lines. See Marquardt et al. ( 2018 ) for details 1 61 Figure 4. 3 Convex polygons of princi pal coordinates analysis of 21 a priori populations of ponderosa pine collected from five sites. Axis1 explained 81.58% of the variance, and Axis2 ex plained 8.53% of the variance. PCoA analysis was based on pairwise genetic distance of chloroplast haplotype frequencies . The three needle types are well resolved. Refer to Table 4.1 for population and site information 1 62 Figure 4. 4 Histogram plot of the fist axis from PCoA with 21 a priori populations of ponderosa pine plotted on the x axis . Y axis is the number of populations counted in each stack. The mixed - needle types are well resolved except for one mixed - needle population and one 5 - needle population 1 63 1 CHAPTER 1 MONITORING F OREST RESPONSES TO CHANGING CLIMATE ABSTRACT Contemporary climate change has unprecedented potential to alter forest structure, productivity, disturbance patterns, and community composition. To meet the challenge of managing forests with climate change, new in terdisciplinary tools need to be developed to monitor forest responses. An important area of research is employed that uses tree - ring widths and population genetic structure to model inter - annual climatic response, which is explained in terms of niche tra deoffs, physiological responses to climate, and adaptation by speciation or hybridization. The distribution of plant species in the Santa Catalina Mountains, Arizona, USA, and significance of the study will also be discussed. The main conclusion of the st udy is that two closely related ponderosa pine species have different ecological requirements for growth at their zone of sympatry. The resource most limiting to growth is seasonal water availability, which also controls the length of the growing season a topic of great interest to land managers in determining range limits for seed planting zones. 2 INTRODUCTION Long - term climatic changes will have profound effects on forest structure, productivity, disturbance patterns, and community composition. Increases in air pollution from the combustion of fossil fuels, deforestation, and land use change elevates the concentrations of greenhouse gases in the atmosphere, which correlates with increases in the mean global temperature. Significant warming has o ccurred over the past fifty years with increased occurrences of extreme weather events such as drought (IPCC, 2007) . Pl ant responses to changing climate can be persistence (if the change is tolerable), range shifts to more favorable areas, adaptation, or extinction (Davis et al., 2005) . The range shifts can be elevational as well as geographic range shifts, and these elevational clines should be predictive of geographic clines. The Pleistocene saw dramatic cha nge in climate 2 Myr (Table 1 .1 ). Plant and animal species retreated to higher elevations to escape warming temperatures and drought of the advancing desert in the American southwest and persisted in isolated mountaintop refuges called the Sky Islands (Re hfeldt 1999). The refuges are scattered throughout the southwestern United States and northern Mexico and are surrounded by deserts that act as barriers to migration. These mountain islands are ideal for studying changes in forest structure as significan t variation in moisture and temperature along steep elevational gradients permits the coexistence of species filling different ecological niches over a highly compressed spatial scale. Elevational gradients are useful for studying constraints on plant gro wth by climate (Adams and Kolb 2004). These studies can be enhanced by considering water availability simultaneously given variability in water supply in high elevation forests with changing climate. For example, as 3 concentrations of atmospheric carbon d ioxide and mean global temperatures rise, snow cover is predicted to decline, and growing seasons lengthen (IPCC, 2007). Potential outcomes of the projected change in environmental conditions due to reduced water availability are drier soils and forests. Furthermore, changes in the water and carbon cycle will impose additional water stress on plants leading to decreased growth. Gymnosperms existing prior to dramatic shifts in climate would have accumulated species during times of climatic stability. The species that survived changes in temperature may have done so because of wide physiological tolerances to climate (Coomes and Grubb, 1996) , biological interacti ons resulting in niche tradeoffs (Tilman, 2004) , or neutral forces such as stochastic dispersal and other population level processes (Kraft et al., 2008) . Alternatively, changing climatic conditions could have been the impetus for speciation. Young pine sp ecies that originated during the Early Tertiary 7 Myr (Table 1 .1 ) through a period of changing climate that brought dramatic changes in vegetation (Millar 1993) were a recently diverged fraction that coexisted with old species after radiation. Subsequentl y, periods of cooling may have repeatedly caused the pines to retreat, followed by independent evolution until drier periods reunited the pines. The repeated warming and cooling cycles could have increased hybridization between closely related species, lea ding to greater diversity . Thus, adaptation by speciation or hybridization would have been governed by harsh climate. Environmental factors can also have a strong influence on the location and maintenance of interspecific hybrid zones (Grant 1981), espec ially on more xeric sites where water availability would be a driving factor. (Swenson et al. 2008). Molecular markers are used to monitor hybrid zones in forest trees (Bennuah et al. 2004), which 4 have been reported to be narrow on southern slopes (Ito et al . 2008), or undetected (Epperson et al. 2009). The Santa Catalina Mountains are located north of Tucson in southeastern Arizona, USA, Rehfeldt (1999) . The plant species assemblage changes abruptly along a steep elevational gradient, 760 m near Tucson, to more than 2740 m at the summit. The vegetation varies from Sonoran Desert Scrub on the lower slopes to Pine - Oak Woodlands at mid - slope, and subalpine forest near the summits. The mountain range is influenced strongly by Mexican and L atin America flora creating a remarkable gradient of plants along its southern slope. Species richness increases along a moisture gradient from higher to lower desert elevations (Whittaker and Neiring, 1965) . Other major ranges in the area that share similar patterns in vegetation are the Pinaleno, Chirricahua, Santa Rita, Huachuca, and Rincon. Although similar in species composition, forest types may be reduced or absent at higher elevations, and ranges further eastward lack the desert communities present at the base of the Santa Catalina mountains (Whittaker and Neiring, 1965) . The southwest high elevation ponderosa pine forests predominately consist of two partially sympatric taxa the well - characterized Pinus arizonica and P. ponderosa var. brachyptera ( Willyard et al. 2017 ) . On the south facing slopes of Mt. Lemmon (Santa Catalina Mountains), the transition between these species is quite dr amatic over 100 meters horizontal distance and is often complete (Epperson et al. 2001; Epperson et al. 2009) . P. arizonica is more successful at the warmer and drier, lower elevations, and extends its range far south into central Mexico. Whereas, P. ponderosa var. brachy ptera survives exclusively on the highest mountaintops and cold air drainages; at 5 the colder and wetter, higher elevations. The sharp elevational transitions observed between these distributions are inconsistent with pine dispersal suggesting there are strong differences in fitness between the two species in response to environmental selection pressures. (Fig. 1 .1 a; Housset et al. 2018). Southern slopes have difficult growing condi tions, which are hot and dry, and P. arizonica is most likely more heat and drought resistant than P. ponderosa var. brachyptera . The goal of this work is to determine how changing climate variables are affecting the growth and genetic makeup of these po nderosa pine forests, as expressed through easy to identify traits used to monitor these changes. The project aims to determine how pine species are generally responding to changing climate factors to predict of how pines will respond genetically and phys iologically to increasing temperatures. The main questions of this study are: 1) can we identify the environmental variable for determining the P. ponderosa var . brachyptera P. arizonica distributions (Chapters 2 - 3), and 2) can we resolve the genetic re lationships between the two species (Chapter 4). To address these questions we will ascertain carbon sequestration by quantifying the growth and genetic responses of natural pine forests to increasing precipitation, temperature, and drought. The main obj ectives of the study are to: 1) correlate annual patterns of growth for two ponderosa pine species with seasonal climate variables (Chapters 2 - 3), 2) examine the temporal stability of the climate - growth correlations (Chapter 3), and 3 ) evaluate the populat ion genetic structure of two pine species (Chapter 4). Five hypotheses tested the growth response of ponderosa pine forests to climate variables: 1) P. ponderosa var. brachyptera growing at its lower - warm elevational limit 6 will be more sensitive to dry con ditions and changes in temperature than P. arizonica growing near its upper cool - moist elevational limit (Chapter 2 ; limits hypothesis ), 2) P. ponderosa var. brachyptera (originating from moist climate) will be more susceptible to moisture stress than P. a rizonica (originating from dry climate; Chapter 2 ; niche tradeoff hypothesis ), 3) the seasonality of limiting factors will change over time because rising temperatures will influence weather patterns and drought in the desert southwest (Chapter 3 ; shifts h ypothesis ), 4) the climate gradient may be a morphological response of P. ponderosa var. brachyptera and P. arizonica to changing environmental conditions (Chapter s 3 - 4 ; morphological hypothesis ), or 5) the climate gradient may reflect the presence of hybridizing populations of P. arizonica and P. ponderosa var. brachyptera (Chapter 4 ; hybridization hypothesis ). RATIONALE AND SIGNIFICANCE A great need exists for collecting fine scale data as part of a forest indicator and monitoring program to validate models that predict climate change effects on individual forest stands (Mar sh et. al. 2009). I ntegrating a physiological and genetics approach (Cavender - Bares and Pahlich 2009) will address lan dscape heterogeneity and uncertainty about growth responses by maximiz ing the common growth signal among intraspecific trees at a site (Babst et al. 2018). Examined will be individuals sampled from southern aspects where the distributions of two species o f ponderosa pine overlap and are potentially hybridizing in the Santa Catalina Mountains. D endrochronology is an ideal approach for monitoring growth trends over decades and centuries because the approach has been used to addresses interspecies relationships for pine species 7 growing along climate gradients where water stress is known to be growth limiting (Marquardt et al. 2018). Also, climate models predict more variable climate patt erns, which should affect the timing of pollen shed and female cone receptivity (Tauer et al. 2012). Consequently, changing environmental conditions are expected to play a role in the timing of reproduction, hybridization , and the formation of population genetic structure. Our discussion about collecting fine scale tree - ring data emphasized the need for field sampling to precisely capture species identification. We used average needle counts per fascicle to rapidly identify species in the field (Epperson et al 2001), a phenotypic measure with high heritability (Rehfeldt 1993). To confirm our initial classifications, more robust genetic analyses are required; although, experimental data is often lacking for studying genetic structure because adequate mol ecular markers are not always available (e.g. we developed 8 new microsatellite markers for ponderosa pine). These knowledge gaps are especially true for non - model conifer species like ponderosa pine with large and complex genomes (Zimin et al. 2014). Cu rrent developments in sequencing technologies are rapidly advancing and improving our ability to study neutral and adaptive genetic variation in wild species (Segelbacher et al. 2010). Such developments allow new possibilities for basic research on plant responses to environmental stress including growth and hybridization. Therefore, basic information was developed because the problem cannot be attacked through straightforward approaches using readily available, clear - cut methods. The lack of sufficient basic knowledge resolving genetic patterns with the growth responses of trees is the primary factor accountable for the knowledge gap about the causes of 8 interspecific differences in annual growth in response to changing climate. Resolving tree - ring growt h with genetic patterns requires a detailed analysis of morphological, tree - ring, and genetic characters to distinguish groups of individuals with different requirements for seasonal precipitation and habitat suitability. Varying site conditions limit the applicability of generalized guidelines; therefore, one goal was to determine how individual species will establish across sites on a local scale to overcome this limitation (Fig. 1 .1 a). Thus, the results of this study that integrates physiology and gene tics will provide a more comprehensive ability to predict the effects of climate change on tree growth and population genetic diversity, which will provide new information for forest management guidelines (Fig. 1 .1 b). Ponderosa pine, a widespread North Ame rican tree species, has high economical, ecological, and historical value for wood production, erosion control and wildlife habitat, and ethnobotanical uses (USDA NRCS National Plants Data Center). Incorporating climate gradients while studying the growth , physiological, and population genetics responses of ponderosa pine will improve our understanding of the role climate tolerances play in driving annual growth and diversity patterns, as temperatures increase and water becomes more limiting. We expand on previous work by Kilgore (2007) and Epperson et al. (2009) by characterizing replicate transition zones of two closely related taxa ( P. arizonica and P. ponderosa var. brachyptera ) for differences in sensitivity to climatic variability and genetic structu re. The project is of value to physiologists and geneticists studying the stress response of closely related plant species. The study results combine estimates of resource use and plant growth, with assessments of genetic diversity, allowing differences i n intra - and inter - specific plant 9 water relations to be coupled to climate. In addition, the conclusions are useful to land managers because the study of quantified differences in tree - ring widths and the genetic diversity of two species of ponderosa pine in the Santa Catalina Mountains will help to determine range limits for seed planting zones. Such an approach allows recommendations to be made that provide viable options for long - term improvement of reforestation efforts and sustainable management of n atural areas. Moreover, a better understanding of the response of these two species to climate change will provide comparisons for P. ponderosa var. scopulorum , a dominant pine in the ponderosa pine ecosystem. METHODS The study sites are located within the Santa Catalina Mountains, near Tucson, Arizona, USA. Climate variables were collected at three contact zones (of similar elevation, slope and southern aspect) that differ in soil water availability, and where sympatr ic morphotypes are potentially hybridizing . The average available water holding capacities are 3.8%, 4,7% and 9.2% for Mt. Bigelow (BIG), Palisades (PAL), and Mt. Lemmon (MTL), respectively (Marquardt et al. 2018) . The small scale of the study minimizes micro - environmental effects from (soil, climate, water, light, elevation) on the analysis of annual growth. The three populations have been undisturbed except for fire, with MTL and PAL more heavily impacted by fire than BIG in recent years. Identificati on of species was based on the average needle number per fascicle (Haller 1965). P. ponderosa var. brachyptera contained two needle types identified as mixed - needle (MN; 3.2 < mean < 4.6 needles per fascicle), and 3 - needle (3N; < 3.2 needles per fascicle; 10 Rehfeldt 1993, 1999). Trees that average > 4.6 needles per fascicle were designated P. arizonica (5N; Peloquin 1984). From each species, 30 continuous mature trees were sampled from MTL and BIG for tree ring s a t each site for genetic analysis. One tree was sample every 10m along a 5x9 grid. The exception was at the center of the transition zone, where one P. ponderosa var. brachyptera and one P. arizonica were sampled at the three innermost transects. The average series length (years + SD) for BIG and MTL combined is 75 + 32 years (Marquardt et al. 2018). Additional trees were sampled for genetic evaluation at six - locus haplotypes (chloroplast) collected from populations located higher in elevation (MTL) and lower in elevation [Green M ountain (GRN), Rose Canyon (ROS)] than the transition zones. Seedlings located near the transition zones, and above the transition zones were also sampled for genetic analysis at MTL and PAL. CONCLUSION S AND FUTURE WORK Because of changing climate , gai ning knowledge about climate - growth relationships of closely related species is vital for forest management. Precipitation is known to be a major factor limiting growth of ponderosa pine in high elevation forests such as the Santa Catalina Sky Islands, bu t the mechanisms have not been thoroughly understood. Therefore, we focused on changes in the annual growth of trees and standardized the growth response in stands (a relative measure) to account for variation introduced by microsite differences in enviro nmental conditions. This allowed us to focus on the growth response to climate stress when comparing sites in contrasting environments. 11 Currently, ponderosa pine is managed as one distinct species under the jurisdiction of the Coronado National Forest. A lthough proper forest management needs accurate knowledge about the structure of genetic variation with respect to environmental conditions, the population genetic structure of ponderosa pine in the Santa Catalina Mountains has yet to be fully determined. Therefore, we also conducted a preliminary study of population genetic structure using distance and frequency based genetic measures to determine the number of operational taxonomic units. In our study, we analyzed the differences in annual growth, climate sensitivity and genetic structure between P. arizonica and P. ponderosa var. brachyptera (both as two species and as three needle - types) at three transition zones in southeast USA, near Tucson, Arizona. Moisture availability is the major factor controlling growth through the length of growing season for these conifers, especially for cool season correlations during winter and spring. The major conclusion from this study are that two closely related ponderosa pine species have different ecologica l requirements for growth and habitat suitability, which resulted in three distinct taxa occupying different elevational niches. Statistical analyses show subtle and significant ecological divergence in climate sensitivities and genetic structure, implyin g that climate change may be a major factor in the differentiation of P. arizonica and P. ponderosa var. brachyptera . The four most important points of the first tree ring study (Chapter 2) are that 1) sympatric pine species growing at high elevation have seasonal differences in precipitation requirements for growth, 2) the most consistent pattern of differences in climatic sensitivity between species is the variability in soil moisture availability, 3) the variation in water availability between sites is further influenced by disturbance especially at the less dry site where 12 growth reductions are larger, and 4) the variability in local site conditions helps to - ring study (Chap ter 3) a dramatic shift of seasonal climate - growth relationships from summer to spring is reported, suggesting 1) water limitations control the length of the growing season , 2) a redirection of resources from maximum biomass production (summer water) to war ds reproductive structures (spring water) , and 3 ) warming trends may be impacting the stability of climate growth relationships, which has implications for conservation management and climate reconstruction . Also, cambial plasticity (a response to variab le winter precipitation) may provide opportunity for hybridization. Lastly, the genetic study of haplotype frequencies (Chapter 4) show 1), genetic variability is reduced for the 5 - needle pine 2), higher divergence was observed for distant (i.e. pure) tha n sympatric populations, and 3) gene flow is high among all populations. Moving forward, we will complete a detailed genetic analysis of the study populations with six - nuclear microsatellite markers. To date w e have developed eight new microsatellite markers for ponderosa pine (two to be used in this study) , verified the motifs by Sanger sequencing, obtained genotypes and complied the allele frequency data. In addition, f ew studies have analyzed tree - ring and genetic data jointly in natural populations (Pluess and Weber 2012; Cole et al. 2010; Babuskina et al. 2016; Heer et al. 2018), and only Heer et al. (2018) found a strong signal between growth and a genetic trait 15 genes linked tree - ring phenotypes to photosynthesis and drought stress. The lack of significant results in the other studies could be attributed to too few 13 genetic data (i.e. lack of power). Therefore, to facilitate the integration of tree - ring and genetic approaches, we propose future work that scales tree data to local conditions by integrating individual chloroplast haplotypes to assess climate response in natural populations. This would be accomplished by rebuilding composite tree - ring chronologies based on haplotypes and conducting correlations with climate to determine which hap lotypes are most sensitive to drought and precipitation. We acknowledge that the tree - ring sample size will be a concern and power of the test would be monitored as the chronologies were being constructed. 14 APPENDICES 15 APPENDIX A Manuscripts 16 This dissertation is the foundation for three peer reviewed papers, which are referred to below by their chapter number in the dissertation: 2. Marquardt, PE, Miranda, BR, Jennings , S, Ginger, T and Telewski, FW (2018) Variable Climate Response Differentiates the Growth of Sky Island Ponderosa Pines. TREES DOI: 10.1007/s00468 - 018 - 1778 - 9 3. Marquardt, PE, Miranda, BR, and Telewski, FW. Shifts in Climate - Growth Relationships of the Sky Island Pines. Manuscript 4. Marquardt, PE, Willyard, A., Miranda, BR, and Telewski, FW. Genetic Differentiation among Sky Island Populations of Ponderosa Pine. Manuscript 17 APPENDIX B Tables and Figures 18 Table 1 . 1 Geologic time scale modified from White et al. (2007; Table 1) Era Period Epoch Myr 1 Cenozoic Quaternary Holocene 0.01 Pleistocene 2.25 Tertiary Pliocene 7 Miocene 26 Paleogene Oligocene 34 Eocene 54 Paleocene 65 1 Million years ago 19 (A ) Costs of site heterogeneity ( B ) Benefits of basic and applied research Fig ure 1 . 1 Establishing the links between tree - ring traits, genotypes, and the climatic constraints limiting to tree growth. (a) Possible consequences of environmental heterogeneity on fitness, population structure, and distribution of taxa. (b) Benefits of linking basic ecology to applied research through the study of population structure 20 REFERENCES 21 REFERENCES Adams HD, Kolb TE (2004) Drought responses of conifers in ecotone forests of 140:217 - 225 Babst F, Bodesheim P, Charney N, Friend AD, Girardin MP, Klesse S, Moore DJ, Seftigen K, Björklund J, Bouriaud O, Dawson A (2018 ) When tree rings go global: Challenges and opportunities for retro - and prospective insight. Quat Sci Rev 197 : 1 - 20 Babushkina EA, Vaganov EA, Grachev AM, Oreshkova NV, Belokopytova L V , Kostyakova TV, Krutovsky KV (2016) The effect of individual genetic heterozygosity on general homeostasis, heterosis and resilience in Siberian larch (Larix sibirica Ledeb ) using dendrochronology and microsatellite loci genotyping. Dendrochronologia 38 : 26 - 37 Bennuah SY, Wang T, Aitken SN (2004) Genetic analysis of the Picea sitchensis x glauca introgression zone in British Columbia. For Ecol Manage 197: 65 - 77 Cavender - Bares J, Pahlich A (2009) Molecular, morphological and ecological niche differentiation of sympatric sister oak sp ecies, Quercus virginiana and Q geminata (Fagace ae). Am J Bot 96: 1690 - 1702 Cole CT, Anderson JE, Lindroth RL, Waller DM (2010). Rising concentrations of atmospheric CO2 have increased growth in natural stands of quaking aspen (Populus tremuloides). Glob Change Biol 16 : 2186 - 2197 Coomes DA, Grubb PJ (1996) Amazonian caatinga and related communities at La Esmeralda, Venezuela: Forest structure, physiognomy and floristics, and con trol by soil factors. Vegetatio 122:167 - 191 Davis MB, Shaw RG, Etterson JR (2005) Evolutionary responses to changing climate. Ecology 86:1704 - 1714 22 Epperson BK, Telewski FW , Plovanich - Jones AE, Grimes JE (2001) Clinal differentiation and putative hybridization in a contact zone of Pinus ponderosa a nd P arizonica (Pinaceae). Am J Bot 88:1052 - 1057 Epper son BK, Telewski FW, Willyard A (2009) Chloroplast diversity in a putative hybrid s warm of Ponderosae (Pinaceae). Am J Bot 96:707 - 712 Grant V (1981) Plant speciation Second edn. New York: Columbia University Press Heer K, Behringer D, Piermattei A, Bässler C , Brandl R, Fady B , Jehl H, Liepelt S, Lorch S, Piotti A, Vendramin GG (2018). Linking dendroecology and association genetics in natural populations: Stress responses archived in tree rings associate with SNP genotypes in silver fir (Abies alba Mill.). Mol Ecol 27 : 1428 - 1438 Housset JM, Nadeau S, Isabel N, Depardieu C, Du chesne I, Lenz P, Girardin MP (2018) Tree rings provide a new class of phenotypes for genetic associations that foster insights into adaptation of conifers to climate change. New Phytol 218:630 - 45 IPCC (2007) Climate change 2007: the Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Pane l on Climate Change -- Solomon S, Daahe Q, Manning M, Avery t K , Marquis M (eds) . Cambridge, Unite d Kingdom and New York, NY, USA Kilgore, J (2007) Distribution and ecophysiology of the Ponderosae in the Santa Catalina Mountains of southern Arizona . Thesis (PhD) , Michigan State University, Michigan , 263 pp Kraft NJB, Valencia R, Ackerl y DD (2008) Functional traits and niche - based tree community ass embly in an amazonian forest. Science 322:580 - 582 Marquardt PE, Miranda BR, Jennings S, Ginger T, Telewski FW (2018) Variable climate response differentiates the g rowt h of sky island ponderosa p ines. Trees 2018:1 - 6 . Marsh A, Mawdsley J, Negra C (2009) Forest observations and indicators needed to re spond to climate change. J For 107: 231 - 232 23 Millar CI (1993) Impact of the Eocene o n the e volution of Pinus - L. Ann Mo Bot Gard 80:471 - 498 Peloquin RL (1984) The identification of three - species hybrids in the ponderosa pine c ompl ex. Southw Nat 29:115 - 122 Pluess AR, Weber P (2012) Drought - adaptation potential in Fagus sylvatica: linking moisture availability with genetic diversity and dendrochronology. PLoS O ne , 7 : e33636 Rehfeldt, GE (1993) Genetic variation in the Ponderosae of the southwest. Am J Bot 80: 330 343 Rehfeldt GE (1999) Systematics and genetic structure of Ponderosae taxa (Pinaceae) inhabiting the mounta in islands of the Southwest. Am J Bot 86:741 - 752 Segelbacher G, Cushman SA, Epperson BK, Fortin MJ , Francois O, Hardy O, Holderegger R, Taberlet P, Waits LP, Manel S (2010) Applications of landscape genetics in conservation biology: concepts and challenges. Conserv Genet 11:375 - 385 Swens on NG, Fair JM, Heikoop J (2008) Water stress and hybridization between Quercus Gambelii a nd Quercus Grisea. West N Am Nat 68: 498 - 507 Tauer CG, S tewart JF, Will RE, Lill y CJ, Guldin JM, Nelson CD (2012) Hybridization leads to loss of genetic integrity in shortleaf pine: unexpected consequences of pine management a nd fire suppression. J For 110: 216 - 224 Tilman D (2004) Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community asse mbly. Proc Natl Acad Sci of the United States of America . 101:10854 - 10861 USDA NRCS National Plants Data Center. https// plants.usda.gov/plantguide/pdf/pg_pipo.pdf/ . Accessed 2 Oct 2018 24 White TW, Adams WT, Neale DB (2007) Evo lutionary genetics. In: Forest G enetics Ca mbridge, MA: CABI publishing Whittaker RH, Neiring WA (1965) Vegetation of the Santa Catalina mountains, Arizona: a gradient analysis of the south slope. Ecology 46:429 - 452 Zimin A, Stevens KA, Crepeau MW, Holtz - Morris A, Koriabine M, Marçais G, Puiu D, Roberts M, Wegrzyn JL, de Jong PJ, Neale DB (2014) Sequencing and assembly of the 22 - Gb loblolly pine genome. Genetics 196: 875 - 890 25 CHAPTER 2 VARIABLE CLIMATE RESPONSE DIFFERENTIATES THE GROWTH OF SKY ISLAND PONDEROSA PINES Paula E. Marquardt 1 * , Brian R. Miranda 1 , Shane Jennings 2 , Ginger Thurston 2 , Frank W. Telewski 2 1 USDA Forest Service, Northern Research Station, Institute for Applied Ecosystem ORCID ID 0000 - 0002 - 3506 - 8297 2 W.J. Beal Bo tanical Garden, Department of Plant Biology, Michigan State University, - 0001 - 7957 - 9841 *Corresponding author: pmarquardt@fs.fed.us ; 715 - 362 - 1121; ORCID ID 0000 - 0001 - 8854 - 0172 Acknowledgments : This paper is part of a dissertation submitted to Michigan State University in partial fulfillment of requirement s for a Doctor of Philosophy degree. The Institute for Applied Ecosystem Studies, USDA Forest Service, and the Paul Taylor Fellowship, Dept. of Plant Biology and Beal Botanical Garden, both of Michigan State University, supported the project. We thank th ree anonymous Reviewers and a Communicating Editor for constructive comments that improved the manuscript, and the following people for help on the project: J. Rayala, J. Grimes, G. Friedlander, J. Kilgore, and B. Epperson assisted with field data collect ion. A. Foss and S. Lietz provided lab, data processing, and GIS support. J. and L. Griffith, and A. and T. Harlan (both 26 deceased) provided logistic support for the field research. The Coronado National Forest provided access to samples and the soils da ta. A. Willyard and B. Sturtevant provided constructive feedback. J. Smith reviewed the paper. L. Rayala edited the paper. 27 ABSTRACT Key message : The seasonally cool and moist conditions of spring improved the growth of two co - occurring ponderosa pine species, which displayed different seasonal climatic responses and length of correlations to drought. Abstract: We examined the climatic sensitivi ty of two partially sympatric pine species growing at their transition zone in the Santa Catalina Mountains, AZ, USA. Pinus arizonica is found at lower elevations compared to P. ponderosa var. brachyptera . Ring widths were measured in trees at two sites and correlated with precipitation, temperature, and Palmer Drought Severity Index to assess the influence of climate on growth. The two species were analyzed within and between sites, which have similar elevation, aspect, and species composition, although soils at the two sites have different water - holding capacities . Response function analyses of P. arizonica [sampled near its upper (and wetter) elevation limit], and P. ponderosa var. brachyptera [sampled near its lower (and drier) elevation limit] indic ated that annual growth correlated positively and strongly with spring precipitation at both study locations. Local site conditions had a major impact on tree growth and variability in site conditions helped resolve the climate. For example, at the less dry site, growth of the lower elevation pine ( P. arizonica ) responded to early - winter precipitation, while P. ponderosa var. brachyptera did not. Also, correlation analysis indicated that P. growth was more sensitive to drought for longer periods than P. ponderosa var . brachyptera . Finally, partial temperature - growth correlations of P. arizonica and P. ponderosa var . brachyptera indicated growth was limited by increased growing season and winter respiration, respectively. Rising night - time temperatures during spring 28 significantly reduced growth of P. arizonica at Mt. Lemmon. These findings demonstrate subtle yet meaningful interspecies differences in sensitivity to seasonal moisture stress and use of carbon resources. Keywords: Dendroecology; drought stress; Pinaceae; response function ; tree ring; Ponderosae INTRODUCTION The Santa Catalina Mountains, near Tucson, Arizona, USA, are among the best - known and well - studied Madrean Sky Islands, which are high elevation mountains defined by pine - oak forests (Shreve 1915, 1917, 1919; Whittaker and Niering 1965, 1975; Whittaker et al. 1968; Bezy 2016). Two partially sympatric pine species grow at the higher elevations: Pinus arizonica Engelm. and P. ponderosa var. brachyptera (Engelm. ) Lemmon, a closely related variety of P. ponderosa [and previously misidentified as P. ponderosa var . scopulorum ( Engelm.)] , also known as Taxon X (Rehfeldt 1999; Epperson et al. 2009; Willyard et al. 201 7). P. arizonica and P. ponderosa var. brachyptera , Diploxylon pines of the subsection Ponderosae , are large trees that grow throughout the mountainous and semiarid regions of the southwestern United States (Perry 1991; Farjon and Styles 1997; Price et a l. 1998). P. ponderosa var. brachyptera , comprising three - needle and mixed - needle morphotypes, survive exclusively at the highest elevations (2300 - 3000 m) and in cold air drainages (Rehfeldt 1999; Epperson et al. 2001), while the five - needle P. arizonica is found at lower elevations (1800 - 2600 m) , and is considered part of the ponderosa pine complex. As mean annual precipitation increases with elevation from 29 orographic processes (Sheppard et al. 2002), P. arizonica is more successful at surviving in warme r and drier habitats, whereas P. ponderosa var. brachyptera survives in cooler and wetter habitats. P. arizonica extends from southeastern Arizona and southwestern New Mexico to the Sierra Madre of northern Mexico (Perry 1991). In comparison, P. ponderos a var. brachyptera ranges from northern Nevada to southern Texas (USDA, NRCS 2017). On south - facing slopes of the Santa Catalina Mountains, shallow soils are mostly homogeneous lithosols of low water - holding capacity that could promote moisture limitation s on annual growth (Shreve 1915, 1919; Whittaker and Niering 1965; Whittaker et al. 1968; Bezy 2016). The transitions in taxa on south slopes are quite dramatic occurring over 100 meters horizontal distance (Epperson et al. 2001). Evaluating geographic va riation in locally adapted populations helped determine 2006; Shinneman et al. 2016; McCullough et al. 2017). Although previous studies indicate North American pine growth is influenced primarily by seasonal precipitation (Dodge 1963; Fritts et al. 1965, Fritts 1976; Norris et al. 2006; Griffin et al. 2013; Dannenberg and Wise 2016; Shinneman et al. 2016; Gonza´lez - Ca´sares et al. 2017 ), few ecological studies have co mpared morphologically distinct taxa to identify interspecific differences (Haller 1965; Peloquin 1984; Rehfeldt 1993; Epperson et al. 2001). Our study design allowed us to identify individual taxa in the field and evaluate differential sensitivity to clim ate and water stress. This was accomplished by correlating seasonal climate data with widths of annual tree rings collected from transition zones where both taxa were present. Taxonomic relationships were determined by estimating 30 (from a small number of branchlets) the frequency of needles per fascicle, an easy way to identify traits in the field and useful for characterizing populations (Epperson et al. 2001). We sampled ring widths to provide information on the climatic factors limiting growth from tre es sampled under environmental stress, which occurs near the lower and upper limits of species distributions (Lamarche and Stockton 1974). The growth of trees in this semiarid region is primarily limited by water stress rather than temperature (Fritts 197 6). Under low moisture conditions, individual trees produce rings that are very narrow. Adapted to different ecological conditions, P. arizonica and P. ponderosa var. brachyptera are high - elevation tree species that express dissimilar tolerances to wate r shortage through their growth responses. Projected climate warming in these semiarid mountainous regions could reduce tree growth as increasing temperatures and decreasing precipitation elevate moisture stress (Barichivich et al. 2014). We expect available for growth is lowest (Sheppard et al. 2002). Thus, at the transition zone where the two species co - occur, our first hypothesis predicts P. ponderosa va r. brachyptera growing at its lower - warm elevational limit to be more sensitive to dry conditions and changes in temperature than P. arizonica growing near its upper cool - moist elevational limit (Lamarche and Stockton 1974; Adams and Kolb 2005). Palmer Dr ought Severity Index (PDSI; a measure of soil moisture availability), was evaluated in a second hypothesis predicting that positive growth - PDSI correlations would differ between taxa. The premise is that trees originating from moist ( P. ponderosa var. brac hyptera ) 31 compared to dry ( P. arizonica ) environments would be more susceptible to moisture stress. METHOD S Study area and Sampling The study areas of Mt. Lemmon (MTL; 32.443°, - 110.788°) and Mt. Bigelow (BIG; 32.414°, - 110.715°) are located within the San ta Catalina Mountains (2,500 m a.s.l.) of the Coronado National Forest, approximately 28 kilometers straight - line distance northeast of Tucson, Arizona, Pima County, USA (Fig. 2. 1). The ponderosa pine species were sampled from steep slopes [average gradie nt of 45% (MTL) and 37% (BIG); data not shown], of Lithic Haplustolls (Buol 1966, Brown 1968) derived from granite and gneiss, gravelly to rocky in texture and shallow in depth ( c. 50 to 140 cm; data not shown); a frigid complex with average available wate r - holding capacity (AWHC) of 9.2% for MTL and 3.8% for BIG ( Suppl. Fig. 2.S 1 ). The mixed conifer forest comprises P. arizonica, P. ponderosa var. brachyptera, P. strobiformis , Pseudotsuga menziesii , Quercus hypoleucoides , Q. gambelii, Q. reticulata (Brown 1968), and rarely Abies concolor . Half of the annual precipitation falls during the summer monsoon (July through September; JAS). Winter storms (November through March: NDJFM) provide an additional third or more of the annual rainfall (Table 2. 2; Sheppard et al. 2002). Ponderosa pine species were identified in the field by conducting average needle counts ( c. > 10 fascicles). Trees with average needle numbers per fascicle < 4.6 and > 4.6 were designated P. ponderosa var. brachyptera and P. arizoni ca , respectively (Haller 1965; Peloquin 1984). While avoiding ridges and cold air drainages, we 32 sampled transition zones within the lower and upper moisture availability limits, respectively, for P. ponderosa var. brachyptera , and P. arizonica. Nine hori zontal transects were established at each study area (MTL and BIG) on southern aspects of similar slope. The nine transects were equivalent to a 5 x 9 grid with 10 - m spacing. The nearest tree was sampled to each point, except for the middle of the transi tion zone (i.e. the three innermost transects), where one P. ponderosa var. brachyptera and one P. arizonica were sampled (Table 2. 1). From 2010 to 2012 all healthy trees > 8 cm diameter at breast height (DBH) were selected and sampled. Individual sample trees were tagged and mapped using North American Datum of 1983 (NAD83) geographic coordinates obtained through global positional system (GPS) data collected at the field plots (n = 120; i.e. 2 sites x 2 species x 30 trees). Two or more tree - core samples were collected from each tree using a 5.5 - mm increment borer, dried and glued in wood core mounts, and sanded to 600 grit to visualize clear ring boundaries. We visually crossdated all samples using the skeleton - plot technique (Stokes and Smiley, 1996). Tree - ring widths were measured to the nearest 0.01 mm using scanned images (2400 dpi resolution) and CooRecorder TM software (Cybis Elektronik, 2010). Accurate calendar dates were assigned to each ring in the time sequence with CDendro TM software (Cybis E lektronik, 2010). Sample dating was statistically verified with the program COFECHA (Holmes et al. 1986). Climate data Specific climate data sets for MTL and BIG were developed with 1 - km resolution using the ANUSPLIN package (Hutchinson and Xu 2013) by creating 1925 - 2009 33 estimates of the temperature and precipitation variables (McKenney et al. 2011). Similarly, the 4 - km gridded PDSI values obtained from the West Wide Drought Tracker website (Abatzoglou et al., 2017) were used to create estimates of the local PDSI variable (model and website development by the University of Idaho and the Desert Research Institute, Reno, Nevada). The composite PDSI index is based on variation in temperature, precipitation, and the local AWHC of the soil, and quantifies l onger - term departures from normal moisture patterns; more negative values indicate drought stress, ranging from < - 4 for extremely dry conditions to > +4 for extremely wet conditions (Palmer 1965). Climatic variables were summarized as monthly values of a verage temperature (TAVG), average minimum temperature (TMIN), average maximum temperature (TMAX), total precipitation (PCP), and average PDSI. To validate the quality of climate data, we conducted correlation analysis using cor function in R (R core Team 2016) with local datasets and displayed a matrix of correlation coefficients created using the corrplot package (Wei and Simko 2017). Correlation was also used to examine the interrelationship between climate variables at both sites (Appendix A) . Trend and Climate - growth analysis Juvenile growth trends were removed from the average ring - width series of individual trees (i.e. de - trended) by fitting a modified negative exponential curve to create standardized indices using the ARSTAN program (Cook 1985). Autocorrelation arises from the persistence of climatic effects on tree growth, which was removed by autoregressive modeling (Cook 1985). The final master ARSTAN chronology was 34 created using a bi - weight robust mean to account for climatic variance and end ogenous disturbance (Cook 1985). We calculated the descriptive statistics of running RBar, expressed population signal (EPS), signal - to - noise ratio (SNR), mean sensitivity (MS), Gini coefficient (G), standard deviation (SD), and autocorrelation (AC) using the dpIR package (Bunn 2008, 2010) for R, which was also used to plot the ARSTAN chronologies. Smoothing splines were applied to the chronologies with defaults in dpIR. All other tree - ring statistics were obtained directly from the ARSTAN DOS output fil es (Cook 1985). The running RBar [i.e. between - tree correlation (Rbt)] is the average correlation coefficient between all possible pairs of indexed series drawn from different trees (Wigley et al. 1984; Briffa and Jones 1990); we chose a 40 - year moving wi ndow with 20 - year overlap. The coefficient of variation (CV) of Rbt was expressed as the SD of the mean for all species and sites combined. EPS measures the common variability of a chronology when series are averaged, ranging from 0 to 1, and a value gre ater than 0.85 indicates a strong common signal of the chronology (Wigley et al. 1984; Briffa and Jones 1990, Speer 2010). SNR in de - trended tree - ring series is the ratio of two variances (common climate signal and random error). Mean sensitivity measure s the relative change in ring widths from year - to - year (Fritts 1976). It theoretically ranges from 0 (all rings have same width) to 2 (locally absent ring every other year), but in practice varies from 0.1 to 0.6 (Biondi and Qeadan 2008). Gini coefficient ranges from 0 to 1, measures diversity in tree - ring chronologies regardless of the degree of autocorrelation, and improves on intra - annual variability comparisons among species and sites (Bunn et al. 2013; DeRose et al. 2015). 35 Because monthly climate vari ables are often highly inter - correlated, tree - ring analyses were conducted using both partial correlations (Meko et al. 2011) and regression involving principal components (Fritts 1974). Partial correlations address collinearity of the primary and secondary climate variables prior to summarizing the seasonal climate signal in the tree - ring data. Partial correlations were used to compare each residual plot chronology to monthly climatic variables to test the null hypothesis of no effects ( r < 0.05). Pearson correlation coefficients were calculated for the primary climate variable (PCP summed; or PDSI averaged), and partial correlations were computed for the secondary climate variable (TAVG), independent of the variance related to the primary variable. Climate variables (PCP, TAVG, and PDSI) were tested over a 14 - month climate window from the preceding August to September of the current growth year, to account for the influence of prior environmental conditions on the current Subsequently, maximum climate - growth correlations were determined across season lengths of 1, 3 and 6 months reflecting seasonal fluctuations of PCP (primary) and TAVG (secondary) climate variables. Season lengths of 1, 3, 6, 12, and 20 months determined the maximum influence of PDSI (primary) and TAVG (secondary) variables on climate - growth correlations over a time window from the previous October to current August. All seasonal correlations were completed with the treeclim package (Zang and Biondi 2015 we quantified the effect size of PDSI using Pearson correlation coefficients ( r ) for single months using cor function in R. We investigated significant variables (p < 0.05) identified by simple correlation analysis (treeclim package; Zang and Biondi 2015), which were consistent in sign and 36 magnitude, as potential predictors to include in response function analysis. Fritts (1974) introduced regression involving principal components to address multicollineari ty by transforming predictor climate variables to produce a set of uncorrelated data points. The calculation of the response function regresses tree - ring data against transformed monthly climatic variables to select the parameters that influence tree grow th (Fritts 1974, 1976). Each principal component variable explains partial variance in the data set. Relationships between climatic variables and ring - width indices were examined using multivariate estimates obtained from the principal component regressi on model. The analysis computed bootstrapped response functions using the treeclim package (Zang and Biondi 2015). To obtain robust parameter estimates, we used bootstrapping to test regression coefficients and the stability of estimates ( Guiot 1991). T wenty - eight monthly climatic variables (i.e. 14 PCP + 14 TMIN) were analyzed from July through December of the previous year (excluding the transition month of October), and January through September of the current year. Climate variables were partitioned into seasons to quantify water balance (i.e. PDSI), precipitation, and temperature effects during key periods in our study plots (winter, spring and summer). Seasonal variables were derived by combining the monthly data into two rainy seasons consisting of three and five months respectively: summer (JAS) and winter (NDJFM) separated by three months of arid spring (AMJ). Summer season spanned the months from previous July to previous September, and current July to current September; winter season spanne d the months from previous November to the current March, and arid spring spanned the months from current April to current June. Climate was also divided into four seasons of three months each, which partitioned winter into fall (OND) and winter (JFM). M ean 37 chronologies and correlation bar graphs were plotted using plotting functions in the dplR package; all other plots were created using standard R plotting functions. Seasonal variables were tested for departure from normality with the Kolmogorov - Smirno = 0.05). Normality tests, and distributions of single variables were analyzed using JMP Pro, v 13.0 (Copyright © 2016 SAS Institute Inc., Cary, NC), or Sigma Plot, v 13.0 (Copyright © 2014 Systat Software, Inc., San Jose, CA). RESULTS Annual and Monthly climate Climate data of the two sites were averaged to obtain a summary record for the Santa Catalina Mountains. Over the 84 - year study period (1925 - 2009), average yearly PCP is 681.1 + 163.9 mm and average yearly TAVG is 10.9 + 0.5 °C (Table 2. 2). Multicollinearity between PCP and TAVG is low, but significant, with negative correlations of r = - 0.3 (Prob > IpI = 0.005; Zar 2010). The average PDSI broadly ranges between + 4.0 with a few values from 5.0 to 7.0 ( Suppl. Fig. 2.S 3 ). We identified four dry periods (PDSI < - - 1939; 1947 - 1948; 2002 - 2003; 2006 - 2007. Seasonal climate PCP was evenly distributed between two rainy seasons: summer (JAS; average 299 + 73.3 mm) and winter (NDJFM; aver age 299.2 + 155.0 mm), with August being the wettest month (Table 2. 2; Suppl. Fig. 2.S4 ). The driest season was spring (AMJ; average 48.3 + 38 median monthly TAVG remained below 6.5 °C , with January being the coldest month (2.7 °C ; Table 2. 2; Suppl. Fig. 2.S4 ). Summer remained above the median monthly TAVG of 18.5 °C, with July being the hottest month (19.5 °C). Spring (the driest season) saw the steepest increase in median TAVG from 9 °C in April, to 18.5 °C in June (a rise of 9.5 °C; Suppl. Fig. 2.S 4B ). Series length, Chronology statistics, and Growth patterns Sawtimber trees > 25 cm diameter DBH were the predominant size class sampled at MTL with some pole - sized trees (10 to 25 cm); mean DBH ( + SD) is 37.1 + 14.4 cm, ranging from 15.5 to 85.3 cm. At BIG the sampled trees were smaller, including pole - sized and two saplings (5 to 10 cm); mean DBH ( + SD) is 27.7 + 13.2 cm, ranging from 7.7 to 61.2 cm. Although there were differenc es in average DBH between sites, the average series length (years + SD) and percentage of locally absent rings were similar between sites. Observed mean series length are 73 + 23 years ranging from 35 to 148 years, and 76 + 41 years ranging from 30 to 162 years, respectively for MTL and BIG. Locally absent rings accounted for 0.32% (MTL) and 0.24% (BIG) of the total rings in all series, combined. The Gini coefficient and mean sensitivity index quantified moderate year - to - year differences that varied littl e between taxa or site (Table 2. 3). In comparison, the indexed ring widths at MTL are significantly narrower for P. arizonica than P. ponderosa var. brachyptera , and P. ponderosa var . brachyptera exhibits the highest correlation in growth among trees (Rbt ), higher SNR than P. arizonica at both sites, and significantly higher SD (14%) than P. arizonica at BIG, (p < 0.05; Table 2. 3). Considering Rbt 2 , the 39 coefficient of variation between species is 18% ( = 0.14; SD = 0.02; Table 2. 3). First - order autocor relation values for the standard chronology are higher for BIG (0.41) on average than for MTL (0.22; Table 2. 3). Growth was reduced from the late 1930s to the late 1940s for P. arizonica and P. ponderosa var. brachyptera at BIG (Fig. 2. 2C - D). A similar re duction in growth for the 1930s - 1940s was not evident at MTL (Fig. 2. 2A - B). The average growth of the two species was calculated to obtain indexed ring widths for BIG (0.70 + 0.12), and MTL (0.91 + 0.08) during seven drought years. The seven years of dec reased growth corresponded to low PDSI values recorded for 1936 - 39 and 1946 - 48 ( Suppl. Fig. 2.S 3 ). Average annual PDSI for the period was greater (wetter) at BIG ( - 2.15) than at MTL ( - 2.45). In contrast, t he combined average index of annual rings measured at MTL (0.76 + 0.23) during the latter eight - year dry period (2002 - 09) is less than for BIG (0.90 + 0.28). Average annual PDSI values were similar between sites for the 2002 - 09 period (MTL, - 2.39; BIG, - 2 .42). Correlations Simple correlation and Response function Correlation analysis ( Suppl. Fig. 2.S 5 ) identified seasonal variables for further response function analysis. Fig. 2. 3A shows significant positive growth responses with spring PCP (AMJ; p < 0 .05) for both taxa , and the correlations were greater at MTL than at BIG. Figs. 2. 3A - B also show significant positive PCP - growth correlations for P. arizonica at MTL during winter (ndJFM; p < 0.05), and early winter (ond; p < 0.05), and 40 growth is negatively correlated with spring TMIN (AMJ; Fig. 2. 3A; p < 0.05). PCP correlation and TAVG partial correlation Tree - ring index regressions with precipitation approach linearity, indicating they were not controlled by outliers ( Suppl. Fig. 2. S 6A - D; p < 0.05 ). Correlations between annual growth and single - month PCP for both species are positive and significant, increasing in strength when summing over multiple months (Fig. 2. 4). The maximum positive PCP - growth correlations for P. arizonica occur for the three - month periods ending in December, May, and June ( 0.4) at MTL (Fig. 2. 4A), and in June and July ( r 0.4) at BIG (Fig. 2. 4B). For P. ponderosa var. brachyptera , maximum correlations occur for the three - month periods ending in December ( r 0.35) and June ( 0.4) at MTL (Fig 2. 4C), and in June and July at BIG ( 0.35, 0.4, respectively; Fig 2. 4D). Partial correlations between annual growth and single - month TAVG at MTL for P. arizonica and P. ponderosa var. brachyptera are negative and significant, increasing in strength when averaging over multiple months (Fig. 2. 4A; 2. 4C). The strongest partial correlations between growth and TAVG for P. arizonica and P. ponderosa var. brachyptera occur for the three - month period ending in July ( - 0 .30, - 0.25, respectively). Although the trend was largely negative, there are no significant negative partial correlations for P. arizonica at BIG (Fig. 2. 4B), and the trend for P. ponderosa var. brachyptera is significant for the three - month period endin g in July ( - 0.25; Fig 2. 4D). 41 PDSI correlation and TAVG partial correlation Correlations between annual growth and single - month PDSI are positive and significant, with maximum values reached in July at both sites for P. arizonica [0.7 (MTL); 0.4 (BIG); Table 2. 4A] and P. ponderosa var. brachyptera [0.5 (MTL); 0.4 (BIG) ; Table 2. BIG. Partial correlations between annual growth and single - month TAVG are nega tive but non - significant for both species at MTL (data not shown), increasing in strength, and becoming significant when averaging over multiple months (Fig. 2. 5A, 2. 5C). There are no significant temperature correlations for either species at BIG (Fig. 2. 5B, 2. 5D). Correlations between annual growth and seasonal PDSI at MTL were of longer duration and stronger for P. arizonica (Table 2. 4A) than P. ponderosa var. brachyptera (Table 2. 4B). Significant positive correlations ( r ) range from r 0.3 to 0.7 for season lengths of 1, 6, 12, and 20 months for P. arizonica , and from r 0.2 to 0.5 for season lengths of 1, 6, and 12 months for P. ponderosa var. brachyptera . In contrast, at BIG there were no observed differences between the two species, the significant correlations only occur for season lengths of 1 month and 6 months, and effect sizes were generally lower ( r 0.3 - 0.4; Table 2. 4 A - B). DISCUSSION Series length, Chronology statistics, and Growth patterns The average age of the sampled trees was c. 74.5 years. This conservative estimate was obtained from the increment cores, not adjusting for pith. Similar aged cohorts (between sites) were sampled to rule out age - related differences. Sympatric 42 areas (wit hin a site) were sampled to rule out disturbance - related differences, such as insect infestation or fire that could influence the ring - width growth for just one species . High mean sensitivity and Gini values indicate variability in annual ring widths kno wn as a sensitive growth response to climate. Variability measured as mean sensitivity (MS) is similar to values reported by Fritts (1974) for P. ponderosa var. brachyptera (average MS = 0.35) and P. arizonica (MS = 0.35) populations sampled near our stud y site at Mt. Bigelow in the early 1960s, which indicate long - term sensitivity to climate and extends the formative work of Fritts. P. ponderosa var. brachyptera expressed the highest common growth signal based on Rbt (Table 2. 3), which suggests the species is responding to strong external climate signals. The dendroclimatic response varies modestly between species and sites, as quantified by the coefficient of variation, which suggests the two closely related ponderosa pine species have different ec ological requirements for growth when sympatric. Autocorrelation is expected to be higher on extreme sites because stressed trees take a year or more to recover following a harsh growing season, thereby conveying persistent physiological effects ( Fritts et al. 1965; Monserud and Marshall 2001). Autocorrelation was significantly higher at Mt. Bigelow than Mt. Lemmon, suggesting the former site is more extreme. Large values of autocorrelation are often explained by the storage of carbohydrates in parenchyma tissue or as the result of temporal autocorrelation found in the precipitation variable (e.g. Matalas 1962; Esper et al. 2015). It is unlikely that annual rainfall is the source of the observed persistence in tree - ring indices because total precipitation between sites is highly correlated. Therefore, the 43 higher value of autocorrelation at Mt. Bigelow is likely caused by biological properties rather than temporal variation in rainfall. One possible explanation for differences in environmental conditions b etween sites affecting autocorrelation is the dissimilarity in elevation between the two mountains. The volume of water released after the melt of snowpack is expected to be greater on Mt. Lemmon than Mt. Bigelow as the former is c. 250 meters higher in el evation with higher snowfall. Also, the AWHC of the soil at Mt. Lemmon (9.2) is more than two times greater than at Mt. Bigelow (3.8). The combined effect of these features is to increase the storage capacity of the winter snowmelt and prolong the moist conditions of the soil into the growing season when the water demand is highest (Barnett et al. 2005). Thus, the growing environment at Mt. Lemmon may be more favorable for tree species adapted to moister growing conditions. The two ponderosa pine speci es experienced recurrent periods of reduced growth. Narrow rings formed during the years 1934, 1936 - 39, and 1947 - 48 correspond to low annual PDSI values, and past reports of reduced streamflow (1932 - 36) and drought (1942 - 64) in Arizona (McNab and Karl 199 1). In both instances, the drought signals of narrow rings were more pronounced at Mt. Bigelow than Mt. Lemmon, a possible result of rockier soils with lower water - holding capacity contributing to lower annual tree - growth at Mt. Bigelow ( Candel - Pérez et al. 2012). More recently, growth was reduced during 2002 - 09 reflecting severe drought impacting the entire Southwest region (Woodhouse et al. 2010). Interestingly, the strongest drought signal for the decade is at Mt. Lemmon. 44 Fire explains the reversal i n site responses observed for the 1930 - 40s drought compared to the 2002 - 09 drought. The change in growth pattern may have resulted from wildfire decreasing the available soil moisture and amount of photosynthetic foliar tissue, thereby reducing annual gro wth at Mt. Lemmon. In fact, throughout the early 21st century, the Santa Catalina Mountains were impacted by several large wildfires. These included the Oracle (2002), Bullock (2002), and Aspen (2003) fires. During the peak growing season of 2003, the A spen wildfire burned the south slope of Mt. Lemmon, including our study site, but left the Mt. Bigelow site relatively undisturbed. The heat of the fire resulted in foliar injury (Telewski, personal observation), which likely decreased cambial growth throu gh a reduction in foliage and photosynthetic capacity. The formation of a hydrophobic soil layer during wildfire may increase runoff. Dyrness (1976) reported that following a large wildfire in the High Cascades of Oregon, USA ( c. 3116 ha), which burned l odgepole pine stands, precipitation runoff increased and recovery (of the soils) was imperceptible until 3 - 5 years after the fire. The phenomenon of increased runoff would occur even during a drought. Therefore, moisture stress induced by fires in 2003 a t Mt. Lemmon would have further reduced observed growth for several more years. Additionally, the PDSI values modeled at Mt. Lemmon would not have captured the effect of increased runoff due to fire because the soil attributes are assumed to be stable. R esponse Function Analysis Response function analysis indicates strong direct positive correlations between precipitation and annual tree growth, which is consistent with previous reports for 45 ponderosa pine species and moisture limitations in the Southwest (Fritts et al. 1965; Brown 1968; Fritts 1974) and Pacific Northwest (Dannenberg and Wise 2016). Our study agrees with summit of Mt. Bigelow, which determined soil moisture is often more lim iting to photosynthesis (hence cambial growth) than low air temperatures for semiarid conifers, even in winter. There are seasonal influences and our results show that cool springtime conditions increased the inter - annual growth for both taxa (April - June) and demonstrate seasonal variation for the lower elevation species ( P. arizonica ). P. arizonica is also sensitive to winter moisture at Mt. Lemmon, which implies growing season use of winter water extracted from deep soil for growth (previous NovembeMa rch). Response function analysis indicates early winter precipitation (previous October - previous December) was influencing the growth of P. arizonica as opposed to late winter (January - March), which suggests that snow melt was not a major factor in the sp ecies response to winter climate . Although responses to temperature were weaker than to precipitation, Fig. 2. 3 shows a growth reduction occurred for P. arizonica when minimum spring temperature (April - June) rose above a critical level at Mt. Lemmon. Th e observed reduction in correlations between tree rings and minimum temperature has been reported in other southwestern ponderosa pine studies and is explained by increased respiration or drought stress (Fritts 1974, 1976; Adams and Kolb, 2005), which nega tively affected the growth of P. arizonica but not P. ponderosa var. brachyptera . 46 Correlations to Climate Analysis of seasonal response functions of ponderosa pine indicated that high precipitation during the cool - moist conditions of winter (November - March), and spring (April - June) are most important for tree growth. Our results broadly agree with González - Cásares (2017) who reported maximum positive correlation for the 9 - month period ending in June (October - June) for the annual grow th of P. arizonica in northwestern Mexico. Further analysis of the climate - growth relationships with correlation analyses indicated that the annual growth of the ponderosa pine species was under significant moisture stress in spring. The strong positive c orrelation to spring precipitation and negative partial correlation to spring temperature indicates that the semiarid local climate of the Santa Catalina Mountains has negatively impacted the growth of P. ponderosa var. brachyptera growing at its lower ele vational limit . In comparison, the growth of P. arizonica growing at its upper elevational limit was only impacted at Mt. Lemmon. A consistent temperature signal is lacking for the tree - ring indices at Mt. Bigelow, demonstrating differences in habitat su itability between sites for ponderosa pine growth. These results correspond to the conceptual models of Fritts (1974, 1976) that describe the relationships leading to moisture stress prior to or during the growing season, resulting in the formation of narr ow rings. High temperature and low precipitation are primarily responsible for decreased growth of semiarid site conifers. With these models in mind, we considered drought conditions in spring, which can lead 47 to high evapotranspirational demand and the f ormation of narrow rings during the growing season. Sensitivity to PDSI P. arizonica and P. ponderosa var. brachyptera are drought - sensitive species (Peltier et al. 2016), which show strong positive PDSI - growth correlations, although the climate - growth relationships vary between species at Mt. Lemmon. The PDSI - growth correlations did not support our hypothesis that P. po nderosa var . brachyptera would be more sensitive to soil moisture availability. P. arizonica has a stronger long - term relationship, while P. ponderosa var correlations were for a shorter duration and weaker. There are differences between s ites; correlation to drought lasts longer at Mt. Lemmon than Mt. Bigelow, regardless of species. A multi - species synthesis of ring - width variance has shown that P. ponderosa has low resistance and slow recovery to drought (Peltier et al. 2016). Similar ly, the observed extended correlation of growth with PDSI in this study implies that P. arizonica and P. ponderosa var. brachyptera may be vulnerable to changing climate (Gonzalez - Casares 2017). Thus, warmer conditions may have a more negative effect on P . arizonica, as periods with low PDSI become longer or more common. Other studies have noted differences in growth sensitivity to climate for the Pinus species of Mexican pine forests (Bickford et al. 2011; Gonzalez - Casares et al. 2017). These dissimilar ities have been partially credited to different drought tolerances of the species under study. Although more work is required to build the appropriate models, the length of correlation to PDSI could be one factor useful in determining tolerance to moistur e stress. 48 It is interesting that tree growth correlated to PDSI for the longest period at Mt. Lemmon, the more favorable site . Although not as extreme as the crevice sites in the Western USA described by Fritts (1976), the ponderosa pine species of Mt. Bigelow also grow in an extreme environment with shallow, rocky soils, and restricted root describes dry soils as constraining root and crown growth for longer periods of time, and trees acclimate to these conditions so that growth responses to precipitation will be less vigorous than trees growing on more favorable sites. Thus, trees growing on more favorable sites will suffer more stress from drought than trees growi ng on drier sites, which have acclimated to low water availability. Therefore, one reason for the longer drought influence at Mt. Lemmon is the greater loss of green needle tissue related to prolonged water stress. This hypothesis of reduction in photosy nthetic capacity is supported by Galiano et al. (2011) who determined that carbon reserves are key to tree recovery and forest resilience following periods of drought; consequently, respiration will influence carbon reserves and tolerance to moisture stres s. Respiration Temperatures on south - facing slopes may exceed the optimum range for plant processes during the growing season or winter months (Fritts 1974). Therefore, slope aspect greatly affects the correlation of ring widths with climate, and annual ring widths are expected to correlate negatively with variations in monthly temperature. The conceptual models of Fritts (1974, 1976) describe the pathways leading to increased respiration, and the formation of narrow rings. These negative temperature effects on 49 gro wth are apparent during the growing season due to high temperature and respiration redu cing the plant photosynthetic processes (Fig. 1 of Fritts 1974). Negative growth effects are also expected from elevated temperatures during the winter months (i.e. p rior to the growing season) due to the direct effects of warm air temperature on respiration resulting in increased carbon consumption and decreased carbon allocation; consequently, cold hardiness is reduced (Fig. 2 of Fritts 1974; Ögren et al. Sundbald 19 winter respiration. The correlation analysis indicated growth of P. arizonica was limited by spring and summer drought stress (+, PDSI) and temperature ( - , TAVG) at Mt. Lemmon, an in dication of increased respiration during the growing season. In comparison, P. ponderosa var. brachyptera was limited by moisture stress (+, PDSI) and temperature ( - different ly to water balance requirements and carbon limitations. The most common cause of increased evapotranspirational demands is high air temperature at low elevations, which may cause increased moisture stress in spring and summer ( P. arizonica ) and decreased cold tolerance in winter ( P. ponderosa var. brachyptera ) . CONCLUSIONS The growth of the two species responded positively to the cool - wet conditions of spring, and winter precipitation also improved the growing conditions for P. arizonica. We did not fin d support for higher climate sensitivity in P. ponderosa var. brachyptera (at its lower elevational limit) above that observed in P. arizonica (at its upper elevational limit). However, we report small differences in ring - width variances between 50 P. ponderosa var. brachyptera and P. arizonica , and that both species are sensitive to the negative effects of respiration on annual growth in different seasons, winter and summer, respectively. Notably, rising night - time temperatures during springtime re duced growth for P. arizonica but not P. ponderosa var . brachyptera at Mt. Lemmon, most likely from increased respiration or moisture stress. Tree rings correlated to PDSI at different temporal scales, which suggests that P. arizonica is sensitive to drou ght for longer periods. Although our data did not support our hypothesis that P. ponderosa var. brachyptera would be more sensitive to soil moisture availability, the temporal differences in drought sensitivities distinguish the two taxa by their varying tolerances to water shortages. Both species exhibit a stronger climatic response at Mt. Lemmon (the less dry site) than Mt. Bigelow, and moisture availability between sites can be further influenced by disturbance. In summary, the growth of P. ponderosa var. brachyptera and P. arizonica at their transition zone demonstrated subtle, biologically meaningful differences in seasonal precipitation requirements for growth, sensitivity to moisture stress and tolerance to cold temperatures. Data availability : The tree - ring datasets generated and analyzed during the current study are available from the International Tree - Ring Data Bank (ITRDB) repository, https://www. ncdc.noaa.gov/data - access/paleoclimatology - data/datasets/tree - ring . These datasets include the raw ring widths for 90 trees (204 radii) and 4 composite chron Methods . Electronic supplementary material: The online v ersion of this article ( https ://doi.org/10.1007/s0046 8 - 018 - 1778 - 9 ) contains supplementary material, which is available to authorized users. The supplementary materials are also included in 51 Appendix D, which comprise Suppl. Fig. 2.S 1 : Figures of soil AWHC, Suppl. Fig. 2.S 2 : Correlations of modeled and locally collected climate data, Suppl. Fig. 2.S 3 : Plots of annual PDSI values, Suppl. Fig. 2.S 4 : Climate graphs, Suppl. Fig. 2.S 5 : Correlation coefficients of seasonal climate variables with tree - ring chronologies, Suppl. Fig. 2.S 6 : Scatterplots of mean tree - ring index against the average primary precipitation variable. Author contributions: P.M. conceived the study; P.M. and F.T. designed the experiments and collected samples in the fie ld; P.M. led and S.J. and G.T. assisted with preparing the samples; P.M. led and B.M assisted with data analysis; P.M. wrote the paper; and B.M, F.T., S.J., and G.T. edited the paper. Conflicts of interest: The authors declare no conflicts of interest. 52 APPENDICES 53 APPENDIX C Validation of C limatic D ata 54 Validation of C limatic D ata Soils data used in the West Wide Drought Tracker PDSI calculation were obtained from Penn State University (Abatzoglou et al. 2017). The applicability of the 4 - km resolution PDSI data was assessed by comparing the AWHC values of the specific soil units in the calculation (Abatzoglou et al. 2017) to those represented in the local soils plots near the research sites (n= 18 0; 60 trees x 3 sites; Top panel; Suppl. Fig. 2.S 1 ). The soils data used in the PDSI calculation indicate AWHC of 6.0% for the area encompassing the study sites. The modeled PDSI data were validated by interpolating AWHC values from soil pedons sampled from the Coronado National Forest (n =14; Lower Panel; Suppl. Fig. 2.S 1 ) and averaging the interpolated values for sampled tree locations (n =120; 60 trees x 2 sites; AWHC = 5.9%). Precipitation and TAVG data were validated by correlating the mod elled site data with locally collected weather data that were accurate but too short for tree - (1965 to 1981) collected at 2425 m a.s.l. within 0.7 km of the study sites ( Vose et al. 2014) , were used to valida te the mo nsoonal patterns ( Suppl. Fig. 2.S2 years of climate data (1960 to 2015; Vose et al. 2014 ) within c . 95 km and similar in elevation (2070m) to the study sites, were used to validate winter PCP and TAVG ( Suppl. Fig. 2.S 2B - C ). Pear son correlations were applied to validate the gridded climate dataset (McKenney et al. 2011). Correlations were consistently highest for the McKenney / NOAA PAL associations. Same month - month correlations are positive and range from 0.8 (June - June) for the TAVG association ( Suppl. Fig. 2.S 2C ), to 1.0 (e.g. Feb. - Feb.) for the PCP associations ( Suppl. Fig. 2.S2 ) . 55 Validation of C limatic D ata (continued) The strength of the collinear correlation is weakest for June - June PCP for the McKenney Mt. Lemmon / Kitt Peak dataset (0.52; Suppl. Fig. 2.S2 ) because summer 56 APPENDIX D Supplementary F igures 57 Supplementary f igures Figure 2.S1 Figures of soil AWHC (%) used to validate PDSI data analyzed in study. 58 Figure 2.S 1 (continued) The upper map depicts an AWHC of 6.0% (the modeled data set). The lower map was generated by values interpolated from 14 soil samples taken from three locations within the Coronado National Forest: MTL (Blue), Palisades (PAL; Green), and BIG (Yellow). The average AWHC 5.9% among study plots (n = 180) was adjusted for depth to bedrock and for rock content (AWHC for MTL, PAL and BIG are 9.2, 4.7, and 3.8% respectively) . The lower map used a conventional weighted average of the soil horizon available water c apacity (AWC) 59 Figure 2.S2 C orrelations of modeled and locally collected climate data are collinear and validate the accuracy of the gridded climate dataset (i.e. same month correlations with high positive r - value, e.g. May - May). The McKenney et al. (2011) modelled variables are on the y - axis and corresponding local NOAA variables are on the x - axis . (A)Total PCP: modeled PAL (y) to local PAL(x) 60 Figure 2.S2 (continued) ( B) Total PCP: modeled MTL (y) to local Kitt Peak (x) 61 Figure 2.S2 (continued) ( C) TAVG: modeled MTL (y) to local Kitt Peak (x) 62 Figure 2.S 3 Plots of average annual PDSI values estimated for two study plots MTL (top) and BIG (bottom), for the reference period 1925 - 2009. Yearly values are depicted along the x - axis and average PDSI index along the y - axis. Negative values denote dry conditions 63 Figure 2.S 4 Climate graphs of mean monthly climate data analyzed for BIG (Mt. Bigelow) research site ( - 110.71495, 32.41378; 2535 m elevation) for the reference period of 1925 - 2009 (McKenney et al., 2011). 64 Figure 2.S4 (continued) Climate variables are total precipitation PC P (A) and mean average temperature TAVG (B). The box shows location of the middle quartile of monthly observations, the median is marked by the horizontal line in the center of the box, and the whiskers extend to the most extreme values not considered outl iers. An outlier is marked by a closed circle and is defined as any data value more or less than 1.5 d beyond the edges of the box, where d is the interquartile range 65 Figure 2.S 5 Correlation analysis for P. arizonica (white) , P. ponderosa var. brachyptera (grey) growing at two transition zones within the Santa Catalina Mountains: MTL (hatched bars) and BIG (non - hatched bars). 66 Figure 2.S5 (continued) Analyzed were climate variables of precipitation (PCP) and minimum temperature (TMIN) with seasons: Previous (Pr) Summer (Jul - Sep), Winter (Nov - Mar), Spring (Apr - Jun), and Summer (Jul - Sep) (A); and climate variables of PCP and maximum temperature (TMAX with 3 - month seasons: Pr Summer, Early Winter (Oct - Dec), Late Winter (Jan - Mar), Spring, and Su mmer (B), which were related to the composite chronologies for the period of 1925 - 2009. The standard index was used in the growth analysis of P. arizonica on MTL, and autoregressive modeling was applied to all other tree - ring time series to remove autoc orrelation from past - 67 Figure 2.S6 Scatterplots of the mean tree - ring index (ARSTAN or Standard chronology) against the average primary PCP variable (summing across months) for the most highly correlated periods in the Seascorr analyses ( Meko et al. 2011). These periods are defined as the three season - lengths with the highest correlation. (A) Plot of Mt. Lemmon P. arizonica tree - ring index ( y) by PCP (x); all p - values are highly significant 68 Figure 2.S6 (continued) (B) Plot of Mt. Lemmon P. ponderosa var. brachyptera tree - ring index (y) by PCP (x) ; all p - values are highly significant (C) Plot of Mt. Bigelow P. arizonica tr ee - ring index (y) by PCP (x) ; p - values are sig. (D) Plot of Mt. Bigelow P. ponderosa var. brachyptera tree - ring index (y) by PCP (x), all p - values are highly significant 69 APPENDI X E Tables and F igures 70 T able 2. 1 Sampling design used for two south - facing slopes, Santa Catalina range (sites MTL and BIG) PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PPB / PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ PAZ Elevation increases from the bottom to top of schematic. Each box represents 1 sample point on a 5 x 9 grid. One tree was sampled c. every 10 m except for the center of the transition zone, where 1 P PB and 1 PAZ each were sampled N = 30 trees collected pe r species PAZ P. arizonica , PPB P. ponderosa var. brachyptera 71 Table 2 . 2 Summaries for climate data obtained from gridded data sets (McKenney et al. 2011) for two study sites: MTL and BIG MTL BIG Geographic Coordinates, WGS84 32.43, - 110.79 32.41, - 110.71 (Latitude, Longitude) Elevation (m) 2577 2534 Years spanned 1925 - 2009 1925 - 2009 Years of record 84 84 Yearly TAVG (mean 0 C + SD) 11.0 + 0.6 10.8 + 0.5 January TAVG (mean 0 C + SD) 2.8 + 1.7 2.6 + 1.8 July TAVG (mean 0 C + SD) 19.5 + 1.0 19.4 + 0.8 Yearly PCP (mean mm + SD) 678.8 + 164.1 685.0 + 163.8 Winter PCP (mean mm + SD) 298.7 + 154.5 299.7 + 155.4 Arid spring PCP (mean mm + SD) 47.8 + 29.2 48.8 + 29.7 Summer PCP (mean mm + SD) 296.9 + 72.9 301.1 + 73.7 Temperature year considers coldest and hottest months: January and July. Precipitation year is divided into three seasons: winter (previous November current March), arid spring (April June), and summer (July September) 72 Table 2. 3 General statistics for the common interval (1925 - 2009) of the de - trended chronologies for P. arizonica and P. ponderosa var. brachyptera from transition zones at two sites (MTL and BIG) MTL BIG P. arizonica brachyptera P. arizonica brachyptera No. trees (No. radii) 25 (58) 24 (64) 24 (40) 20 (40) Ave No. + yrs. SD 67 + 10 63 + 12 63 + 20 59 + 22 RWI 0.96 a 1.00 b 0.96 0.97 SD 0.38 0.36 0.40 a 0.46 b MS 0.35 0.34 0.36 0.36 Gini 0.21 0.20 0.22 0.25 Rbt Rbt 2 0.35 0.12 0.39 0.15 0.33 0.11 0.40 0.16 EPS 0.94 0.94 0.90 0.92 SNR 14.93 17.09 9.44 12.08 Autoregressive model 0 2 2 2 Variance due to Autoregression n.a. 16.4 23.4 28.7 AC1 a 0.13 0.32 0.38 0.45 Statistically significant differences are indicated by lowercase letters (p<0.05) RWI ring width index; SD standard deviation; MS mean sensitivity; Gini Gini coefficient; Rbt between tree correlations; EPS expressed population signal; SNR signal - to - noise ratio; AC1 first - order autocorrelation a Standard Chronology 73 Table 2. 4 Correlation coefficients of seasonal P DSI with tree - ring chronologies (4A) PAZ SIT E MTH LEN COE F SIG LEN COE F SIG LEN COE F SIG LEN COE F SIG MTL J 1 0.44 * 6 0.33 * 12 0.29 * 20 0.26 MTL F 1 0.45 * 6 0.38 * 12 0.31 * 20 0.29 * MTL M 1 0.49 * 6 0.44 * 12 0.34 * 20 0.31 * MTL A 1 0.53 * 6 0.48 * 12 0.37 * 20 0.33 * MTL M 1 0.59 * 6 0.52 * 12 0.41 * 20 0.36 * MTL J 1 0.61 * 6 0.55 * 12 0.45 * 20 0.38 * MTL J 1 0.67 * 6 0.58 * 12 0.49 * 20 0.42 * MTL A 1 0.58 * 6 0.61 * 12 0.53 * 20 0.44 * MTL S 1 0.53 * 6 0.62 * 12 0.56 * 20 0.46 * BIG M 1 0.27 * 6 0.21 12 0.15 BIG J 1 0.29 * 6 0.23 12 0.17 BIG J 1 0.36 * 6 0.26 * 12 0.19 BIG A 1 0.31 * 6 0.29 * 12 0.22 BIG S 1 0.29 * 6 0.31 * 12 0.24 74 Table 2. 4 ( continued ) (4B) PPB Site M LEN COEF SIG LEN COEF SIG LEN COEF SIG LEN COEF SIG MTL J 1 0.23 * 6 0.06 12 - 0.01 20 - 0.06 MTL F 1 0.24 * 6 0.13 12 0.02 20 - 0.03 MTL M 1 0.28 * 6 0.20 12 0.05 20 - 0.01 MTL A 1 0.33 * 6 0.25 * 12 0.09 20 0.02 MTL M 1 0.40 * 6 0.30 * 12 0.14 20 0.05 MTL J 1 0.43 * 6 0.33 * 12 0.19 20 0.08 MTL J 1 0.48 * 6 0.38 * 12 0.24 * 20 0.12 MTL A 1 0.40 * 6 0.40 * 12 0.29 * 20 0.15 MTL S 1 0.35 * 6 0.42 * 12 0.32 * 20 0.17 BIG M 1 0.31 * 6 0.25 12 0.22 BIG J 1 0.32 * 6 0.27 12 0.24 BIG J 1 0.39 * 6 0.30 * 12 0.26 BIG A 1 0.35 * 6 0.32 * 12 0.28 BIG S 1 0.33 * 6 0.34 * 12 0.30 (A) PAZ P. arizonica (B) PPB P. ponderosa var. brachyptera, MTH month, LEN length of season in months. One - , 6 - , and 12 - month seasons were analyzed at both sites; 20 - month increment is also analyzed at MTL, significance level of 0.05. Partial correlation coefficients for the secondary climate variable (TAVG) are not shown 75 Table 2. 5 Effect sizes ( r ) of PDSI as a measure of drought sensitivity ( Strongest Month Strength Site Taxa Score highest r r effect MTL P. arizonica Single month July 0.67 Large MTL P. pond. var. brachyptera Single month July 0.45 Medium BIG P. arizonica Single month July 0.36 Medium BIG P. pond. var. brachyptera Single month July 0.39 Medium All p - values are highly significant (p < 0.001) 76 Inset shows the study area located in southeast AZ, USA (red circle) Fig. 2.1 Locations of two study sites sampled for tree - ring analysis: MTL (Mt. Lemmon) and BIG (Mt. Bigelow). Black triangles mark the study sites northeast of Tucson, AZ USA 77 Fig. 2. 2 ARSTAN or Standard chronologies of P. arizonica (A,C) and P. ponderosa var. brachyptera (B,D) when Expressed Population Signal (EPS) is greater than 0.85, and sample depth is > 14, for the reference period of 1925 2009. Two sites were evaluated: MTL (A,B) and BIG (C,D). The corresponding sample depth (number of tree cores) is indicated by shading 78 Figure 2.3 Climate - growth relationships determined by principal components multiple regression for two species of ponderosa pine [ P. arizonica = white , P. ponderosa var. brachyptera = grey] growing at two transition zones within the Santa Catalina Mountains: MTL (hatched bars) and BIG (non - hatched bars). 79 Figure 2.3 (continued) Climate variables of precipitation (PCP) and minimum temperature (TMIN) w ith seasons: Previous (Pr) Summer (July - S e ptember), Winter (N ovember - March), Spring (April - June), and Summer (July - September ) (A); and climate variables of PCP and maximum temperature (TMAX with 3 - month seasons: Pr Summer, Early Winter (October - December), Late Winter (January - March), Spring, and Summer (B), which were related to the composite chro nologies for the period of 1925 - 2009. The standard index was used in the growth analysis of P. arizonica on MTL, and autoregressive modeling was applied to all other tree - ring time series to remove autocorrelation from past - 80 Figure 2. 4 Climate - growth associations obtained by relating total precipitation (PCP) as primary climate variable (t op in each panel) and mean average temperature (TAVG) as secondary climate variable (bottom in each panel) to the ARSTAN ring - width chronologies of P. arizonica (A,B) and P. ponderosa var. brachyptera (C,D) at two sites MTL (A,C) and BIG (B,D) [Note: Stand ard chronology for (A)]. PCP was summed and TAVG was averaged for a period of three months whose ending months are shown on the x - axes (previous August to current September). Significant correlations and partial correlations (p < 0.05) are shown by dark bars 81 Fig ure 2. 5 Climate - growth associations obtained by relating average PDSI as primary climate variable (top in each panel) and mean average temperature (TAVG) as secondary climate variable (bottom in each panel) to the ARSTA N ring - width chronologies of P. arizonica (A,B) and P. ponderosa var. brachyptera (C,D) at two sites MTL (A,C) and BIG (B,D) [Note: Standard chronology for (A)]. PDSI was summed and TAVG was averaged for a period of three month whose ending months are sho wn on the x - axes (previous August to current September). Significant correlations and partial correlations (p < 0.05) are shown by dark bars 82 REFERENCES 83 REFERENCES Abatzoglou JT, McEvoy DJ, Redmond KT (2017) The west wide drought tracker: drought monitoring at fine spatial scales. Bull Am Meteorol Soc. https://doi.org/10.1175/BAMS - D - 16 - 0193.1 Adams HD, Kolb TE (2005) Tree growth response to drought and temperature in a mountain landscape in northern Arizona, USA. J Biogeogr 32:1629 1640. https://doi.org/10.1111/j.1365 - 2699.2005.01292.x Barichivich J, Briffa K, Myneni R, Schrier G, Dorigo W, T ucker C, Osborn T, Melvin T (2014) Temperature and snow - mediated moisture controls of summer photosynthetic activity in northern terrestrial ecosystems between 1982 and 2011. Remote Sens 6:1390 1431. https://doi.org/10.3390/rs6021390 Barnett TP, Adam JC, Lettenmaier DP (2005) Potential impacts of a warming climate on water availability in snow - dominated regions. Nature 438:303 309. https://doi.org/10.1038/nature0414.1 Bezy J (2016) A guide to the geology of the Santa Catalina Mountains, Arizona: the geology and life zones of a Madrean s ky i sland. Arizona Geological Survey, Tucson. ISBN:978 - 0 - 9854798 - 2 - 4 Bickford IN, Fulé PZ, Kolb TE (2011) Growth sensitivity to drought of co - occurring Pinus spp. along an elevation gradient in northern Mexico. West N Am Nat 71:338 348 Biondi F, Qeadan F (2008) Inequality in paleorecords. Ecology 89:1056 1067 Briffa K, Jones P (1990) Basic chronology statistics and assessment. In: Cook ER, Kairiukstis LA (eds) Methods of dendrochronology. Springer, Dordrecht, pp 137 152. ISBN:978 - 90 - 481 - 4060 - 2 Brown JM (1968) The photosynthetic regime of some southern Arizona ponderosa pine. Thesis ( PhD ) , University of Arizona, Arizona 84 Bunn AG (2008) A dendrochronology program library in R (dplR). Dendrochronologia 26:115 124. https://doi.org/10.1016/j.dendro.2008.01.002 Bunn AG (2010) Statistical and visual crossdating in R using the dplR library. Dendrochronologia 28:251 258. https://doi.org/10.1016/j.dendro.2009.12.001 Bunn AG, Jansma E, Korpela M, Westfall RD, Baldwin J (2013) Using simulations and D endrochronologia 31:250 254 Buol SW (1966) Soils of Arizona. Univ. of Ariz. Ag ric. Exp. Sta. Tech. Bull., Ari zona, p 171 Candel - Pérez D, Linares JC, Viñegla B, Lucas - Borja ME (2012) Assessing climate growth relationships under contrasting stands of co - occurring Iberian pines along an altitudinal gradient. For Ecol Manage 274:48 57. https://doi.org/10.1016/j.foreco.2012.02.010 Cook ER (1985) A time series analysis approach to tree ring standardization. Thesis ( PhD ) , University of Arizona, Arizona Cybis Elektronik (2010) CDendro & CooRecorder program package for dendrochronology ring width and blue channel measurements and for crossdating. http://www.cybis.se/forfun/dendro. Accessed 7 Nov 2013 Dannenberg MP, Wise EK (2016) Seasonal climate signals from multiple tree ring metrics: a case study of Pinus ponderosa in the upper Columbia River Basin. J Geophys Res Biogeosci 121:1178 1189. https://doi.org/10.1002/2015JG003155 DeRose RJ, Bekker MF, Wang SY, Buckley BM, Kjelgren RK, Bardsley T, Rittenour TM, Allen EB (2015) A millennium - length reconstruction of Bear River stream flow, Utah. J Hydrol 529:524 534. Plants, Soils, and Climate Faculty Publications. Paper 741. https://digitalcommons.usu.edu/psc_facpub/741. Accessed 11 June 2018 Dodge RA (1963) Investigations into the ecological relationships of ponderosa pine in southeast Ariz ona. Thesis ( PhD ) , University of Arizona, Arizona 85 Dyrness C (1976) Effect of wildfire on soil wettability in the high Cascades of Oregon; USFS Res. Pap. PNW - 202, Portland Epperson BK, Telewski FW, Plovanich - Jones AE, Grimes JE (2001) Clinal differentia tion and putative hybridization in a contact zone of Pinus ponderosa and P. arizonica (Pinaceae). Am J Bot 88:1052 1057. https://doi.org/10.2307/2657087 Epperson BK, Telewski FW, Willyard A (20 09) Chloroplast diversity in a putative hybrid swarm of Ponderosae (Pinaceae). Am J Bot 96:707 712. https://doi.org/10.3732/ajb.0800005 Esper J, Schneider L, Smerdon JE, Schöne BR, Büntgen U (2015) Signals and memory in tree - ring width and density data. Dendrochronologia 35:62 70. https://doi.org/10.1016/j.dendro.2015.07.001 Farjon A, Brian T (1997) Styles. Flora neotropica monograph. Pinus (Pinaceae). New York Botanical Garden, New York City, p 75. ISBN:0893274119, 9780893274115$4 Fritts HC (1974) Relationships of ring widths in arid - site conifers to variations in monthly temperature and precipitation. Ecol Mon ogr 44:411 440. https://doi.org/10.2307/1942448 Fritts HC (1976) Tree rings and climate. Academic, London. ISBN:0 12 268450 - 8 Fritts HC, Smith DG, Budelsky CA, Cardis JW (1965) The variabili ty of ring characteristics within trees as shown by a reanalysis of four ponderosa pine. Tree Ring Bull 27:3 18 Galiano L, Martínez - Vilalta J, Lloret F (2011) Carbon reserves and canopy defoliation determine the recovery of Scots pine 4 year after a dro ught episode. New Phytol 190:750 759. https://doi.org/10.1111/j.1469 - 8137.2010.03628.x González - Cásares M, Pompa - García M, Camarero JJ (2017) Differences in climate growth rela tionship indicate diverse drought tolerances among five pine species 86 coexiting in northwestern Mexico. TREES 31:531 544. https://https//doi.org/10.1007 /s00468 - 016 - 1488 - 0 Griffin D, Woodhouse CA, Meko DM, Stahle DW, Faulstich HL, Carrillo C, Touchan R, Castro CL, Leavitt SW (2013) North American monsoon precipitation reconstructed from tree - ring latewood. Geophys Res Lett 40(5):954 958. https:// doi.org/10.1002/grl.50184 Guiot J (1991) The bootstrapped response function. Tree Ring Bull 51:39 41 Haller JR (1965) The role of 2 - needle fascicles in the adaptation and evolution of ponderosa pine. Brittonia 17:354. https://doi.org/10.2307/2805029 Holmes RL, Adams RK, Fritts HC (1986) Tree - ring chronologies of western North America, California, eastern Oregon and northern Great Basin with procedures used in the chronology development work including user manuals for computer programs COFECHA and ARSTAN. Chronol Ser VI:182. http://hdl.handle.net/10150/304672. Accessed 12 May 2016 Hutchinson M, Xu T (2013) ANUSPLIN version 4.4 user guide. Australia National University, Canberra LaMarche VC, Stockton CW (1974) Chronologies from temperature - sensitive bristlecone pines at upper treeline in western United States. Tree Ring Bull 34:21 45 Matalas NC (1962) Statistical properties of tree ring data. Hydrol Sci J 7:39 47 McC ullough IM, Frank W, Davis P, Williams A (2017) A range of possibilities: assessing geographic variation in climate sensitivity of ponderosa pine using tree rings. For Ecol Manage 402:223 233. https://doi.org/10.1016/J.FORECO.2017.07.025 McKenney DW, Hutchinson MF, Papadopol P, Lawrence K, Pedlar J, Campbell K, Milewska E, Hopkinson RF, Price D, Owen T (2011) Custo mized spatial climate models for North America. Bull Am Meteorol Soc 92:1611 1622. https://doi.org/10.1175/2011BAMS3132.1 87 McNab A, Karl T (1991) Climate and droughts. In: Paulson RW, Cha se EB, Roberts RS, Moody DW (eds) N ational water summary 1988 89 H ydrologic events and floods and droughts. US Geological Survey, Reston, pp 89 95 Meko DM, Touchan R, Anchukaitis KJ, Seascorr (2011) A MATLAB program for identifying the seasonal climate signal in an annual tree - ring time series. Comput Geosci 37:1234 1241. https://doi.org/10.1016/j.cageo.2011.01.013 Monserud RA, Marshall JD (2001) Time - tree rings. I. Time trends and autocorrelation. Tree Physiol 21:1087 1102 Norris JR, Jackson ST, Betancourt JL (2006) Classification tree and minimum - volume ellipsoid analyses of the distribution of ponderosa pine in the western USA. J Biogeogr 33:342 3 60. https://doi.org/10.1111/j.1365 - 2699.2005.01396.x Ögren E, Nilsson T, Sunblad LG (1997) Relationship between respiratory depletion of sugars and loss of cold hardiness in coniferous seedlings over - wintering at raised temperatures: indications of different sensitivities of spruce and pine. Plant Cell Environ 20:247 253. https://doi.org/10.1046/j.1 365 - 3040.1997.d01 - 56.x Palmer WC (1965) Meteorological drought. U.S. Weather Bureau Research Paper, 30 Peloquin RL (1984) The identification of three - species hybrids in the ponderosa pine complex. Southw Nat 29:115. https://doi.org/10.2307/3670776 Peltier DMP, Fell M, Ogle K (2016) Legacy effects of drought in the southwestern United States: a multi - species synthesis. Ecol Monogr 86:312 326. https://doi.org/10.1002/ecm.1219 Perry JP (1991) The pines of Mexico and Central America. Timber, Portland. ISBN:10 1604691107 Price R, Liston A, Strauss S (1998) Phylogeny and systematics of Pinus . In: Richardson D (ed) Ecology and bioge ography of Pinus . Cambridge University Press, Cambridge, p 527, ISBN:0 521 55176 5 88 R Core Team (2016) R: A language and environment for statistical computing. R foundation for Statistical Computing, Vienna, Austria. https://www.R - project.org/ . Accessed 1 Feb 2017 Rehfeldt GE (1993) Genetic variation in the Ponderosae of the southwest. Am J Bot 80:330 343. https://doi.org/10.2307/2445357 Rehfeldt GE (1999) Systematics and genetic structure of Ponderosae taxa (Pinaceae) inhabiting the Mountain Islands of the Southwest. Am J Bot 86:741. https://doi.org/10.2307/2656584 Sheppard PR, Comrie AC, Packin GD, Angersbach K, Hughes MK (2002) The climate of the US Southwest. Clim Res 21:219 238. https://doi.org/10.3354/cr021219 Shinneman DJ, Means RE, Potter KM, Hipkins VD (2 016) Exploring climate niches of ponderosa pine ( Pinus ponderosa Douglas ex Lawson) haplotypes in the western United States: I mplications for evolutionary history and conservation. PLoS One 11:e0151811. https://doi.org/10.1371/journal.pone.0151811 Shreve F (1915) The vegetation of a desert mountain range as conditioned by climatic factors. Carnegie Institution of Washington, Washington DC. https://doi.org/10.5962/bhl.title.45648 Shreve F (1917) The density of stand and rate of growth of Arizona yellow pine as influenced by climatic conditions. J For 15(6):695 707 Shreve F (1919) A comparison of the vegetational features of two desert mountain ranges. Plant World 22:291 307 Speer J (2010) Fundamentals of tree - ring research. University of Arizona Press, Tucson. ISBN:978 - 0 - 8165 - 2684 - 0. https: //doi.org/10.3959/1536 - 1098 - 70.2.161 Stokes M, Smiley T (1996) An introduction to tree - ring dating. University of Arizona Press, Tucson. ISBN:978 - 0 - 8165 - 1680 - 3 89 USDA NRCS, The PLANTS, Database (2017) National plant data team, Greensboro, NC 27401 - 4901 USA. http://plants.usda.gov . Accessed 8 Oct 2017 Vose RS, Applequist S, Squires M, Durre I, Menne MJ, Williams C N, Fenimore C, (CLIMDIV7). NOAA Nat Clim Data Cen. https://doi.org/10.7289/V5M32STR Wei T, Simko VR (2017) package atrix (Version 0.84). https://github.com/taiyun/corrplot. Accessed 1 May 2017 Whittaker RH, Niering WA (1965) Vegetation of the San ta Catalina Mountains, Arizona: a gradient analysis of the south slope. Ecology 46:429 452. https://doi.org/10.2307/1934875 Whittaker RH, Niering WA (1975) Vegetation of the Santa Catalina Moun tains, Arizona. V. Biomass, production, and diversity along the elevation gradient. Ecology 56:771 790. https://doi.org/10.2307/1936291 Whittaker RH, Buol SW, Niering WA, Havens YH (1968) A soi l and vegetation pattern in the Santa Catalina Mountains, Arizona. Soil Sci 105(6):440 450 Wigley TML; Briffa, K.R.; Jones, P.D. (1984) On the average value of correlated time series, with applications in dendroclimatology and hydrometeorology. J Clim Appl Meteorol 23:201 213. https://doi.org/10.1175/1520 - 0450(1984)023%3C0201:OTAVOC%3E2.0.CO;2 Willyard A, Gernandt DS, Potter K, Hipkins V, Marquardt P, Mahalovich MF, Langer SK, Telewski FW, Cooper B, Douglas C, Finch K (2017) Pinus ponderosa : a checkered past obscured four species. Am J Bot 104:161 181. https://doi.org/10.3732/ajb.1600336 Woodhouse CA, Meko DM, MacDonald GM, Stahle DW, Cook ER (2010) A 1,200 - year perspective of 21st century drought in southwestern North America. Proc Natl Acad Sci 107:21283 21288. https:/ /doi.org/10.1073/pnas.0911197107 Zang C, Biondi F, Treeclim (2015) An R package for the numerical calibration of proxy climate relationships. Ecography 38:431 436. https://doi.org/10.1111/ec og.01335 90 Zar JH (2010) Biostatistical analysis: Pearson new international edition. Pearson Higher Ed, London 91 C HAPTER 3 SHIFTS IN CLIMATE - GROWTH RELATIONSHIPS OF THE SKY ISLAND PONDERSOA PINES Paula E. Marquardt 1 * 1 USDA Forest Service, Nort hern Research Station, 5985 Highway K, Rhinelander, WI 54501 USA *Corresponding author: pmarquardt@fs.fed.us ; 715 - 362 - 1121; ORCID ID 0000 - 000 - 8854 - 0172 Acknowledgments : This paper is part of a dissertation submitted to Michigan State University in partial fulfillment of requirements for a Doctor of Philosophy degree. The USDA Forest Service, and the Dept. o f Plant Biology, Michigan State University provided financial support for the project. We thank the following people for help on the F. Telewski helped with experimental design . B. Miranda provided processing support. J. Stanovick provided statistical support. H. Jenson provided GIS support. D. Donner reviewed the paper. We also thank the USDA Coronado National Forest for providing permission to sample trees in the Santa Catalina Mountains. 92 ABSTRACT Abstract: Steep climate gradients permit the coexistence of species filling different ecological niches. We compared the sensitivity and temporal stability of climate - growth relationships across two contact zones in southeas tern, Arizona, USA to quantify the growth responses of sympatric morphotypes (3 - , mixed - , 5 - needle) of ponderosa pine to changing climate. Positive climate - growth correlations in semiarid forests indicated a seasonal shift from summer - to spring - dominant precipitation since 1950 , impacting tree growth and reproduct ion . Mixed - and 5 - needle responded to winter precipitation , a possible hybrid condition . Furthermore, opposing trends in response to temperature were observed at the less dry site, which can dilute the climate signal when data are combined into a single chronology . Our results highlight important functional differences among morphotypes in response to climate. Tree phenology and sensitivity of growth functions were impacted by climatic shifts ; therefore, a climate trend that increases local moisture stress may impact hybridity and the stability of climate - growth relationships. Keywords: dendrochronology; ecology; Pinaceae; Pinus; Ponderosae ; response function ; tree rings 93 INTRODUCTION The Madrean Sky Islands are rugged mountain ranges isolated by desert that signify natural ecological laboratories heavily influenced by the last glacial period 10 - 20,000 years ago (Warshall 1995). The steep rocky terrain offers opportunities to study biologi cal populations across varied microclimate gradients and landscapes. Tree growth at the elevational limits of gradients are particularly susceptible to varying climate because small changes in the environment may have a large impact on annual growth (Jump et al. 2006; Tegel et al. 2014). In the Western United States, the consequences of rising temperatures and seasonal shifts toward earlier onset of spring and reduced snowpack are being observed at high elevation, and these trends in shifting climate are expected to increase the length of hydrological drought by the end of the century (USGCRP 2017). The Sky Islands of Southeastern Arizona include the well - known Santa Catalina Mountains, which contain habitat suitable for stands of ponderosa pine that include three morphological variants representing two distinct species (Rehfeldt 1999; Epperson et al. 2009) . At high elevation two species co - occur, Pinus arizonica Engelm. and the closely related P. ponderosa Lawson & Lawson var . brachyptera (Engelm.) L emmon . The latter species is also known as Taxon X and was previously misidentified as P. ponderosa Lawson & Lawson var. scopulorum Engelm. (Rehfeldt 1999; Epperson et al. 2009; Willyard et al. 2017), even though the two are clearly distinguishable by nee dle traits with high heritability (Rehfeldt 1992). P. ponderosa var. brachyptera exists in two forms, a nearly pure 3 - needle type that survives at the highest altitudes on southern slopes (2300 - 3000 m), northern slopes and cold air drainages, and a mixe d - needle 94 type that is interspersed with the 3 - needle type at transition zones. P. arizonica is a 5 - needle type found at lower elevations (below 2600 m) . Thus, the 5 - needle type is more successful at warmer and drier, lower elevations; whereas, the 3 - need le type survives at colder and wetter, higher elevations . On steep southern slopes the sharp transition of species occurs over just c. 130 meters slope distance (Epperson et al. 2001, 2009). Fritts 1963, 1974), part of an earlier meta - analysis which consisted of the three co - occurring needle types (3 - , mixed - and 5 - needle) sampled across a narrow climate gradient . We sampled the same population c. 50 - years later for a comparative analysis of similar aged cohorts, and also sampled the less dry Mt. Lemmon location to determine the effect of site conditions on growth responses. Our goal was to determine whether the cl imate - growth relationships have changed over the last century. Shifting seasonality or relationships with limiting factors would have implications for reliably predicting the vulnerability of tree species and needle types to climate change, conservation m anagement, and climate reconstruction. Previously, we reported on the two species of ponderosa pine displaying different growth responses to moisture stress that varied based on the microsite environment (Marquardt et al. 2018b). growth wa s reduced for longer periods by drought than P. ponderosa var. brachyptera , and the climatic response was greater at the site with higher soil moisture content. Since 1906, average global temperatures have been steadily increasing (0.74 °C + 0.18 °C; IPCC 2007). Because rising temperatures and predicted shifts in the monsoon season will influence weather patterns and drought in the Southwestern USA (Pascale et al. 2017), we hypothesized 95 that the seasonality of limiting factors will change over time. The objective of the study was to compare the correlations of ring widths with temperature and precipitation to assess shifts in the seasonality of climate - growth responses. METHODS Study area and Sampling Tree - ring cores were sampled from two rugged southern slopes in the Santa Catalina Mountains located 28 km northeast of Tucson, Arizona, USA: Mt. Lemmon (MTL) and Mt. Bigelow (BIG; Fig. 3. 1). BIG ( 32.41, - 110.71, 2534 m a.s.l. ) is drier than MTL ( 32.43, - 110.79 , 2577 m a.s.l. ), with average available water holding capacities of 3.8% and 9.2%, respectively. The climate of the desert southwest is semiarid and warm with two rainy seasons, the summer monsoons (July through September) and winter cyclones (Nove mber through March) each deliver up to half of the annual precipitation to the region (Sheppard et al. 2002). P. arizonica , and P. pond. var. brachyptera are dominant species in the mixed conifer forest. We sampled transition zones where all needle types grew together while avoiding cold - air drainages, which included the lower and upper moisture availability limits for P. ponderosa var. brachyptera and P. arizonica , respectively. Because the objective of this study is to compare our tree - ring data with t he analysis conducted c . 50 years ago by Fritts (1963, 1974), we analyzed a minimum of six trees of each needle type selected from a larger population of trees sampled from our earlier study (Marquardt et al. 2018b), which provides full site and sampling d etails. Trees that average > 4.6 needles per fascicle were designated P. arizonica (Peloquin 1984), and P. ponderosa var. brachyptera contained two needle types identified as a 96 mixed - needle tree (3.2 < mean < 4.6 needles per fascicle), and a 3 - needle tree (< 3.2 needles per fascicle; Rehfeldt 1993, 1999; Epperson et al. 2001). Climate data The local meteorological stations provided only short but accurate climate records; so, specific high - resolution gridded climate (1km) datasets were developed using the ANUSPLIN package (McKenney et al. 2011; Hutchinson and Xu 2013). Climate dataset devel opment and validation were described earlier (Marquardt et al. 2018b). Climatic variables were summarized to monthly values of average temperature (TAVG), average minimum temperature (TMIN), and total precipitation (PCP) because they were previously found to be significant drivers of growth responses (Marquardt et al. 2018b) . Tree growth chronologies We obtained tree - ring widths for P. arizonica [chronology (crn) 651000], P . ponderosa var. brachyptera (3 - needle type; crn 631000, and mixed - needl e type; crn 641000), located at 2500 meters in elevation on south facing slopes , from the Laboratory of Tree - Ring Research (LTRR, Tucson, AZ; Suppl. Table S 2; Suppl. Fig 3.S 1). These early period (E) data ( 1881 - 1960 ) were acquired for BIG as three composite chronologies that had been standardized with a negative exponential curve and converted to unit - less ring - width indices (RWI; Dodge 1963; Fritts 1963; 1974). In addition, we developed six new chronologies (3 from BIG, 3 from MTL) from raw ring - widths archived with the International Tree - Ring Data Bank (ITRDB; Suppl. 97 Table 3.S1 , Suppl. Figs 3.S 1 - S 2; Marquardt et al. 2018a) from trees characterized by averaged needle number per fascicle (Table 3. 1). Consistent with the early period d ata, the recent period (R; 1950 - 2007) raw chronologies were de - trended using a modified negative exponential curve to remove juvenile and geometric age trends, then converted to standardized ring - width indices using the ARSTAN program (Cook 19 85 ). The mea n chronologies (computed by ARSTAN or acquired from the LTRR) were visualized by using plotting functions in the dplR package (Bunn 2008; 2010) for R (R Core team 2016). Climate - growth relationships Our first analysis of climatic influence on tree - ring growth examined the correlations between climatic variables and ring - width index values using multivariate estimates obtained from the principal component regression model (Fritts 1974). We used the re lationship between local climatic variables and inter - annual growth to determine the differences in seasonal growth responses. Full methods of the response function correlations can be found in Marquardt et al. (2018b), with modifications: standard chron ologies were analyzed rather than ARSTAN chronologies; fewer climate variables were analyzed in this study; needle type was the subject of analysis rather than species; and two time periods were considered: early period and recent period. Briefly, the ana lysis computed bootstrapped response functions using the TREECLIM package (Zang and Biondi 2015), at three site x period combinations (BIG - E, BIG - R, MTL - R). Eight seasonal climatic variables (4 PCP and 4 TMIN) were partitioned to quantify precipitation an d temperature effects during peak times in our study areas. 98 Seasonal variables were derived by combining the monthly data into two rainy seasons consisting of three and five months respectively: summer (JAS) and winter (NDJFM) separated by three months o f arid spring (AMJ). For specific seasons with significant response function values (p < 0.05), those values lying above the axis show a positive response between tree - ring width and climatic variables, and values lying below the axis show a negative resp onse. To examine the impact (on growth) of the climate reversal in 1961 from dry winter precipitation to normal precipitation in 1962, we used the percent difference to quantify the growth variability between years with contrasting moisture conditions. To evaluate the temporal stability of climate - growth relationships, we compared climate data and composite chronologies using a univariate moving window correlation 2015). The 18 - month dendroclimatic window was set from September (current year) to April (previous year) with 30 - year moving interval and 5 - year offset. Two climate variables were considered separately for analysis. Total precipitation (PCP) and average tempe rature (TAVG) were divided into four seasons of three months each: winter (JFM), spring (AMJ), summer (JAS), and fall (OND). The previous summer and fall seasons were also considered, which increased the number of seasons analyzed to six. Temporal instab ility of the moving correlation functions was tested with G - test to determine which time series fluctuations were significantly different from those of a random time series (Zang and Biondi 2015). The TREECLIM results of the moving window analyses were pl 2017) for R. 99 RESULTS Seasonal climate - growth relationships Winter precipitation was averaged between sites for two years with contrasting moisture patterns. The winter of 1961 was one of the d riest on record (137.2mm), but the following year (1962) was of normal moisture conditions (319.43mm). Considering mean ring - width indices, growth for 5 - needle type is 28% greater on average in 1962 than 1961 (Table 2). A smaller increase in growth occur red for the 3 - needle type (20%) and mixed - needle type (16%). Current spring PCP was the strongest predictor of tree growth for all needle types for BIG - R (average r = 0.27 + 0.04) and MTL - R (average r = 0.36 + 0.04) after 1950. The strongest growth resp onses were positive, producing significant correlations with spring PCP for all needle types at both sites for the recent period (Fig 3. 2). In contrast, significant positive PCP - growth correlations were observed in summer rather than spring for BIG - E for all needle - types (aver r = 0.30 + 0.02; Fig 3. 2). The response function analysis (independent of period and location) indicated that 44% (4 out of 9) of mixed - and 5 - needle type populations recorded a significant climate signal during winter prior to the growing season (Fig. 3. 2). For the two periods, significant positive PCP - growth correlations in winter were observed only for BIG - E (aver age r = 0.25 + 0.03) and at MTL - R (average r = 0.33 + 0.03) for mixed - and 5 - needle types. Temporal stability of climate - growth relationships The PCP - growth relationships were stable; the G - tests for PCP are non - significant (p > 0.05) for all needle types analyzed across sites and periods (BIG - E; BIG - 100 R; MTL - R; data not shown). Because these correlation results supported the response function analysis reported above, we have only shown data for the PCP respons e function analysis. The TAVG - growth relationship was stable from 1896 - 1960 for BIG - E across needle types ( i.e., G - test p values > 0.05; Fig. 3. 3A) . The season with the largest number of individual significant negative correlations ( c . < - PrevSummer with 5 - 7 correlations (out of 8), followed by Summer with 4 - 5 correlations (out of 8). In contrast, the strong stable relationship between TAVG and growth became unstable during the recent period from 1958 to 2007 (p < 0.05 ; Fig. 3. 3B - C). This instability was indicated by significant G - tests, and changes in sign for seasonal TAVG - growth correlations, from negative to positive for PrevSummer for BIG - R (Fig. 3. 3B) and from positive to negative for PrevFall at MTL - R (Fig. 3. 3C ) across all three needle types. For BIG - R, the Fall season had the greatest frequency of significant positive correlations ( c - 5 correlations (out of 5). In comparison, the growth of MTL - R trees experienced mostly negative correl ations with TAVG. The season with the largest number of significant negative correlations ( c . < - was Winter with 2 - 3 correlations (out of 5) for all needle types, followed by Spring with 2 correlations (out of 5) for mixed - and 5 - needle typ es. DISCUSSION Seasonal climate - growth relationships Our chronology data highlighted a dramatic shift in annual growth in response to 2 levels began rising steadily (Kienast 101 and Luxmoore 1988). Mt. Bigelow ring widths (early period) correlated most strongly with summer precipitation, indicating moisture conditions were sufficient for stomatal conductance and annual growth to occur during the entire summer season (Levesque et al. 2014) ; t hus, maximizing biomass production . In comparison, recent climate - growth relationships were cool season correlations that occur when trees restrict their wood production to early spring while conditions are favorable for growth at dry sites (Fritts 1974; Levesque et al. 2014) . Thus, our recent season data indicated s pring precipitation was the most important variable for tree growth , a crucial time for the formation of male and female reproductive structures . Therefore, these results support our hypothes is that the seasonality of climate variables important to annual growth has shifted over time in response to changing climate. We also found that winter precipitation improved the growth of mixed - and 5 - needle types but not the 3 - needle type, suggesting h ybridization occurred between the pure needle types on southern slopes (Dodge 1963; Fritts 1963; Rehfeldt 1999; Epperson et al. 2009). We report important fundamental differences within and between species in seasonal response to climate suggesting that adaptation of tree species to length of growing season may have occurred during the last glacial period. These results support the influence of historical precipitation patterns and climate restrictions for different taxa (Norris et al. 2006; Kilgore 2007 ). Ziaco et al (2018) also determined the length of growing season of P. ponderosa in Nevada was controlled entirely by moisture conditions and not temperature by monitoring the cambial activity for two consecutive years with contrasting moisture, one we t winter and one dry spring. S ummer water was not used by the trees unless the start of the growing season was delayed by dry winter 102 and/ or spring conditions . Because moisture availability determines the length of growing season for all needle types of ponderosa pine studied in the Santa Catalina Mountains (Marquardt et al. 2018b), we hypothesize that one needle type may tolerate a shorter or longer growing season (5 - needle, 3 - needle, respectively; Dodge 1963), and a third needle type may grow best under a variable growing season (mixed - needle or hybrid; unpublished data; i.e., dry cool - season hypothesis). In support of this dry cool - season hypothesis, Dodge (1963) measured cambial activity during the dry 1961 growing season (average PDSI c. - 3.0; Marquar dt et al. 2018b), and found that P. ponderosa var. brachyptera (the 3 - and mixed - needle types) started cambial expansion under dry winter conditions in February, followed by the 5 - needle type ( P. arizonica ) six weeks later in April. Cambial growth occurred for all needle types during the summer, and cessation of growth coincided in September. All needle types co - existed at the Mt. Bigelow transition zone (2500 m a.s.l.) near our study site including the tre es that Fritts (1963) analyzed for climatic response. Although Dodge (1963) did not determine the initiation of radial growth for the following normal moisture year (average PDSI c. - 0.5; Marquardt et al. 2018b), we can predict the growth response of 5 - ne edle type to wet winter conditions from the recent literature on moisture driven cambial activity described above (Ziaco et al. 2018). Because 5 - needle type responds positively to winter precipitation, the taxon should initiate growth earlier. If the mix ed - needle tree is a hybrid possessing new adaptive gene complexes, then the mixed trees in this system may have inherited flexibility in using winter moisture (from 5 - needle type) and cold tolerance (from 3 - needle type; Kilgore 2007), providing increased f itness under variable growing season conditions. Further support was provided by our 103 chronology data, which shows that growth for 5 - needle type is robust when winter moisture conditions are normal (compared to dry winter). Also, an 8 - 12% increase in grow th for 5 - needle type was observed above that for the 3 - and mixed - needle types, respectively, for a normal moisture winter. This increase in growth is suggestive that 5 - needle starts cambial expansion earlier in wet springs than dry springs. Therefore, o ur study results reinforce that moisture availability determines the length of growing season of ponderosa pine, and that some needle types may display variability in the timing of xylogenesis (Ziaco et al. 2018). Temporal stability of climate - growth rel ationships We found opposing temperature trends in response to temperature after 1950 in the tree - ring data collected from the less dry Mt. Lemmon site, which can dilute the climate signal when data are combined into a single chronology . Wilmking et al. ( 2004) first reported opposing temperature - growth trends in white spruce stands occurring at tree - line in two mountain ranges sampled across Alaska, and temperature explained more variability in annual tree growth post - 1950. These regional forests are comp osed of individual trees growing in heterogeneous environments where temperature thresholds operate, influencing tree biology and the final chronology response to limiting factors. Our report extends these tree - line studies to semiarid regions primarily c ontrolled by precipitation rather than temperature as the factor most limiting to tree annual growth. We also observed a significant shift from positive to negative growth correlations with temperature under less dry conditions for 5 - and mixed - needle typ es during cool seasons, which is indicative of temperature - induced drought stress (Wilmking et al. 104 2004), or growing season respiration (Marquardt et al. 2018b) limiting the ponderosa pine annual growth. In comparison, at the drier site we observed a shif t from negative to positive temperature correlations during Fall season post - 1950, indicating reduced water stress. By extending our analysis to include the Fall season, we observed that growth in individual trees is positively and significantly correlate d to temperature. These results suggest that all needle types may respond positively to a range of temperature increases late in the season at the contact zone of dry habitats . However, an extended growing season most likely will be determined by moisture availability and not rising temperatures, as recently described for P. ponderosa, a closely related species (Ziaco et al. 2018). C ONC LUSIONS A temporal shift in limiting factors correlated with a shift in climate, indicating spring is the most important growing season for recent period; thus , the three taxa that established in summer wet habitats pre - 1950 are now responding to spring wet habitats post - 1950 , influencing reproduction . In comparison, 5 - and mixed - needle types have similar winter responses s uggesting soil moisture was controlling the length of growing season. The 3 - needle type did not respond to winter moisture, which may be suggestive of a hybrid condition . Although phenological plasticity plays a role in the formation of annual rings, our data suggest a genetic component is also acting on wood formation. In addition to the hybrid index developed using leaf morphology defined by the average needle number per fascicle, a full genetic analysis of the tree populations is required to provide f urther insight about hybrid formation and wood phenology under 105 increased evaporative demand, which determines length of growing season. Finally, we found significant fluctuations in temperature - growth correlations during recent period for all needle types , with the effects more variable at the less dry Mt. Lemmon site. Thus, shifts in climate that impact the growth sensitivity of trees, such as warming trends that increase local moisture stress, may also impact the stability of climate - growth relationship s. Supplementary materials : Supplementary data include Suppl. Table 3.S 1 : Ponderosa pine tree - ring data archives, Suppl. Fig. 3.S 1: Figures of average yearly chronologies for 3 - , mixed - , and 5 - needle types, Mt. Bigelow (Fritts 1963), Suppl. Fig. 3.S 2: Figures of average yearly chronologies for 3 - , mixed - , and 5 - needle types, Mt. Bigelow (new data), Suppl. Fig. 3.S 3: Figures of average yearly chronologies for 3 - , mixed - , and 5 - needle types, Mt. Lemmon (new data). Author contributions: P.M. conceive d the study, designed the experiments, and analyzed data; B.M assisted with data analysis; P.M and F.T. collected samples in the field; F.T. contacted the LTRR to obtain historical tree - ring data files; P.M. wrote the paper; F.T. and B.M. edited the paper. Conflicts of interest: The authors declare no conflicts of interest . 106 APPENDICES 107 APPENDIX F Supplementary Tables and Figures 108 Table 3.S 1 Ponderosa pine tree - ring data obtained from the ITRDB, except where noted. Three needle types were analyzed [5 - from an earlier period (pre - 1950) t han Marquardt et al. ( xxxxx chronologies will be deposited with ITRDB within 1 yr.) Site Type ID Lat Long Elev (m) N Scholar Source BIG 5 65100 110° 43´ 2500 6 Fritts LTRR, Tucson, AZ; citation doi.org/10.2307/1942448 BIG 3 63100 110° 43´ 2500 6 Fritts LTRR, Tucson, AZ; citation doi.org/10.2307/1942448 BIG M 64100 110° 43´ 2500 6 Fritts LTRR, Tucson, AZ; citation doi.org/10.2307/1942448 BIG 5 xxxxx 32.414 - 110.715 2534 7 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx BIG 3 xxxxx 32.414 - 110.715 2534 7 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx BIG M xxxxx 32.414 - 110.715 2534 6 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx MTL 5 xxxxx 32.434 - 110.790 2577 6 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx MTL 3 xxxxx 32.434 - 110.790 2577 7 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx MTL M xxxxx 32.434 - 110.790 2577 7 Marquardt et al. www.ncdc.noaa.gov/paleo/study / xxxxx 109 Fig ure 3.S1 Average yearly standard chronology for 3 - , mixed - , and 5 - needle types in the early BIG tree - ring study (Dodge 1963; Fritts 1963). Each chronology represents the mean of 12 cores from 6 trees for the reference period of 1881 - 1960. A = 3 - needle pine, B = mixed - needle pine, and C = 5 - needle pine . 110 Fig ure 3.S2 Average yearly standard chronologies for 3 - , mixed - , and 5 - needle types in the recent BIG tree - ring study. Each chronology represents an EPS > 0.82 and sample depth > 12, for the reference period of 1950 2007. A = 3 - needle type, B = mixed - needle type, and C = 5 - needle type 111 Figure 3.S3 Average yearly standard chronologies for 3 - , mixed - , and 5 - needle types in the recent MTL tree - ring study. Each chronology represents an EPS > 0.84 and sample depth > 15, for the reference period of 1950 2007. A = 3 - needle type, B = mixed - needle type, a nd C = 5 - needle type 112 APPENDIX G Tables and Figures 113 Table 3.1 Average needle counts for individual trees sampled for tree - ring correlation analyses at two sites: MTL and BIG Site Type Ave (#) Range (#) SD (#) N BIG 3 3.1 3.0 - 3.1 0.1 7 MTL 3 3.1 3.0 - 3.2 0.1 7 BIG M 3.7 3.2 - 4.6 0.6 6 MTL M 3.8 3.4 - 4.6 0.5 7 BIG 5 4.9 4.8 5.0 0.1 7 MTL 5 4.9 4.9 5.0 0.0 6 The needle counting method is described in Marquardt et al. (2018b) Type needle type, Ave (#) grand mean of average needle number per fascicle, SD standard deviation, N number of trees sampled. Three needle types were sampled: P. ponderosa var. brachyptera (3 - needle type), P. ponderosa var. brachyptera [mixed (M) needle type], and P. arizonica (5 - needle type) 114 Table 3.2 Percent difference (% DIFF) in growth between one dry winter (1961) and one normal winter (1962) for 3 - , mixed - , and 5 - needle types 3 - needle type Year BIG - R MTL - R AVE SD % DIFF 1961 730 769 749.5 027.6 20 1962 998 886 942.0 079.2 Mixed - needle type Year BIG - R MTL - R AVE SD % DIFF 1961 832 673 752.5 112.4 16 1962 967 828 897.5 098.3 5 - needle type Year BIG - R MTL - R AVE SD % DIFF 1961 716 627 671.5 062.9 28 1962 966 898 932.0 048.1 Growth data for each needle type were obtained from standardized chronologies developed for recent period tree - ring correlation analy ses at two sites: BIG - R and MTL - R 115 Fig ure 3. 1 Locations of the two populations sampled for growth analysis: Mt. Lemmon (MTL) and Mt. Bigelow (BIG). Black symbols locate study plots northeast of Tucson, AZ where 3 - , mixed - , and 5 - needle pines are sympatric 116 Fig ure 3. 2 Total PCP - growth relationships were determined by principal components multiple regression for three needle types [3 - needle = white , mixed - needle = blue, 5 - needle = grey] growing at two transition zones within the Santa Catalina Mountains. Within a plotted color, standard index compos ite chronologies analyzed were BIG - E (double hatched bars) for the early period of 18 81 1960 , and BIG - R (hatched bars) and MTL - R (non - hatched bars) for the recent period 1950 2007 . Climate variables analyzed were seasonal PCP: PrSummer (Jul - Sep), Winter ( Nov - Mar), Spring (Apr - Jun), Summer (Jul - Sep)] 117 (A) BIG - E (B) BIG - R (C) MTL - R Figure 3.3 TAVG - growth relationships were determined by moving window correlations. From top to bottom are correlations for 3 - , mixed and 5 - needle types. A = BIG - E data, B = BIG - R data, and C = MTL - R data. Color scale = correlation coefficients; significant correlations (p < 0.05) are denoted by *. Intervals during which a variable was not significantly correlated to tree growth are shaded pale red or pale blue ( c. + 0.2). 118 Figure 3.3 (continued) X axis = time periods and a significant G - test denoted by #, which means TAVG displays a change in correlation sign with tree growth over time (i.e. a sign of temporal instability). Y axis = seasons, which are defined as PrevSummer = (April - June), PrevFall = (July - September), Winter = (October - December), Spring = (January - March), Summer = (April - June), and Fall = (July - September). The standard index was used in the growth analyses 119 REFERENCES 120 REFERENCES Bunn AG (2008) A dendrochronology program library in R (dplR). Dendrochronologia 26 : 115 124 . https:// doi.org.10.1016/j.dendro.2008.01.002 Bunn AG (2010) Statistical and visual crossdating in R using the dplR library. Dendrochronologia 28 : 251 - 258 . https://doi.org.10.1016/j.dendro.2009.12.001 Cook ERA (19 85 ) Time s eries a nalysis a pproach to t ree r ing s tandardization. Thesis ( PhD ) , University of Arizona, Arizona Dodge RA (1963) Investigations into the ecological relationships of ponderosa pine in southeast Arizona. Thesis ( PhD ) , University of Arizona, Arizona , 117 pp Epperson B K , Telewski FW, Plovanich - Jones AE, Grimes JE ( 2001) Clinal differentiation and putative hybridization in a contact zone of Pinus ponderosa and P. arizonica (Pinaceae). Am J Bot 88 : 1052 1057 . https://doi.org.10.2307/2657087 Epperson BK, Telewski FW, Willyard A (2009) Chloroplast diversity in a putative hybrid swarm of Ponderosae (Pinaceae). Am J Bot 96 : 7 07 712 . https://doi.org.10.3732/ajb.0800005 Fritts HC (1963) Computer programs for tree - ring research. Tree Ring Bull 25 : 2 - 7 Fritts HC ( 1974) Relationships of r ing w idths in a rid - s ite c onifers to v ariations in m onthly t emperature and p recipitation. Ecol Monogr 44 : 411 440 . https://doi.org.10.2307/1942448 Hutchinson M, Xu T ( 2013 ) ANUSPLIN version 4.4 u ser g uide. Australia National University, Canberra, Australia IPCC (2007) Climate change 2007: the Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Pane l on 121 Climate Change -- Solomon S, Daahe Q, Manning M, Averyt K, Marquis M (eds) . Cambridge, Unite d Kingdom and New York, NY, USA Jump AS, Hunt JM, Penuelas J (2006) Rapid climate change related growth decline at the southern range edge of Fagus sylvatica. Glob Chang Biol 12 : 2163 - 2174 . https:// doi.org.10.1111/j.1365 - 2486.2006.01250.x Kienast F, Luxmoore RJ (1988) Tree - ring analysis and conifer growth responses to increased atmospheric CO2 levels . Oecologia 76 : 487 - 495 Kilgore J (2007) Distribution and ecophysiology of the Ponderosae in the Santa Catalina Mountains of southern Arizona , Thesis (PhD) , Michigan State University , Michigan , 263 pp Lévesque M, Rigling A, Bugmann H, Weber P, Bran, P (2014) Growth response of five co - occurring conifers to drought across a wide climatic gradient in Central Europe. Agr F orest M eteorol 197 : 1 - 12 . https://dx.doi.org/10.1016/j.agrformet.2014.06.00 1 Marquardt PE, Miranda BR, Telewski FW (2018a) Raw ring widths and crossdated series from: Variable c limate r esponse d ifferentiates the growth of s ky i sland p onderosa pines - International Tree Ring Data Bank c hronology series AZ595 - AZ598; corresponding study numbers 23650 - 23653. https://www.ncdc.noaa.gov/paleo/study/xxxxx Marquardt PE, Miranda BR, Jennings S , Thurston , G, Telewski FW (2018b) Variable c limate r esponse d ifferentiates the g rowth of s ky i sland p onderosa p ines, TREES . https://doi.org/10.1007/s00468 - 018 - 1778 - 9 McKenney D W , Hutchinson MF, Papadopol P, Lawrence K, Pedlar J, Campbell K, Milewska E, Hopkinson RF, Price D, Owen T ( 2011) Customized spatial climate models for North America. Bull Am Meteorol Soc 92 : 1611 1 622 . https://doi:10.1175/2011BAMS3132.1 Norris JR , Jacks on ST , Betancourt JL (2006) Classification tree and minimum volume ellipsoid analyses of the distribution of ponderosa pine in the western USA. J Biogeogr 33 : 342 - 360 . https:// doi.org/10.1111/j.1365 - 2699.2005.01 396.x 122 Pascale S, Boos WR, Bordoni S, Delworth TL, Kapnick SB, Murakami H , Vecchi GA, Zhang W (2017) Weakening of the North American monsoon with global warming. Nat Clim Chang 7 : 806 . https:// doi:10.1038/NCLIMATE3412 Peloquin RL ( 1984 ) The identification of three - species hybrids in the ponderosa pine complex. Southwest Nat 29 : 115 . https// doi:10.2307/3670776 R Core Team. (2016) R: A language and environment for statistical computing. R F oundation for Statistical Computing, Vienna, Austria. https://www.R - project.org/ Rehfeldt G E (1992) Early selection in Pinus ponderosa : compromises between growth potential and growth rhythm in developing breeding strategies. Forest Sci 38 : 661 - 677 . https:// doi.org/10.1093/forestscience/38.3.661 Rehfeldt GE (1993) Genetic variation in the Ponderosae of the southwest. Am J Bot 8 0 : 330 343 . https://doi:10.2307/2445357 Rehfeldt G E ( 1999) Systematics and g enetic s tructure of Ponderosae t axa (Pinaceae) i nhabiting the m ountain i slands of the Southwest. Am J Bot 86 : 741 . https://doi:10.2307/2656584 Sheppard PR, Comrie AC, Packin GD, Angersbach K, Hughes MK (2002) The climate of the US Southwest . Clim Res 21: 219 238 . https// doi:10. 3354/cr021219 T egel W, Seim A, Hakelberg D, Hoffmann S, Panev M, Westphal T, Büntgen U (2014) A recent growth increase of European beech (Fagus sylvatica L.) at its Mediterranean distribution limit contradicts drought stress. Eur J For Res 133 : 61 - 71 . https:// doi:10.1007/s10342 - 013 - 0737 USGCRP (2017) Climate s cience s pecial r eport: Fourth n ational c limate a ssessment -- Wuebbles D, Fahey K , Hibbard D , Dokken B , Stewart, Maycock T (eds). U.S. Global Change Research Program Washington, DC, USA , volume 1 , 470 pp . https// doi:10: 10 .7930/J0J964J6 W arshall P (1995) The Madrean sky island archipelago: a planetary overview. In: DeBano LF, Folliott PF , Ortega - Rubio A, Gottfried GJ, Hamre RH, Edminster CB (Eds.). Biodiversity and management of the Madrean a rchipelago: T he sky 123 islands of southwestern United States and northwestern Mexico , pp 6 - 18 . https:// doi.org/10.2737/RM - GTR - 264 Wei T, Simko V (2017) R package "corrplot": Visualization of a c orrelation m atrix (Version 0.84). https://github.com/taiyun/corrplot Willyard A, Gernandt DS, Potter K, Hipkins V, Marquardt P, Mahalovich MF, Langer SK, Telewski FW, Cooper B, Douglas C, Finch K (2017) Pinus ponderosa : a checkered past obscured four species. Am J Bot 104:161 181. https://doi.org/10.3732/ajb.1600336 Wilmking M, Juday GP, Barber VA, Zald HS (2004) Recent climate warming forces contrasting growth responses of white spruce at treeline in Alaska through temperature t hresholds. Glob Chang Biol , 10 : 1724 - 1736 . https:// doi.org/10.1111/j.1365 - 2486.2004.00826.x Zang C, Biondi F ( 2015) T reeclim: an R package for the numerical calibration of proxy climate relationships. Ecography 38 : 431 - 436 . https:// doi:10.1111/ecog.01335 Ziaco E, Truettner C , Biondi F, Bullock S (2018) Moisture driven xylogenesis in Pinus ponderosa from a Mojave d esert mountain reveals high phenological plasticity. Plant, C e ll Environ 41 : 823 - 836 . https:// doi.org/10.1111/pce.13152 124 C HAPTER 4 GENETIC DIFFERENTIATON AMONG SKY ISLAND POPULATIONS OF PONDERSOA PINE Paula E. Marquardt 1 * 1 USDA Forest Service, Nort hern Research Station, 5985 Highway K, Rhinelander, WI 54501 USA *Corresponding author: pmarquardt@fs.fed.us ; 715 - 362 - 1121; ORCID ID 0000 - 000 - 8854 - 0172 Acknowledgments : This paper is part of a dissertation submitted to Michigan State University in partial fulfillment of requirements for a Doctor of Philosophy degree. The Northern Research Station, USDA Forest Service, and the Dept. of Plant Biology, Michigan State Univer sity provided financial support for the project. We thank the following people for help on the project: F. Telewski helped conceptualize research questions . F. Telewski, J. Rayala, J. Grimes, G. Friedlander, and B. Epperson assisted with field data collection. B. Miranda, A. Willyard, H. Stricker, A. Foss, H. Jensen, S. Lietz, J. Lund, B. Prom and M. Meisenheimer provided lab, data processing, or GIS support. J. and L. G riffith, and A. and T. Harlan (both deceased) provided logistic support for the field research. K. Scribner reviewed the paper. We also thank the Coronado National Forest for providing permission to sample trees in the Santa Catalina Mountains . 125 ABSTRACT Discontinuous distributions limit gene flow and promote population differentiation through selective and stochastic processes. Ponderosa pine individuals were identified based on needle type (3 - , 5 - , and mixed), which easily distinguished between ecotypes of two closely related ponderosa pine species sampled from three contact zones in southeastern Arizona, USA . We collected similar data for P. arizonica and P. ponderosa var. brachyptera at lower and higher elevations; respectively, to describe pure - needl e types. Subsequently, the needle type identities were confirmed with haplotype frequencies representative of each phenotypic population sampled. Mixed needle is more closely related to P. ponderosa var. brachyptera than P. arizonica based on paternal lineage . Genetic analyses of chloroplast microsatellites indicated that variability is lower for the 5 - needle P. arizonica compared to 3 - and mixed - needle type. Higher divergence in haplotypes was observed in areas with distant groups (3 - and 5 - needle ty pes) than sympatry (3 - , mixed, and 5 - needle types) . Finally, integrating morphologic characters and genetic data is an effective way of identifying ponderosa pine taxa and quantifying population genetic structure, which can be used to classify geneticall y differentiated populations that could be considered separately in management plans. Keywords: edge effect; genetic structure; microsatellites; plastids; Pinaceae; Ponderosae 126 INTRODUCTION Genetic divergence among populations is often described as caused by geographic barriers or isolation by distance (Wright 1946; Slatkin 1987 ) despite high dispersal. The southwestern USA high elevation ponderosa pine forest is one system where species tra nsitions occur over short distances. Marquardt et al. (2018) recently showed that two closely related ponderosa pine species growing sympatrically have different ecological requirements for growth , primarily based on water availability . Ponderosa pine ( P inus subsection Ponderosae ) are large trees that grow throughout the mountainous and semi - arid regions of western North America, broadly distributed from British Columbia to northern Mexico. From the wide distribution of the Ponderosae (Conkle and Critchf ield 1988) evolved regionally adapted ecotypes that are collectively referred to as subspecies or varieties within a ponderosa pine complex, or as separate species (see Appendix Table 4. S 1 for glossary of terms). The coastal P. ponderosa var . ponderosa , a nd the interior P. ponderosa var . scopulorum , are two commonly accepted varieties with a sharp transition at their Montana contact zone (Potter et al. 2013; 2015; Shinneman et al 2016), corresponding to a climatic boundary (Norris et al. 2016). A recent biogeographic study of ponderosa pine describe d four distinct species ( P. benthamiana Hartw. , P. brachyptera Engelm. , P. ponderosa Laws. a nd P. scopulorum Engelm. (Willyard et al. 2017). Two of these species also co - occur in the southwestern United States ( P. brachyptera and P. scopulorum; refer to Appendix J f or taxonomy of P. brachyptera ). In the Santa Catalina Mountains, P. ponderosa var. brachyptera is a spatial mix of two needle types (3 - and mixed - needle), which is synonymous with Taxon X (Table 127 4.A 1 ). P. ponderosa var. brachyptera grows along with the well characterized P. arizonica at several transition zones on the mountain range (Epperson et al 2001, 2009), a possible glacial refuge for ponderosa pine during the Wisconsin epoch (Betancourt et al. 1990). P. arizonica is formally known as P. arizonica Engelm. A nd sometimes has been recognized as a subspecies or variety of P. ponderosa (Table 4.A2). These high elevation forests of the Santa Catalina Mountains are of interest because they are dom inated by the two ponderosa pine species that vary with elevation, and the differences in elevation are repeated throughout their ranges (Peloquin 1984). Not only do ecological differences exist between the two species from climate limitations (Marquardt et al. 2018), but P. arizonica is also often confused with P. ponderosa var. brachyptera south of the Mogollon Rim, where the ranges of the two species overlap , with possible hybridization and introgression (Farjon and Styles 1997), as a mixed - needle phenotype. Moreover, at transition zones in the Santa Catalina Mountains, similarity in climate - growth responses between P. arizonica and the mixed - needle phenotype is suggestive of a hybrid condition ( Marquardt et al. 2018 ) . An elevational gradient in average needle number per fascicle is shown to decrease from the lower to upper elevational limits of two ponderosa pine species sampled on Mt. Lemmon (Dodge 1963; Epperson et al. 2001). This gradient may be a morphological response of ponderosa pine to changing environmental conditions (Morphology hypothesis), or it may reflect the presence of hybridizing populations of P. arizonica and P. ponderosa var. brachypter a (Hybridization hypothesis) . The use of microsatellite DNA markers should allow discrimination between these two hypotheses presented. 128 In this study we analyzed the genetic structure of P. arizonica , P. ponderosa var. brachyptera , and a mixed - needle ponderosa pine with chlor oplast microsatellite markers, which are also useful for detecting admixture as fewer haplotype frequency changes are required for differentiation of the admixed and reference populat ions (Haasl and Payseur 201 1). We r eport preliminary results of tests for chloroplast microsatellites sampled from three transition zones (Mt. Lemmon, Mt. Bigelow, Palisades) and adjacent areas above and below the zone of contact. We examined evidence for genetic differentiation a mong ponderosa pine by considering habitat suitability for moisture availably. Our goal is to determine whether genetic differentiation can be attributed to the transition zone between P. arizonica and P. ponderosa var. brachyptera . Future work for hypot hesis testing of the chloroplast haplotypes will consider shifts in haplotype nucleotide frequencies. We will also analy ze nuclear microsatellite markers to examine hybridization and gene flow in the two species . METHODS Study s ite Gene frequencies were estimated from trees on three southern slopes in the Santa Catalina Mountains located 28 km northeast of Tucson, Arizona, USA: Mt. Lemmon (MTL), Mt. Bigelow, (BIG), and Palisades (PAL; Fig 4. 1 ; Table 4.1). A buildup of fuel from wildfire suppression in th e Santa Catalina Mountains resulted in destructively hot wildfires: Bullock (2002) and Aspen (2003) fires burned nearly 114,000 impacted the Mt. Lemmon and Palisades st udy areas, allowing regeneration to occur. 129 The semiarid climate of the desert southwest is warm with the summer monsoons (July through September) and winter cyclones (November through March) each delivering up to half of the annual precipitation to the reg ion (Sheppard et al. 2002 ; Fig 4.1 ). P. arizonica and P. pond. V ar. brachyptera dominate the canopy of the mixed conifer forest. P. ponderosa var. brachyptera contained two needle types, a mixed - needle tree (3.2 < mean < 4.6 needles per fascicle), and a 3 - needle tree (< 3.2 needles per fascicle; Rehfeldt 1993, 1999; Epperson et al. 2001). Trees that average > 4.6 needles per fascicle were designated P. arizonica (Peloquin 1984). We sampled from transition zones of steep altitudinal gradients where the distributions of taxa overlap and are potentially hybridizing, while avoiding cold - air drainages, which include the upper and lower moisture availability limits, respectively, for P. arizonica and P. ponderosa var. brachypt era . We also sampled pure populations in adjacent areas above and below the zone of contact. Sampled were three x 3 - needled sites above the elevation limit for P. arizonica ( c. 2700 m) near the summit of MTL (MLA, MLN, MLS; Fig. 4. 1 ; Table 4.1) , and two x 5 - needled populations below the elevation limit for P. ponderosa var. brachyptera ( c. 2200 meters) near Rose Canyon (ROS) and Green Mountain (GRN). Full site details for MTL and BIG are provided elsewhere (Marquardt et al. 2018). Sampling and M arker analysis Needle tissue was collected for DNA analysis on a 5 x 9 grid with 10 - m spacing. F rom 60 to 75 over story ponderosa pine ( P. arizonica, P.ponderosa var. brachyptera ) at each of three sites using a pole pruner during May and June, 2010 2012 (Mt. Lemmon, Mt. Bigelow, Palisades). In total, 195 adult trees were sampled at transition 130 zones. Eighty regenerated seedlings c. 10 60 cm in height were harvested on a 2 x 10 grid with 10 - m spacing at two transition zones during May 2011 2012 (Mt. Lemmon, Palisades). During 2012 - 13 , 20 - 40 seedlings per site were harvested at higher elevations near pure 3 - needle adults (Mt. Lemmon, Palisades). In total, 140 entire seedlings includ ing roots were harvested by hand. In addition, sixty pure 3 - needle trees ( P. ponderosa var. brachyptera ) were sampled on a 6 x 10 grid from three high - elevation Mt. Lemmon sites (MLA, MLN, MTS) and fifty - four pure 5 - needle trees ( P. arizonica ) were sample d from two low - elevation sites (Rose Canyon, Green Mountain). [Note: MLA was sampled along a horizontal transect, 1 tree every 30 meters]. In total, 114 pure - needle trees were sampled. All sample locations ( n = 449) were mapped using North American Dat um of 1983 (UTM - NAD83) through global positioning system (GPS) data gathered at the field sites. Community base station (Trimble Navigation Ltd., Sunnyvale, California) was used to differentially correct the GPS data for accuracy to 10 m (Adults only for M TL and BIG). Individual adult trees were tagged with aluminum tags. Fresh needle tissue was transported on dry ice, then freeze dried and stored at - 20C (2010), or dried and transported on silica desiccant, then stored at room temperature (2012 - 2013) unt il DNA extractions were completed. Total DNA was extracted from 20 mg dried needle tissue using Plant D n easy Mini Kits (QIAGEN Inc., Germantown, Maryland, USA) with following modifications to manufacturers protocol: needles were broken then ground dry for 1 minute using the Mini Beadbeader (BioSpec Products, Bartlesville, OK) with 0.25 - inch ceramic bead and 0.14g garnet sand in 2 - ml Fastprep extraction tubes. Dry ground needles were homogenized in AP1 buffer and R n ase for 1 minute. Eluted DNA eluted was quantified 131 using the Qubit fluorimeter (Life Techonologies, Carlsbad, CA, USA). Stock DNA solutions were diluted to 4.0 ng/ul in T10E1 (10 mmol/L Tris - Cl pH 8.0, 1 mmol/L EDTA) for amplification by the polymerase chain reaction (PCR). We surveyed vari ation in 6 - chloroplast microsatellite markers [ Pt71936 (cp13) , Pt100783 (cp15) , Pt87268 (cp14) , Pc10 (cp10) , PcG2R1 (cp9) , Pc2T1 (cp8) ; Vendramin et al. 1996, Stoehr and Newton 2002 ] that were multiplexed as a set of six mononucleotide loci with four color genotyping (Wofford et al. 2014) using three - primer fluorescent labeling during PCR (Culley et al. 2008, 2013). The reverse primer was primers, and the third primer contained the same unique tagged sequence with a - Fam, NED, VIC or PET ). Forward primers with tags ( Suppl. Table 4.S 2 ) were purchased from Integrated DNA Technologies, Coralville, Iowa, and the fluorescent primers from Lif e Technologies, Carlsbad, with the following modification; a 1um primer master mix with six forward, six reverse, and four labeled primers in a 1:4:4 (forward: reverse: labeled) volume ratio to limit forward primers (0.25um) to ¼ of the reverse and labelled primers (1um; Culley et al. 2008, 2013). Reaction mixtures contained 1ng/ul DNA, 1uM primers, and 3Mm MgCl in 10 - ul reaction volumes consisting of 5 ul (2x) Multiplex Master Mix, 2 ul (5x) primer mix, and 3 ul dH20. DNA was dried in wells prior to amplification. Cycling conditions were modified from Wofford (2014) by decreasing the number of cycles to thirty using Biorad Icyclers (Biorad Inc, Hercules, CA ): 95C for 15 min, followed by 30 cycles each of 94C for 30 sec., 58C for 90 sec, 72C for 90 sec, and final extension at 72C for 10 min. 132 A subsample of each PCR was evaluated for correct size product and quantified by agarose gel electrophoresis (2 ul of reaction mix ture on 2% gel). PCR products were diluted 1:20 for fragment analysis . S amples were run on an ABI 3730xl DNA analyzer (Life Technologies), with the GS600 LIZ as the internal size standard at the Purdue Genomics Core Facility, and fragments visualized and recorded using GeneMarker genotype analysis software (Softgenetics, State College, PA). Haplotypes were confirmed by independent scoring by two trained lab oratory personal . Thirteen - percent of multiplexed haplotypes were randomly selected for retesting in pcr and fragment analysis to assess scoring error (< 0.07 %) . Diversity and Population s tructure We defined a population as a discrete group of individuals of the same species, and samples collected from a population can be readily identified by a common morphological characteristic that is easily recorded and analyzed. Genetic profiles were develope d from nearly pure 3 - and 5 - needle pine populations selected from sampled trees removed from the contact zone. H aplotype diversity w as calculated separately for each population and locus. F irstly, we analyzed the plastid genetic structure of 21 populati ons of ponderosa pine (Table 4.1) by principal coordinate analysis (P C oA) to resolve clusters based on a distance matrix among all pairs of populations using GenAlEx6 (Huff et al. 1993; Peakall and Smouse 2012 ) . Genetic PC oA (an ordination technique) visually characterizes population structure by delineating species and populations wit hin contact zones (Mallet 1995). S amples clustered closer together have smaller genetic distance values than samples clustered further apart. Secondly, convex 133 polygons were applied to more easily visualize genetic structure. Thirdly, we evaluated whether haplotype frequencies of mixed - needle population s disproportionately reflect frequencies of one pure population ( suggesting genetic relatedness) by plotting a histogram of the first axis of P C oA . Fourthly , we assessed whether genetic differentiation between species in distant populations exceeds that in sympatry, as expected with hybridization and introgression (pairwise population PhiPT values, Peakal l and Smouse 2012). Lastly, genetic differentiation ( using a distance matrix among haplotypes) was analyzed by analysis of molecular variance (AMOVA) using Poppr function ( Kamvar et al. 2014 ) in R ( ade4 package ; Bougeard and Dray 2018 ; R Core Team 2016 ) . W ithi n and among region phiRT and pairwise phiRT was calculated with GenAlEx6 . Populations were divided into three regions based on the results of the P C oA (i.e. the three needle types: 3 - , 5 - and mixed - needle). Histograms, polygons, and P C Wickham 20 10 ) for R . RESULTS Identification of groups In our study, all tree populations had average needle counts that ranged from 3.0 to 4.9 (Table 4.1). The P C oA based on chloroplast haplotype frequencie s showed that all 21 populations grouped with the three needle types (Figs. 4.3 , 4.4) except for one mixed - needle and one 5 - needle population near 0.05 on the x axis (Fig. 4. 4 ). There is a gap between the mixed - needle populations in the center of the distribution (between 0 and 0.05 on the x axis) . Mixed needle populations var ied in number of needles fo rmed from year to year, allowing them to be distinguished from pure needle populations based on 134 the SD of average needle number per fascicle ( NF; Table 4.1). Pines retain needles, and growth is counted and averaged over 3 - 5 years that needles are retained . We set SD > 3.0 as the limit for mixed - type; subsequently, five populations were selected as mixed - needle (BGM, LSI2, MLM, PLM, PSI2) . The three sigma rule allows for 99.7% of values to fall within three standard deviations of the mean ( for n ormal dist ributions ) , ind icating a significant result . Two more population s were selected as mixed - needle ( bringing the total to 7) based on > 3.2 average needle count per fascicle (LSI2; NF=3.32 ; Rehfeldt 1999 ) or SD > 2 (PSI1; SD=.21 ). T wo sigma rule allows for 95% of the population to lie within two standard deviations of the mean . Genetic variation among haplotypes and populations The chloroplast locus was polymorphic in all populations and diversity was high (Table 4.2). Among the populations, a total of 42 multi - locus haplotypes were detected in 427 samples based on a 6 - locus chloroplast haplotype , the most common haplotype is 12 (Table 4. S 3). Referring to Table 4.2, the number of haplotypes ranged from 4 to 1 6 per population . T he adult 3 - needled populations (11.5 + 2.8 ; n=6 ) had more haplotypes on average + SD than the adult 5 - needled populations (6.2 + 1.6 ; n=5 ). Haplotype 12 is present in all popu lations but MLN, ranging in frequency from 0.053 - 0.741 (Table 4.S4) . All popula tion samples showed moderate to high unbiased haplotype diversity, which ranged from 0.486 to 0.94. In comparison , 5 - needle (n=5) had the lowest diversity (0.453 - 0.600), contrasted to 3 - needle (n=9) with the highest diversity (0.914 - 0.943). Diversity of m ixed needle was more varied, ranging from 0.638 - 0.911 (n=7). The number of private haplotypes was low and ranged from 0 to 2 per population , the 135 majority occurring in reference populations: 2 private haplotypes occurred in GRN (5 - needle), 2 in MLA (3 - nee dle), and 2 in MLS (3 - needle). Few private haplotypes occurred at the transition zones: 1 each in 5 - needle populations BGZ and PLZ, 1 in 3 - needle population MLB, and 1 in mixed needle population PLM. Population s tructure The AMOVA results indicated that 13% and 87% of the haplotype variation can be attributed to variation among regions (needle types) and among individuals within populations, respectively (p < 0.01; Table 4.3). There was no variation among populations with in needle types. The average pairwise PhiPT values (Table 4. 4 ) were large between distant populations and lowest between populations in sympatry. The distant populations represent the pure needle reference populations (3 - needle vs 5 - needle) and the popul ations in sympatry represent the populations in contact at the transition zone (3 - , 5 - and mixed - needle) . The PhiPT pairwise values are consistently highest between populations of 5 - and mixed - needle pine s (e.g. MLZ - BGB ; Table 4.5 ) . In comparison, 3 - and mixed - needle pines are most similar in genetic structure ( e.g. BGB - MLB ; Table 4.5 ). P C oA revealed three groups of ponderosa pine populations (Fig. 4. 3 ). Axes 1 and 2 accounted for 81.58% and 8.53% of the total variation. The 3 - needle populations group i n the left quadrants of the first principle coordinate (e.g. MLN, ML S , PSU1 ) , and the 5 - needle populations group in the right upper quadrant (e.g. PLZ, ROS, MLZ ) . 136 Mixed needle populations group between 3 - and 5 - needle populations of the first principle c oordinate (e.g. MLM , PSI2, LSI2) . DISCUSSION We analyzed the spatial population genetic structure of P. arizonica and P. ponderosa var. brachyptera on the Santa Catalina Mountains and found evidence for three distinct genetic groups ; g enetic variability is reduced for the low elevation speciees compared to the high elevation specie s . Below we discuss the relevance of our results with respect to the management of the Sky Island pine forest. Population structure and D ifferentiation The ponderosa pine populations from the Santa Catalina Mountains cluster genetically into three distinct groups based on chloroplast haplotype frequency . The mixed - needle populations clustered between the 3 - and 5 - needle populations (Figs 4.3; 4.4). The gaps in the center of the histogram plot (between 0 and 0.05 on the x axis) could be indicative of distant populations being differentiated from transition zone populations (Adams 1982). The groups corresponded to the environmental boundary of high and low elevation slopes, which separates cool - moist habitat from dry - hot habitat (Fig. 4.2). Geographically, P. arizonica is located at lower elevation s , P. ponderos a var. brachyptera is located at higher elevation s , and mixed - needle type is located between populations of the two species. Consistent with a hypothesis about climatic barriers to gene flow, the mixed - needle populations were closest in genetic structure to 3 - needle 137 populations. D ivergence was lowest in the transition zone indicating higher levels of gene flow than between distant populations. An earlier study of P. arizonica and P. ponderosa var. brachyptera on Mt. Lemmon used two of the six chloroplast loci in our study (Epperson et al. 2009). The four loci used in the 2009 study revealed similar haplotype diversity (mean h = 0.85) to our study and approximately based on two groups ( P. arizonica and P. ponderosa var . brachyptera ), our measure (Phi) was based on three groups (5 - , 3 - , mixed). Both studies found significant differences in genetic structure for the two ponderosa pine species . A third study using chloroplast haplotypes to analyze ponderosa pine genetic structure did not resolve P. ponderosa var. brachyptera from P. arizonica (Willyard et al. 2017) , which suggests these two species are more closely related to each other than t hey are to western ponderosa pine species or to P. ponderosa var. scopulorum . Like our study, the authors report ed reduced diversity for P. arizonica as compared to other ponderosa pine species sampled in this range - wide study. For our study , we found th at 5 - needle type encompassed c. 40% lower levels of haplotype diversity than did 3 - needle type , which is consistent with P. arizonica being near its northern range boundary and having experienc ed a bottleneck during range expansion (Eckert et al. 2008). In comparison , a range wide study of ponderosa pine with nuclear microsatellite markers and isozymes did not find reduced genetic diversity near edges of species ranges (Potter et al. 2015). Epperson et al. (2009) conducted spatial autocorrelation analysis of 3 - , 3 - plus mixed - , and 5 - needle types separately, indicating randomly mating populations. A 138 fourth analysis was of the entire transect (all needle types were grouped together), which i ndicated strong autocorrelation at short distances. Although the significance of this result was not discussed by the authors, the pattern of fine scale genetic structure within populations (of the needle types combined) may decrease the chances of hybrid ization thru reduced gene flow . Our study results differ from Epperson et al. (2009) who reported no correlation between needle number and haplotypes for P. ponderosa var. brachyptera (3 - plus mixed - needle types). Our research found a strong relationshi p between the two factors in our principle coordinate analysis, whether we grouped mixed - needle with 3 - needle or analyzed them separately. One main difference between the two studies is the type of analyses conducted . Epperson et al. (2009) analyzed chlo roplast alleles as independent observations in a spatial analysis . Our study used PCoA look ed for a relationship between haplotype genetic distance and needle type , an easy to measure trait useful for identifying genetic structure. CONC LUSIONS Two competing hypotheses explain the climate gradient in average needle number observed in southeast Arizona, the greatest expression of which is found in the Santa Catalina Mountains (Dodge 1963). This gradient may be a response to environmental conditio ns or it may reflect the presence of hybridizing populations of P. arizonica and P. ponderosa var. brachyptera . The morphology hypothesis states that the gradient in average needle number is a habitat response, indicated by a reduction in needle number th at is expressed in trees more suited to cold, moist environments, i.e. the 3 - needle type is more tolerant of cold 139 than 5 - needle type (Kilgore 2007). In comparison, an increase in average needle - number will be expressed in trees that show greater effects o f cold damage on cold sites (Haller 1962, 1965) , or adaptation with water loss on dry sites ( this thesis ) . T wo separate taxa occur on the Santa Catalina Mountains; P. arizonica, a 5 - needle pine and P ponderosa var . brachyptera , a putative spatial mix of 3 - and mixed - needle types (Epperson et al. 2009). Mixed needle type contains 3 - 4 - 5 needles per fascicle, which occurs only at the contact zones of the two species, and the gradient in average needle number increases from high to low elevation. If the mix ed - needle trait is a result of habitat response to arid conditions the previous fall (unpublished data), then variance in needle number possibly modulates inter - annual water stress, increasing carbon gain. The gradient could also be the result of hybridizi ng populations of P. arizonica and P. ponderosa var . brachyptera that were established in a colder and drier environment, where pure individuals outperformed the hybrids (or there was no opportunity to hybridize). Then when changes in phenology occurred b ecause growing seasons became longer and/ or more variable based on water availability (Ziaco et al. 2018), interbreeding also had opportunity to occur. Hybrids that formed may have outcompeted in the warmer climate with increased evaporative demand, and if they were fertile, could have interbred with parents. The ponderosa pine populations of the Santa Catalina Mountains may be in jeopardy because high elevation species not only occupy a small but important niche, they are also more vulnerable to changing environmental conditions. We found high levels of genetic variation for 3 - and mixed - needle types; levels are reduced by c. 40% for P. arizonica compared to P. ponderosa var. brachyptera . These results highlight the 140 dissimilar genetic structure of the two species, suggesting they should be considered differently in management plans in areas where they coexist. P. arizon ica populations are near their n orthern range boundaries, have reduced diversity and are most at risk from climate change because they may be less fit due to environmental stress; therefore, conservation efforts may in want to represent these peripheral po pulations. The three taxa form three distinct genetic groups reflective of elevational niches (3 - , mixed, and 5 - needle) ; t he upper and lower elevation divisions contribute most to the observed genetic variation with divergence greatest among distant populations. 141 APPENDICES 142 APPENDIX H S upplementary Tables and Figures 143 Table 4. S1 Glossary of genetic terms Admixture: Two or more previously isolated populations begin interbreeding, introducing new genetic lineages into a population. Conspecific: Belonging to the same species. For example, it describes the interactions between two or more individuals of the same species. Ecotype: a genetically distinct group that is adapted to its local environment. An ecotype is a variant in which the phenotypic differences are too few or too subtle to warrant being classified as a subspecies. Hybridization: the movement of genes from one species into the gene pool of another resulting in an even mix of two parental species in the first generation. Hybridity: The offspring of genetically dissimilar parents produced by breeding plants or animals of different varieties, sp ecies, or races. Infraspecific: variation or taxonomic division within a species, e.g. subspecies or variety. Introgression: the incorporation of alleles from one species into the gene pool of a second species, usually through hybridization and repeated backcrossing resulting in a complex mixture of parental genes. Morphology: the study of the form and structure of organisms and their evolutionary history. Morphotype: in taxonomy, an alternate form or phenotype to illustrate variation within a populatio n (of a single species). Population: a locally discrete group of individuals of the same species where individuals can be readily identified by a common morphological characteristic. Population structure: populations are subdivided in some way, which allows them to diversify and evolve. Such populations have deviations from Hardy - Weinberg proportions, for example when there is inbreeding, selection, migration, or genetic drift. Species: a genetically distinct group that forms a single cluster from a genetic sample if all individuals belong to a single species, and multiple clusters if there are two or more species present. Subspecies: a taxonomic group that is a division of a species, and usually occurs from migration barrier s causing isolation within a species. 144 Table 4. S1 ( continued ) Variety: a unique group that displays genetic divergence that is not considered high enough to be a different subspecies. The origin of a new variety may reflect neutral mutation without changes in fitness and would not be a different ecotype. Also, the variety may be the result of adaptation by natural selection producing an ecotype . 145 Table 4.S 2 F luorescently labeled unique sequences used in 3 - primer PCR Primers amplify m ononucleotide microsatellite repeats in the chloroplast genome (Vendra m in et al. 1996 ; Stoehr and Newton 2002 ). Refer to Wofford et al. 2014 (Table 3) for repeat information and locus characteristics. In 3 - primer PCR (Culley et al. 2008) , the forward pri mer incorporates a unique tailed sequence, the reverse primer is unlabeled, and the third primer incorporates the same tailed sequence that is attached to a fluorescent label 146 Table 4.S3 Summary information for the 42 chloroplast haplotypes analyzed Haplotype Count 1 2 3 4 5 6 Ex. Sample Ex. Pop 1 1 277 104 206 163 176 125 PALS001 PSI1 2 1 277 104 206 163 177 125 GRN008 GRN 3 1 277 104 207 163 177 125 PALS051 PSU1 4 1 277 104 207 163 178 125 PAL006 PLZ 5 1 277 104 207 164 178 124 BIG043 BGZ 6 1 277 104 207 164 179 124 MTLA011 MLA 7 1 277 105 206 161 177 125 MTLS034 LSI2 8 3 277 105 206 162 176 125 ROS079 ROS 9 3 277 105 206 162 177 125 PAL004 PLZ 10 8 277 105 206 162 177 126 MTL277 MLZ 11 18 277 105 206 163 176 125 ROS094 ROS 12 166 277 105 206 163 177 125 ROS093 ROS 13 3 277 105 206 163 178 125 PAL002 PLZ 14 1 277 105 206 163 179 126 GRN009 GRN 15 12 277 105 206 164 177 125 PALS049 PSI2 16 17 277 105 206 164 178 124 PALS090 PSU2 17 15 277 105 206 164 179 124 PALS018 PSI1 18 4 277 105 206 164 180 124 PALS087 PSU2 19 2 277 105 206 164 181 124 MTL050 MLM 20 1 277 105 206 165 179 124 MTLA004 MLA 21 1 277 105 207 162 177 125 PAL071 PLM 22 5 277 105 207 162 179 124 PALS069 PSU1 23 28 277 105 207 163 177 125 ROS083 ROS 24 5 277 105 207 163 177 126 MTLA027 MLA 25 9 277 105 207 163 178 124 PALS091 PSU2 26 2 277 105 207 163 178 125 PAL039 PLB 27 3 277 105 207 164 177 124 PALS015 PSI1 28 1 277 105 207 164 177 125 MTL239 MLS 29 32 277 105 207 164 178 124 PALS089 PSU2 30 33 277 105 207 164 179 124 PALS088 PSU2 31 1 277 105 208 162 179 124 MTL236 MLS 32 3 277 105 208 163 177 125 PALS084 PSU2 33 11 277 105 208 164 177 124 PALS016 PSI1 34 1 277 105 208 164 177 125 PALS035 PSI2 35 15 277 105 208 164 178 124 PALS083 PSU2 36 3 277 105 208 164 179 124 PALS057 PSU1 37 6 277 105 208 165 178 124 PALS063 PSU1 147 Table 4.S3 ( continued ) Haplotype Count 1 2 3 4 5 6 Ex. Sample Ex. Pop 38 1 277 105 209 164 179 124 MTLS065 LSU 39 2 277 105 209 164 180 124 MTLS050 LSI2 40 1 277 106 207 164 178 124 MTL249 MLB 41 1 278 105 207 164 178 124 PALS075 PSU2 42 3 278 105 207 164 179 124 MTL123 MLN Count = number of haplotypes, 1 - 6 = mononucleotide loci, Ex. = Example 148 Table 4. S 4 Hapl otype frequencies (42 haplotypes; 1 locus; 427 samples; N, number of individuals per 21 population s ) BGZ GRN MLZ ROS PLZ BGM LSI1 LSI2 MLM PLM PSI1 PSI2 N 28 30 27 24 30 10 15 22 11 9 23 14 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 2 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 4 0.000 0.000 0.000 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 5 0.036 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 6 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 7 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.045 0.000 0.000 0.000 0.000 8 0.000 0.000 0.000 0.042 0.000 0.000 0.133 0.000 0.000 0.000 0.000 0.000 9 0.000 0.000 0.000 0.000 0.033 0.100 0.067 0.000 0.000 0.000 0.000 0.000 10 0.000 0.000 0.074 0.000 0.000 0.000 0.133 0.182 0.000 0.000 0.000 0.000 11 0.036 0.067 0.074 0.125 0.033 0.000 0.067 0.000 0.182 0.000 0.000 0.286 12 0.679 0.700 0.741 0.708 0.600 0.300 0.600 0.500 0.273 0.444 0.348 0.429 13 0.036 0.000 0.000 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 14 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 15 0.107 0.067 0.037 0.000 0.033 0.000 0.000 0.045 0.000 0.000 0.043 0.071 16 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 17 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.091 0.000 0.087 0.000 18 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 19 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.091 0.000 0.000 0.000 20 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 21 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.111 0.000 0.000 22 0.000 0.000 0.000 0.000 0.000 0.100 0.000 0.000 0.000 0.000 0.043 0.000 23 0.071 0.033 0.037 0.125 0.233 0.100 0.000 0.045 0.091 0.222 0.000 0.143 24 0.036 0.033 0.000 0.000 0.000 0.100 0.000 0.000 0.000 0.000 0.000 0.000 25 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.111 0.000 0.000 26 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 149 Table 4.S 4 ( continued ) BGZ GRN MLZ ROS PLZ BGM LSI1 LSI2 MLM PLM PSI1 PSI2 N 28 30 27 24 30 10 15 22 11 9 23 14 27 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 28 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 29 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.182 0.000 0.087 0.000 30 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.091 0.000 0.000 0.174 0.000 31 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 32 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 33 0.000 0.000 0.000 0.000 0.000 0.100 0.000 0.000 0.091 0.000 0.043 0.000 34 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.071 35 0.000 0.000 0.037 0.000 0.000 0.200 0.000 0.000 0.000 0.000 0.043 0.000 36 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 37 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.111 0.000 0.000 38 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 39 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.091 0.000 0.000 0.000 0.000 40 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 41 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 42 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 150 Table 4.S 4 ( continued ) BGB LSU MLA MLN MLS MLB PLB PSU1 PSU2 N 21 20 30 14 15 21 23 21 19 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.048 0.000 4 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 5 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 6 0.000 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 7 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 11 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 0.053 12 0.143 0.200 0.100 0.000 0.200 0.095 0.217 0.286 0.053 13 0.000 0.000 0.000 0.000 0.000 0.048 0.000 0.000 0.000 14 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 15 0.048 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 16 0.000 0.050 0.100 0.286 0.200 0.095 0.043 0.000 0.105 17 0.095 0.150 0.067 0.214 0.000 0.048 0.043 0.000 0.000 18 0.095 0.000 0.000 0.000 0.000 0.000 0.043 0.000 0.053 19 0.048 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 20 0.000 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 21 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 22 0.000 0.000 0.033 0.000 0.000 0.000 0.043 0.048 0.000 23 0.000 0.050 0.033 0.071 0.067 0.143 0.000 0.000 0.000 24 0.048 0.000 0.033 0.000 0.000 0.000 0.000 0.000 0.000 25 0.000 0.050 0.067 0.071 0.000 0.000 0.000 0.095 0.105 26 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 0.000 151 Table 4.S 4 ( continued ) BGB LSU MLA MLN MLS MLB PLB PSU1 PSU2 N 21 20 30 14 15 21 23 21 19 27 0.000 0.000 0.000 0.071 0.000 0.048 0.000 0.000 0.000 28 0.000 0.000 0.000 0.000 0.067 0.000 0.000 0.000 0.000 29 0.095 0.300 0.200 0.143 0.200 0.048 0.174 0.048 0.158 30 0.190 0.050 0.067 0.071 0.067 0.143 0.217 0.286 0.211 31 0.000 0.000 0.000 0.000 0.067 0.000 0.000 0.000 0.000 32 0.000 0.000 0.000 0.000 0.000 0.000 0.043 0.000 0.105 33 0.048 0.050 0.067 0.000 0.067 0.143 0.000 0.000 0.000 34 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 35 0.095 0.050 0.067 0.000 0.000 0.095 0.043 0.048 0.105 36 0.048 0.000 0.000 0.000 0.000 0.000 0.043 0.048 0.000 37 0.000 0.000 0.033 0.000 0.067 0.048 0.000 0.095 0.000 38 0.000 0.050 0.000 0.000 0.000 0.000 0.000 0.000 0.000 39 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 40 0.000 0.000 0.000 0.000 0.000 0.048 0.000 0.000 0.000 41 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.053 42 0.048 0.000 0.033 0.071 0.000 0.000 0.000 0.000 0.000 Populations are grouped by needle - t ype: five 5 - needle (e.g. BGZ, GRN, MLZ, ROS, PLZ) , seven mixed - needle populations, and nine 3 - needle populations 152 APPENDIX I Tables and Figures 153 Table 4.1 Summary information for the 21 pine populations analyzed in this study SITE POP N TYPE NF SD SPECIES COHORT ELEV PAL PLB 23 3N 3.05 0.07 brachyptera adult 2501 BIG BGB 21 3N 3.05 0.06 brachyptera adult 2510 MTL MLB 21 3N 3.08 0.08 brachyptera adult 2568 MTL MLA 30 3N 3.11 0.17 brachyptera adult 2756 MTL MLS 15 3N 3.03 0.06 brachyptera adult 2761 MTL MLN 14 3N 3.01 0.02 brachyptera adult 2768 PAL PSU2 19 3N 3 .00 0.09 brachyptera seedling 2553 PAL PSU1 21 3N 3.03 0.16 brachyptera seedling 2566 MTL LSU 20 3N 2.94 0.08 brachyptera seedling 2741 GRN GRN 30 5N 4.87 0.11 P. arizonica adult 2155 ROS ROS 24 5N 4.8 0 0.22 P. arizonica adult 2172 PAL PLZ 30 5N 4.91 0.08 P. arizonica adult 2457 BIG BGZ 28 5N 4.9 0 0.08 P. arizonica adult 2498 MTL MLZ 27 5N 4.88 0.09 P. arizonica adult 2529 PAL PLM 9 MN 3.82 0.59 brachyptera adult 2488 BIG BGM 10 MN 3.85 0.6 0 brachyptera adult 2514 MTL MLM 11 MN 3.98 0.49 brachyptera adult 2557 MTL LSI1 15 MN 3.32 0.13 brachyptera seedling 2453 MTL LSI2 22 MN 4.14 0.36 brachyptera seedling 2454 PAL PSI2 14 MN 4.12 0.29 brachyptera seedling 2470 PAL PSI1 23 MN 3.14 0.21 brachyptera seedling 2471 Five study sites are Mt. Lemmon (MTL), Mt. Bigelow (BIG), Palisades (PAL), Rose Canyon (ROS), and Green Mountain (GRN) which comprise 21 ponderosa pine populations that were sampled for genetic analysis. POP is population, NF is average number of needles per fascicle, brachyptera is P. pondersoa var. brachyptera , and mean elevation of the area sampled for each population is reported in meters 154 Table 4.1 ( continued ) The study areas included three transition zones (where two species of ponderosa pine grow together) located on Mt. Lemmon (MLB, MLZ), Palisades (PLB, PLZ), and Mt. P. ponderosa v ar. brachyptera and P. arizonica , respectively. Also sampled were reference groups of three pure 3 - needle sites located at high elevation on North Slope (MLN), MTLA (MLA), and South slope (MLS), and two pure 5 - needle sites sampled at low elevation on Rose Canyon (ROS) and Green Mountain (GRN). In addition to these adult populations, seven groups of ponderosa pine seedlings were sampled: three at elevation 3 - needle group (LSU) a nd two intermediate elevation mixed - needle groups - needle groups (PSU1; PSU2) and two intermediate elevation mixed - needle groups (PSI1; PSI2) 155 Table 4.2 Sites, needle type, populations, sa mple sizes, and haplotype polymorphism characterized for 21 populations of ponderosa pine at 5 sites Site T ype Pop n ph Na h uh BIG 5 BGZ 28 1 7 0.512 0.537 GRN 5 GRN 30 2 8 0.496 0.513 MTL 5 MLZ 27 0 6 0.436 0.453 PAL 5 PLZ 30 1 7 0.580 0.600 ROS 5 ROS 24 0 4 0.465 0.486 BIG M BGM 10 0 7 0.820 0.911 MTL M LSI1 15 0 5 0.596 0.638 MTL M LSI2 22 2 7 0.694 0.727 MTL M MLM 11 0 7 0.826 0.909 PAL M PLM 9 1 5 0.716 0.806 PAL M PSI1 23 1 11 0.820 0.858 PAL M PSI2 14 1 5 0.704 0.758 BIG 3 BGB 21 0 12 0.893 0.938 MTL 3 LSU 20 1 10 0.830 0.874 MTL 3 MLA 30 2 16 0.909 0.940 MTL 3 MLN 14 0 8 0.827 0.890 MTL 3 MLS 15 2 9 0.853 0.914 MTL 3 MLB 21 1 12 0.898 0.943 PAL 3 PLB 23 0 12 0.858 0.897 PAL 3 PSU1 21 1 9 0.807 0.848 PAL 3 PSU2 19 1 10 0.875 0.924 mean 20.3 0.8 8.4 0.734 0.779 SD 6.6 0.7 3.0 0.157 0.170 n, number of individuals; ph , number of private haplotypes; Na, number of different haplotypes; h, haplotype diversity = 1 - Sum pi^2; uh, unbiased diversity = (N / (N - 1)) * h; where pi is the frequency of the ith allele for the population and Sum pi^2 is the sum of the squared popul ation haplotype frequencies 156 Table 4.3 Analysis of molecular variance of plastid SSRs for 9 populations of ponderosa pin e at contact zones in southwest Arizona, USA Hierarchical level* df SS MS Est Variance Var (%) Stat Phi p Among needle types (3, 5, M) 2 13.404 6.702 0.1069 12.8 PhiRT 0.129 0.0001 Among populations within needle types 6 4.480 0.747 0.0008 0.1 PhiPR 0.001 0.5243 Within populations 171 124.771 0.730 0.7296 87.1 PhiPT 0.128 0.0001 Total 179 142.655 0.797 0.8374 100.0 *3 regions (needle types), 9 populations, 180 trees, 28 haplotypes; Hierarchical level, source; df, degrees of freedom; SS, sum of squared observations; MS, mean of squared observations; Est. Variance, Estimated Variance; Var(%), percentage variance; Phi, value of variation; p, p - value; PhRT, proportion of the total genetic variation that are between region; PhiPR, proportion of the total variation that are among populations within a region; PhiPT, proportion of the total genetic variance that are among individuals within a population 157 Table 4.4 Pairwise PhiPT values between sympatric and distant popula tions of ponderosa pine Sympatry BLB BGZ BGM Distant GRN ROS MLB 0.000 0.208 0.000 MLA 0.215 0.221 MLZ 0.228 0.000 0.148 MLN 0.325 0.329 MLM 0.002 0.123 0.000 MLS 0.190 0.194 Average 0.077 0.110 0.049 Average 0.244 0.248 SD 0.131 0.104 0.086 SD 0.07 0.07 Twenty - one populations were evaluated for genetic distance. Shown are comparisons between populations sampled at two zones of contact (Sympatry) and populations located far from zones of contact (Distant) 158 Table 4.5. Pairwise PhiPT values be tween populations of ponderosa pine from contact zone s in the Santa Catalina Mountains BGB MLB PLB BGZ MLZ PLZ BGM MLM PLM BGB 0.000 0.510 0.380 0.010 0.010 0.010 0.380 0.490 0.070 BGB MLB 0.000 0.000 0.180 0.010 0.010 0.010 0.460 0.360 0.120 MLB PLB 0.000 0.010 0.000 0.010 0.010 0.010 0.270 0.460 0.070 PLB BGZ 0.185 0.208 0.162 0.000 0.360 0.280 0.030 0.040 0.220 BGZ MLZ 0.228 0.253 0.196 0.000 0.000 0.100 0.030 0.010 0.060 MLZ PLZ 0.165 0.158 0.143 0.003 0.030 0.000 0.070 0.060 0.380 PLZ BGM 0.003 0.000 0.020 0.105 0.148 0.067 0.000 0.450 0.410 BGM MLM 0.002 0.009 0.000 0.123 0.163 0.086 0.000 0.000 0.360 MLM PLM 0.063 0.044 0.053 0.035 0.080 0.000 0.000 0.000 0.000 PLM BGB MLB PLB BGZ MLZ PLZ BGM MLM PLM Nine populations were evaluated for genetic distance. Shown are comparisons between populations sampled for regional comparisons of needle types 159 Figure 4. 1 Black circles mark the locations of five study sites comprising 21 populations that were sampled for genetic analysis. 160 Figure 4. 1 ( continued ) These include three transition zones (where two species of ponderosa pine grow together) located on Mt. Lemmon, Palisades, and Mt. Bigelow. Highlighted are Insets mark the location of individual trees and seedlings on Mt. Lemmon and Palisades, and mature trees only on Mt. Bigelow, for three transition zones (blue oval) . These data reflect the locations of individuals from 13 populations: Mt. Bigelow (BGM, BGB, BGZ); Mt. Lemmon (MLM, ML B, MLZ, LSI1, LSI2); and Palisad es (PLM, PLB, PLZ, PSI1, PSI2). Not shown on the map are pure 3 - needle individuals from Mt. Lemmon (MLA, MLN, MLS ; red oval ); and pure 5 - needle individuals from Green Mountain (GRN ; green oval) and Rose Canyon (ROS ; green oval ) , which range in elevation from 2685 - 2780m (3 - needle) and 2140 - 2171m (5 - needle) . Locations for individuals from two additional high elevation seedling sites (3 - needle type) are also not shown ( PSU1, PSU2 ; 2542 - 2576m elevation ) . Refer to Table 4.1 for full population and site information 161 Figure 4. 2 Climate diagram for the Santa Catalina Mountains, for the 84 year reference period spanning 1925 - 2009 for (a) average annual temperature (°C; red line) and (b) total annual precipitation (mm; blue line), respectively. The average temperature increased by 1.3°C over the course of the study (from 10.6°C to 11.9°C). Site specific climate data sets (MTL, BIG) were generated using the methods of McKenney et al. (2011) and averaged together to construct this diagram. The dotted lines are the linear trend lines. See Marquardt et al. ( 2018 ) for details 162 Figure 4. 3 Convex polygons of p rincipal coordinate s analysis of 21 a priori populations of ponderosa pine collected from five sites . Axis1 explained 81.58% of the var iance, and Axis2 explained 8.53 % of the variance. PCoA analysis was based on pairwise genetic distance estimated using haplotype frequencies . The three needle types are well resolved . Refer to Table 4.1 for population and site information 163 Figure 4. 4 Histogram plot of the fi r st axis from PCoA with 21 a priori populations of ponderosa pine plotted on the X axis . Y axis is the number of populations counted in each stack . The mixed - needle types are well resolved except for one mixed - needle population and one 5 - needle population 164 APPENDIX J Taxonomy of p onderosa p ine 165 Taxonomy of p onderosa p ine P. brachyptera (Southwestern ponderosa pine) has been regarded as a species, a single taxon within P. ponderosa ( P. ponderosa var. brachyptera ), or the accepted name P. ponderosa var. scopulorum (Rocky Mountain ponderosa pine; Table 4.A1). Thus, P. ponderosa var. scopulorum refers to two major varieties, the Southwestern ponderosa pine ( P. ponderosa var. brachyptera ) and the Rocky Mountain ponderosa pine (P. ponderosa var. scopulorum) , causing taxonomic confusion because the two taxa are clearly distinct. Characte ristics that distinguish the Southwestern form are its open foliage and lack of 2 - needle fascicles, which contrasts with the typical short and stiff, 2 - needle fascicles of the Rocky Mountain form (Conkle and Critchfield 1998). In contrast to the narrow co ntact zone of coastal and interior varieties, the transition zone Basin into central Colorado (Conkel and Critchfield 1998), possibly indicating different stages of speciation (Kindler et al. 2017). We will refer to the Southwestern form as P. ponderosa var. brachyptera 16 6 Table 4.A1 Authorities for three taxa of the genus Pinus growing on Mountain Islands in southwestern USA: Sky Island 7 ; P. arizonica and P. brachyptera were obtained from The Plant List (2013) Species Infraspecific Authority Year Pub. Citation The Plant List Status Ref. 1,2,4,5 arizonica Engelm. 1879 i. accepted 4 ponderosa subsp. arizonica (Engelm.) A.E.Murray 1982 ii. synonym ponderosa var. arizonica (Engelm.) Shaw 1909 iii. synonym brachyptera Engelm. 1848 iv synonym 7 ponderosa subsp. brachyptera (Engelm.) Silba 2011 v. synonym ponderosa var. brachyptera (Engelm.) Lemmon 1888 vi. synonym ponderosa var. scopulorum Engelm. 1880 vii. accepted ponderosa Sky Island 2017 synonym 7 ponderosa Taxon X 1999 synonym 3,6 167 Table 4.A1 ( continued ) Citation refers to the place and date of original publication of name of taxa. Provided also are the accepted name for the taxonomic unit, and the synonyms, which are alternative names used to refer to the taxa that are not currently accepted. Ref. are ad ditional references which review accepted species names and synonyms (1,2,4,5), P. arizonica and P. brachyptera as species unique from P. ponderosa (5,7), and a proposed new species (3,6,7) i. Rep. U.S. Geogr. Surv., Wheeler 6(Bot):260.1878. (1879). 1. Callaham 2013 ii. Kalmia 12: 23. (1982). 2. Conkle and Critchfield 1988 iii. Publ. Arnold Arbor. 1: 24. (1909). 3. Epperson et al. 2001, 2009 iv. Wisl. Mem.Tour N Mexico [Wislizenus] 89. (1848). 4. Farjon and Styles 1997 v. J. Int. Conifer Preserv. Soc. 18: 16. (2011). 5. Lauria 1996 vi. Bienn. Rep. Calif. State Board Forest. 2: 73. (1888). 6. Rehfeldt et al. 1996, Rehfeldt 1999 vii. S.Watson [W.H.Brewer] Bot. California 2: 126 (1880). 7. Willyard et al. 2017 168 REFERENCES 169 REFERENCES Adams RP (1982) A comparison of multivariate methods for the detection of hybridization . Taxon 31: 646 - 661 Betancourt JL, Van Devender TR , Martin PS (eds) (1990) Packrat middens: the last 40,000 years of biotic change. University of Arizona Press Bezy, JV (2016) A g uide to the g eology of the Santa Catalina Mountains, Arizona: t he g eology and l ife z ones of a Madrean s ky i sland. Arizona Geological Survey Down - to - Earth # 22, 83 pp Bougeard S, Dray S (2 018) Supervised m ultiblock a nalysis in R with the ade4 p ackage J Stat Softw 86 : 1 17. https// doi .org. 10.18637/jss.v086.i01 Callaham, RZ (2013) Pinus ponderosa: a taxonomic review with five subspecies in the United States. Research Paper PSW - RP - 264. Albany, CA: USDA Forest Service, Pacific Southwest Research Station Conkle MT, Critchfield WB (1988) Genetic variation and hybridization of ponderosa pine . In: Ponderosa Pine: the species and its managemen t, Washington State University Cooperative Extension, p p . 27 - 43 Culley TM, Stamper TI, Stokes RL, Brzyski JR, Hardiman NA, Klooster MR, Merritt BJ (2013) An efficient technique for primer development and application that integrates fluorescent labeling a nd multiplex PCR . Applications in plant sciences 1:1300027 Culley TM, Weller SG, Sakai AK, Putnam KA . (2008) Characterization of microsatellite loci in the Hawaiian endemic shrub Schiedea adamantis (Caryophyllaceae) and amplification in related species and genera . Molecular Ecology Resources 8: 1081 - 1084 Dodge RA (1963) Investigations into the ecological relationships of ponderosa pine in southeast Arizona. Thesis (PhD), University of Arizona, Ari zona, 117 pp 170 Eckert geographical ranges: the central - marginal hypothesis and beyond. Mol Ecol 17:1170 1188 Epper son BK, Telewski FW, Plovanich - Jones AE Grimes JE (2001) Clinal differentiation and putative hybridization in a contact zone of Pinus ponderosa and P. arizonica (Pinaceae). Am J Bot 88: 1052 - 1057 Epperson BK, Telewsk FW , Willyard A (2009) Chloroplast diversity in a putative hybrid swarm of Ponderosae (Pinaceae ). Am J Bot 96: 707 - 712 Farjon A, Styles BT (1997) . Flora N eotropica. Pinus (Pinaceae). New York Botanical Garden Haasl RJ, Payseur BA (2011) Multi - locus inference of population structure: a comparison between single nucleotide polymporphisms and microsatellites. Heredity 106: 158 - 171 Haller, JR (1962) Variation and hybridization in p onderosa and j effrey p ines. University of California Publications in Botany 34: 123 - 65 Haller, JR (1965) The role of 2 - needle fascicles in the adaptation and evolu tion of ponderosa pine. Brittonia 17: 354 - 382 Huff DR, Peakall R, Smouse PE (1993) RAPD variation within and among natural populations of outcrossing buffalograss Buchloe dactyloides (Nutt) Engelm. Theor Appl Genet 86 : 927 - 934 Kamvar ZN, Tabima JF, Grünwald NJ. (2014) Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2:e281. https:// doi .org. 10.7717/peerj.281 Kilgore, J (2007) Distribution and ecophysiology of the Ponderosae in the Santa Catalina Mountains of southern Arizona . Thesis (PhD) , Michigan State University, 263 pp 171 Kindler C, Chèvre M, Ursenbacher S, Böhme W, Hille A, Jablonski D, Vamberger M, Fritz U (2017) Hybridization patterns in two contact zones of grass snakes reveal a new Central European snake species. Scientific Reports 7: 7378 Lauria, F (1996) The identity of Pinus ponderosa Douglas ex C. Lawson (Pinaceae) . Linzer Biologische Betraege 28: 9 99 - 1052 The Plant List (2013) Version 1.1. http://www.theplantlist.org. A ccessed 27 March 2017 Mallet, J (1995) A species definition for the modern synthesis. Trends Ecol Evolut 10: 294 - 299 Marquardt PE, Miranda BR, Jennings S, Ginger T and Telewski FW ( 2018 ) Variable c limate r esponse d ifferentiates the g rowth of s ky i sland p onderosa p ines. TREES . https://doi.org.10.1007/s00468 - 018 - 1778 - 9 Norris JR, Jackson ST , Betancourt, JL (2006) Classification tree and minimum volume ellipsoid analyses of the distribution of ponderosa pine in the western USA. J Biogeogr 33 : 342 - 360 Peakall R , Smouse PE (201 2 ) GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and re search - an update. Bioinformatics 28: 2537 - 39 Peloquin, RL (1984) The identification of three - species hybrids in the ponderosa pine complex. Southwest Nat 29: 115 - 122 Po tter KM, Hipkins VD, Mahalovich MF, Means RE (2013) Mitochondrial DNA haplotype distribution patterns in Pinus ponderosa (Pinaceae ): range - wide evolutionary history and implications for conservation. Am J Bot 100: 1562 - 1579 Potter KM, Hipkins VD, Mahalovich MF, Means RE (2015) Nuclear genetic varia tion across the range of ponderosa pine ( Pinus ponderosa ): Phylogeographic, taxonomic and conservation implications. Tree Genet Genomes 11: 1 - 23 172 R Core Team. (2016) R: A language and environment for statistical computing. R foundation for Statistical Com puting, Vienna, Austria. https://www.R - project.org/ Rehfeldt, GE (1993). Genetic variation in the Ponderosae of the southwest. Am J Bot 80 : 330 343 . https://doi.org.10.2307/2445357 Rehfeldt GE, Wilson BC, Wells SP, Jeffers RM (1996) Phytogeographic, taxonomic, and genetic implications of phenotypic variation in the Ponderosae of the Southwest. Southw Nat 41 : 409 - 418 Rehfeldt, GE (1999) Systematics and genetic structure of Ponderosae taxa (Pinaceae) inhabiting the mountain islands of the Southwest. Am J Bot 86: 741 - 752 Sheppard PR, Comrie AC, Packin GD, Angersbach K, Hughes MK (2002). The climate of the US Southwest. Clim Res 21: 219 238 . https://doi.org/10.3354/cr021219 Shinneman DJ, Means RE, Potter KM, Hipkins VD (2016) Exploring climate niches of ponderosa pine (Pinus ponderosa Douglas ex Lawson) haplotypes in the western United States: implications for evolutionary history and conservation. PloS one 11 : e0151811 Slatkin M (1987) Gene flow and the geographic structure of natural populations. Science 236: 787 - 793 Stoehr MU, Newton CH (2002) Evaluation of mating dynamics in a lodgepole pine seed orchard using chloroplast DNA markers. Can J For Res 32 : 469 - 476 Vendramin GG, Lelli L, Rossi P, Morgante M (1996) A set of primers for the amplification of 20 chloroplast microsatellites in Pinaceae. Mol Ecol 5: 595 - 598 Wickham, H (2010) g gplot2: elegant graphics for data analysis. J Stat Softw 35 : 65 - 88 Willyard A, Gernandt DS, Potter K, Hipkins V, Marquardt P, Mahalovich MF, Langer SK, Telewski FW, Cooper B, Douglas C, Finch K (2017) Pinus ponderosa : A checkered past obscured four species. Am J Bot 104: 161 - 181 173 Wofford AM, Finch K, Bigott A, Willyard A (2014) A set of plastid loci for use in multiplex fragment length genotyping for intraspecific variation in Pinus (Pinaceae). Appl Plant Sci 2: https://dx.doi.org/10.3732/apps.1400002 Wright, S (194 6) Isolation by distance under diverse systems of mating. Genetics 31: 39 Ziaco E, Truettner C, Biondi F, Bullock S (2018). Moisture driven xylogenesis in Pinus ponderosa from a Mojave Desert mountain reveals high phenological plasticity. Plant C ell E nviron 41 : 823 - 836 . https:// doi.org/10.1111/pce.13152