Fine-scale trait variation of five native forbs along environmental gradients

Site Selection and Plant Species

Our survey occurred on MPG Ranch (46.68000 N, 114.02778 W), which is a 6000 ha (14,826 ac) conservation property in western Montana. The landscape habitat types include riparian corridors, former and current agricultural areas, sagebrush steppe and intermountain grasslands, and montane mixed-conifer forests (Montana Natural Heritage Program 2017). This area of Montana is characterized by short, dry summers, cold winters, and cool, wet springs. Mean annual precipitation ranges from 300 to 350 mm (11.8–13.8 in; valley to mountain summits) and mean annual temperature of 7.6 °C (45.7 °F), ranging from –4.7 °C (23.5 °F, January) to 19.4 °C (66.9 °F, July; 30-y trend from a weather station in nearby Missoula, Montana; https://ncdc.noaa.gov). All survey sites occurred in the grassland habitat, which had never been cultivated and was dominated by cool-season bunchgrasses including bluebunch wheatgrass (Pseudoroegneria spicata (Pursh) Á. Löve [Poaceae]), Sandberg bluegrass (Poa secunda J. Presl [Poaceae]), Idaho fescue (Festuca idahoensis Elmer [Poaceae]), and the perennial forbs silky lupine (Lupinus sericeus Pursh [Fabaceae]) and arrowleaf balsamroot (Balsamorhiza sagittata (Pursh) Nutt. [Asteraceae]) (Mueggler and Stewart 1980).

View looking up at the general location of our sites to show the range of habitats.

We selected 5 native forbs based on their widespread occurrence and their relevance for restoration; blue-eyed Mary (annual; Collinsia parviflora Lindl. [Scrophulariaceae]), thread-leaf phacelia (annual; Phacelia linearis (Pursh) Holz. [Hydrophyllaceae]), cutleaf daisy (perennial; Erigeron compositus Pursh [Asteraceae]), blanketflower (perennial; Gaillardia aristata Pursh [Asteraceae]), and prairie smoke (perennial; Geum triflorum Pursh [Rosaceae]). We used elevation and insolation data to divide the area into 11 strata using a k-means clustering algorithm (Hartigan and Wong 1979; R Core Team 2020, stats package). Insolation was calculated as the megajoules of solar radiation each site received on the 2016 summer solstice (Doherty and others 2021); 2021 values showed minimum deviation relative to 2016 and were strongly correlated, (r = 0.99). This technique for modeling insolation considers both direct and diffuse irradiance for solar zenith angles throughout the day, as well as shadows cast by surrounding terrain and site slope (Corripio 2019).

Within this selected area, elevation ranged from 1130 to 1670 m (3707–5479 ft), and insolation ranged from 27.43 to 31.86 MJ m⁻² (2.55 MJ ft⁻²–2.96 MJ ft⁻²). Within each stratum we surveyed the focal forb species starting in mid-March 2021. Two 50 m² (538 ft²) sites were designated within each stratum where the greatest number of focal species were found, except in one stratum where only 1 suitable site was located (Figure 1). Sites within a specific stratum were typically located on separate ridgelines at least 300 m (984 ft) apart to increase the likelihood of independence. Sites in different strata along the same ridgeline were at least 145 m (475 ft) from other stratums. All 5 plant species occurred along the elevation gradient, but thread-leaf phacelia and cutleaf daisy were not found in the lowest insolation sites. Our study sites spanned 2 Generalized Provisional Seed Zones (Bower and others 2014): 1) 2 sites in minimum temperature class –9.4 to –6.7 °C (15–20 °F), and annual heat moisture index 2 to 3, and 2) 19 sites in minimum temperature class –9.4 to –6.7 °C (15–20 °F), and annual heat moisture index 3 to 6.

Measurements of Abiotic and Biotic Site Conditions

At each site, we measured soil temperature, soil moisture and nutrients, and site productivity. Soil temperature was monitored every 90 min from May (5th, 6th, 11th, or 17th depending on site) to 20 or 22 July with iButton temperature loggers (Thermocron, Sydney, Australia) placed 5 cm into the soil. Minimum, maximum, and mean monthly temperatures were calculated for each site. On 6 to 7 June, we collected 5 soil cores (0–10 cm [0–0.39 in]) within each site, which were pooled, sieved through a 2 mm (0.08 in) sieve, and placed in resealable plastic bags in a cooler to bring to the laboratory for further processing. We subsampled 100 g (0.03 oz) that was air-dried and then sent for analyses of soil pH, percent organic matter, available nitrogen (NO₃), phosphorus (Merlich), potassium, calcium, sulfur, magnesium, manganese, iron, sodium, zinc, and copper (Ward Laboratories, Kearney, Nebraska, USA). Ten g was dried at 105 °C (221 °F) for 48 h to determine gravimetric water content, and 10 g (0.35 oz) was freeze-dried for molecular analyses of fungal communities (see below). Gravimetric water content was also measured 12 to 13 July following the same protocols. During the July sampling, we also collected 2 samples of aboveground biomass using a 0.5 × 0.5 m (3.3 × 3.3 ft) frame as a measure of site productivity. Biomass was clipped at the soil and dried at 65 °C (149 °F) for 48 h before being weighed.

