Provenance Variation in Early Survival, Growth, and Carbon Isotope Discrimination of Southwestern Ponderosa Pine Growing in Three Common Gardens across an Elevational Gradient

Abstract

We investigated early survival, growth, and carbon isotope discrimination of ponderosa pine (Pinus ponderosa Lawson & C. Lawson var. scopulorum Engelm.) seedlings from different provenances using common gardens across an elevational gradient in order to examine the potential for adaptation to extreme environments and constraints to artificial regeneration. Twenty-one provenances from a range of elevations across Arizona and New Mexico were planted in three common gardens: a high-elevation meadow in aspen-mixed conifer forest, a mid-elevation ponderosa pine forest, and a low-elevation pinyon juniper woodland. Two years after planting in 2018, survival was highest at the mid-elevation site (54%), low at the high-elevation site (1.5%), and 0% at the low-elevation site. At the hot and dry low-elevation site, provenances from low-elevations survived longer than provenances from mid- and high-elevations, which suggests greater drought tolerance of low-elevation provenances. Mortality agents changed from abiotic (drought) to biotic (herbivory) with an increase in elevation across sites. High mortality of seedlings planted at high-elevation sites from biotic agents, such as rodents, may challenge efforts to establish ponderosa pine in assisted migration projects. Seedlings had significantly higher growth rate and carbon isotope discrimination (∆¹³C) at the mid-elevation site than the high-elevation site. Provenances differed significantly in diameter, and ∆¹³C, but not in height growth rate for the first year after planting. Provenance variation in ∆¹³C suggests genetic variation in water use efficiency that may be useful for future evaluation of southwestern ponderosa pine seed sources for reforestation.

1. Introduction

Climate warming and associated high-severity fires and extended droughts have caused large-scale tree mortality in the southwestern United States [1,2]. In southwestern ponderosa pine (Pinus ponderosa Lawson & C. Lawson var. scopulorum Engelm.) forests, severe wildfires and droughts have also reduced natural regeneration [3,4]. A continuation of high tree mortality and meager regeneration threatens to shrink the geographical distribution of these forests [5]. Range contraction of southwestern ponderosa pine forests after severe drought has already been observed [6,7]. Climatic constraints on natural regeneration of ponderosa pine by drought and high temperature are expected to increase in the future [8].

In the absence of sufficient natural regeneration, it is important to explore the use of tree planting as a possible tool to compensate for recent deforestation [9,10]. Future plantings should include seedlings that are genetically and physiologically pre-adapted to future arid conditions through enhancing drought tolerance [11,12]. Non-local populations may have more pre-adaptation to the future climate of a site than local populations because populations that currently occur at a site are adapted to past local conditions [13,14] that, in the American Southwest, were cooler and wetter [15]. Assisted migration of populations has been proposed to mitigate climate change impacts on forests by preemptively moving drought tolerant populations to sites predicted to increase in aridity [11]. In addition to climate, assisted migration strategies and tree planting in general should consider constraints to success of planted seedlings from biotic agents [16,17], that have been previously reported to negatively impact ponderosa pine seedlings [18,19,20].

Previous common garden investigations of ponderosa pine in the southwestern United States set the stage for our work by showing clinal patterns of genetic differences among provenances in seedling traits, such as growth rate, spring budburst timing, and frost tolerance [21,22]. These previous investigations, however, did not specifically include low-elevation provenances [21], which may be arid adapted [23], did not cover a wide environmental and geographic range [24], or only studied seedlings growing under resource-rich nursery or greenhouse conditions [22]. More studies of ponderosa pine provenance performance in field conditions over environmental gradients are needed in order to understand impacts of abiotic and biotic factors on planted seedlings and to refine seed transfer guidelines based on a changing climate.

Our study builds on earlier investigations of provenance variation in ponderosa pine growth rate, drought tolerance and carbon isotope discrimination by specifically focusing on fine-grained spatial variation in the southwestern United States. Kerr et al. [25] reported a tradeoff between seedling growth rate and drought tolerance for ponderosa pine provenances from Oregon. Carbon isotope discrimination (∆¹³C) is a proxy for water use efficiency, a potentially important mechanism of drought tolerance [26]. Carbon isotope discrimination is determined by leaf internal CO₂ concentration, which is controlled by the ratio of net photosynthesis to stomatal conductance [27]. Carbon isotope discrimination differed among widely separated ponderosa pine provenances from the interior United States [28,29], but these previous studies did not specifically focus on southwestern or low-elevation provenances within a region.

