A history of species co-occurrence in plant communities is hypothesized to lead to greater niche differentiation, more efficient resource partitioning, and more productive, resistant communities as a result of evolution in response to biotic interactions. A similar question can be asked of co-occurring populations: do individual species or community responses differ when communities are founded with plants sharing a history of population co-occurrence (sympatric) or originating from different locations (allopatric)? Using shrub, grass, and forb species from six locations in the western Great Basin, USA, we compared establishment, productivity, reproduction, phenology, and resistance to invaders for experimental communities with either sympatric or allopatric population associations. Each community type was planted with six taxa in outdoor mesocosms, measured over three growing seasons, and invaded with the annual grass Bromus tectorum in the final season. For most populations, the allopatric or sympatric status of neighbors was not important. However, in some cases, it was beneficial for some species from some locations to be planted with allopatric neighbors, while others benefited from sympatric neighbors, and some of these responses had large effects. For instance, the Elymus population that benefited the most from allopatry grew 50% larger with allopatric neighbors than in single origin mesocosms. This response affected invasion resistance, as B. tectorum biomass was strongly affected by productivity and phenology of Elymus spp., as well as Poa secunda. Our results demonstrate that while community composition can affect plant performance in semi-arid plant communities, assembling communities from sympatric populations is not sufficient to ensure high productivity and invasion resistance. Instead, we observed an idiosyncratic interaction between sampling effects and evolutionary history, with the potential for seed source of individual populations to have community-level effects.

Target species were selected after conducting plant surveys in 21 locations in western Nevada, eastern California, and southeastern Oregon, with the goal of identifying the most commonly co-occurring plant community dominants. We selected six of the most common native plant species representing a variety of life forms. Species included three perennial grasses (Elymus spp. L., Achnatherum thurberianum (Piper) Barkworth, Poa secunda J. Presl), two deep-rooted, long-lived shrubs (Artemisia tridentata Nutt. and Ericameria nauseosa (Pall. ex Pursh) G.L. Nesom & Baird), and one annual/biennial forb (Chaenactis douglasii (Hook.) Hook. & Arn.). We included sites where Elymus elymoides (Raf.) Swezey and E. multisetus M.E. Jones co-occur and likely hybridize, given their shared distribution in the western Great Basin (Barkworth et al., 2007). For simplicity, we refer to Elymus taxa jointly as a “species”, even though Elymus collections may include two hybridizing species; the Austin site was the only location where either no E. multisetus was observed in the field and all seed-grown plants were confirmed as E. elymoides, which is consistent with the more western Great Basin distribution of E. multisetus.

We identified six collection sites within the Cold Deserts ecoregion where all focal species occurred (Omernik & Griffith, 2014). Sites were primarily mid-elevation, sagebrush steppe communities, with average annual precipitation between 237 and 384 mm and elevation ranging from 1395-2055 m (PRISM Climate Group, 2004). We bulk-collected mature seeds for each species in each location from a minimum of 50 individual plants, with a minimum of 1 m separating individuals, between 1 June and 15 December of 2015. We prepared 90 mesocosms (200L) outside of the Valley Road Greenhouse Complex at the University of Nevada, Reno (39.537924, -119.804757). Each mesocosm was filled with locally sourced topsoil to 0.9m depth (soil surface area = 0.25m², soil volume =0.3m³). Seeds were examined on a light table to ensure seed fill, and planted in November 2016, following the same orientation and planting design in each mesocosm. Each mesocosm included two replicates of Elymus spp. and A. thurberianum, six replicates of P. secunda, and one replicate of C. douglasii, E. nauseosa, and A. tridentata to approximate a dense natural community, and two seeds were sown into each planting position to increase the chances of seedling emergence. In cases where seedlings from both seeds established, each planting position was thinned to one seedling based on a random coin flip. To augment direct seeding, plants of each species and site were also grown in the greenhouse to serve as replacements for plants failing to establish directly from seed. By species, this was 45.5% of Elymus spp., 76.7% of A. thurberianum, 73.3% of P. secunda, 100% of C. douglasii (an annual or biennial, so all were transplanted at least once over the three-year period), 83.3% of E. nauseosa, and 100% of A. tridentata (a small-seeded shrub that necessitated transplanting to ensure controlled placement) were transplanted.

We created two treatment types, sympatric or allopatric, with six different combinations for each treatment. As discussed above, the terms sympatric and allopatric are often used to describe species, rather than populations, but we follow the terminology of Grady et al. (2017) and use these terms to describe population history. In the sympatric treatments, seeds of all six species from one of the six locations were sown together. In the allopatric treatments, seeds from each species and location were randomly assigned to one of six allopatric mixes, with one representative species from each location in each allopatric mix. This led to a total of 12 unique communities, which were each planted with seven or eight replicates, for a total of 90 mesocosms.

