Cytotype and local adaptation drive phenotypic variation in two subspecies of big sagebrush (Artemisia tridentata)

Abstract

Big sagebrush (Artemisia tridentata) is a widespread and locally dominant shrub throughout many ecosystems in western North America. A. tridentata ssps. tridentata and wyomingensis are two subspecies whose populations occupy the warm‐arid regions of the species range and whose trailing edge is threatened by climate change. Previous studies have presented conflicting results in relation to the genetic control of physiological variation in A. tridentata. Understanding how different genetic factors contribute to physiological variation can provide insight into how these two subspecies may respond to future climate change. To explore possible variation among and within two subspecies of A. tridentata, we measured physiological and morphological traits in A. t. tridentata and A. t. wyomingensis during mid‐summer (July), seven years after establishment in a common garden. Contributions to trait variation were quantified for both genetic (subspecies and cytotype) and environmental (climate‐of‐origin) factors. Measurements revealed an unequal contribution to phenotypic variation by subspecies, cytotype, and climate‐of‐origin. Ploidy and climate‐of‐origin were more important than subspecies in driving phenotypic variation in A. tridentata. These findings suggest that A. tridentata has a highly plastic drought response, or that culling (mortality over time due to environmental factors) in the common garden over seven years has led to a lack of genetic diversity within the garden. Understanding what factors drive phenotypic expression in big sagebrush can provide better insight into how climate change may affect migration and extirpation and may aid in the effectiveness of restoration efforts.

INTRODUCTION

Recent climate models predict geographic differences in precipitation regime change across the Intermountain West for the mid-21st century (Abatzoglou & Kolden, 2011). Specifically, they suggest that along with a significant warming, the southwest region will experience a significant drop in summer and winter precipitation, while the northern region of the Great Basin will experience an increase in winter precipitation, but with a larger percentage occurring as rain instead of snow (Abatzoglou & Kolden, 2011). Changes in climatic regimes such as these can have a profound impact on both populations of individual species and entire ecosystems. Artemisia tridentata (big sagebrush) is a locally dominant species of the Intermountain West and is one of the most widespread shrubs in western North America (Beetle, 1960; Connelly et al., Unpublished report; Knick et al., 2003; McArthur et al., 1978). It is a xerophytic, drought-deciduous shrub whose productivity is mostly regulated by wintertime recharge of soil moisture by snow and rain (Beetle, 1960; McArthur et al., 1978; Miller & Shultz, 1987; Schlaepfer et al., 2012). A. tridentata can be found in a variety of environments, ranging from relatively mesic mountainsides to xeric bajadas (McArthur et al., 1978), where it serves as an integral part of the ecosystem that it occupies (Bechtold & Inouye, 2007; Connelly et al., Unpublished report; Richards & Caldwell, 1987; Welch & McArthur, 1979). It drives carbon, water, and nitrogen cycling in its communities (Bechtold & Inouye, 2007; McAbee et al., 2017; Richards & Caldwell, 1987) and is high in crude protein compared with other Great Basin woody shrubs, and thus provides nutritious forage to species that find it palatable, including mule deer (Odocoileus hemionus) and the greater sage grouse (Centrocercus urophasianus; Connelly et al., Unpublished report; Welch & McArthur, 1979). Therefore, any future negative impacts to A. tridentata populations driven by climate change could result in large-scale consequences to Great Basin ecosystems.

There are three recognized site-specific subspecies of A. tridentata, two of which exist at the driest regions of the species' range: A. tridentata ssp. tridentata (basin big sagebrush) and A. tridentata ssp. wyomingensis (Wyoming big sagebrush). Both A. t. tridentata and A. t. wyomingensis can be found in foothills and valley floors, but A. t. tridentata occurs across a wider range of elevations (610–2140 m) than A. t. wyomingensis (1530–2150 m; McArthur et al., 1978). However, despite being found in similar environments and elevations, the two subspecies appear to exist on a soil-moisture gradient (Bonham et al., 1991; Shumar & Anderson, 1986), and often grow in contiguous populations with low hybridization (Barker & McKell, 1983; McArthur et al., 1988). In addition, the two subspecies vary morphologically. A. tridentata ssp. tridentata is much larger and denser than A. t. wyomingensis, with longer and thinner leaves, while A. t. wyomingensis exhibits more dieback and flowers that subtend the crown of the shrub (Beetle, 1960; McArthur et al., 1978; Rosentreter, 2005).

In addition to habitat and morphological differences, these two subspecies vary physiologically. Several papers have reported differences in survivorship and physiology between mature A. t. tridentata and A. t. wyomingensis, both in and outside their native habitats (Brabec et al., 2017; Frank et al., 1984; Kolb & Sperry, 1999; Lazarus et al., 2019; McArthur & Welch, 1982). A common garden study performed by Frank et al. (1984) reported higher water use efficiency and growth rates in A. t. tridentata compared with A. t. wyomingensis. Kolb and Sperry (1999) reported Ψ₅₀ (the amount of tension required for xylem to lose 50% of its water transport efficiency) varied significantly between A. t. wyomingensis (Ψ₅₀ = −4.9 MPa) and A. t. tridentata (Ψ₅₀ = −3.9 MPa). These differences were observed both in native populations and in common garden shrubs, suggesting an underlying genetic basis. However, previous studies have not considered the effect of other genetic variables when considering phenotypic expression in A. tridentata subspecies.

