Premise of the study: Understanding edaphic specialization is crucial for conserving rare plants that may need relocation due to habitat loss. Focusing on Eriogonum crosbyae, a rare soil specialist in the Great Basin, US, we asked how site-level variation among volcanic soil outcrops affected plant growth and population distribution.

Methods: We measured emergence, survival, size, and biomass allocation of E. crosbyae seedlings planted into soils collected from forty-two outcrops of actual and potential habitat. We also measured phenotypic variation in the wild, documented abiotic and biotic components of E. crosbyae habitat, re-surveyed Nevada populations, and evaluated occupancy changes over time.

Key results: Plants responded plastically to edaphic variation, growing larger and allocating relatively more tissue above-ground in soils with greater nutrient availability, and growing smaller in soils higher in copper in the field and the greenhouse. However, the chemical and physical soil properties we measured did not predict site occupancy, nor was plant phenotype in the greenhouse different when plants were grown in soils from sites with different occupation status. We observed occupation status reversals at five locations.

Conclusions: E. crosbyae performed well in soils formed on hydrothermally altered rocks that are inhospitable to many other plants. Extirpation/colonization events observed were consistent with metapopulation dynamics, which may partially explain E. crosbyae’s patchy distribution among outcrops of potential habitat. While soil properties did not predict site occupancy, early life stages showed sensitivity to soil variation, indicating that seedling dynamics may be important to consider for the conservation of this soil specialist.

Sampling locations