Mesocosm soil was initially watered to maximum water holding capacity and lightly watered once weekly if no natural precipitation had occurred. Each mesocosm was monitored for seedling emergence from November 2016 to May 2017, and seedlings were thinned as needed. In May 2017, for any plants that failed to recruit from seed, individuals were transplanted from seeds sown in the greenhouse. At the end of August 2017, when many of the perennial grasses had become senescent, every plant was assessed for height, crown size (length ✕ width), senescence index (a visual estimate using a continuous integer between 0-3; 0 = no live tissue and 3 = >75% green tissue), and the number of inflorescences. Starting in late September 2017, the second growing season, a new round of monitoring began following several rainstorms. We tracked new leaf growth and survival through May 2018. As in the previous year, in May 2018, we replaced any dead individuals. From May to October 2018, we noted whether plants had green tissue (0/1) on a weekly basis and used this information to estimate the number of days a given plant had potentially photosynthetically active tissue during the growing season. At the end of August 2018, each plant was assessed for crown size, height, senescence index, and the number of inflorescences.

In October 2018, each mesocosm was invaded with 130 B. tectorum seeds to test invasion resistance, with seeding density based on a field survey of seed production in a moderately invaded site. Seeds were lightly raked into the soil surface, in interspaces between plants. All plants were grown through August 2019 and watered as needed as in previous years. In August 2019, the aboveground biomass of all plants was harvested, oven-dried, and weighed.

Annual productivity was estimated from measurements of native plant volume in years 1 and 2 (calculated as crown size multiplied by plant height) and from aboveground biomass in 2019, and reproductive output was estimated as the number of inflorescences per plant. Overall survival was measured from the number of transplants used to replace dead plants required at any planting location in a mesocosm over all three growing seasons. The total mass of B. tectorum per mesocosm was used as the indicator of invasion resistance. For species that have senescence as part of their life history strategy (all but A. tridentata), phenology was evaluated in two different ways: 1) the total number of days that a given plant had green tissue from fall through spring (September through May), or 2) the senescence index measured in August.

We first used community composition models to analyze overall community composition and plant survival in each year, using linear models implemented in R version 4.1.0 (R Core Team, 2021). These models included species, community, and a species x community interaction as fixed effects with survival, plant volume, or biomass as response variables. Response variables were transformed as needed to better fit model assumptions. F statistics and p-values were calculated using the R package ‘car’ (Fox & Weisberg, 2019). The goodness of fit (R2-values) was calculated with the R package ‘MuMIn’ (Barton, 2020) for this and all glm models and post-hoc contrasts were generated using the Tukey honest significant differences method with the ‘agricolae’ R package (de Mendiburu & Yaseen, 2020). Residual histograms and scatterplots of predictors and response variables were created to verify model assumptions of normality, homogeneity of variances, and linearity of relationships.

Q1: Do plants perform differently when grown with sympatric or allopatric neighbors, and do responses differ among species or populations?

We created a main effects model asking whether individual plants had responses to sympatric or allopatric neighbors, using linear models with species, collection site, and treatment (allopatric or sympatric) as fixed effects. For grass species, we used linear mixed models that included mesocosm as a random (intercept) effect, as they had multiple replicates in a mesocosm (Bates et al., 2015). Response variables included survival, size, reproduction, and phenological variables. Based on our previous work, we expected considerable variation among species, thus we also built species-specific models asking whether there were species- or location-specific differences in response to sympatric or allopatric neighbors. First, we first created linear mixed models for each species that included survival as response and collection site, treatment, and the interaction between site and treatment as fixed effects. Then, we asked whether surviving plants of each species differed in size, phenology, or reproduction in allopatric or sympatric communities using the same model structure described above, but also included plant age (days since transplanting) as a covariate; response variables were transformed as needed.

Q2: Do community-level responses differ among allopatric or sympatric treatments?

We then quantified whether community-level responses (total productivity, total survival, inflorescence production, or B. tectorum suppression) differed among unique communities or by allopatric/sympatric treatment using linear models on the aggregated community level values per mesocosm. Total responses were calculated by summing response values across all individual plants in each mesocosm. Models included treatment (allopatric/sympatric) and unique community (one of the 12 allopatric or sympatric combinations, nested within treatment) as fixed effects, and separate models were created for each response variable and each year of measurement, with transformations as needed. Means were again compared using the Tukey honest significant differences method for multiple comparisons.

Q3: Can we identify factors that can predict differences in invasion resistance?

After assessing whether communities differed in key performance metrics, we sought to identify factors that might explain the observed differences in B. tectorum invasion resistance. For this, we asked whether any individual species was having an outsized predictive effect on B. tectorum biomass. We took a tiered model building approach. First, we considered potential explanatory variables (survival, inflorescence production, volume, senescence level of native plants prior to invasion, and the number of live green days experienced during the fall-winter of invasion) separately for each species. For grass species, we summed the values across the replicates within a community to represent the total performance of a species. We ran a set of generalized linear models separately for each metric of performance and each species, considering whether any individual response for an individual species predicted B. tectorum biomass, log transformed for analysis. We then retained a subset of variables for a multiple regression model, which included the most impactful explanatory variables (p < 0.05) from the previous analyses, further reducing as needed to ensure no variable in the model had a VIF greater than 3 (Fox & Weisberg, 2019). Finally, we retained significant variables to create a final multiple regression model of species-level responses that best explained B. tectorum biomass.