Genetic variation can strongly affect phenotypic expression. In widespread plant species such as A. tridentata, it is common to observe populations occupying a range of environmental conditions and possessing a divergent suite of phenotypic traits resulting from habitat-specific selection (Chaney et al., 2017; Richardson & Chaney, 2018). More recent studies considering genetic variables such as autoploidy and local adaptation have found these factors to be significantly correlated with survivorship, phenology, and physiology of A. tridentata shrubs (Chaney et al., 2017; Germino et al., 2019; Lazarus et al., 2019; Richardson et al., 2021). A few common garden studies report that survival (Chaney et al., 2017) and freeze resistance (Lazarus et al., 2019) in A. tridentata shrubs were correlated with parental source populations, and thus variation in traits could be explained by the climate of origin. In addition, it is generally held that populations with larger cytotypes possess wider ecological distributions, potentially due to a higher capacity for plasticity and local adaptation (Münzbergová, 2007). Evidence of this phenomenon can be assessed in A. tridentata, which is an autoploidy species, with diploid and tetraploid populations of A. t. tridentata and tetraploid populations of A. t. wyomingensis (McArthur et al., 1981). Recently, Richardson et al. (2021) found that the cytotype of A. tridentata was a strong predictor for growth and seed yield. Specifically, 15% and 35% of variation in growth and seed yield were attributed to ploidy, respectively, with diploid individuals exhibiting higher growth and seed yield than tetraploid shrubs, regardless of subspecies or environment (Richardson et al., 2021). These results suggest that genome size may play an important role in the physiological processes of A. tridentata and highlight potential gaps in the A. tridentata physiological literature.

Despite the vast literature highlighting A. tridentata subspecies differentiation, little is still known about whether ecophysiological traits underlie those differences. Although some recent work (Chaney et al., 2017; Richardson et al., 2021) has called into question consistent genetic differences in A. tridentata subspecies, we consider subspecies designations (based on phenotype) to be an encompassing genetic component here. Therefore, understanding how traits relate to genetic variables (e.g., cytotype) could provide insight into how this shrub might respond to future climate change. In this study, we explore the contribution of multiple genetic variables to variation of ecophysiological traits in A. t. ssp. tridentata and A. t. ssp. wyomingensis. We examined a wide range of physiological traits relating to morphology, gas exchange, chlorophyll fluorescence, and water relations (hereafter referred to as traits) to provide physiological snapshots. We utilized an A. tridentata common garden that housed individuals from the two subspecies, grown from numerous parent populations across A. tridentata's range. This allowed quantification of morphological and physiological traits (e.g., gas exchange, chlorophyll fluorescence, and water relations) while controlling for environmental effects. Thus, this setting allowed for isolation of subspecies, cytotype, and local adaptation (climate-of-origin) effects on A. tridentata's physiology. Because subspecies, cytotype, and local adaptation can produce climate-specific stress responses, we hypothesized that all three genetic factors would influence morphological and physiological responses. However, we predicted that, based on recent findings (Chaney et al., 2017; Germino et al., 2019; Lazarus et al., 2019), subspecies classification would have a relatively small effect on variation in these measured traits, and that climate-of-origin (local adaptation) and cytotype would be the main drivers for phenotypic variation in mature A. tridentata.

METHODS

Study location and design

This study was conducted at the A. tridentata “Orchard” common garden site, established in April of 2010, at a location east of Boise, Idaho (43.32° N, −115.9° W, 974 m). The common garden site has been characterized as typical of an A. t. wyomingensis native environment: warm intermountain steppe with a dry climate, and was confirmed by the identification of shrubs from surrounding populations as A. t. wyomingensis (Chaney et al., 2017). Shrubs within the common garden originated from 55 different populations of A. tridentata from 11 different states in the western United States (see fig. 3 in Chaney et al., 2017), occupying a wide breadth of climates (see Appendix S1: Table S1). Subspecies and cytotypes were determined through a combination of morphology, UV fluorescence, flow cytometry, and genetic markers (Richardson et al., 2012, 2017). All three subspecies were represented in the common garden, but not equally. Of the 458 seedlings originally planted in the garden, 227 were A. t. tridentata, 114 were A. t. wyomingensis, and 117 were A. t. vaseyana. This study focused on only A. t. tridentata and A. t. wyomingensis and therefore, A. t. vaseyana shrubs were not measured.