In this study, we compared early survival, growth, and ∆¹³C among seedlings from 21 ponderosa pine provenances from Arizona and New Mexico at three planting sites that varied in elevation, temperature, and precipitation. The planting sites were selected to bracket the current climatic niche of ponderosa pine in the southwestern U.S. by including a site in the middle of the current species distribution; a hot, low-elevation site with a temperature expected to occur within the core range of the species in the next 100 years [15], and a cool high-elevation site where assisted migration of the species could occur in the future. We address the following questions: (1) How does common garden environment, provenance, and their interaction influence early survival, growth, and carbon isotope discrimination of planted ponderosa pine seedlings? (2) How does seedling performance relate to provenance environmental characteristics?

2. Materials and Methods

2.1. Common Garden Sites

The study was conducted in three common gardens across an elevational gradient. The common gardens were planted in the core of the species range (mid-elevation site; elevation 2200 m, Mean Annual Temperature (MAT) 7.6 °C) and outside the current range at high-elevation (elevation 2780 m, MAT 4.9 °C) and low-elevation (elevation 1930 m, MAT 9.7 °C) sites (Figure 1). Mean annual temperature of the low-elevation site is similar to the projected temperature of sites in the core of the species range by the end of the 21st century. Each site was equipped with a weather station and a datalogger (Campbell Scientific, Logan, UT, USA) to measure air temperature and precipitation. The planted areas at all sites were protected against ungulate herbivory by fences that prevented access by large mammals such as cattle, deer, and elk, but not small animals.

