Phenological Variation in Bluebunch Wheatgrass ( Pseudoroegneria spicata ): Implications for Seed Sourcing, Harvest, and Restoration

: To reduce maladaptation in cultivated seed lots, seed transfer zones (STZs) have been developed for grasslands and other habitats using morphological traits and phenological measurements that only capture the ﬁrst day of events such as ﬂowering and seed ripening. Phenology is closely linked to plant ﬁtness and may affect genetic loss during harvests of seed raised for ecological restoration. Here, we measured the detailed phenologies of 27 populations from six STZs of bluebunch wheatgrass ( Pseudoroegneria spicata ) (Pursh) Á . Löve (Poaceae) raised in a common garden to test whether existing STZs created using a combination of plant morphology and “ﬁrst-day” phenological measurements adequately capture population-level variation in season-long, detailed phenologies. We also used detailed phenologies to test whether genetic losses may occur during single-pass harvests of commercial seed. Mixed and random effect models revealed differences in detailed reproductive phenology among populations within two of six STZs. The number of individual plants within an STZ not producing harvestable seed during peak harvest levels indicated that 10–27% of individuals from a seed lot could be excluded from a single-pass harvest. Although our ﬁndings generally support current STZ delineations for P. spicata , they point to the possible precautionary importance of sourcing from multiple populations and harvesting with multiple passes when resources permit. Author Contributions: Conceptualization, K.P., M.R.O., F.F.K. and H.R.P.; methodology, K.P., M.R.O., F.F.K., R.J.R. and H.R.P.; formal analysis, K.P.; investigation, K.P. and H.R.P.; data curation, K.P.; writing—original draft preparation, K.P.; writing—review and editing, K.P., M.R.O., F.F.K., R.J.R. and H.R.P.; project administration, F.F.K. and R.J.R.; funding acquisition, M.R.O., K.P., F.F.K. and H.R.P.


Introduction
Ecological restoration is needed over large areas of the Intermountain West in grasslands disturbed by wildfire, over-grazing, and invasive species [1,2]. Native local seed sources adapted to project site conditions are in need of restoration because their use may improve population establishment and persistence [3][4][5][6]. While genetic diversity has been successfully conserved in some propagated lines of Intermountain West plants [7], commercially available sources often represent a fraction of the genetic diversity that exists across the areas in which they are planted [8,9].
Genetic variation plays a fundamental role in determining whether species can adapt to variable environmental conditions and thus requires consideration in restoration [10][11][12]. Plants tend to be adapted to local environmental conditions such as aridity and temperature [13,14]. Conserving genetic variation also can be crucial for sexual reproduction [15], disease resistance [16], and natural recovery after disturbance [17]. Genetically diverse seed sources are at times more successful at colonizing new habitats than those with less variation, regardless of seeding density [18]. Consequently, meeting commercial demands for native seed requires knowledge of adaptation and genetic variation across landscapes.
Maintaining genetic diversity while producing commercially viable quantities of native plant material can be difficult [11,19]. To address this challenge, land managers have developed seed transfer zones (hereafter referred to as 'STZs') to guide the selection of plant material. STZs identify geographic areas across which seed from a particular native species can be collected and planted with reduced risk of maladaptation at the planting site [20,21]. In STZ studies on grassland species in the Intermountain West [22,23], strong trait-by-environment relationships have been found for phenology, suggesting local adaptation [13].
In species where flowering, seed ripening, and dispersal patterns consist of recurring flushes throughout the growing season, STZs ideally should consider patterns of reproductive phenology through the entire season. Methods exist for quantifying season-long patterns in phenology [24], but until now, variation in only the first day of anthesis and ripening, along with morphological traits, have been used to develop STZs for restoration species in the Intermountain West [23,25]. It is therefore of interest to test whether STZs designed using a combination of simpler phenological metrics and morphology capture the range of phenologies that may exist among populations and STZs assessed throughout the growing season.
A better understanding of season-long phenologies has direct implications for successful restoration for two reasons. First, the timing of reproductive events in seasonal environments is generally correlated to plant fitness [26,27]. For instance, two years post-fire, restoration seeding in bottlebrush squirreltail (Elymus elymoides) experienced directional selection for earlier flowering phenology [28]. Likewise, climate-mediated flowering phenology can impact pollination success [29], alter pollination dynamics and cause assortative mating [30], influence seed number [31], and alter seed predation dynamics [29].
A second reason to study season-long phenologies in STZs is that restoration seed diversity can be reduced not only during sourcing but also during harvest. In typical agronomic practice, seeds are harvested in a single pass on a single day [12]. Single-pass harvests of native seed crops have been shown to reduce within-population variance in the timing of reproductive stages and cause inadvertent selection for populations with later phenology [32]. This and other forms of unintentional selection during seed increase can lower the probability that diverse genotypes are included in the resulting seed lots and jeopardize restoration success [19].
In this study, we focus on bluebunch wheatgrass (Pseudoroegneria spicata), a bunchgrass commonly used for ecological restoration in the Intermountain West [33,34], that has high phenotypic variance among individuals and populations [23], and which enters flowering, ripening, and dispersal in recurring flushes. We examined the following two null hypotheses with regard to commercial sourcing and harvest of seed for restoration: (1) single populations are sufficient to capture variation in season-long phenologies within an STZ. We tested this hypothesis by examining phenological variation in anthesis, ripening, and dispersal within and between zones. (2) A single-pass harvest is sufficient to capture seed from most or all individuals within a zone. We tested this hypothesis by examining the proportion of individuals that lacked seed in the dispersal stage on the date of peak dispersal across a zone.