Measurements were taken on three separate dates in mid-summer of 2019: 1 July, 11 July, and 25 July, at a time when we assumed that shrubs would be experiencing higher stress responses, resulting in potentially strong adaptive differences among genotypes. Our measurements were based on samples of convenience, with sampling starting at the far southeast corner of the garden and working northward up and down the rows. However, the randomized planting of the shrubs during initial establishment and the environmentally homogenous nature of the garden provided a general setting for nonbiased sample collection. Final measurements were taken on a total of 90 shrubs originating from 10 different states and 28 different populations. Fifty-one shrubs were diploid (2n) A. t. tridentata, 24 were tetraploid (4n) A. t. tridentata, and 15 were tetraploid A. t. wyomingensis.

Plant morphology

To determine the maximum height, shrubs were measured from ground level to the tip of the highest point of green tissue on the plant. To determine specific leaf area (SLA, leaf area/leaf dry mass), 8 cm branch-tip samples were collected from each shrub. No collections were made for shrubs whose health/survival would have been compromised by collection. Branch samples were stored in sealed plastic bags in a cooler and taken back to the Reinhardt plant physiology lab at Idaho State University (Pocatello, ID). All leaves at the distal end of collected branches were removed. From these leaves, a subset of 25 leaves from each branch was randomly selected and their one-sided leaf area was measured using Image J software (Schneider et al., 2012). The fresh mass of leaves was measured; leaves were then dried at 78°C in a drying oven for 24 h, and the dry mass was measured. To determine the sapwood area, thin sections were cut from each branch sample (one from the middle, and one each from the distal and proximal branch ends) and dyed using Toluidine blue stain (Carolina Biology Supply). The thin sections were then mounted onto microscope slides and photographed at 10× resolution using a digital camera software (AMScope, Irvine, CA, USA). As A. tridentata only has one active ring of xylem (the distal ring), only the area of the outermost ring was measured using ImageJ software (Ferguson, 1960; Schneider et al., 2012).

Gas exchange, chlorophyll fluorescence, and water relations

Measurements of dark respiration (R) and dark-adapted photosystem II efficiency (Fv/Fm) were taken predawn (between 3 and 7 am MST), while photosynthesis (A), stomatal conductance (gsw), transpiration (E), and photosystem II efficiency (PSII) were taken midday (between 10 am and 1 pm MST) to avoid total stomatal shutdown. Measurements were collected at the shoot-tip level using a Li-Cor LI-6800 portable photosynthetic system equipped with an LED chamber (Li-Cor Biosciences, Lincoln, NE). These measurements were then used to calculate the values of the photosynthesis-to-respiration ratio (A:R), water use efficiency (WUE; A/E), and intrinsic water use efficiency (WUEi; A/gsw), All gas exchange measurements were taken on randomly chosen shrub branches. Due to the complex nature of the leaf structure of A. tridentata, measurements were normalized by silhouette leaf area (Smith et al., 1991). Before measurements, selected branch tips were placed on a piece of paper within a drawn outline of the exact size and shape of the Li-Cor LI-6800 chamber. The natural arrangement of the leaves was maintained, and a picture was taken of the selected leaves to determine the true silhouette leaf area of the branch tip within the chamber. Silhouette leaf area was determined using ImageJ software (Schneider et al., 2012). The orientation of the shoot tip was replicated within the chamber of the Li-Cor for measurements. T_air and percent relative humidity inside the leaf chamber were kept at ambient conditions (1 July: 12–33°C, 20%–55%; 11 July: 17–35°C, 25%–54%; 25 July: 14–34°C, 14%–35%). During measurements, air flow in the chamber was set at a fixed rate of 500 mmol m⁻² s⁻¹.

Water potentials (Ψ) of branch tips were measured at predawn (Ψ_p) and midday (Ψ_m) using a pressure chamber (PMS Instrumental, model-1000, Albion, OR) on the same day and time that gas exchange measurements were obtained. Leaf water content was calculated using the formula (González & González-Vilar, 2001): $$$ \mathrm{Leaf}\kern0.3em \mathrm{Water}\kern0.3em \mathrm{Content}=\frac{\mathrm{Fresh}\kern0.3em \mathrm{Leafmass}(g)-\mathrm{Dry}\kern0.3em \mathrm{Leafmass}(g)}{\mathrm{Fresh}\kern0.3em \mathrm{Leafmass}(g)}, $$$ where leaves were weighed fresh, dried at 78°C in a drying oven for 24 h, and then weighed again once dry. Lastly, branch maximum hydraulic conductivity (K_max) was determined using a gravimetric “Sperry apparatus” (Sperry et al., 1988). Measurements were obtained from mature (>1-year-old and bark-covered) branch samples using a straight branch segment of ≥8 cm in length (larger than A. tridentata vessel length) to prevent the occurrence of open vessels (Kolb & Sperry, 1999). Excised branch segments were fully hydrated in an acidified (pH = 2), filtered (<20 μm) water bath in a vacuum chamber overnight to eradicate any existing embolisms in the xylem. The water was acidified to prevent microbial/fungal growth while the branches sat in the water overnight. Rubber tubing was then attached at each end of the branch and secured with Parafilm (Parafilm M). The branch lay horizontally, while secondary tubing was attached to the proximal end of the branch and ran vertically (to apply a gravimetric pressure head; ~0.009 MPa) to a flask of acidified, filtered water. At the distal end of the branch, a pipette was attached and secured to the workbench where the movement of water could be measured over time. All conductivity measurements were normalized to sapwood area and leaf area.