The low-elevation site was planted on 21 July 2018, in an open area within a pinyon–juniper woodland in northern Arizona. The soil at the site is a clay loam derived from basalt and limestone parent material (https://sega.nau.edu/gardens/blue-chute, accessed on 19 September 2020). The year of planting (2018) was 1.4 °C warmer and received only 33 percent of average annual precipitation compared to 30-year normals (1981–2010; https://prism.oregonstate.edu, accessed on 19 September 2020). The mid-elevation site was planted on 20 July 2018, in an open grass-dominated meadow in a ponderosa pine forest in northern Arizona. The soil type at this site is clay loam derived from basalt parent material (https://sega.nau.edu/gardens/arboretum-forest, accessed on 19 September 2020). At this site, the planting year was 0.2 °C warmer and received only 66 percent of mean annual precipitation compared to 30-year normals. Conditions at this site were particularly dry in the year 2020 with only 33 percent of mean annual precipitation. The high-elevation site was planted on 24 July 2018, in a grassy meadow within an aspen and mixed conifer forest in southern Utah. The soil type at this site is a silty loam [30] derived from volcanic parent material. The planting year was 1.5 °C warmer and received about 85 percent of mean annual precipitation compared to the 30-year mean. Conditions at this site were particularly dry in the year 2020 with about 40 percent of 30-year mean annual precipitation. Overall, the study period (2018–2020) was drier compared to the long-term average at all three sites (Supplementary Figures S1 and S2). During this period, temperatures were warmer at all sites except for the year 2019 when the mean annual temperatures were 0.8 °C and 0.6 °C lower than long term average at the mid- and high-elevation sites, respectively (Supplementary Figure S1). The role of mycorrhizae specific to ponderosa pine on early seedling survival and performance could not be addressed in our study.

2.2. Experimental Design and Plant Material

At each of the three sites, 1008 seedlings were planted in a randomized complete block design. At each site, seedlings were planted in 4 blocks, each consisting of one 12-tree linear plot for each provenance. The arrangement of provenance plots was randomized within each block. Spacing among seedlings was 1.2 m. Each of the 21 provenances was represented by 48 seedlings from between 3 and 6 mother trees or a pooled collection of seeds in cases where mother tree level information was not available (5 provenances). We planted a buffer row of ponderosa pine seedlings spaced at 1.2 m outside of the experimental trees to reduce potential edge effects. Existing vegetation within the fenced exclosures at the mid-and high-elevation sites was reduced using herbicide (Ranger Pro, Bayer Crop Science, Chesterfield, MO, USA) prior to planting. Herbicide was not needed at the low-elevation site due to sparse grass cover.

To produce the seedlings, seeds were sown in March 2018 in 164 mL containers (Ray Leach Cone-tainers SC10 Super, Stuewe & Sons, Inc., Tangent, OR, USA) at the John T. Harrington Forestry Research Center with New Mexico State University located in Mora, New Mexico. Seedlings were grown in the greenhouse nursery under standard operational protocols for approximately four months until planting in July 2018. Detailed greenhouse growing conditions, and seed and provenance information are described in [31].

The provenances were located over a gradient of elevation, temperature, and precipitation throughout Arizona and New Mexico (

Table 1. Name (code), latitude, longitude, elevation (m a.s.l.), mean annual temperature (MAT), and mean annual precipitation (MAP) for provenances and common garden sites, ordered by increasing elevation. Climate data (30 year normal, 1981–2010) are from PRISM (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 19 September 2020). Values are means over mother tree locations at each provenance.

	Latitude/Longitude	Elevation (m)	MAT (°C)	MAP (mm)
Provenance (code)
Cherry Road (CR)	34.586/–112.057	1592	14.2	408
Blue River (BR)	33.555/–109.193	1674	11.8	581
Mesa Del Medio (MDM)	35.116/–105.217	1714	12.3	397
Prescott—Iron Springs Road (PIS)	34.585/–112.559	1846	11.0	546
Townsend Winona (WIN)	35.254/–111.415	1934	10.2	398
Hualapai Mountains (HM)	35.084/–113.875	1969	11.5	403
Ruidoso Service Office (RSO)	33.350/–105.583	1976	11.1	506
Sapello Rt. 3 (SAP)	35.700/–105.250	2050	9.8	456
South Kaibab Tusayan Dist. (SKT)	35.939/–112.084	2067	8.9	410
Northern Arizona University (NAU)	35.182/–111.655	2104	8.4	540
Rim District (RD)	34.487/–111.343	2244	9.4	767
HH Ranch (HHR)	34.183/–107.525	2270	9.2	364
Mineral Hill (MH)	35.633/–105.461	2277	8.5	528
Mud Springs (MUD)	36.463/–106.859	2277	6.9	443
Manzano Mountains (MZN)	34.623/–106.400	2366	8.6	638
Vallecitos-Jemez Springs (VJS)	35.809/–106.589	2436	6.8	571
Hartman Ridge (HR)	35.550/–105.533	2500	8.8	526
Borrego Mesa (BM)	35.990/–105.794	2560	6.3	470
Magdalena Mountains (MAG)	34.006/–107.215	2565	8.9	512
Green’s Peak (GP)	34.126/–109.535	2760	5.7	671
Mount Taylor (TAY)	35.266/–107.633	2814	6.1	722
Common Garden Site
Low-elevation	35.586/–111.969	1930	9.7	482
Mid-elevation	35.160/–111.730	2200	7.6	556
High-elevation	37.620/–113.025	2780	4.9	685

). Provenances ranged from ~1600 m to ~2800 m in elevation, 364 mm to 767 mm mean annual precipitation, and 5.7 °C to 14.2 °C mean annual temperature. Provenance elevation had a strong negative correlation with mean annual temperature (MAT) (r = −0.902, p < 0.0001) and a moderate positive correlation with mean annual precipitation (MAP) (r = 0.500, p = 0.020); provenance MAT and MAP had a moderate negative correlation (r = −0.434, p = 0.048). Latitude and MAT had a moderate negative trend (r = −0.406, p = 0.067) and longitude and elevation had a moderate positive trend (r = 0.370, p = 0.098).

2.3. Survival and Growth Measurements

At the mid- and high-elevation sites, survival was measured once every fall and summer between 2018 and 2020. Height and diameter at the soil surface (ground line diameter) were measured every fall (October or November) between 2018 and 2020. At the low-elevation site, survival was assessed once per week between August and September 2018. We did not measure growth at the low-elevation site due to complete lack of seedling survival by the end of fall 2018. Seedlings were classified as dead when they were 100 percent brown or were removed from planting locations by herbivores. We also observed and recorded potential causes of mortality and signs of mortality agents such as scats, gopher holes, and ant nesting.

2.4. Carbon Isotopic Discrimination

In November 2019 at the end of the second full growing season after planting, we collected current year true needles from between 7 and 13 seedlings from each provenance at each of the mid- and high-elevation sites. At each site, seedlings were selected for sampling from each provenance–block combination (row plot) using stratified random sampling. We used mother trees from each row plot as strata to ensure balanced sampling over all surviving mother trees. Provenance rows in some of the blocks had no survival, and thus, samples were not taken. For the 5 provenances without mother tree information, we randomly sampled 3 seedlings from each row plot per block. The needles were oven-dried (Sheldon Manufacturing, INC, Cronelius, OR) at 65 °C for 72 h and ground to homogenous powder using a ball mill grinder (Mixer Mill MM200, Retsch, Haan, Germany). The samples were analyzed for carbon stable isotopic composition at the Colorado Plateau Stable Isotope Laboratory at Northern Arizona University, Flagstaff, Arizona using a DELTA V Advantage isotope ratio mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA), which is configured through a Finnigan ConFlo III for automated continuous-flow analysis of δ¹³C using a Carlo-Erba NC2100 elemental analyzer for combustion and separation of carbon and nitrogen. Carbon isotope discrimination (∆¹³C) was calculated using the following equation [32] where δ¹³C is the isotopic ratio of ¹³C to ¹²C in needle tissues, and −0.008 is the approximate δ¹³C of atmospheric CO₂, compared with the Pee Dee Belemnite standard.