Figure 1

Overview map of the 21 sites included in this study.

Plant Traits Measurements

Every week we visited sites to assess phenology. We recorded the time of flowering and seed ripening based on the prominent growth stage of the entire site (not individual plants). At each site, we randomly selected 2 individuals of each species to measure plant traits. We counted the number of flowers and then cut the shoots for measurements of specific leaf area (SLA) and aboveground biomass. While this involved destructive sampling that precluded seed collections on the same individuals, it was necessary to capture peak biomass before shoots started to senesce. Shoots were wrapped in wet paper towels, placed in a cooler, stored at 4 °C (39.2 °F), and processed within 2 d. Two leaves were removed per plant at a consistent stem location for SLA analysis to reduce the variation of leaf morphology within a plant (Shipley 1995). Leaves were flattened onto a sheet of white printer paper and scanned at 200 dpi using a scanner (Dell H825cdw Printer/Scanner) according to Pérez-Harguindeguy and others (2013). Leaf area was calculated in R using k-means algorithm (Hartigan and Wong 1979) and the HDBSCAN algorithm (McInnes and others 2017) for image processing. Specifically, we assigned each pixel to a white background group or leaf foreground group with k-means. We removed the group corresponding to white background pixels from further analysis. To isolate each of the 2 leaves in the scanned images, we applied HDBSCAN to the remaining leaf pixels. This method resulted in 2 clusters, one corresponding to each leaf. We counted the pixels in each leaf cluster and multiplied the counts by the area of a pixel (1.6e-8 meters-squared in a 200 DPI scanner) to derive area for each leaf. SLA was then calculated by dividing leaf area by leaf biomass after drying leaves at 65 °C (149 °F) for 48 h. The remainder of the shoot was also dried at 65 °C (149 °F) for 48 h. Values from the 2 plants per species were averaged within sites. SLA was not measured on threadleaf phacelia because plants began to dry down due to a heat wave before collection and we could not get accurate leaf area scans.

We collected ripe seeds from each species to measure seed mass. We collected from a minimum of 10 individuals at each site to capture seed weight diversity and to ensure that we had enough seeds for future reciprocal transplant studies. At a few sites, blue-eyed Mary seeds were collected farther than 10 m (32.8 ft; but not more than 30 m [98.4 ft]) away from the site center given the low seed numbers per plant. Seeds were stored in paper envelopes at room temperature to after-ripen, and we measured seed mass in September 2021 (Sartorius scale, precision = 0.1 mg [3.5 e⁻⁶ oz]). For larger-seeded species (blue-eyed Mary, blanket flower, and prairie smoke), we weighed 5 seeds individually from 5 plants for a total of 25 seeds per species per site. Smaller-seeded species (cutleaf daisy and threadleaf phacelia) had seeds that were too small to be weighed individually, so 5 seeds were weighed together and the average individual seed mass was calculated. In total, 5 replicates of 5 seeds from 5 plants were weighed for 25 replicates (125 seeds) per species per site.

Amplification of Fungal Communities and Assessment of Fungal Guilds