Quantification of impacts of climatic variables on shrub morphology and physiology

To determine the contribution of climate-of-origin to morphological and physiological variation, climate data from parent populations spatial locations (Rehfeldt et al., 2006; ) were reduced in dimensionality using principal components analysis (PCA). Separate PCAs were applied to a 13-variable temperature dataset and a 7-variable precipitation dataset to generate a smaller number of independent (albeit synthetic) predictors for temperature and precipitation, respectively (see Appendix S1: Table S1).

Parent-population climate variables were reduced to two principal components in the separate PCAs for temperature and precipitation. Ninety-six % of variability in the data was explained by the first two principal components from the PCA of temperature variables, whereas 97% of the variability in the data was explained by the first two principal components in the PCA of precipitation variables (Figure 1A,B, respectively). Temperature principal component one (TPC1) was strongly correlated with all 11 raw temperature variables. Variables with positive loadings with TPC1 included those positively associated with the length of winter, whereas variables with negative loadings included those positively associated with average temperature (Figure 1A). Temperature principal component two (TPC2) was strongly correlated with nine temperature variables. Variables with positive loading included those positively associated with the length of summer, whereas variables with negative loadings included those positively associated with colder temperatures (Figure 1A). Precipitation principal component one (PPC1) was strongly correlated with five precipitation variables. Variables with positive loadings included those that were positively associated with winter precipitation, whereas variables with negative loadings were those positively associated with summer, spring, and growing season precipitation (Figure 1B). Precipitation principal component two (PPC2) was strongly correlated with three precipitation variables. For this axis, variables with negative loading were those positively associated with growing season length, and mean annual precipitation increased (Figure 1B).

[IMAGE OMITTED. SEE PDF]

The distinctiveness of parent population environments from the Orchard Common Garden environment was measured by calculating the Euclidean distance of their respective climates, based on the raw climatic data used in PCAs. This variable will simply be referred to as “Euclidian distance” for the remainder of this paper.

Hypothesis testing

General linear models were used to test the association of response variables (physiology and morphology) to predictor variables (subspecies, cytotype, precipitation and temperature principal components, and Euclidian distance). In preliminary tests, no significant difference was found between the three dates of mid-summer measurements for responses under consideration; therefore, all three dates were grouped into a single dataset to allow for a larger sample size. The subspecies and cytotype predictors consistently demonstrated high levels of multicollinearity in preliminary models as measured by variance inflation factors (see Kutner et al., 2004). To address this, subspecies and cytotype were combined into a single categorical predictor with three designations: A. t. tridentata 2n (T2), A. t. tridentata 4n (T4), and A. t. wyomingensis 4n (W4). For models in which subspecies:cytotype was a significant predictor variable in omnibus models, post hoc pairwise tests were used to compare the three designations using Tukey's honestly significant difference (HSD) test, which optimally controls for familywise type I error in pairwise testing (Aho et al., 2014).

Model building

To identify useful models for response variables based on particular sets of predictors, backward and forward stepwise model selection was used to determine the most parsimonious model for each measured physiological and morphometric character response variable (Table 2). Model parsimony was gauged using the Akaike information criterion (AIC; Akaike, 1998; Aho et al., 2014).

Multivariate analyses

Multivariate similarities in the traits of shrubs were visualized using a two-dimensional nonmetric multidimensional scaling ordination (NMDS; Kruskal, 1964). Bray–Curtis dissimilarity was used as the underlying resemblance measure (Bray & Curtis, 1957). Vector fitting (Oksanen et al., 2013) of climatic principal components and Euclidean distance was used to visualize the contribution of these variables to trait variation in the summarized context of NMDS. This analysis allows us to better visualize how trait characteristics from individual shrubs varied simultaneously with respect to subspecies:cytotype and climate of origin.

Software

The R computational environment (R Core Development Team, 2022) was used for all statistical analyses with heavy reliance on the packages MASS (Venables & Ripley, 2002) and asbio (Aho, 2023) for model selection and simultaneous inference procedures, the package car (Fox & Weisberg, 2019) for general linear model diagnostics, the package vegan (Oksanen et al., 2013) for NMDS analyses, and the package ggplot2 (Wickham, 2016) for graphics.

RESULTS