$${\Delta }^{13}\mathrm{C}=\frac{-0.008+\left[\left({\mathsf{\delta }}^{13}\mathrm{C}\right)/\left(-1000\right)\right]}{1-\left[\left({\mathsf{\delta }}^{13}\mathrm{C}\right)/\left(-1000\right)\right]}\times 1000$$

2.5. Data Analyses

For seedling growth, we used a mixed-effects model with provenance, common garden site, and the interaction between provenance and site as fixed effects and blocks nested within sites as a random effect. We conducted the analyses on block-level provenance means and regarded the row plot of seedlings from same provenance within a block as the experimental unit. We used post hoc Tukey’s Honestly Significant Difference (HSD) comparison of means to evaluate significant common garden site and provenance differences (α = 0.05). Provenance effect on growth was not tested for the year 2020 due to very low survival at the high-elevation site at the end of the year. We used relative height growth to assess growth since planting calculated as: (Fall height—Height at planting)/Height at planting. We used survival analysis to assess the impact of common garden site on seedling survival between July 2018 and October 2020. All data approximated a normal distribution. The relationships between provenance traits and provenance environmental characteristics were evaluated on provenance means with correlation and regression analyses. We interpreted the strength of relationship based on the value of the correlation coefficient (r) as: weak between 0 and 0.3, moderate between 0.3 and 0.7, and strong between 0.7 and 1.0 [33]. JMP Pro version 14 (SAS Institute Inc., Cary, NC, USA) was used to perform all analyses.

3. Results

3.1. Survival Analyses

Early seedling survival differed significantly (p < 0.0001) among the three common garden sites (Figure 2). By fall 2020, survival at the mid-elevation site was much higher than the other two sites; 54 percent of the seedlings survived at the mid-elevation site compared to 1.5 and 0 percent at the high and low-elevation sites, respectively. Survival duration since outplanting at the low-elevation site differed significantly among provenances (p = 0.0006) and had a moderate negative correlation with provenance elevation (r = −0.43; p = 0.049; Figure 3). Most seedlings (80 percent) at the low-elevation site died in the first month after planting, and survival was 0 percent at the end of fall 2018. At the mid- and high-elevation sites, survival at the end of 2020 did not differ significantly among provenances (p = 0.600 and 0.058, respectively). Two major mortality events occurred at the mid-elevation site; summer 2019 (30 percent mortality) and fall 2020 (10 percent mortality). At the high-elevation site, survival decreased continuously between 2018 and 2020, with the largest decrease (40 percent) in summer 2020.

