Bumble bees exhibit body size clines across an urban gradient despite low genetic differentiation

Environmental heterogeneity resulting from human-modified landscapes can increase intraspecific trait variation. However, less known is whether such phenotypic variation is driven by plastic or adaptive responses to local environments. Here, we study five bumble bee (Apidae: Bombus) species across an urban gradient in the greater Saint Louis, Missouri region in the North American Midwest and ask: (1) Can urban environments induce intraspecific spatial structuring of body size, an ecologically consequential functional trait? And, if so, (2) is this body size structure the result of plasticity or adaptation? We additionally estimate genetic diversity, inbreeding, and colony density of these species—three factors that affect extinction risk. Using ≥ 10 polymorphic microsatellite loci per species and measurements of body size, we find that two of these species (Bombus impatiens, Bombus pensylvanicus) exhibit body size clines across the urban gradient, despite a lack of population genetic structure. We also reaffirm reports of low genetic diversity in B. pensylvanicus and find evidence that Bombus griseocollis, a species thought to be thriving in North America, is inbred in the greater Saint Louis region. Collectively, our results have implications for conservation in urban environments and suggest that plasticity can cause phenotypic clines across human-modified landscapes.


Site Descriptions
Calvary Cemetery (CC) is a Catholic cemetery located in the city of Saint Louis, Missouri (MO) (human population density = 1,889 people km -2 ), which contains 25 acres of prairie managed by the Missouri Department of Conservation along the cemetery's northwestern edge, for which a conservation plan was implemented in 2005 [1]. EarthDance Farms (ED) is an organic farm located in Ferguson, MO (human population density = 1,293 people km -2 ), comprising 14 acres and a variety of native and agricultural plants, which has been a location of organic food production since 1883 [2]. Castlewood State Park (CW) is a state park adjacent to the Meramac River in Ballwin, MO (human population density = 1,297 people km -2 ), comprising 1,818 acres of land and was established in 1974 [3]. Shaw Nature Reserve (SNR) is a private nature reserve located on the edge of the Missouri Ozarks in Gray Summit, MO -an unincorporated community near Pacific, MO (human population density = 472 people km -2 ) -comprising 2,500 acres of land and upwards of eight biomes [4], which was established in 1925 [5].

Power Analysis
To ensure that our data had sufficient statistical power to detect true genetic differentiation, we performed a power simulation per species with the program POWSIM 4.1 [6]. POWSIM tests the null hypothesis of no genetic differentiation between subpopulations, given different combinations of sample size, loci, and alleles [6]. Each simulation estimates power via chi-square and Fisher's exact tests, while sampling from populations that diverge following a Wright-Fisher model [6].
For all simulations, we set the expected differentiation between subpopulations to FST=0.05, which is an appropriate minimum value for true genetic structure [7]. This FST is equivalent to each subpopulation having Ne=100 after 10 generations of drift [8]. We parameterized each simulation with its respective species' observed sample size, loci number, allele number, and allele frequencies. We ran 1,000 iterations of each simulation with default parameters for dememorizations, batches, and iterations per batch. These simulations indicate the power of our sampling protocol to detect an FST=0.05 and do not represent the true evolutionary history of our study populations.

Weinberg Equilibrium
Prior to performing population genetic analyses, we removed loci from our microsatellite data following various quality control measures. Specifically, we removed loci that had >20% amplification failure, noisy amplification (making a locus unreliable to score), >25% null allele frequency following Chakraborty et al. (1992) [9], or showed significant linkage disequilibrium (LD) with one or more loci. Analyses were performed in R Statistics. The package PopGenReport 2.0 [10] identified null alleles. The package Genepop '007 1.1.4 [11] identified loci in LD. In the following, we describe the results of these quality control measures per species. See Table S1 for loci retained per species for analyses.
Following these quality measures, 10 loci remained for the population genetic analyses of B. griseocollis.
Accordingly, we removed BTMS0062, BTern02, BTMS0044, and BTern01 from B. pensylvanicus. Following these quality measures, 14 loci remained for the population genetic analyses of B. pensylvanicus.

Bumble Bee Sampling
In addition to the sampling performed in 2018, in the summer of 2017, we sampled worker bumble bees at Shaw Nature Reserve (SNR). From late-June through mid-August, we sampled foraging workers of B. impatiens, B. griseocollis, B. auricomus, and B. pensylvanicus by hand-netting 3-4 days per week. After capture, we immediately transferred bees to individual vials containing 100% ethanol. We did not sample B. bimaculatus in 2017, as the onset of our sampling corresponded with the latter half of their seasonal period of foraging activity. Sample sizes collected per species can be found in Table S5.

Microsatellite Genotyping, Colony Density, and Allelic Richness
We  Table S5.     Figure   Fig S1. Global allelic richness at each locus after sample size rarefaction (n=176 alleles/species).

Supplemental Tables
Loci with zero values were either unamplified or dropped from analyses. Means (+SE) are computed for intraspecific loci with nonzero values.