Our initial prediction that subspecies would have minimal contribution to variation in phenotypic traits compared with cytotype and climate-of-origin was partially supported. Hypothesis testing revealed all measured responses to be correlated with different combinations of the explanatory variables under consideration, but cytotype and climate-of-origin (save for Euclidian distance) were responsible for the most significant variation in traits (p < 0.1; Table 1). However, overall phenotypic variation in our measurements was low, with only 18 of 102 null hypothesis tests of model parameters being statistically significant (Table 1). In addition, all significant general linear models had R² values between 0.04 and 0.40, with models possessing an R² value of ≥0.20 containing significant correlations with cytotype and climate-of-origin only (Table 1). Optimal approximating models (i.e., minimum AIC models) also revealed all measured responses corresponded with different combinations of the explanatory variables, with 6 of 14 models containing subspecies:cytotype and 12 of 14 models containing climate-of-origin predictor variables (e.g., TPC1 and TPC2, PPC1 and PPC2, and Euclidian distance). Lastly, NMDS showed a distinct lack of grouping corresponding to subspecies:cytotype, while vector fitting of the climate variables, emphasized temperature principal component two (TPC2), precipitation principal component two (PPC2), and Euclidian distance of climate-of-origin to the climate of the Orchard common garden.

TABLE 1 Results of multiple general linear regressions between response variables (measurements) and predictor variables (subspecies:cytotype, temperature principal components [TPCs] 1 and 2, precipitation principal components [PPCs] 1 and 2, Euclidean distance).

Trait	Subspecies:cytotype	TPC1	TPC2	PPC1	PPC2	Euclidian distance	R2
E (mmol m−2 s−1)							0.04
$$$ \hat{\upbeta} $$$		0.04	−0.15	0.14	−0.09	<0.01
t74	0.26	0.49	−0.89	0.87	−0.64	−0.01
p	0.77	0.63	0.38	0.39	0.52	0.99
A (mmol m−2 s−1)							0.18
$$$ \hat{\upbeta} $$$		0.29	−0.17	0.13	−0.14	<−0.01
t74	4.40	1.54	−0.45	0.37	−0.41	−1.65
p	0.016*	0.13	0.66	0.71	0.68	0.10
R (mmol m−2 s−1)							0.22
$$$ \hat{\upbeta} $$$		0.06	0.24	−0.10	0.05	−0.01
t74	1.76	1.21	2.14	−1.02	0.52	−1.51
p	0.18	0.23	0.04*	0.31	0.61	0.14
A:R							0.28
$$$ \hat{\upbeta} $$$		0.07	−0.03	0.12	0.02	<0.01
t74	1.37	2.00	−0.34	1.84	0.34	−0.90
p	0.27	0.05***	0.74	0.07***	0.74	0.38
gsw (mol m−2 s−1)							0.05
$$$ \hat{\upbeta} $$$		<0.01	<0.01	<0.01	<0.01	<−0.01
t74	0.05	0.71	0.02	0.33	0.50	−0.75
p	0.95	0.48	0.98	0.75	0.62	0.46
PSII							0.24
$$$ \hat{\upbeta} $$$		0.01	<0.01	−0.07	0.04	<0.01
t74	8.94	0.44	0.08	−2	1.08	−1.82
p	0.0003**	0.66	0.94	0.05*	0.29	0.07***
Fv/Fm							0.08
$$$ \hat{\upbeta} $$$		<0.01	−0.01	<0.01	<−0.01	<0.01
t74	0.20	0.11	−1.22	0.89	−0.76	1.48
p	0.82	0.91	0.23	0.38	0.45	0.14
WUE (mmol/mmol)							0.14
$$$ \hat{\upbeta} $$$		−0.06	−0.05	−0.06	0.43	<−0.01
t74	0.60	−0.33	−0.15	−0.19	1.41	−1.11
p	0.55	0.74	0.88	0.85	0.17	0.27
WUEi (mmol/mol)							0.08
$$$ \hat{\upbeta} $$$		−5.77	7.32	−2.12	16.97	<0.01
t74	0.04	−1.33	0.82	−0.29	2.15	0.06
p	0.96	0.19	0.42	0.80	0.04*	0.96
LWC (%)							0.03
$$$ \hat{\upbeta} $$$		0.01	−0.01	<0.01	−0.01	<−0.01
t74	0.23	1.31	−0.79	0.13	−0.90	−0.50
p	0.79	0.19	0.44	0.90	0.37	0.62
SLA (g/cm2)							0.08
$$$ \hat{\upbeta} $$$		−1.09	0.47	1.90	−0.58	0.01
t74	0.53	−0.74	0.17	0.73	−0.23	1.33
p	0.59	0.46	0.87	0.47	0.82	0.19
LA:SA (cm2/mm2)							0.10
$$$ \hat{\upbeta} $$$		−4.10	−2.53	−1.90	8.34	0.01
t74	2.68	−0.87	−0.27	−0.22	1.01	0.45
p	0.075***	0.39	0.79	0.82	0.31	0.65
Kmax (mol m−1 s−1 mPa−1)							0.08
$$$ \hat{\upbeta} $$$		13.62	−101.60	26.19	53.69	−0.27
t74	2.14	0.30	−1.15	0.32	0.69	−0.95
p	0.13	0.76	0.25	0.75	0.50	0.35
Kmax_leaf (mol m−1 s−1 mPa−1)							0.10
$$$ \hat{\upbeta} $$$		<0.01	<0.01	<0.01	0.01	<−0.01
t74	1.97	0.29	0.30	0.11	0.82	−0.95
p	0.15	0.80	0.77	0.92	0.42	0.35
Ψp (MPa)							0.04
$$$ \hat{\upbeta} $$$		−0.03	0.06	−0.03	−0.02	<0.01
t74	0.68	−0.79	0.77	−0.42	−0.27	0.37
p	0.51	0.43	0.44	0.68	0.79	0.72
Ψm (MPa)							0.11
$$$ \hat{\upbeta} $$$		−0.08	0.14	−0.03	0.13	<0.01
t74	1.81	−2.14	1.94	−0.46	2.08	0.55
p	0.17	0.04*	0.06***	0.65	0.04*	0.59
Maximum height (cm)							0.40
$$$ \hat{\upbeta} $$$		1.04	−2.2	1.20	−4.34	−0.01
t74	20.77	0.87	−0.95	0.56	−2.08	−0.18
p	<0.00001**	0.39	0.35	0.58	0.04*	0.86