Mortality agents changed from abiotic to biotic factors with an increase in elevation over the site elevational gradient. At the low-elevation site, 99.7 percent of seedling mortality was attributed to desiccation (visually observed dry seedlings with intact roots). The only exception was herbivory-induced mortality of three seedlings that were planted close to ant colonies. At the mid-elevation site, important mortality agents were rabbit herbivory (~30 percent seedlings; obvious scat and feeding patterns) and desiccation (~10 percent seedlings). Rabbit herbivory occurred over an 8-week period in the first spring after planting (2019). At the high-elevation site, most seedling mortality (90 percent) was attributed to pocket gopher (Thomomys bottae) herbivory as evidenced by gopher holes, complete absence of seedlings, and seedlings with missing roots, with a smaller amount (10 percent) attributed to desiccation.

3.2. Growth

Provenances did not vary significantly in relative height growth between planting and fall 2019 (

Table 2. Provenance relative height growth ((Fall 2019 height—Height at planting)/Height at planting) and ground line diameter (GLD) at mid- and high-elevation sites (ordered by increasing elevation). Values are means (1 SE); NA indicates SE not calculated due to low survival in most blocks. Means in the same column followed by the same letters do not differ significantly (Tukey’s HSD tests; α = 0.05).

Provenance	Fall 2019 Relative Height Growth (mm/mm)	Fall 2019 GLD (mm)
Mid-elevation site
CR	1.0 (0.1) a	7.66 (0.41) a
BR	1.3 (0.2) a	7.96 (0.45) a
MDM	1.2 (0.0) a	6.24 (0.06) a
PIS	1.4 (0.0) a	7.27 (0.23) a
WIN	1.3 (0.0) a	7.53 (0.45) a
HM	1.5 (0.1) a	7.00 (0.31) a
RSO	1.4 (0.2) a	7.38 (0.33) a
SAP	1.4 (0.1) a	7.21 (0.21) a
SKT	1.0 (0.0) a	7.80 (0.23) a
NAU	1.1 (0.1) a	7.75 (0.34) a
RD	1.6 (0.1) a	8.46 (0.45) a
HHR	1.2 (0.0) a	7.49 (0.20) a
MH	1.5 (0.1) a	7.55 (0.46) a
MUD	1.4 (0.2) a	6.50 (0.25) a
MZN	1.2 (0.2) a	7.86 (0.36) a
VJS	1.3 (0.1) a	7.90 (0.24) a
HR	1.5 (0.2) a	7.80 (0.22) a
BM	1.6 (0.2) a	7.41 (0.32) a
MAG	1.1 (0.1) a	6.40 (0.52) a
GP	1.7 (0.3) a	7.45 (0.58) a
TAY	1.7 (0.1) a	7.19 (0.19) a
High-elevation site
CR	0.08 (NA) a	5.15 (NA) ab
BR	0.3 (0.0) a	5.32 (0.01) ab
MDM	0.3 (0.1) a	4.46 (0.20) ab
PIS	0.4 (0.0) a	5.39 (0.26) ab
WIN	0.5 (0.0) a	5.58 (0.59) ab
HM	0.3 (0.1) a	5.10 (0.09) ab
RSO	0.5 (0.2) a	6.01 (0.07) a
SAP	0.4 (0.0) a	5.05 (0.84) ab
SKT	0.4 (NA) a	5.32 (NA) ab
NAU	0.3 (0.0) a	5.75 (0.22) ab
RD	0.2 (NA) a	6.42 (NA) a
HHR	0.3 (NA) a	3.80 (NA) ab
MH	0.3 (NA) a	4.58 (NA) ab
MUD	0.2 (NA) a	3.90 (NA) ab
MZN	0.3 (0.0) a	5.70 (0.12) ab
VJS	0.6 (NA) a	4.72 (NA) ab
HR	0.2 (0.0) a	4.14 (0.56) b
BM	0.4 (0.0) a	4.50 (0.05) ab
MAG	0.0 (NA) a	5.80 (NA) ab
GP	0.6 (NA) a	6.45 (NA) a
TAY	0.3 (0.2) a	5.60 (0.40) ab

and

Table 3. p-value of seedling relative height growth, ground line diameter (GLD), and carbon isotope discrimination (∆¹³C) for the provenance effect, planting site effect, and provenance × site interaction. Bolded values significant at p < 0.05.

	Provenance	Site	Provenance × Site
2019 relative height growth	0.926	<0.0001	0.982
2019 GLD	0.008	0.001	0.143
∆13C	0.006	0.0001	0.239

). However, relative height growth differed among planting sites for years 2019 and 2020 and was significantly higher at the mid-elevation site than the high-elevation site for both years (Figure 4). The provenance effect on growth was not tested for the year 2020 because only 1.5 percent of seedlings survived at the high-elevation site during that year. Site differences were consistent over provenances as indicated by the non-significant provenance x site interaction (Table 3).

