Study System
The study took place at the University of California, Donald & Sylvia McLaughlin Reserve in California’s inner North Coast Range. The study site consisted of serpentine grasslands, located in Morgan Valley (Lake County, CA 38.861°N, 122.408°W). Prior to the experimental treatments, the site was heavily invaded Aegilops triuncialis but had few additional non-native species (Aigner & Woerly 2011).
Study Design: The study was carried out within the context of a larger Aegilops triuncialis removal experiment (see Aigner & Woerly 2011). Focal 2m x 2m plots were arranged within 10 experimental blocks based on A. triuncialis density and randomly positioned within each block (Fig. 1). The experiment included five treatments: control, hand-pulling of A. triuncialis, mowing, fluazifop applied mid-season, and clethodim applied mid-season. Experimental treatments were carried out during the spring of 2008–2010. In 2011, after three years of treatment, A. triuncialis density was very low with a mean of 4.50% cover following restoration (Aigner & Woerly 2011) in the 40 pollinator plots with a consistent reduction in A. triuncialis for each of the A. triuncialis removal methods. For pollinator data, we pooled all A. triuncialis removal treatments together and compared them to the control.
Aigner and Woerly (2011) previously found that hand-pulling, grass-specific herbicides, and mowing restoration treatments relative to control treatments effectively decreased Aegilops triuncialis abundance as well as decreased the abundance of other species of non-native annual grasses that were present in small amounts. A. triuncialis removal treatments increased the frequency of native grasses and forbs (Aigner & Woerly 2011).
All plant cover data were obtained in 2011 and 2012 based on visual estimates using a 1 m x 1 m frame placed in the center of each plot. We estimated cover using square pieces of cardboard representing 0.1%, 1%, 5%, 10%, and 25% cover as references (see Aigner & Woerly 2011). Cover was estimated by two trained observers per plot in 2011 and one trained observer per plot in 2012. Within each treatment, we divided plots evenly between observers.
Floral Visitor Observations: For the pollinator observations, we used 50 plots with ten plots for the control treatment and 40 plots for the restoration treatments. Trained observers visited each plot and recorded plant-floral visitor interactions and floral visitor abundance for a gi ven time period. All observations were done during morning and afternoon on sunny days when pollinators were active. Each plot was visited for one-minute observation periods in 2011 and for two-minute observation periods in 2012 & 2013. We recorded all flower visits by flying insects during the observation period, recording the taxon of the floral visitor and the flower. Pollinators were collected as voucher specimens and later identified outside the observation period by Prof. Robbin Thorp. We identified pollinators to the level of morphospecies. We further classified floral visitors into the following functional groups based on existing literature: butterflies/moths, syrphids, bombyliids, non-bee mimic flies, short-tongued bees, medium to long-tongued bees, and kleptoparasitic bees. We conducted 17 observation sessions over a three-year period. In 2011, observation dates were April 16, 22, 27, and June 30. In 2012, observation dates were April 21, 22, May 2, and September 7, 8, and 11. In 2013, observation dates were April 12, May 22, 24, 31, and August 14, 15, and 28. This made for a total of 25 observer hours. A floral visit was recorded if an insect made contact with the reproductive parts of a flower. If a plot had no flowering plants present, pollinator visitation data were not collected.
Statistical Analysis
For all linear and generalized linear models, we scaled the variable year using the scale() function in R that centers a numeric matrix for each of the years, and we calculated plant and pollinator community variables at the plot level for each year. To determine whether Aegilops triuncialis removal would affect wildflower cover, we used linear models to examine if A. triuncialis cover, total cover for all annual grass species present, Lasthenia californica cover, and total cover for all native forb species visited by pollinators varied by treatment, year, and their interaction (e.g. Cover ~ Treatment*Year). For the wildflower cover analysis, we used all native forb species whose flowers were visited by pollinators (Table 1). We found no significant differences in pollinator community response between A. triuncialis removal treatments, so we pooled all A. triuncialis removal treatments together in our analyses.
Table 1
Native forb species observed being visited by pollinators. Taxonomic names follow Baldwin et al. (2012). Abbreviation used network diagram noted in second column. Navarretia pubescens does not appear in network diagram, as plant was visited by an unidentified bee, whose morphospecies could not be determined.
Species Name | Abbreviation | Family | Naming Authority |
Acmispon wrangelianus | LOWR | Fabaceae | (Fisch. & C.A. Mey.) D.D. Sokoloff |
Agoseris heterophylla | AGHE | Asteraceae | (Nutt.) Greene |
Calycadenia pauciflora | CAPA | Asteraceae | A. Gray |
Delphinium variegatum | DEVA | Ranunculaceae | Torr. & A. Gray |
Dichelostemma capitatum | DICA | Themidaceae | (Benth.) Alph. Wood |
Eriogonum nudum | ERGU | Polygonaceae | Benth. |
Gilia tricolor | GITR | Polemoniaceae | Benth. |
Grindelia camporum | GRCA | Asteraceae | Greene |
Hemizonia congesta | HECO | Asteraceae | D.C. |
Holocarpha virgate | HOVI | Asteraceae | A. Gray |
Lasthenia californica | LACA | Asteraceae | Lindl. |
Leptosiphon bicolor | LIBI | Polemoniaceae | Nutt. |
Lomatium hooveri | LOHO | Apiaceae | (Mathias & Constance) Constance & Ertter |
Microseris douglasii | MIDO | Asteraceae | (D.C.) Sch. Bip. |
Navarretia jepsonii | NAJE | Polemoniaceae | Jeps. |
Navarretia pubescens | NAPU | Polemoniaceae | Benth. (Hook & Arn.) |
Ranunculus californicus | RACA | Ranunculaceae | Benth. |
Trifolium albopurpureum | TRAL | Fabaceae | Torr. & A. Gray |
Trifolium fucatum | TRFU | Fabaceae | Lindl. |
Viola douglasii | VIDO | Violaceae | Steud. |
To test whether Aegilops triuncialis removal would affect pollinator diversity and abundance, we used a linear model to further examine if Shannon diversity for pollinators varied by treatment, year, and their interaction. After testing for overdispersion, we used a generalized linear model with a Poisson distribution to test whether pollinator morphospecies richness varied by treatment, year, and their interaction. To examine whether removing A. triuncialis improved habitat for ground-nesting bees, we used a generalized linear model with a negative binomial distribution to test whether ground-nesting bee abundance varied by treatment, year, and their interaction. We further tested whether the proportion of ground nested bees (ground nested bee abundance/total abundance) varied by treatment, year and their interaction using a generalized linear model with a gaussian distribution. We classified bees as ground-nesters based on information obtained from existing literature (LeBuhn 2013). We used a generalized linear model with a negative binomial distribution to test whether total pollinator abundance across all morphospecies varied by treatment and year, and their interaction.
To test whether changes in floral cover had an additional effect on pollinator diversity within restored plots, we used generalized linear models with Poisson distributions to test whether pollinator morphospecies diversity varied by preferred forb diversity, total preferred forb cover and Lasthenia californica cover respectively. Within restored plots, we further used generalized linear models with negative binomial distributions to test whether pollinator abundance varied by preferred forb diversity, total preferred forb cover and L. californica cover. All models were run using the ‘lme4’ package in R (Bates et al. 2015).
We used the adonis2 function in the ‘vegan’ package in R to compare pollinator community composition between the restored and control treatments using a permanova with a Bray-Curtis dissimilarity index with 999 permutations (Dixon 2003). To test for differences in community dispersion, we used the betadisper(_) function followed by permutest to run beta dispersion tests.
To test whether Aegilops triuncialis removal would affect plant-pollinator network structure, we calculated the observed NODF (nestedness) of the networks using the bipartite package (Dormann et al. 2009). To address the sensitivity of NODF to network size, we calculated a measure of NODF normalized to network size (NODFc) using the maxnodf package (Hoeppke & Simmons 2021). However, the small size of the control network warrants caution in interpreting these values for normalized nestedness as the number of links in the control network was below the minimum threshold.
To account for size and sampling effects, we further used Bayesian inference on a model that accounts for the relative abundance of species in the network, the effect of observation error, and the probability a true interaction exists (Young, Valdovinos & Newman 2021). This model included all interactions across the restored and control treatments since experimental plots were located in the same meadow and were therefore populated by species of the same community. We sampled the model 10000 times. For each sample we calculated two NODFc values for the interaction matrix among species observed in the control plots and the restored plots, respectively. Together, this procedure accounts for both size and sampling effects to more accurately compare nestedness between the control and restored treatments.
We then compared observed network niche properties to the ‘swap web’ null model to test for whether plant-pollinator network structure differed between restoration treatments using the ‘bipartite’ package in R (Dormann et al. 2009). The ‘swap web’ null model generates marginal totals identical to those observed and the same connectance as observed. We compared mean number of shared partners for both plants and pollinators, and niche overlap, as calculated by the Morisita-Horn Index, for both pollinators between observed and null models using z-scores. The Morisita-Horn Index calculates dissimilarity in plant species visited for pollinators and dissimilarity in pollinator visitors for plant species (Horn 1966). To account for temporal turnover in plant-pollinator interactions, we then repeated these analyses for data subsetted to include only the spring plant community (plant-pollinator interactions recorded in April and May) (see Supplementary Methods for further detail).
To test our fourth hypothesis that the most abundant forb species would strongly influence network structure, we further assessed individual contribution of each plant species to nestedness using the nestednesscontribution() function in the R bipartite package (Dormann et al., 2009). This function estimates the degree to which the interactions of each plant species increase or decrease community nestedness by comparing to a random null model that is designed to control for the effect of differences in degree. Networks included plant-pollinator interaction data from both spring and summer across three years. To account for differences in network size, we compared networks to null models via z-scores. To account for potential interaction turnover, we then repeated these analyses for only the spring plant-pollinator interactions of species temporally co-occurring with the abundant, generalist forb Lasthenia californica. For each network type across all seasons, we then tested whether mean plant species cover, as collected independently of the network data, explained plant species’ individual contribution to nestedness using a generalized linear model with a gaussian distribution.