While subspecies:cytotype was a significant predictor in some models (Table 1), post hoc pairwise analyses revealed only A, LA:SA, leaf area-specific K_max, and maximum shrub height varied significantly between the designations T2, T4, and W4, and in most cases, this variation was attributable to cytotype, not subspecies (Figure 2). For example, photosynthesis in diploid A. t. tridentata was significantly higher than that in tetraploid A. t. tridentata (p = 0.007) and tetraploid A. t. wyomingensis (p = 0.058), but no difference was evident between tetraploid A. t. tridentata and A. t. wyomingensis shrubs (Figure 2A). In addition, LA:SA was significantly higher in diploid A. t. tridentata shrubs than in A. t. wyomingensis (p = 0.059), and diploid A. t. tridentata had significantly lower values of leaf area-specific K_max (Figure 2B,C). The only trait to significantly vary by subspecies and cytotype was maximum shrub height, in which diploid A. t. tridentata was the tallest, followed by tetraploid A. t. tridentata and A. t. wyomingensis (Figure 2D).

[IMAGE OMITTED. SEE PDF]

Cytotype and climate-of-origin were responsible for the most significant variation in traits (coinciding with our initial prediction), with climate-of-origin responsible for 9 of 18 significant correlations observed in hypothesis testing. Climate-of-origin was a significant predictor in linear models for multiple traits (p $$$ \le $$$ 0.1; Table 1). Temperature principal components (TPC1 and TPC2) varied significantly with traits relating to respiration and water status (Table 1). Shrubs from climates with longer winters and colder temperatures exhibited higher values of A:R and lower values of R and Ψ_m than shrubs originating from climates with longer summers and hotter temperatures (Figure 3). Precipitation principal components (PPC1 and PPC2) varied significantly with traits associated with physiology and morphology (Table 1). Shrubs that originated from climates with more winter precipitation, but lower summer and spring precipitation exhibited higher values of A:R and lower values of ΦPSII than shrubs from climates with less winter precipitation but more summer and spring precipitation (Figure 4A,C). Additionally, traits associated with water status and efficiency were significantly correlated with PPC2 (Table 1). WUEi increased and Ψ_m became less negative as growing season shortened and mean annual precipitation decreased. Maximum shrub height decreased as growing season shortened and mean annual precipitation decreased (Figure 4E). Lastly, a simple linear regression showed ΦPSII decreased as Euclidian distance of the climate-of-origin increased from the Orchard Common Garden (Figure 5).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

In addition to displaying the highest explained variation in traits than all other predictor variables, optimal approximating models included climate-of-origin in most models (12 of 14; Table 2). Temperature principal components appeared in 8 of 14 models, while precipitation principal components appeared in 9 of 14 models, and despite a lack of statistical significance (Table 1), Euclidean distance appeared in 6 of 14 of the chosen models (Table 2). Optimal approximating models also generally included subspecies:cytotype in 6 of 14 models, but most of those models corresponded with the same traits that showed cytotype-specific variation (e.g., A, LA:SA, K_max, and maximum height; Table 2). Therefore, except for models for maximum height, in all approximating models in which subspecies:cytotype was selected, cytotype was most likely the driving genetic factor. Lastly, model selection resulted in intercept-only models for E, LWC, K_max, and Ψ_p, indicating the predictor variables under consideration had little to no effect on these response variables.

TABLE 2 Results of Akaike information criterion (AIC) model selection (chosen model and its corresponding Akaike weight).