Correlation analysis provided evidence of relationships between provenance environmental characteristics and site-specific growth traits (

Table 4. Correlation coefficients between provenance environmental characteristics and trait means for mid and high elevation sites. Values in parentheses are p values. Boldface type indicates significance (p < 0.05). Ground line diameter (GLD); carbon isotope discrimination (∆¹³C).

	Latitude	Longitude	Elevation	MAT	MAP
Mid-elevation site
2019 relative height growth (mm/mm)	0.09 (0.668)	0.13 (0.561)	0.56 (0.007)	−0.51 (0.015)	0.54 (0.010)
2019 GLD (mm)	−0.12 (0.586)	−0.26 (0.244)	0.001 (0.993)	−0.04 (0.843)	0.43 (0.051)
∆13C (‰)	0.32 (0.154)	−0.05 (0.806)	−0.001 (0.996)	0.03 (0.864)	−0.09 (0.683)
High-elevation site
2019 relative height growth (mm/mm)	0.09 (0.708)	0.007 (0.973)	0.10 (0.657)	−0.30 (0.175)	0.04 (0.850)
2019 GLD (mm)	−0.53 (0.014)	−0.40 (0.066)	0.06 (0.789)	0.01 (0.960)	0.60 (0.004)
∆13C (‰)	−0.42 (0.054)	0.01 (0.942)	0.02 (0.936)	0.09 (0.699)	0.31 (0.161)

). At the mid-elevation site, 2019 relative height growth had a positive correlation with provenance elevation (r = 0.56; p = 0.007; Supplementary Figure S3) and provenance MAP (r = 0.54; p = 0.010), and a negative correlation with MAT (r = −0.51; p = 0.015). At the high-elevation site, 2019 GLD was correlated negatively with provenance latitude (r = −0.53; p = 0.014) and positively with provenance MAP (r = 0.60; p = 0.004).

3.3. Carbon Isotope Discrimination

Provenances (p = 0.006; Figure 5) and common garden sites (p = 0.0001; Figure 6) varied significantly in ∆¹³C, whereas the provenance x site interaction was not significant (p = 0.239) (Table 3). ∆¹³C was about 1‰ higher in seedlings growing at the mid-elevation site compared with the high-elevation site (Figure 6). ∆¹³C ranged among provenances from 18.4 ‰ for SKT to 20.2 ‰ for PIS (Figure 5). Correlations between mean ∆¹³C and provenance environmental characteristics were weak and non-significant at the mid-elevation site. However, at the high-elevation site, provenance latitude had a moderate negative trend with ∆¹³C (r = −0.42; p = 0.054; Table 4).

4. Discussion

Our results show that early seedling survival was highest at the mid-elevation site in the core of the species current range and decreased substantially at sites outside the current range with no survival at the low-elevation site and almost none at the high-elevation site. Our finding of differences in survival among sites is similar to a recent report of significant site effects for survival of Jeffrey (Pinus jeffreyi Grev. & Balf.) and ponderosa pine (Pinus ponderosa Lawson & C. Lawson) seedlings planted at three different elevations in California [34]. In our study, the lowest survival occurred at the low-elevation site where almost all seedlings died within the first two months after planting due to desiccation. The ~2 °C warmer MAT of the low-elevation site than the mid-elevation site is an example of the future warmer climate in the current range of ponderosa pine in the southwestern United States by the end of this century [15]. Our result of 0 percent survival at this low-elevation site implies future challenges in successful planting of ponderosa pine within the core of its current geographical range due to high temperature and water stress. In the future, reforestation via planting might need to be timed with particularly wet years, if possible.

