The fragment sites exhibited marked reductions in bee species richness and assemblage evenness relative to reserve sites , making this an excellent system in which to examine the extent to which ecological filtering causes the restructuring of assemblages. Our dataset also possesses a number of other desirable properties. First, since the original intent of data collection was to examine the effect of habitat fragmentation in isolation from other effects of urbanization-induced landscape change , we selected study plots representative of intact CSS flora to the extent feasible. Accordingly, reserves and fragments did not differ with respect to the plot-level species richness or composition of native insect-pollinated plants . Second, study plots exhibited no spatial autocorrelation with respect to the species composition of native bee assemblages, minimizing the potential for patterns in the spatial arrangement of study plots to drive patterns of bee taxonomic or functional distribution. Lastly, reserves and fragments did not differ with respect to overall bee abundance ; minimizing the potential for sampling effects to contribute to any differences that we detect between reserve and fragment plots with respect to bee taxonomic or functional composition. Functional trait assignment and analyses of functional diversity: Every identified bee species was assigned a category with respect to each of the following natural history traits: lecty , large plastic pots for plants nest location, nest building behavior, sociality, body size, and flight season.
Table 2-1 lists each trait, how it is analyzed in functional diversity models , and its method of assignment; Table 2-S1 lists the bee species and their associated traits. With the exception of body size and flight season, all traits were assigned to individual species using literature syntheses for the species or species group and subgenus , as well as revisionary publications on lower taxa and our own field observations. Due to the lack of data for many species in our region, we also relied on phylogenetic inference when such data are available and appropriate: e.g., all Lasioglossum species were scored as polylectic, ground-nesting, actively constructing nests, and eusocial . Cleptoparasitic bees are included in all analyses, although excluding them does not qualitatively alter our results. Cleptoparasites are classified to a unique lecty category , the same nest location as their presumed hosts, and always as nest renters. Although some species of Sphecodes in our system may be social parasites of eusocialHalictini species , we classified all cleptoparasites as solitary because they do not exhibit reproductive division of labor . Two traits, body size and flight season, were not assigned based on published data. To estimate body size, we measured the intertegular lengths of three haphazardly selected females of each solitary species or four haphazardly selected females of each eusocial species, when possible. For species for which no females were collected, we measured the intertegular lengths of males. To assign flight season, we assembled an individual-level database of collection dates from our own field data and the database of the University of California, Riverside Entomology Museum.
Since climatic conditions may drive intraspecific variation in the timing and duration of bee flight seasons, we including only specimen records from south of 36.00° latitude, within 80 km of the Pacific coast, and below an altitude of 1000 m. From this database, we performed 1,000 random subsamples of 30 specimens per species and calculated the tenth and ninetieth percentile collection dates . We then scored each bee species with respect to whether or not they are active in the early, middle, and late part of our study period using these percentiles. Rather than assigning quantitative measures of flight season duration and median flight date , our approach of assigning flight season as presence-absence in coarser-grained season categories minimizes biases inherently present in most databases, such as non-uniform distribution of sampling dates and non-random captures of different taxa by different collectors. This approach also minimizes the impact of the sampling effect, in which rarer taxa in the database are likely to be recorded as having shorter flight seasons . For our metric of functional diversity, we chose functional dispersion , a widely used metric that can provide insight into how native bee functional diversity responds to anthropogenic alterations of natural habitat . FDis is calculated as the mean distance of each species from its site level, multivariate centroid of functional traits . Thus, FDis is mathematically independent of species richness and can take into account differences in the relative abundances of functional trait combinations . These attributes make FDis relatively insensitive to rare species and functionally equivalent species; thus FDis particularly well-suited for our dataset, in which we have detected strong differences between reserves and fragments in both species richness and assemblage evenness Statistical analyses: Except where noted otherwise, we analyzed data from 2011 and 2012 separately because of differences in sampling location and frequency.
In order to avoid biases associated with variation in aerial netting efficacy in different terrains, bees collected by aerial netting were excluded from our main analyses; however, inclusion of these netted specimens in our analyses yielded qualitatively similar results . We also repeated all analyses at the genus level using genus-level mean or modal averages for functional traits to ensure that particularly species-rich genera did not disproportionately influence our findings. All analyses were conducted in R version 3.3.1 . Taxonomic and functional alpha diversity: To assess plot-level alpha diversity, we calculated Shannon-Weiner diversity H and FDis, where each diversity metric was weighted by relative abundances of each species. To account for variation in the number of bees sampled per plot, we calculated both diversity metrics after rarefying our data to the lowest plot-level bee abundance recorded each year . Diversity metrics in reserves and fragments were compared with two-sample t-tests. Additionally, we constructed a linear mixed-effects model , lmerTest , and MuMIn to examine the relationship between taxonomic and functional diversity. In this linear mixed-effects model, data from the two years were combined; functional diversity was the dependent variable, Shannon-Weiner diversity was the independent variable, and study year, study plot identity, and habitat category were included as random effects. While a direct comparison of diversity metrics provides information on how habitat categories differ from one another, assessing whether observed differences are driven by stochastic or deterministic processes requires testing null models . Here, we generated random bee communities for each study plot that have species richness and Shannon-Weiner diversity equivalent to their respective observed communities in order to test whether observed differences in functional diversity are simply due to underlying differences in species richness . In this analysis, we generated random communities for each study plot by first permuting observed species-level individual abundances across all species within the species pool for the study year in question , and then rarefying each randomly permuted community to the lowest plot level bee abundance recorded each year . This permutation procedure resulted in 10,000 random communities for each study plot in each year. We then assembled 100,000 datasets by randomly selecting one permuted community from each study plot, compared the FDis scores of reserve and fragment plots via two-sample t-tests, and extracted the test statistic of each comparison. Finally, we compared the test statistic of our empirical dataset against the null distribution of test statistics to assess the frequency with which null datasets yielded FDis differences between reserves and fragments that equal or exceed those observed in our empirical dataset. Taxonomic and functional beta diversity and assemblage composition: To assess spatial beta diversity among plots, we used analyses of multivariate dispersion , plastic pot plant containers which compare habitat categories with respect to the degree of compositional similarity among their constituent study plots. Multivariate dispersion is calculated by first performing non-metric multidimensional scaling ordinations based on all combinations of pairwise among-plot dissimilarity in bee taxonomic or functional composition , and then comparing the non-metric distances of plots from the centroids of their respective habitat categories via a permutation test .
In our analysis of taxonomic beta diversity, pairwise among-plot dissimilarity was calculated as the abundance-weighted Bray-Curtis distance between each pair of plots. In our analysis of functional beta diversity, we first calculated the coordinates of each plot’s abundance-weighted functional centroid in multivariate trait space , and then calculated pairwise among-plot dissimilarity as the non-metric distances between these functional centroids. We performed 10,000 permutations for calculations of both functional and taxonomic beta diversity. In addition to examining beta diversity among plots within each habitat category, we also assessed whether reserves and fragments differed from each other with respect to the taxonomic and functional assemblage composition of their bee faunas. To accomplish this comparison, we performed permutational multivariate analyses of variance on the same pairwise among-plot dissimilarity scores used to calculate beta diversity, described above, with 10,000 permutations . Lastly, we also performed a Mantel test to examine the relationship between taxonomic and functional composition, based on the same pairwise distance matrices discussed above. Unbalanced designs such as that used in year 2012 of our study are known to introduce bias into PERMANOVA tests when within-group heterogeneity is unequal among habitat categories . Thus, to aid in the interpretation of our results, we also performed tests of beta diversity and assemblage composition for the 2012 dataset on random subsamples of 6 fragment plots and examined the proportion of results in which the findings of the subsamples agreed with those of the full dataset with all 11 fragment plots included.To assess the drivers underlying differences between reserves and fragments with respect to taxonomic and functional diversity and composition, we performed two additional analyses. First, we used two-sample t-tests to compare bee assemblages in reserve and fragment plots with respect to the relative representation of each functional trait. In this analysis, we used the plot-level mean average for intertegular length , and the plot-level proportional representation by each categorical or binary state for all other traits. Proportion data were logit-transformed prior to analysis as recommended by Warton and Hui , and all calculations were weighted by the relative abundance of each species. Second, we performed an indicator species analysis to identify bee species or functional groups associated with each habitat category. We used the Indval.g association index in the indicator analysis to account for the unbalanced sampling design in 2012. To assign bee species to functional groups, we constructed a dendrogram of all bee species collected in the study using hierarchical clustering based on functional trait data . We used Ward’s algorithm to perform hierarchical clustering , and assigned bees into 25 functional groups based on their positions in the dendrogram. Functional group membership of each species is given in Table 2-S1. Geographical range sizes: For bee taxa identified to described species, we calculated their geographical range size based on our field data and the database of the Bee Research Laboratory of the United States Department of Agriculture. This database includes specimens collected from throughout the United States, Canada, and Mexico; and represents one of the most comprehensive and unbiased databases of bees in our study region. These data enabled us to calculate range size for 171 bee species ; range size calculations were not possible for species with too few specimen records or taxa not identified to described species . Using geographical information systems analyses available via ArcGIS and QGIS , we mapped individual records of each bee species and constructed concave polygons bounding the set of location data points for each species. Range size for each species was calculated as the internal area of each species’ concave polygon, which yields the smallest area that encompasses the set of location data points for each species and allows for more accurate determination of species range size compared to convex polygons. We then calculated the average range size of all bee individuals at each study plot , and compared average range sizes between reserve and fragment plots using two-sample t-tests.Functional trait distribution among species is often influenced by phylogeny , including in bees . Thus, to aid in the interpretation of our results in view of phylogenetic signals present in our data, we quantified the variation in each trait attributed to each major taxonomic rank by constructing nested generalized linear mixed-effects models . In these models, the value of each trait is the dependent variable, and taxonomic ranks were included as random effects .