Trait	Best (low AIC) model	Akaike weight
A (μmol m−2 s−1)	= subspecies:cytotype + TPC1 − Euclidian distance	0.62
R (μmol m−2 s−1)	= TPC2	0.33
A:R	= TPC1 + PPC1	0.47
gsw (mol m−2 s−1)	= TPC1 − Euclidian distance	0.54
PSII	= subspecies:cytotype + TPC1 − PPC1 − Euclidian distance	0.61
Fv/Fm	= PPC1 + Euclidian distance	0.46
WUE (μmol m−2 s−1/mmol m−2 s−1)	= PPC2 − Euclidian distance	0.59
WUEi (μmol m−2 s−1/mol m−2 s−1)	= TPC1 + PPC2	0.55
SLA (g/cm2)	= TPC1 + PPC1 + Euclidian distance	0.61
LA:SA (cm2/mm2)	= subspecies:cytotype	0.42
Kmax (mol m−1 s−1 mPa−1)	= subspecies:cytotype + PPC1	0.54
Kmax_leaf (mol m−1 s−1 mPa−1)	= subspecies:cytotype	0.49
Ψm (mPa)	= TPC1 + TPC2 + PPC2	0.40
Maximum height (cm)	= subspecies:cytotype − PPC2	0.54

Our prediction that subspecies would have a lesser contribution to phenotypic variation in A. tridentata when compared with cytotype and climate-of-origin was further supported by results from our nonmetric multidimensional scaling (NMDS) projection. NMDS revealed that the three subspecies–cytotype designations (T2, T4, and W4) were physiologically and morphologically variable within designations and indistinct from each other. Notably, two outlier shrubs occurred in the ordination at NMDS1 ~0.75 (tetraploid A. t. wyomingensis) and NMDS1 ~1.6 (diploid A. t. tridentata). When environmental vectors were placed over the plot, TPC2 and PPC2 were the strongest drivers of physiological variation, followed by Euclidean distance (Figure 6).

[IMAGE OMITTED. SEE PDF]

DISCUSSION

In this study, we measured the physiological and morphological traits of A. tridentata during mid-summer (July), seven years after individual shrubs were planted in a common garden to better understand how genetic (subspecies:cytotype) and environmental variables (climate-of-origin) contribute to phenotypes in a widespread, drought-adapted shrub species. Based on earlier findings in this garden (Chaney et al., 2017; Lazarus et al., 2019; Richardson et al., 2021), we predicted that local adaptation and cytotype would be primary drivers for phenotypic variation in A. tridentata and that subspecies would be less influential. We evaluated these predictions using formal hypothesis testing, model selection, and multivariate graphical summaries.

While overall phenotypic variation within our study was low (Table 1), most of our findings lent support to the idea that climate of origin was the predominant driver of A. tridentata phenotypic variation. The results from our general linear models, AIC model selection, and NMDS ordination were consistent with several recent studies considering A. tridentata physiology (Chaney et al., 2017; Germino et al., 2019; Lazarus et al., 2019). Regression models with the highest explanatory power in this study were those containing climate-of-origin predictor variables (i.e., PCAs and Euclidian distance; Table 1). Models with the highest explanatory power were those concerning photosynthetic efficiency (A:R and PSII) and carbon metabolism (R; R² = 0.28, 0.24, and 0.22 respectively). These R² values are comparable to Richardson et al. (2021) who linked ploidy and climate of origin to growth and seed production in A. tridentata.

The trends seen between measured traits and climate-of-origin predictor variables (e.g., PCAs and Euclidian distance) were also what we would expect to see from locally adapted shrubs in a common garden setting. For example, respiration (R) increased with TPC2 and A:R increased with TPC1 (Figure 3A,B), meaning shrubs from populations that experienced longer summers, shorter winters, and overall warmer temperatures exhibited higher R and A:R than shrubs whose originating populations experienced more mesic conditions. This response was expected, as shrubs from seed sources adapted to climatic conditions more xeric than the Orchard Common Garden (e.g., longer summers, shorter winters, and warmer temperatures) would most likely show a greater photosynthetic efficiency and carbon relations than those from climatic conditions more mesic than Orchard. In addition, Ψ_m was negatively correlated with TPC1, but positively correlated with TPC2 (Figure 3C,D), meaning shrubs that originated from hotter climates with longer summers and shorter winters had less midday water stress—an unsurprising trend, given trends for R and A/R. Higher winter precipitation and low summer, spring, and growing season precipitation imply long periods of drought throughout the year. Thus, shrubs from these environments are more likely to be adapted to use water more efficiently than shrubs from populations that experienced little drought.

While climate-of-origin was the predominant predictor of phenotypic variation in this study, cytotype also played a significant role. Post hoc pairwise analyses showed significant differences in A, LA:SA, leaf-area specific K_max, and maximum shrub height between diploid and tetraploid individuals both within and between subspecies (Table 1, Figure 2). Fewer traits were correlated with climate-of-origin variables, and fewer still were significantly correlated with subspecies. In addition, subspecies:cytotype was selected in 6 of 14 AIC-driven models. NMDS ordination, however, showed a lack of subspecies:cytotype grouping of shrubs, indicating an absence of distinct morphological and physiological patterns among cytotypes. This outcome could potentially be attributed to the lack of phenotypic variation within the study system.