Experimental Species and Design
Pseudoroegneria spicata (Pursh) Á. Löve (bluebunch wheatgrass) is an outcrossing, wind-pollinated, cool-season, perennial bunchgrass native to western North America and common in grasslands, shrub steppe habitats, and woodlands [33]. Phenotypic variation among P. spicata populations is correlated with seed source climates, which is a signature of local adaptation [13,23]. To help land managers limit maladaptation in restoration, the Intermountain West was delineated into 11 STZs for P. spicata, with lower numeric classifications signifying warmer and drier climates and earlier phenology and higher numeric classifications signifying cooler and wetter climates and later phenology [23]. Despite the wide geographic range of some STZs (Figure 1), plants within a zone share similar morphologies and first-day phenological patterns [23].
To examine genetic variation in phenology within an STZ, in 2013, wild-source seeds were collected from 4-5 populations each in STZs 3a, 4, 6a, 6b, 7a, and 7b (29 populations total; Figure 1, Table S1). Our populations were sourced both within seed zones and within level III ecoregions, which is an approach recommended in the original zone designations for P. spicata for land managers who wish to be conservative about maintaining adaptation and local genetic structure [23]. Seeds were planted in Super Cone-tainers (Stuewe and Sons, Inc, Tangent, Oregon, USA) and grown for four months in a greenhouse at the USDA Forest Service Rocky Mountain Research Station nursery facility in Moscow, ID. Seedlings were transplanted in early November 2014 to a common garden on private land located near the Crooked River Grasslands in Central Oregon within zone 7b, with a mean annual temperature of 7.8 • C (46 • F), 356 mm (14 inches) of annual precipitation, and a plant community characterized by Festuca idahoensis (Idaho fescue), Artemisia tridentata (big sagebrush), Juniperus occidentalis (western juniper), and P. spicata. Common gardens are a standard approach to testing for genetic variation because by minimizing environmental variation it can be inferred that phenotypic variation is explained by genes [35]. Prior to planting, we used a walk-behind tiller to remove existing vegetation from the garden and promote plant establishment.

Data Collection
Reproductive phenology was monitored in all surviving individuals from late May to early August 2016 approximately every four days, except for one visit, which was skipped in mid-July. Monitoring captured the first instance of anthesis, as confirmed by two visits pre-anthesis in mid-May. Many plants produced more than 200 flowering We divided the garden into four replicate blocks. Each block contained one randomly placed plot per STZ for each of the 6 STZs. Each plot contained plants from five or more individuals from each population within its zone. Individuals were interspersed randomly by population within each plot. In total, 580 P. spicata were planted (4 replicate blocks x 6 plots/block x 4-5 populations per plot x 5 individuals per population), for a total of 20 individuals per population in the entire garden (5 individuals per population per block x 4 blocks). Plants were spaced 0.6 m apart within rows and 0.8 m apart within columns. The garden was surrounded by two rows of P. spicata to reduce edge effects and was hand-weeded twice during the growing season in both 2015 and 2016. Plants were allowed to grow for one year prior to collecting data to reduce any greenhouse effects. After one year, 461 (79%) of plants survived (Table S1).