We assessed current occupancy status and collected and compared soil samples from a subset of 71 previously occupied and 90 unoccupied outcrops across the species’ range in Nevada. The light color and barren nature of E. crosbyae’s preferred volcanic habitat allows it to stand out in comparison with the surrounding landscape, such that potentially suitable habitat can be identified from a distance or using aerial photography. As a result, many potentially suitable sites have been surveyed and have had their occupation status described (Kaye, 1990; Morefield, 2003). Sampling locations for this experiment were chosen primarily from locations compiled in the Conservation Status Report for E. crosbyae (Morefield, 2003), and were last surveyed during (two sites) or prior to (thirty-three sites) 1993. Additional sampling locations were included based on survey locations documented in the Consortium of Intermountain Herbaria (USA) online repository (http://intermountainbiota.org). The four sites chosen for inclusion from this source by the methods below were surveyed between 1995-2010.

Potential survey locations were located on topographic map layers in Esri ArcMap (Redlands, California, USA) and a 150x180km (East to West x North to South) fishnet was overlaid, with the resolution set to ensure that twenty-six cells contained at least one occupied site. An occupied site was randomly selected from each cell for sampling. Where possible, the nearest known unoccupied site to each occupied location along the same aspect and landform was also chosen for sampling, with a 2km maximum separation distance, chosen based on the range of wind dispersal for comparable plumed seeds (Vittoz and Engler, 2007). Not every occupied site could be paired with a previously-surveyed unoccupied site that met these criteria, so two additional unoccupied sites were chosen using aerial imagery, along the same aspect and landform as occupied sites. Ultimately, accessibility issues, reversals of occupation status in several sampling locations, and discovery of one new population resulted in sampling of 25 occupied sites and 17 unoccupied sites (Appendix S3).

Soil properties and occupancy status (questions 1-3)—

We visited all 42 of our sampling locations during October of 2017 to collect soil and document their occupancy status. We collected twenty 10cm diameter x 10cm deep cylindrical soil cores from each site. We estimated that these cores would represent the rooting zone of E. crosbyae seedlings in the relatively shallow soil present at our sampling locations. In occupied locations, soil samples were taken from within 165cm of E. crosbyae individuals to account for potentially abrupt edaphic boundaries between suitable and unsuitable habitat that varied on a fine scale. Samples were bulked and homogenized for each site. From each composited sample, one 0.5L soil sample was taken and sieved to 2mm, then sent to A&L Western Laboratories (Modesto, California, USA) (http://www.al-labs-west.com/) for analysis. Soil properties quantified included extractable  concentrations of major soil cations, selected macro and micro-nutrients, saturation percent, salinity, organic matter content, pH, and soil texture (Table 1). Soil analysis methods followed A&L Western Laboratories (http://www.al-labs-west.com/) protocols.

Plant responses to soil variation (questions 4-5)—

Plant phenotypes in the field—

Occupied sites were visited at peak bloom during the summer of 2018 for vegetation data collection, including a 15-minute meandering walk to record presence of all associated plant species. We measured E. crosbyae size and growth characteristics thought to be potentially adaptive and observed to be variable in the field. These traits included: mat area (length of the longest side of the mat multiplied by the length of the mat perpendicular to the first direction), mat depth (measured from the top of the root crown to the estimated average height of the leaves), height of the tallest flower stalk (from the base of the flower stalk to the base of the inflorescence), and number of flower stalks. Twenty plants per site were measured, beginning at the lowest elevation in the occupied area and proceeding in a systematic zig-zag pattern wide enough to cover the entire sampling area, measuring the closest plant to the sampler after walking 10 paces (~4.5m) from the previous plant.

Seedling phenotypes in the greenhouse—

Seeds were collected in bulk from roughly 100 individuals of the Hycroft population on June 21, 2017 and cleaned by hand. Seeds were stored in a paper bag in the dark at room temperature (18°C) until planting. To determine how soil characteristics affected plant growth, we set up a soil preference experiment in the University of Nevada, Reno (Reno, Nevada, USA) greenhouses using soils collected at every sampled field location, plus a native field soil from Washoe Valley, Nevada (39.322771, -119.812121, 1549m elevation). The Washoe Valley soil is characterized as a Surpass gravelly sandy loam (coarse-loamy, mixed, superactive, mesic Aridic Haploxeroll) according to the NRCS Web Soil Survey (https://websoilsurvey.nrcs.usda.gov/app/WebSoilSurvey.aspx). This field soil was included as a standard for comparison, rather than a potting soil or other commercial soil mix, because it is the substrate we routinely use for growing a variety of Nevada native plants, including other species in the genus Eriogonum. Thus, we would expect any issues with seed quality or growing methodology to have been apparent if the seeds failed to grow in this standard soil.

We employed a randomized complete block design with 12 replicate containers per soil type. Rocks larger than 2.5cm were removed and replaced with an equal volume of field soil, and a 50/50 mix of field soil and washed decomposed granite was used to improve drainage in 8.9cm x 8.9cm x 7.6cm pots. Two seeds were planted into each pot on March 1, 2018 to increase likelihood of seedling establishment in all pots. Pots were watered lightly every day during germination, then 2-3 times per week as needed throughout the experiment. All pots were watered evenly when a majority appeared dry on the surface, to prevent desiccation of small seedlings and encourage germination. By watering a constant amount, rather than to a constant weight or field capacity, we allowed plant performance to be affected by the differences in soil texture and water holding capacity among our different soil types. Germination, survival, and phenology were monitored weekly.

Temperatures in the glass-ceilinged greenhouse were set to approximate minimum and maximum natural temperature limits, with no supplemental lighting. From March 1 to March 28, minimum and maximum temperatures were maintained between 1.5-15.5°C; from March 28 to April 8, between 4.5-18.5°C, and from April 9 to harvest, between 4.5- 21°C. Pots containing more than one seedling after several weeks of growth were thinned using a coin toss to determine which seedling to remove. Above and belowground biomass was harvested during the week of July 22, 2018, after approximately 4.5 months of growth. Roots were washed clean, dried to a constant weight in a drying oven, and weighed.

Data Analysis—

Soil properties and occupancy status (questions 1-3)—

Soil properties characteristic of typical E. crosbyae habitat were attained by calculating the mean of individual soil properties from sites that produced average plant growth in the greenhouse (defined as sites ± 1 SD from the overall means for biomass and root mass ratio; Table 1; Appendix S4). This excluded sites where the plants performed either unusually poorly or unusually well. We also tabulated changes in site occupancy status over time, comparing the results from our surveys to previous ones.

To quantify variation in soil properties between occupied and unoccupied sites, soil variation across all sites was summarized and visualized using Principal Components Analysis (PCA). In addition, 2-sided t-tests were used to test for differences in individual soil variables between occupied and unoccupied sites. For these analyses, all sites with extant populations of E. crosbyae at the time of sampling were considered occupied, while all sites without populations were considered unoccupied, regardless of whether occupancy status had changed since previous surveys. One-way analysis of variance was used to test whether emergence, total biomass, and root mass ratio (RMR; root weight/total biomass) of seedlings grown in the greenhouse varied based on field occupation status of the soil they were grown in; for these tests, we also included categories for recently colonized or abandoned sites. Tukey’s HSD post-hoc tests were performed on significant models (α < 0.05) for multiple comparisons.

Plant responses to soil variation (questions 4-5)—

An exploratory approach to variable selection and model generation was used to discern which soil characteristics had the greatest effects on plant characteristics in the field and total biomass, root mass ratio, emergence, and survival in the greenhouse. Soils from the more fertile Washoe Valley site were excluded during variable selection and model creation to improve comparability of greenhouse and field data. To reduce multicollinearity in measurements of plant characteristics in the field, we used PCA to identify uncorrelated responses with the strongest loadings on orthogonal axes. We chose two variables, mat area and mat depth, with Pearson correlation coefficients < 0.30. In the greenhouse, emergence and survival were recorded as the percentage of planted seeds per pot that emerged, and the percentage of un-thinned seedlings that survived from emergence until the date of harvest, respectively. For these analyses, site occupation was coded as a binary variable, with both occupied and colonized sites considered occupied, and all unoccupied and abandoned sites considered unoccupied.

Thirty-two possible soil variables were initially considered during variable selection (Table 1; Appendix S4). To reduce multicollinearity, we eliminated variables with Pearson correlation coefficients ≥ 0.70. To do this, out of each pair of highly correlated soil variables, the variable with highest deviance with the response of interest in univariate generalized linear models was removed, with this selection process proceeding separately for each response (mat area, and mat depth, and greenhouse total biomass, RMR, emergence, and survival). Additional thinning of soil variables was performed using the R package ‘randomForestSRC’, to facilitate final model selection (Breiman, 2001; Ishwaran et al., 2008; Ishwaran and Kogalur, 2019). Random Forest is a nonparametric machine learning technique that uses iterative decision trees to estimate the predictive capability of variables, and is an ideal method for variable selection. The top 16 variables with the highest RandomForest VIMP (variable importance) were chosen for inclusion. In all RandomForest runs for greenhouse response variables, block was shown to be unpredictive, and so this factor was omitted from further analyses.

Univariate and final multivariate models for the effects of soil variation on plant characteristics were constructed with gaussian distributions for measures of plant characteristics in the field, as well as for total biomass and RMR in the greenhouse, and with binomial distributions with a logit link function for percentages of emergence and survival in the greenhouse. For models of total biomass and RMR, plant age (days since emergence) was included as a potential covariate. Sites with four or fewer surviving seedlings were removed from analyses of the effects of soil variation on total biomass and RMR due to high variation in responses; sites excluded were two unoccupied sites, one colonized site, and one abandoned site. A Tukey Ladder of Powers transformation was performed on total biomass and RMR to help meet the assumption of normally distributed residuals, and transformed values were scaled to a mean of zero and standard deviation of one. All potential explanatory variables (except block) were similarly scaled.

Final model selection and multi-model averaging were performed using a genetic algorithm in the R package “glmulti” using AICc for model comparison. Coefficient estimates, in units of standard deviations from the mean for both soil variables and plant responses (Appendix S5), were calculated based on the subset of models which yielded 95% of the total evidence weight (generally within 2 AIC units of the top model). Comparisons of R², Adjusted R², and Allen’s (1971) predicted R²  (calculated using leave-one-out cross-validation) for the top models were made to test for model over-fitting, and model assumptions were checked using Normal Q-Q and Residuals vs. Fitted plots of top models.

Finally, we constructed linear mixed models to ask whether emergence, survival, total biomass, or root mass ratio of plants grown in the greenhouse differed based on site occupation status of the soils they were grown in, here with occupation status coded as one of four categories: occupied, unoccupied, colonized, or abandoned. Site was included in each model as a random variable nested within occupation status. Linear mixed models were constructed with gaussian distributions for total biomass and RMR, and generalized linear mixed models were constructed with binomial distributions and a logit link function for percentages of emergence and survival. Models were fit using maximum likelihood.

Please see 'ReadMe.txt' for descriptions of the contents of each worksheet within the dataset.