The lack of overall phenotypic variation evident in this study could be due to a number of factors, including the age of the common garden and phenotypic plasticity. Our initial prediction was that strong adaptive differences would occur in mid-summer when environmental stress was high, coinciding with the timing of drought stress/tolerance measurements. However, our results (e.g., relatively small genetic effects in models)s (Table 2) and a lack of grouping of shrubs in our NMDS ordination (Figure 5) suggest that drought stress mitigation via plasticity may override genetic effects. This hypothesis is supported by evidence from other studies that emphasize the importance of seasonality for detecting subspecies and cytotype-level differences in A. tridentata (Booth et al., 1990; Kolb & Sperry, 1999; Welch & Jacobson, 1988; Zaiats et al., 2020). For example, Kolb and Sperry (1999) showed variation in hydraulic functional traits was greatest between subspecies of A. tridentata in spring/early summer, with significant differences decreasing into mid- late-summer, and Zaiats et al. (2020) indicated that differences in root water uptake timing were most evident in the growing season, as opposed to the summer months. Unfortunately, the use of a single common garden prevented the assessment of the contribution of phenotypic plasticity to population phenotypic variation.

While common garden studies can be used to effectively consider the effects of genotypic variation in plant species, the utilization of a single garden in a single environment prevents the implementation of a full genecological study (sensu Turesson, 1923). Genecological studies are the primary tool for determining how landscape and spatial variation lead to genetic differences among plant populations and also allow consideration of the effects of phenotypic plasticity (expression of different phenotypes from a single genotype in response to environmental change) on physiological variation (Chaney et al., 2017; Matesanz et al., 2010). Evidence of phenotypic plasticity has been previously observed in A. tridentata in genecological studies (Richardson et al., 2017, 2021); however, the range of measured traits in Richardson et al.'s studies was limited, and further exploration of the capacity for phenotypic plasticity in A. tridentata would be beneficial.

In addition, the lack of phenotypic variation within our study could be due to the age of the common garden. A paper by Germino et al. (2019) suggests that a significant amount of time may be required to detect evidence of local adaptation when working in a common garden of perennial shrubs. While the number of measured traits in Germino et al. (2019) were limited (size, survival, and fecundity; time since garden establishment, concentration of measurements into a short month/year, etc.), the authors found that >20 years were required to see the evidence of local adaptation in common gardens of A. tridentata, likely due to the long lifespan conservative growth strategy of the perennial shrub (Germino et al., 2019). The Orchard Common Garden was established in 2012, making the age of the shrubs only seven at the time of measurements, and decreasing potential variations in physiology due to parent population climate. Alternatively, it is possible that less fit genotypes had already been culled out in the common garden by the time of our study, leaving less potential for physiological variability. The Orchard Common Garden experienced high rates of mortality within the first year of establishment (Chaney et al., 2017). While mortality was highest in A. t. vaseyana, culling also occurred for A. t. tridentata and A. t. wyomingensis.

CONCLUSION

Measurements of physiological and morphological variation of mature A. tridentata during mid-summer in a common garden revealed an imbalanced contribution to phenotypic variation by subspecies, cytotype, and climate of origin to phenotypic variation. While post hoc pairwise analysis revealed subspecies-level differences in maximum height, formal linear model hypothesis testing and AIC model selection emphasized the importance of ploidy level and climate of origin to phenotypic variation, and nonmetric multidimensional scaling revealed an overall lack of correspondence between subspecies and physiological and morphological variation. These findings suggest A. tridentata potentially has highly plastic physiology or that culling in Orchard over the seven years since its establishment has led to a lack of genetic diversity within the garden.

A. tridentata populations currently face several anthropogenic stressors, including habitat fragmentation and invasive species that increase steppe fire severity and decrease fire return intervals. Further, arid/semiarid environments are predicted to undergo major transformations in temperature and precipitation regimes in the coming decades (Abatzoglou & Kolden, 2011). Understanding drivers of physiological variation could help properly manage steppe ecosystems and predict climate change-induced events. However, the link between phenotypic variation and genotype is still widely unexplored for arid/semiarid perennial forbs and shrubs. Thus, understanding the drivers of phenotypic expression in A. tridentata could clarify how climate change may affect survival and extirpation and help increase the effectiveness of steppe management and restoration efforts.

ACKNOWLEDGMENTS

This project described was supported by NSF award number OIA-1757324 from the NSF Idaho EPSCoR Program and by the National Science Foundation. The findings and conclusions in this publication are those of the authors and should not be construed to represent any official USDA determination or policy. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. This material is based upon work supported by (while serving at) the National Science Foundation. We thank Adler Patch, Harrison Seitz, and Tierin Osterfeld for their assistance in the field/laboratory and the two anonymous reviewers who helped improve this manuscript.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

Data and code (Roop, 2025) are available from Zenodo: .

References

Abatzoglou, J. T., and C. A. Kolden. 2011. “Climate Change in Western US Deserts: Potential for Increased Wildfire and Invasive Annual Grasses.” Rangeland Ecology & Management 64: 471–478. [DOI: https://dx.doi.org/10.2111/REM-D-09-00151.1].