A general decline in survival duration with an increase in provenance elevation at the hot, low-elevation site (Figure 3) suggest better adaptation of low-elevation provenances to warmer and drier conditions. This result is consistent with a previous report of longer survival of ponderosa pine seedlings from provenances from arid locations under experimental drought in the greenhouse [24]. While we found a negative correlation between survival duration and provenance elevation, under extreme field conditions, it is important to note that no provenances successfully established at the low-elevation site. Ponderosa pine provenances from low-elevation and drier locations have been shown to have traits related to drought tolerance such as lower specific leaf area and higher growth allocation to roots [22,24]. These traits are one of the reasons that lower elevation or southern sources have been recommended for future planting to mitigate impacts of climate warming [11,12]. However, mechanisms of drought tolerance in such provenances are not fully understood and should be investigated in future studies.

At the mid- and high-elevation sites, the most important source of seedling mortality was herbivory followed by desiccation. Herbivory and drought have been recognized as major challenges in the outplanting success of nursery grown seedlings for reforestation [10,35]. Since the mid-elevation site was located in a ponderosa pine-dominated forest and in the core of current range in northern Arizona, high survival of planted seedlings was expected. However, at this site, ~30 percent of seedlings died due to rabbit herbivory. At the high-elevation site, gophers were the most important mortality agent, accounting for ~90 percent seedling mortality. Our findings of seedling mortality due to rabbits and gophers is consistent with ponderosa pine mortality agents previously reported [19]. However, the use of fences at our sites to prevent ungulate herbivory may have indirectly increased herbivory by small mammals due to restricted access of the fenced areas by predators [36]. As of fall 2020, about 10 percent of seedlings died due to desiccation at the mid- and high-elevation sites each. This could be attributed to the dry conditions in summer 2020 at both sites (Supplementary Figures S1 and S2). Overall, our results of high seedling mortality due to herbivory suggests severe biotic constraints to planting success at some high-elevation sites [37]. Future research should investigate the effectiveness of protective measures against herbivory such as tree shelters or tubes and if herbivory patterns can be spatially predicted [38,39,40].

Relative height growth differed significantly across sites with seedlings at the mid-elevation site growing faster than seedlings at high-elevation site (Figure 4). Seedling growth at the high-elevation site may have been limited by a shorter frost-free season and cooler temperatures compared with the mid-elevation site. At the mid-elevation site, provenances from high-elevation, cooler, and wetter areas had a higher relative height growth rate than low-elevation, warmer, and drier provenances. This result shows that high-elevation provenances can take advantage of warmer conditions, which in our study were created by planting these high-elevation provenances at the mid-elevation common garden site.

Results for leaf ∆¹³C suggest an increase in water use efficiency with increasing elevation of the planting site. Our results for ∆¹³C are consistent with Hultine and Marshall [41], which suggested that a decrease in ∆¹³C with increasing elevation may result from a decrease in stomatal conductance with elevation due to lower temperature. Our findings of a decrease in ∆¹³C with elevation of the planting site are similar to an earlier study of ten provenances of ponderosa pine growing at three sites in Nebraska [42]. We also found a significant provenance effect on ∆¹³C, which was mainly due to provenance SKT, which had the lowest ∆¹³C among all provenances across both sites; this provenance needs to be further investigated as it might be useful for future reforestation due to its high water use efficiency. Provenances from more southern latitudes had higher ∆¹³C and by inference lower water use efficiency than northern provenances at the high-elevation site. This latitudinal pattern is similar to a greenhouse study involving 21 provenances of Populus balsamifera L. [43], but differs from a common garden study of Pinus strobiformis seedlings [44]. This latitudinal pattern could result from local adaptation of provenances to abiotic or biotic factors that vary with geographic location, or perhaps are related to genetic variation among putative subspecies of ponderosa pine [45].

In summary, our results provide evidence of site and provenance effects on early survival and performance of planted southwestern ponderosa pine seedlings growing at different elevations. Our finding of longer survival of low-elevation provenances under dry, hot field conditions support previous recommendations for planting ponderosa pine provenances from lower elevations in a warming climate [12]. In addition to taking site environment and provenance origin into consideration, our study highlights the importance of accounting for mortality from biotic agents in reforestation projects [10], especially at high-elevation meadow sites [37]. Additional studies should focus on predicting the influence of biotic agents over a longer period of time and the interactions between abiotic and biotic stresses on seedling survival and performance [46]. Overall, our findings enhance understanding of constraints on the success of planted seedlings and thus inform strategies to maintain ponderosa pine during changing climate by active reforestation.