We chose to focus on fungi because they constitute the most important group of soil mutualists and pathogens (Smith and Read 2010; Raaijmakers and others 2009). Briefly, we extracted DNA from 250 mg (0.009 oz) soil per site using the DNeasy PowerSoil Pro DNA isolation kit (Qiagen, Redwood City, California) as well as a blank to monitor potential background or cross-contamination among samples. Amplification, sequencing, and bioinformatics followed Bullington and others (2020). Briefly, we amplified the internal transcribed spacer 2 (ITS2) region using the fungal-specific forward primer pairs fITS7 and ITS7o and the general eukaryotic reverse primer ITS4 and attached barcodes in a second PCR. Amplicons were pooled based on band intensity and sequenced using 2 × 300 paired-end run on an Illumina MiSeq sequencing platform (Illumina Inc, San Diego, California, USA). Raw sequence data were processed using Quantitative Insights into Molecular Microbial Ecology 2 (QIIME2 version 2018.4; https://qiime2.org/), and forward and reverse reads were trimmed to where median quality score fell below 30. We used the q2-dada2 plugin to produce sequence clusters with 100% similarity. We then clustered the sequences at 97% similarity using the q2-vsearch plugin to produce operational taxonomic units (OTUs), and taxonomy was assigned using the database UNITE (Nilsson and others 2019). Taxa were placed into fungal guilds based on potential function (plant pathogens and plant mutualists) with the database FUNGuild (Nguyen and others 2016). To represent plant pathogen communities, we used OTUs identified solely as “plant pathogens” within “highly probable” and “probable” confidence rankings, whereas mutualists were identified as any sequence of Glomeromycota origin. Samples were rarified to 3500 sequences before statistical analyses, in which plant pathogens made up 17% and mutualists 4% of total fungal sequences. Sequences have been submitted to the National Center for Biotechnology Information’s Sequence Read Archive under accession number PRJNA983821.

Statistics

To determine if abiotic and biotic site characteristics or plant traits varied along elevation and insolation gradients and if a significant interaction occurred between elevation and insolation, we conducted two-way ANOVAs (type II) on multiple regression models with insolation and elevation as predictors. We anticipated opposing effects of elevation and insolation on drought conditions, with lower elevation and higher insolation corresponding to the strongest drought conditions. For this reason, we specified linear models with the inverse of insolation to test for multiplicative interactions among drivers of drought. However, we discuss results using the original insolation definition wherein increased insolation equals increased solar radiation and decreased topographical shade. Flower count and phenology were also analyzed using a two-way ANOVAs (type II) on multiple linear regression models but required a generalized linear model with a Poisson distribution. Normality for all site and plant trait data was determined using Shapiro-Wilk tests and visually using QQ-plots. Variables were log-transformed when necessary. Statistics for each species were run independently. Given strong pairwise correlations (>0.70, assessed using Pearson’s correlation), we removed gravimetric water content for July (correlated with gravimetric water content for June), and soil temperatures for May and July (correlated with June temperatures, which had a more complete data set for the entire month). There were no significant differences in phenology within species among sites, potentially attributable to measurements being too coarse (only once per week), and these data are therefore not presented. To assess if certain species exhibited greater trait variation and if certain traits were more variable than others, we relativized absolute trait values around the mean for each species and performed one-way Levene’s tests. We conducted Tukey adjusted post-hoc comparisons to reveal the nature of group differences. Statistical analyses were carried out using R.4.0.2 (R Core Team 2020), the car (v3.0.8; Fox and others 2012), and the hmisc (v4.4.0; Harrell 2020) packages.

What Is the Variation in Soil Properties along Gradients in Elevation and Insolation?

Soil moisture and temperature, both known to affect plant phenology and growth (Reid 1973; Sundqvist and others 2013; Kunkel and others 2016), varied by orders of magnitude across our site. For example, June soil moisture ranged tenfold, minimum soil temperature ranged nearly threefold, and maximum temperature nearly twofold (Table 1). This corresponded with variation in other soil variables, including a sixfold difference in soil organic matter, a tenfold difference in available N, a threefold difference in S, Z, Mn, Ca, and a twofold difference in K, Fe, Mg, Na, and P (Table 1). Finally, site productivity differed eightfold among sites, and we documented more than a sixfold difference in the relative abundance of putative fungal mutualists and pathogens (Table 1).

We then determined if and how this substantial variation changed predictably along elevation and insolation gradients and found that both were suitable proxies for at least some of the abiotic and biotic properties measured. Soil moisture decreased with increased insolation in accordance with previous findings (Reid 1973; Sundqvist and others 2013) but did not change along the elevational gradient (Table 1). The availability of several nutrients (Ca, K, Mg, Fe) also decreased with increasing insolation, although the effect of insolation depended on elevation for three of them (Ca, K, Fe) (Table 1). We also observed a significant interaction between insolation and elevation for productivity, most likely because the effect of insolation was greater at lower elevation sites. Soil moisture has been shown to drive differences in soil fertility and productivity (Daly and Porporato 2005; Meisner and others 2013), which may explain the responses to insolation we observed. Elevation and insolation also affected soil temperature in predictable ways; minimum June temperature decreased with increased elevation (P = 0.002), and maximum temperature increased with insolation (P = 0.004). At one high insolation site with extensive rock cover and low productivity, soil temperature reached 66 °C (150 °F). Even with these extreme conditions, we were able to locate 4 of 5 target species at this site. Similarly, high soil temperatures have been recorded in eroded patches in the Alps with plants exhibiting heat resistance exceeding 60 °C (140 °F) (Larcher and others 2010), suggesting that high-elevation plants may be adapted to large temperature ranges. High-elevation plants may also associate less with AM fungi at our site, because the relative sequence abundance belonging to Glomeromycota decreased with increasing elevation (P = 0.008, although this was primarily driven by high abundances at the 2 lowest elevation sites (Table 1). Previous work has observed greater proportions of non-mycorrhizal plants in the Arctic, possibly because host benefits from AM fungi decrease in colder soils (Kytöviita 2005). Putative soil pathogens, by comparison, did not change predictably along elevation or insolation gradients, even though their abundance varied greatly within this area (Table 1).

Figure 2

Trait variation by species. Absolute values for each trait were relativized around the mean for each species (panel A). The median is represented by the line across the box and dots represent outliers. After relativizing values, we tested for homogeneity of variance across species (panel B) and trait (panel C). Groups with larger values on the x-axis had greater variance. Post hoc testing (letters left of boxes in B and C) showed blue-eyed Mary had greater trait variance than all other species except threadleaf phacelia, which was among the most widely distributed in our study. Flower count and shoot biomass appear to have the largest amount of variation among traits. We did not measure SLA on threadleaf phacelia.

Overall, elevation and insolation were joint—possibly multiplicative—drivers of the soil environment in our grassland study area. However, where only main effects were significant predictors, insolation effects were more abundant, notably for water content. At this time, these results are limited to this system and additional studies would be necessary to determine if this is a local, regional, or more widespread pattern. Furthermore, we did not test for nonlinear relationships among soil properties and environmental gradients, though we encourage investigating this possibility with a larger data set capable of resolving these complex patterns. For future work studying the influence of soil properties on ecological phenomena, such as trait-by-environment interactions, we advise including both insolation and elevation in the sampling strategy and analysis.

Do Plant Traits Vary along Fine-Scale Elevation and Insolation Gradients and Is Variation Trait Species-Specific?

Despite the large variation in environmental conditions outlined above, we found few consistent patterns in trait differentiation. For example, no measured traits in threadleaf phacelia or prairie smoke varied along insolation and elevation gradients, although for threadleaf phacelia this could be caused by its more restricted distribution to only half of the sites surveyed. Absence of trait relationships with gradients may be attributable to extensive seed dispersal across the landscape obscuring adaptive responses and (or) an inability to respond to variations in environmental conditions (that is, low plasticity). Also possible is that even greater environmental variation is required to prompt a trait response, or that relationships with the environments are nonlinear, which we could not assess here given sample size constraints. Additional fitness parameters important to plant establishment could be analyzed, including seed viability and germination rates.

Cutleaf daisy (Erigeron compositus) at a high insolation and high-elevation site.

For the other 3 species, we observed responses in selected traits that may be explained by the preferred microhabitats of each species within the heterogeneous grassland landscape. Blanketflower SLA was smaller in high insolation areas (Figure 2; Table 2). Smaller SLA is associated with slower growing leaves that have increased thickness and density with less transpiring leaf surface and therefore reduced water loss in drier soils (Shipley and others 2006; Poorter and others 2009). A similar trend was observed in cutleaf daisy but was primarily driven by the 2 low-elevation sites that had lower insolation. SLA in these 2 Asteraceae species also changed with elevation, but responses differed and for cutleaf daisy depended on insolation (Table 1); SLA increased with elevation in blanketflower and decreased in cutleaf daisy. The underlying reason for this difference is uncertain but may be related to their preferred habitats. Cutleaf daisy tended to grow on exposed, high insolation sites including high-elevation ridgelines, while blanketflower preferred the more north-facing, lower insolation, high-elevation slopes. Cutleaf daisy also produced larger seeds at higher elevation sites that were drier and hotter (Table 2). This increase in seed size may be a response to stress because larger seeds contain more resources that can produce faster-growing, larger seedlings with higher colonizing ability (Westoby and others 1992; Coomes and Grubb 2003). Blue-eyed Mary had more flowers at higher elevation sites that tended to have increased shrub cover.

Blanketflower (Gaillardia aristata) at a high insolation and low-elevation site.

Overall, traits differed in variability and so did plant species (Figure 2), but the drivers of this variability remain unknown. SLA and seed weight, both heritable traits (Westoby and others 1992; Leishman 2001; Shipley 1995; Morecroft and Woodward 1996; Scheepens and others 2010) correlated with gradients in elevation and insolation in some species in ways that support previous findings (Qi and others 2015; Poorter and others 2009). We do not know if this trait variation is related to local adaptation or plasticity as it requires reciprocal transplant experiments (Gibson and others 2016) and longer-term monitoring of survival and fitness of seeds collected along elevational and insolation gradients.