Data Collection
Reproductive phenology was monitored in all surviving individuals from late May to early August 2016 approximately every four days, except for one visit, which was skipped in mid-July. Monitoring captured the first instance of anthesis, as confirmed by two visits pre-anthesis in mid-May. Many plants produced more than 200 flowering spikes. Because it was impossible to monitor every spike in the garden, and because we wished to partition monitoring effort across the entire growing season, we labeled and subsequently monitored a random sub-sample of up to five spikes entering anthesis per individual per visit (but fewer if fewer new spikes were present) over each of the five distinct 12-14 day time periods (hereafter cohorts), resulting in a maximum of 5 spikes x 5 cohorts = 25 spikes per individual. We labeled each spike loosely with tape on the stem on the date of first anthesis and tracked it until the study ended on 4 August 2016.
Plants that lived but did not produce flower spikes were excluded, resulting in 389 reproductive plants (i.e., plants that entered anthesis) in the dataset (65.0 ± 12.2 sd plants per STZ; Table S1). Each STZ was represented by at least four populations, and each population by 13.4 ± 4.2 sd individuals per population (range 6-19, Table S1). One population (#133 from STZ 3a, Table S1) was removed from the dataset because it had fewer than six flowering individuals, which we deemed the minimum to represent populationlevel variation.
Flower spikes were scored on their progression through growth stages according the extended BBCH (Biologische Bundesantalt, Bundessortenamt and CHemische Industrie) scale [24], composed of 2-digit codes. The first digit corresponds to the primary growth stage (e.g., anthesis, seed ripening, dispersal) and the second describes intermediate stages within the primary stage (e.g., half of the florets have produced anthers, seeds have reached early milk stage, half of the seeds have dispersed, Table S2). We focused on anthesis, seed ripening, and dispersal because the timing of anthesis and seed ripening are commonly measured in STZ studies [23,25] and because these traits correlate to plant fitness [26,27]. The BBCH scores of flower spikes within a cohort on the same plant were averaged to produce multiple cohort-level BBCH scores for each plant for each site visit. For a heatmap analysis, we also used a plant-level BBCH score, which averaged all cohorts for a plant.
Cohort-level BBCH scores were linked to the cumulative growing degree day values for each site visit. Anthesis was defined as a BBCH score of 6.1-6.9, ripening as 7.5-8.0, and dispersal as 8.1-8.9 (Table S2). BBCH scores of 7.0-7.4 reflect the fruit development stage, marked by an absence of anthers and preceding viable seed, and were thus excluded from our statistical analysis. Cumulative growing degree days are often used instead of calendar days to detect and/or predict phenological events in plants [36] and they have been used to predict growth stage in various perennial grass species [37]. Cumulative growing degree day values were calculated by averaging the daily minimum and maximum temperatures starting on 1 January 2016 from the nearest weather station (Grizzly Weather Station, High Plains Regional Climate Center) with a base of 2 • C, which was used because it is a minimum for vegetative growth for P. spicata [38].

Statistical Analysis
All analyses were carried out in R statistical software version 3.5.1 (R Core Team, 2015). We used generalized linear mixed effects models (GLMMs) and likelihood ratio tests to examine differences in the timing of anthesis, ripening, and dispersal among populations within the same STZ. A GLMM was developed for each response variable: cohort-level anthesis, cohort-level ripening, and cohort-level dispersal timing. In each model, the population was a fixed effect (explanatory variable), and block, plant, and the interactions (due to crossing) between blocks and populations were random effects. A negative binomial error distribution provided the best fit for anthesis and a Poisson error distribution with a log link best fit for ripening and dispersal. The glmer.nb or glmer function from the 'lme4' package (Bates et al., 2015) was used to construct all GLMMs with the negative binomial and Poisson error distribution, respectively. All GLMMs were visually inspected for violations to model assumptions using the simulateResiduals function in the package 'DHARMa' [39]. To test for differences in cohort-level anthesis, ripening, and dispersal timing among populations from the same STZ, a likelihood ratio test was performed for each stage using the anova function of the 'lme4' package [40].
To examine broader patterns of variation among STZs, GLMMs were used to calculate intra-class correlation coefficients (ICCs) for STZ, population, plant and block, and the interactions block: zone, and block: population. ICCs are a measure of the proportion of variance that can be attributed to particular objects of measurement [41]. The same error distributions were used as stated above. ICC calculations for all random effects models were executed with the function icc from the 'sjstats' package in R statistical software version 3.5.1 [42,43].
To quantify the potential for unintentional selection during seed harvest, we used heatmaps of BBCH scores for each individual during each site visit to determine (a) the optimum date to collect seed from the largest number of individuals during a hypothetical single-pass seed harvest and (b) the percentage of plants that lacked ripe seeds (and would therefore be excluded from the resulting seed lot) on that optimum date. BBCH scores for each individual were the average of its cohort BBCH scores (described above). The optimum harvest date was defined as the date when the greatest number of plants from within each STZ exhibited ripe seeds (individual BBCH scores 7.5-8.1). Defining the optimum harvest date in this way is consistent with seed harvest approaches for other grass species with progressive ripening (pp. 191-212, [38]) [39]. All plants that were in pre-(anthesis) or post (dispersal) ripening stages on the optimum harvest date were tallied and counted towards the percentage of "lost" plants. This analysis of seed loss during a singlepass harvest was performed for each STZ as well as all populations without accounting for STZ.

Results
By the end of the growing season, 390 of the 461 surviving plants (84%) had at least one cohort that reached reproductive stage (i.e., BBCH anthesis); of those, 94% of plants had at least one cohort that reached BBCH dispersal. Differences in the timing, as measured by cumulative degree days, of plant-level BBCH scores for anthesis, ripening, and dispersal were uncommon among populations sourced from the same STZ (Table 1). The timing of anthesis was only statistically different among populations within STZ 7b, a generally cooler and wetter climate that supports plants with later phenology and larger stature [23]. For the ripening and dispersal stages, statistical differences among populations were only found within STZ 6b (Table 1). For each of the three phenological traits, variance partitioned differently among STZs, populations, and individual plants, with much of the variance unexplained (Table 2). We found that a hypothetical single-pass harvest at peak ripeness could reduce seed collection by 10%-27% of individual plants, depending on the STZ (Table 3, Table 4, and Table S3a-e).  Table 2. Intra-class correlation coefficient (ICC) scores based on random effect models built for each reproductive stage. ICC scores within a column indicate that the structuring of the data into each random effect accounted for the given percentage of variation in that stage. The "unexplained" column is the total unexplained variance for each phenological stage.   Table 4. Heatmap of BBCH phenology scores for P. spicata plants sourced from STZ 3a. 'n' is the total number of plants in the zone that reached reproductive stage. Cell values are the sum of all plants in the zone exhibiting a given BBCH score (rows) on a given calendar date (CD) or cumulative growing degree day (GDD) (columns). Color scale darkens with more plants observed in each stage per date of observation. Since P. spicata produces multiple flowering flushes, the progression from anthesis to dispersal is not always linear. The dark box surrounding the July 18 visit date was added for emphasis to indicate the optimum date of seed harvest (i.e., highest number of plants with harvestable seeds) in the zone. The percentage of individuals in an STZ not in a ripening stage at the time of peak harvest (Table 3) was calculated as the "sum" cell beneath the dark lines divided by n. Table S3a-e show heatmaps for remaining zones.

Discussion
Long-term prospects of restored plant communities are improved by present and future adaptation, which requires genetic diversity [10,44]. To maintain landscape-level adaptive genetic variation within a given STZ, native plant materials ideally would be sourced to contain a complete range of flowering, ripening and seed dispersal phenologies because these traits are associated with plant fitness [26,27]. In four of the six STZs, we found no statistically significant difference in season-long phenology among populations within an STZ (Table 1). A lack of phenotypic variation within STZs has been used as evidence for proper boundary delineation [45,46] and indicates that a single population may be sourced to represent season-long phenologies throughout much of the zone. The general applicability of STZs developed for P. spicata is further supported by three differences between the populations in our study and the original populations used to develop the zones: our populations were collected eight years later (2013 vs. 2005); they were sourced from different locations; and they were raised in a different STZ (7a vs. 4, 6a and 6b) [23].
It is possible that had we collected across multiple level III ecoregions in an STZ, we would have discovered more variation within zones. However, our approach was consistent with recommendations of the original seed zone designations for P. spicata, which state that a conservative approach to seed sourcing is to collect both within STZs and within level III ecoregions [23]. We allocated our experimental effort to test for variation consistent with this recommendation.
We did, however, find exceptions to the general applicability of the STZs for capturing variation in season-long phenologies. Anthesis differed among populations from one STZ (7a), and both ripening and dispersal differed among populations from another STZ (6b). Other species exhibited genetic variation in reproductive phenology across the topographically and environmentally complex Intermountain West. For instance, Achnatherum hymenoides (Indian ricegrass) varied by population in three traits relating to reproductive phenology [25]. In a meta-analysis of species native to the Great Basin, population variation and correlations among traits and environmental variables existed in over half of 76 studies that evaluated the date of flowering and 14 studies that examined seed shattering, and pointed toward a higher prevalence of local adaptation in the Great Basin than other sites globally [13].
The only STZ that exhibited significant population-level variation in ripening, 6b (Table 1), also had the largest number of individuals excluded from a single pass harvest (Table 3). Moreover, the STZ with the lowest chi square values for population-level variation in ripening, 7a (Table 1), lost the fewest individuals in a single-pass harvest (Table 3). These relationships make sense, since both were derived from the same ripening data; they highlight the fact that STZs most at risk for genetic loss of phenological variation during collection may remain at risk for loss during harvest, even if multiple populations are sourced. Loss during harvest was shown to be a risk in two native perennial grass species, Elymus glaucus and Nassella pulchra, which changed in genetic variation after two years of single-pass mechanical harvesting [32]. In another study, awned slender wheatgrass (Elymus trachycaulus) exhibited a potential loss of up to 8% of the overall genetic variation after two generations of typical seed increase [47]. Multiple pass harvests have been recommended to avoid such losses [32].
Our study is too limited to identify general characteristics, if any, of P. spicata STZs most at risk for genetic losses caused by single population collections and/or single-pass harvests. Neither the STZ with the greatest distances between populations, 3a (Table S1), nor the STZ with the greatest variance in elevation among populations, 4 (Table S1), exhibited significant variation in season-long phenological traits among populations. The best advice we can offer restoration practitioners is that single collections of P. spicata within a zone, as long as they occur in the same level III ecoregion, are likely to capture most important genetic variation in phenology, but extra populations might be sourced for seed when feasible. If extra populations are sourced, then extra attention should be given to multiple-pass harvests in order not to lose variation in phenology during production.
Another limitation of our study was the use of a single common garden location, which made it impossible to assess possible influences of phenotypic plasticity across multiple rearing environments. In contrast, the study that delineated STZs for P. spicata utilized three garden sites and found a site effect on the majority of morphological traits and all "first-day" phenological traits [23]. Thus, the results reported here must be considered sitespecific, which, due to gene-by-environment (GxE) interactions, is a limitation of any study that infers genetic effects from phenotypic variation [35]. By fostering similar phenotypes among plants in a common environment, phenotypic plasticity could have contributed to the large amount of unexplained variance in phenology among all populations in our study ( Table 2). In terms of future studies, it would be useful to know whether GxE interactions in agricultural environments providing irrigation or fertilizer can magnify or diminish phenotypic variation in ripening, which could either exacerbate or mitigate, respectively, the risks of genetic loss during single-pass harvests.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/land10101064/s1. Table S1: seed transfer zones (STZs) and populations used in our analysis, Table S2: criteria used to quantify reproductive phenology, Table S3a-e: heatmaps of reproductive phenologies for all STZs not shown in Table 4.