Pollinator assemblage and pollen load differences on sympatric diploid and tetraploid cytotypes of the desert-dominant Larrea tridentata.

PREMISE: Whole-genome duplication (polyploidy) is an important force shaping flowering-plant evolution. Ploidy-specific plant-pollinator interactions represent important community-level biotic interactions that can lead to nonrandom mating and the persistence of mixed-ploidy populations.
METHODS: At a naturally occurring diploid-tetraploid contact zone of the autopolyploid desert shrub Larrea tridentata, we combined flower phenology analyses, collections of bees on plants of known cytotype, and flow cytometry analyses of bee-collected pollen loads to investigate whether (1) diploid and tetraploid plants have unique bee pollinator assemblages, (2) bee taxa exhibit ploidy-specific visitation and pollen collection biases, and (3) specialist and generalist bee taxa have ploidy-specific visitation and pollen collection biases.
RESULTS: Although bee assemblages overlapped, we found significant differences in bee visitation to co-occurring diploids and tetraploids, with the introduced honeybee (Apis mellifera) and one native species (Andrena species 12) more frequently visiting tetraploids. Consistent with bee assemblage differences, we found that diploid pollen was overrepresented among pollen loads on native bees, while pollen loads on A. mellifera did not deviate from the random expectation. However, mismatches between the ploidy of pollen loads and plants were common, consistent with ongoing intercytotype gene flow.
CONCLUSIONS: Our data are consistent with cytotype-specific bee visitation and suggest that pollinator behavior contributes to reduced diploid-tetraploid mating. Differences in bee visitation and pollen movement potentially contribute to an easing of minority cytotype exclusion and the facilitation of cytotype co-occurrence.

Polyploidy—whole‐genome duplication—is a pervasive and important force shaping angiosperm evolution (Barker et al., 2015) and contemporary interspecific, community‐level, and ecosystem interactions (Segraves, 2017; Gaynor et al., 2018). Yet, persistent questions remain over the potential cascading ecological consequences and biodiversity implications of polyploidy (Ramsey and Ramsey, 2014). Chromosome number differences arising from genome duplication typically lead to strong intrinsic reproductive isolation between populations differing in ploidy (e.g., Coyne and Orr, 2004; Madlung, 2013; Ramsey and Ramsey, 2014; Gross and Schiestl, 2015; Sutherland and Galloway, 2017), and rapid selection against the minority cytotype resulting from a greater proportion of low fitness intercytotype matings (i.e., minority cytotype exclusion; Levin, 1975). Studies over the last several decades, however, have repeatedly documented phenotypic and ecological differences between plants differing in ploidy (e.g., Lumaret, 1988; Segraves and Thompson, 1999; Husband and Schemske, 2000; Levin, 2004; Maherali et al., 2009; Ramsey, 2011; Madlung, 2013; Gross and Schiestl, 2015; McCarthy et al., 2016; McIntyre and Strauss, 2017). Even if only slight, differences in traits such as cell size, secondary compound production, water use, flowering time, and flower color can facilitate the exploitation of novel ecological niches, ease competition between ploidies, and result in the long‐term maintenance of multiple intraspecific cytotypes. What remains less clear is the degree to which the suite of phenotypic alterations typically accompanying shifts in ploidy might contribute to cascading biotic interactions at the community level and the origins of new biodiversity (Segraves and Anneberg, 2016).

Polyploid plant–pollinator interactions remain relatively understudied despite representing ecological differentiation that simultaneously influences fitness in species reliant upon animal pollinators (Segraves and Thompson, 1999; Coyne and Orr, 2004; Kennedy et al., 2006; Halverson et al., 2008; Segraves, 2017). Differing animal‐mediated reproductive consequences between co‐occurring diploids and polyploids have only been carefully documented within the last few decades (Segraves and Thompson, 1999), and recent investigations exploring ploidy‐specific pollinator interactions suggest visitation differences may influence the population genetics of populations where cytotypes co‐occur (i.e., mixed cytotype populations; Segraves and Thompson, 1999; Husband and Schemske, 2000; Nuismer and Thompson, 2001; Husband and Sabara, 2003; Thompson et al., 2004; Kennedy et al., 2006; Thompson and Merg, 2008; Nghiem et al., 2011; Borges et al., 2012; Gross and Schiestl, 2015; Roccaforte et al., 2015; Barringer and Galloway, 2017). For example, the bee assemblages of Erythronium mesochoreum and Erythronium albidum (diploid and autotetraploid, respectively) overlap, but differ significantly in the frequency of visits where the two species co‐occur (Roccaforte et al., 2015). Similarly, the major pollinators of Heuchera grossulariifolia (Segraves and Thompson, 1999; Thompson and Merg, 2008), Chamerion angustifolium (Husband and Schemske, 2000; Kennedy et al., 2006), and Libidibia ferrea (Borges et al., 2012) differ in their visitation frequency to sympatric diploids and autotetraploids and in their pollination effectiveness on the two cytotypes. While influencing the frequency of intercytotype gene flow and potentially playing a role in maintaining the contemporary co‐occurrence of multiple ploidies by easing minority cytotype exclusion, the bee assemblage and visitation differences documented on these species may also represent cytotype‐specific specialization and the exploitation of novel ecological niches.

Plant species with very large and diverse pollinator assemblages present a unique opportunity to gain insight into the role that polyploidy plays in altering plant–animal interactions, as different pollinator species may have distinct ploidy‐specific interactions. For example, large pollinator assemblages comprising pollen specialists (species that consistently collect pollen from a single plant species or group of related species in the presence of alternative pollen sources) and pollen generalists (species that are not identifiably limited to particular pollen sources; Hurd and Linsley, 1975) may reveal unique interspecific interactions with differing consequences for the maintenance of cytotype co‐occurrence. Specialist pollinators, in particular, may be finely attuned to subtle phenotypic variation of their floral host (Waser, 1986; Minckley et al., 1999; Vaudo et al., 2016) and may thus be most likely to cause assortative mating where cytotypes occur in sympatry. In contrast, generalist species may move pollen indiscriminately between cytotypes resulting in minority cytotype exclusion and/or intercytotype gene flow. Yet, the degree to which pollinators facilitate or prevent intercytotype pollen movement is difficult to study and remains relatively undercharacterized in natural mixed ploidy populations.

The North American creosote bush [Larrea tridentata (DC.) Coville; Zygophyllaceae], a classic autopolyploid complex (Hunziker et al., 1977; Lewis, 1980), represents a unique opportunity to investigate pollinator cytotype specialization and to better understand how ecological processes such as plant–pollinator interactions may be contributing to nonrandom pollen movement and the maintenance of mixed‐ploidy populations. Larrea tridentata is a characteristic arid‐adapted, predominantly outcrossing, but self‐compatible, multiflorous shrub (Simpson 1977), that comprises three cytotypes distributed throughout the Chihuahuan (diploids, 2n = 2x = 26), Sonoran (predominantly tetraploids, 2n = 4x = 52), and Mojave Deserts (hexaploids, 2n = 6x = 78) of the southwestern United States and northern Mexico. The ploidies naturally occur sympatrically in geographically restricted areas at their distributional boundaries (Hunter et al., 2001; Yang, 1970; Laport et al., 2012) where rare triploid and pentaploid intercytotype hybrids have also been documented (Laport and Ramsey, 2015). Analyses of DNA molecular markers also suggest occasional ongoing gene flow restricted to areas of cytotype sympatry and parapatry (Laport et al., 2016).

Ecologically dominant across thousands of hectares of desert biome, L. tridentata responds rapidly to modest rainfalls, typically initiating spring flowering before co‐occurring species (Barbour et al., 1977; Benson and Darrow, 1981; Turner et al., 1995; Whitford et al., 1996), and represents a major reliable pollen and nectar resource for the hyperdiverse bee communities of the North American warm deserts (estimated to be nearly 900 species across the deserts, with ~500 species documented in the San Bernardino Valley spanning the United States–Mexico border between Arizona and Chihuahua; Moldenke, 1979; Minckley et al., 2008; Danforth et al., 2019, Meiners et al., 2019; R. L. Minckley and W. R. Radke, unpublished data). A diverse assemblage of 120 native pollen specialist and pollen generalist bees visit L. tridentata throughout its range, including 20 native pollen specialists and the non‐native and recently naturalized generalist, Apis mellifera (Hurd and Linsley, 1975; Minckley et al., 2000). The native pollinator species are primarily small, ground nesting, active for a month or less per year, and are solitary (Danforth et al., 2019), whereas the introduced A. mellifera is social with perennial colonies active year‐round.

The large and diverse assemblage of generalist and specialist bees on L. tridentata presents multiple opportunities for cytotype‐specific specialization and cryptic assortative mating that may facilitate the persistence of mixed ploidy populations (Segraves and Thompson, 1999; Husband and Schemske, 2000; Coyne and Orr, 2004; Gross and Schiestl, 2015; Roccaforte et al., 2015). We investigated bee assemblage differences and pollen movement biases in a naturally occurring sympatric population of diploid and tetraploid L. tridentata, combining field observations of flower production, collections of bees on plants of known cytotype, and flow cytometry analyses of individual bee‐collected pollen loads. These data allowed us to determine whether (1) bee assemblages differ between diploid and tetraploid plants, (2) native bees and non‐native A. mellifera differ in visitation to, and pollen collection from, diploid and tetraploid plants, and (3) native pollen specialist and pollen generalist bees differ in visitation to, and pollen collection from, diploid and tetraploid plants. We found that sympatric diploids and tetraploids had modestly differentiated bee assemblages, and flow cytometry analysis of bee‐collected pollen loads revealed pollen load composition differed significantly from that expected from a model of random mating, consistent with pollinator‐mediated assortative mating in excess of that expected from pollinator assemblage overlap alone.

MATERIALS AND METHODS

Sampling site and flowering phenology

We sampled 10 permanently marked diploid and 10 permanently marked tetraploid plants in the spring of 2014, 2015, and 2016 at a previously identified sympatric site in the San Pedro River valley of southeastern Arizona, United States (Fig. 1; Laport and Ramsey, 2015, “San Pedro 3”, 32°35.800′N, 110°32.300′W). Larrea tridentata is dominant at the study site and occurs nearly continuously along the river valley terraces with other characteristic desert perennial species [e.g., Carnegiea gigantea (Engelm.) Britt. & Rose, Prosopis velutina Wooton, Parkinsonia L. spp., Opuntia Mill. spp., Acacia Martius spp., Fouqueria splendens Engelm., Calliandra eriophylla Benth., and Ferocactus wislizeni (Engelm.) Britt. & Rose; Laport and Minckley, 2013; Laport and Ramsey, 2015]. Prior sampling and cytotype screening suggest diploids comprise ~69% and tetraploids comprise ~31% of the plants at the site, with the cytotypes often intermingling within ~2–5 m of each other (Laport and Ramsey, 2015), but it is unclear whether these estimates reflect the true cytoytpe frequencies because subsequent analyses suggested the cytotypes might be similarly abundant at the site (i.e., ~54% diploids, ~46% tetraploids; R. G. Laport, unpublished data). Prior analyses at this site indicate flower phenology differs slightly, but overlaps between diploids and tetraploids (Laport and Ramsey, 2015; Laport et al., 2016). We confirmed that flowering on the marked diploid and tetraploid plants overlapped by counting flowers on marked branches on each plant and estimated total flower production during the periods we collected bees. Flowers were considered “open” when buds had opened enough for bees to access the stigmas and anthers. We calculated total and mean flower production for each cytotype over the three collection seasons. Flower counts were transformed by adding 0.5 to all values before square‐root transformation for analysis with repeated measures ANOVA models in the JMP statistical package (version 13; SAS Institute, Cary, NC, USA) that included ploidy, date, and ploidy × date as effects.

Figure 1

The range of Larrea tridentata, indicated by the gray shading in (A), extends throughout the Chihuahuan, Sonoran, and Mojave Deserts of the southwestern United States and northern Mexico. (A) Diploids and tetraploids naturally co‐occur near the boundary between the Chihuahuan and Sonoran Deserts, denoted with a star, where we sampled bees visiting sympatric plants. (B) The cytotype boundary between diploids and tetraploids was previously characterized by sampling 50 plants every 2–5 km along a transect from pure diploid sites (white squares) to pure tetraploid sites (gray circles) in the San Pedro River Valley. The sympatric site (San Pedro 3) is denoted with a pie chart showing the apparent relative proportions of diploids (~69%) and tetraploids (~31%) determined from a previous study (Laport and Ramsey, 2015). (C) Detail of the relative spatial distribution of the 10 diploid and 10 tetraploid plants monitored for flower production and bee visitation at the sympatric site. At the sympatric site, L. tridentata is dominant, and the 20 focal plants are embedded within a nearly continuous matrix of L. tridentata for which cytotype information is not known.

Bee collections

We collected and identified bees from the 20 permanently marked plants of known cytotype during the spring blooms of 2014, 2015, and 2016 (Fig. 1; Laport and Ramsey, 2015). Sampled plants were ≥5 m apart to avoid sampling clonal ramets. Bees were netted while foraging on flowers of each marked focal plant. Because many of the native bees were small and difficult to see until flying, we additionally collected bees flying within ~0.5 m of focal plants (≤50% of the canopy span), assuming they already had, or would have, visited flowers on the focal plant. Sampling effort on each plant was standardized by netting bees for 5 min/plant/sampling bout. We conducted one to four sampling bouts per day at approximately 08:00, 10:00, 13:00, and 15:00 hours. We randomized the order of plants to be netted upon for each sampling bout. All collected bees were pinned, preserved over silica desiccant, and identified to the lowest possible taxonomic level (all vouchers held in the collection of R. L. Minckley).

We calculated species richness (S) and Shannon–Wiener diversity (H′) indices from the specimens collected on diploids and tetraploids over the three collection seasons and tested whether bee assemblages differed between co‐occurring cytotypes. First, we tested whether entire bee assemblages differed between cytotypes with permutational MANOVA (PERMANOVA) implemented with the adonis function in the R (v3.3.2; R Core Team, 2016) package vegan (Oksanen et al., 2017). Individual species differences were analyzed with ANOVA. We additionally tested for differences in bees caught on each ploidy with a general linear model (GLM) in JMP that included ploidy and year as fixed effects, mean flower number as a random effect, as well as ploidy × year and ploidy × mean flower number interactions.

The overlap in total bee assemblages on diploids and tetraploids was calculated using Pianka’s niche overlap index, O (Pianka, 1974), where j = bees collected on diploid plants and k = bees collected on tetraploid plants over the three surveyed seasons. O values near 0 indicate low bee assemblage (“niche”) overlap, while those close to 1 indicate high bee assemblage (“niche”) overlap. To assess the significance of differences in bee assemblage overlap, we simulated expected overlap 1000 times by randomizing bee taxa identities to diploids and tetraploids, preserving the observed number of bees on each cytotype in the R (v3.3.2) base package, and compared the observed value to the distribution of simulated expectations. We considered an observed value of overlap falling outside the 95% confidence interval of the simulated data as a significant difference. O and randomizations were calculated with a custom script (Appendix S1) with A. mellifera included and with A. mellifera excluded.

Pollen‐load composition and pollen–plant ploidy mismatches

We investigated whether bee visitation differed between diploids and tetraploids in two ways. First, we removed and determined the cytotype of whole pollen loads collected by bees (i.e., pollen rinsed from entire bee) and tested whether the proportion of pure diploid (containing only haploid pollen), pure tetraploid (containing only diploid pollen), and mixed‐pollen loads (containing both haploid and diploid pollen) removed from the bees differed from a random‐mating expectation derived from the total mean estimated flower production for the 20 focal diploids and tetraploids pooled across the three seasons: expected diploid–diploid matings = (proportion diploid flowers)2, expected tetraploid–tetraploid matings = (proportion tetraploid flowers)2, and expected diploid–tetraploid matings = 2(proportion diploid flowers)(proportion tetraploid flowers). Second, we combined pollen load composition and plant cytoytpe information to test whether the number of pollen load–plant ploidy matches and pollen load–plant ploidy mismatches (including mixed pollen loads) differed from the random mating expectation derived from total mean estimated flower production of the 20 focal diploids and tetraploids. Both analyses were conducted for A. mellifera, all native bees combined, and native pollen specialist and generalist bees separately for all 3 years combined to obtain adequate pollen load sample sizes. The first analysis leverages ploidal determinations of pollen present on bees and assumes that pollen collection and deposition is concordant with floral visitation, but may underestimate realized intercytotype pollen movement and deposition because of unsampled foraging bees. The second approach may overestimate intercytotype pollen movement by counting all pollen–plant ploidy mismatches and mixed–cytotype pollen loads as intercytotype matings. Both models assume that all bee taxa have an equal probability of visiting diploid and tetraploid plants in sympatry, that diploid and tetraploid pollen collection is equally likely, and that pollination is equally likely on both cytotypes. Because of uncertainty in diploid and tetraploid frequencies and distributions at the sympatric site, our models also assume that the cytotypes are equally abundant and randomly distributed.

We used flow cytometry to determine the cytotype of pollen loads from collected bees. Flow cytometry has recently been optimized for ploidal analysis of pollen grains (Kron and Husband, 2012) and to estimate pollen movement among intraspecific polyploids (Kron et al., 2014), and we followed the procedure of Kron et al. (2014) to determine pollen cytotype. Briefly, we rinsed pollen loads from silica‐preserved bees with 1 mL of LB01 buffer (Doležel et al., 2007). Resuspended pollen grains were passed through a 100‐µm pre‐filter (Partec CellTrics, Görlitz, Germany) to remove large debris. Nuclei were extracted by gently rubbing pollen grains against a 10‐µm “bursting” filter (Partec CellTrics, Görlitz, Germany) with a glass stir rod, and then rinsed through the filter and stained with 500 µL of LB01 buffer containing 50 µL of propidium iodide at 1 mg/mL and 25 µL of RNase at 1 mg/mL. All samples were run on a FACSCalibur flow cytometer (B‐D Biosciences, San Jose, CA, USA) at the University of Nebraska‐Lincoln Flow Cytometry Service Center. Using CellQuest Pro Software (version 5.2.1; B‐D Biosciences), we inferred ploidy from the relative fluorescence (FL2A) of each sample compared to L. tridentata tissue of previously determined DNA content (Laport et al., 2012), or plant tissue recommended by Doležel et al. (2007) as external standards run at the beginning of each session (Raphanus sativus cv. Saxa, 2C DNA content = 1.11 pg; Glycine max cv. Polanka, 2C DNA content = 2.50 pg). The standards allowed us to determine the approximate range of expected fluorescence for diploid‐ and tetraploid‐derived pollen (± approximately 10–15%), and we scored the presence or absence of DNA fluorescence modality in these expected ranges to infer whether each pollen load comprised only diploid‐derived, only tetraploid‐derived, or both diploid‐ and tetraploid‐derived pollen (additional flow cytometry details in Appendix S2).

RESULTS

Flowering phenology

Diploid and tetraploid flowering phenologies overlapped in all years (Fig. 2). The combined mean flower production from 2014 to 2016 was not different for diploids and tetraploids (2014: F 1, 16 = 3.439, P = 0.169; 2015: F 1, 10 = 3.070, P = 0.053; 2016: F 1, 8 = 1.888, P = 0.162), though peak flowering times were not concordant. Tetraploids tended to produce more flowers earlier in the season (before mid‐March) than diploids, which tended to produce more flowers later in the season (after mid‐March; Fig. 2). Tetraploids (82.2 mean flowers/plant) produced slightly more flowers than diploids (61.4 mean flowers/plant) over the three collection years, but the difference was not significant.

Figure 2

Spring flower production by, and bee visitation to, sympatric diploid (squares) and tetraploid (circles) Larrea tridentata from (A) 2014, (B) 2015, and (C) 2016. Mean flower counts are shown for the 10 diploid and 10 tetraploid plants on which bees were collected each spring as solid black and gray lines, respectively. Mean bee visitation (measured as the number of bees collected per plant) for the 10 diploid and 10 tetraploid plants is shown as black and gray dashed lines, respectively. Flower counts were made on days bees were collected. Error bars indicate ±1 SE.

Bee assemblages

We sampled for 135 h over 37 days between spring 2014 and spring 2016, collecting 1272 bees, representing 61 taxa, foraging on or flying near marked focal L. tridentata. We also observed an additional 463 A. mellifera visiting focal plants during sampling in 2016 that were intentionally not collected to increase sampling of native bees. Thus, a total of 1735 bees were caught or observed. Of the 1272 collected bees, 19 individual specimens on diploids and 25 individual specimens on tetraploids could not be identified and ~55% of collected bees were non‐native A. mellifera (Appendix S3). Forty‐one bee taxa were collected on diploids (S 2 = 41, = 1.62) and 50 taxa were collected on tetraploids (S 4 = 50, = 1.45).

The number of bees collected on individual plants ranged from 43 to 161 bees/plant (mean = 86.8 bees/plant), with the largest number of bees being caught on both cytotypes within a few days of the peak bloom of diploids and tetraploids (Fig. 2). The number of bees collected differed between years for the two co‐occurring ploidies (ploidy × year; F 2, 37 = 3.524, P = 0.040). Significantly fewer bees were collected on diploids (69.4 bees/plant) than on tetraploids (104.2 bees/plant; F 1, 18 = 7.631, P = 0.013). However, the number of bees collected/flower on diploids (1.86) was not significantly different from the number of bees collected/flower on tetraploids (2.01; F 1, 18 = 0.029, P = 0.867), and there was not an effect of flower number on the number of bees caught (ploidy × flower number; F 16, 36 = 1.009, P = 0.470).

The bee assemblages visiting sympatric diploids and tetraploids were not identical (Fig. 3; Appendix S4; F 1, 18 = 2.687, R 2 = 0.130, P = 0.037). Apis mellifera (F 1, 18 = 9.382, P = 0.007) and a species of Andrena (species 12; F 1, 18 = 8.707, P = 0.009) were more commonly collected on tetraploids (Fig. 3). Native bee taxa comprised a slightly larger proportion of all bees collected on diploids (42.4% of bees) than on tetraploids (32.6% of bees), and several taxa appeared to be more common in collections on diploids (Hoplitis biscutellae, Colletes clypeonitens, Lasioglossum microlepoides, Halictus tripartitus; Fig. 3; Appendix S4), though the proportions of native bees on diploids and tetraploids were not significantly different (F 1, 18 = 2.219, P = 0.154).

Figure 3

Frequency of native bee taxa collected on sympatric diploid (2x) and tetraploid (4x) Larrea tridentata. Most of the native bee taxa visiting L. tridentata were rare, but 14 taxa were represented by 10 or more occurrences on diploids and tetraploids combined: Hesperapis larreae (126 specimens), Ancylandrena larreae (65 specimens), Colletes salicicola (45 specimens), Megachile xerophila (44 specimens), Lasioglossum microlepoides (30 specimens), Trachusa larreae (26 specimens), Ashmeadiella breviceps (21 specimens), Colletes clypeonitens (20 specimens), Hoplitis biscutellae (19 specimens), Andrena species 12 (18 specimens), Halictus species 1 (18 specimens), Halictus tripartitus (16 specimens), Andrena species 9 (10 specimens), and Perdita exclamans (10 specimens). Differences in individual bee species abundance between diploids and tetraploids were slight, however, Andrena species 12 was more commonly collected on tetraploids than on diploids.

The observed overlap in bee assemblages on diploids and tetraploids was O = 0.985 (O = 0.946 excluding A. mellifera), indicating very high overlap. However, this value was significantly lower than the simulated distribution of overlap values (range 0.991–0.999), suggesting that the observed intercytotype pollinator assemblage overlap was lower than expected if pollinator visitation was random (P < 0.01). With A. mellifera excluded, overlap values were shifted lower (range 0.905–0.989), and the observed overlap did not differ from the simulated expectation (P > 0.05).

Pollen load composition and pollen–plant ploidy mismatches

Initial flow cytometry analysis of known diploid and tetraploid pollen produced bimodal FL2A histograms of relative fluorescence, consistent with L. tridentata having binucleate pollen (Brewbaker, 1967). Flow cytometry analyses of pollen loads removed from bees containing a mix of diploid and tetraploid pollen produced trimodal fluorescence histograms (Appendix S2). Approximately 18% of the collected bees had visible pollen loads, though we attempted to remove and analyze pollen from all collected native bees. For 115 bees, pollen loads were of sufficient size to produce acceptable fluorescence histograms via flow cytometry (105 samples were rejected for producing poor quality histograms or containing non‐L. tridentata pollen). Fluorescence peaks for these 115 histograms had coefficients of variation averaging 4.68% (range 1.32–10.45%) with an average of 766 events (range 20–3739; Appendix S2). Of the 115 pollen loads, 23.5% comprised diploid pollen, 29.6% comprised tetraploid pollen, and 46.9% comprised both diploid and tetraploid pollen. Mixed pollen loads, and plant ploidy–pollen ploidy mismatches occurred on diploid and tetraploid plants throughout the site.

Of the 86 A. mellifera pollen loads, 17.4% comprised diploid pollen, 37.2% comprised tetraploid pollen, and 45.4% comprised both diploid and tetraploid pollen. However, of the 29 native bee pollen loads, 41.4% comprised diploid pollen, 6.9% comprised tetraploid pollen, and 51.7% comprised both diploid and tetraploid pollen (Table 1). From the observed total mean flower production from 2014 to 2016, we expected 18.3% of matings to be diploid–diploid, 32.8% to be tetraploid–tetraploid, and 49.0% to be diploid–tetraploid. Assuming the pollen carried by bees could be deposited on receptive stigmas, the proportions of diploid, tetraploid, and mixed pollen loads from A. mellifera were not significantly different from a random‐mating expectation ( = 0.778, P = 0.678; Table 1). In contrast, the pollen load proportions from the native bees were significantly different from random mating with a bias toward diploid pollen and fewer tetraploid pollen loads than expected = 14.417, P < 0.001; Table 1).

Table 1

Expected proportion of diploid and tetraploid matings of Larrea tridentata from a model of random mating and mean ploidy compositions of pollen loads from 2014 to 2016. The expected proportions were derived from the mean flower production for the 10 focal diploid and 10 focal tetraploid plants over the three surveyed years of the study. The percentage (and number) of bee‐collected pollen loads comprising diploid only, mixed, and tetraploid only pollen by Apis mellifera and native bees was inferred from flow cytometry analyses. Apis mellifera pollen load compositions were not significantly different from random mating while native bees exhibited a bias toward diploid pollen. Native pollen specialist and generalist bees (according to Hurd and Linsley [1975]) exhibited similar biases toward diploid pollen.

	2x–2x	2x–4x	4x–4x
Expected matings	18.3%	49.0%	32.8%

Proportions for category are significantly different from random mating.

Categorizing native bees as pollen specialists and generalists following the classifications of Hurd and Linsley (1975) revealed similar pollen collection patterns. Of the 15 specialist‐bee pollen loads, 40.0% comprised diploid pollen, 6.7% comprised tetraploid pollen, and 53.3% comprised both diploid and tetraploid pollen. Of the 14 generalist‐bee pollen loads, 42.9% comprised diploid pollen, 7.1% comprised tetraploid pollen, and 50.0% comprised both diploid and tetraploid pollen. Pollen load compositions for both pollen specialist bees ( = 7.047, P = 0.030; Table 1) and generalist bees differed significantly from the random mating expectation ( = 7.433, P = 0.024; Table 1).

We also analyzed the proportion of ploidy matches and mismatches between pollen loads and the plant on which bees were collected (Table 2). Among A. mellifera, 16.3% of pollen loads removed from diploid visitors comprised diploid pollen, 34.9% of pollen loads removed from tetraploid visitors comprised tetraploid pollen, and 48.8% of pollen loads represented a mismatch (i.e., diploid pollen on tetraploid plants, tetraploid pollen on diploid plants, or mixed pollen loads on either cytotype). These match/mismatch proportions were not significantly different from the random expectation ( = 0.306, P = 0.858; Table 2). Among native bees, 37.9% of pollen loads removed from diploid visitors comprised diploid pollen, 6.9% of pollen loads removed from tetraploid visitors comprised tetraploid pollen, and 55.2% of pollen loads represented a mismatch. These proportions indicated a native bee bias toward diploid pollen, but also slightly more pollen–plant ploidy mismatches than expected if mating was random ( = 12.265, P = 0.002). Notably, no native visitors to diploid plants had tetraploid pollen loads (Table 2).

Table 2

Expected proportion of diploid and tetraploid matings of Larrea tridentata from a model of random mating and percentage (and number) of bees exhibiting plant ploidy–pollen load ploidy matches and mismatches from 2014 to 2016. The expected proportions were derived from the mean flower production for the 10 focal diploid and 10 focal tetraploid plants over the three surveyed years of the study. The most common mismatches were mixed pollen loads on both diploid and (especially) tetraploid plants for both Apis mellifera and native bees. Native pollen specialist and generalist bees (according to Hurd and Linsley [1975]) exhibited similar plant ploidy–pollen ploidy matches and mismatches.

	2x‐2x	2x‐4x	4x‐4x
Expected matings	18.3%	49.0%	32.8%

Proportions for category are significantly different from random mating.

Classifying the native bees as pollen specialists and generalists revealed similar pollen–plant ploidy mismatches (Table 2). Among pollen specialists, 40.0% of pollen loads removed from diploid visitors comprised diploid pollen, 6.7% of pollen loads removed from tetraploid visitors comprised tetraploid pollen, and 53.4% of pollen loads represented a mismatch. Plant–pollen ploidy mismatches among specialist bees only involved mixed pollen loads on both diploid and tetraploid plants, with most of the mixed pollen loads occurring on tetraploids rather than diploids. Among pollen generalists, 35.7% of pollen loads removed from diploid visitors comprised diploid pollen, 7.1% of pollen loads removed from tetraploid visitors comprised tetraploid pollen, and 57.1% of pollen loads represented a mismatch. These match/mismatch proportions represented a significant departure from the random expectation for pollen specialists ( = 7.047, P = 0.030), but not for pollen generalists ( = 5.324, P = 0.070; Table 2).

DISCUSSION

Our in‐depth investigation into pollinator assemblage differences, and pollen movement within and between co‐occurring diploid and tetraploid L. tridentata indicates that bee assemblage and foraging behavior differences may be playing a subtle but important role in facilitating the persistence of mixed‐cytotype populations. While bee assemblages were different on diploids and tetraploids, we also found that native specialist and generalist taxa exhibited pollen load biases toward diploid pollen relative to tetraploid pollen. In contrast, the non‐native, recently naturalized A. mellifera was collected more frequently on tetraploids, but pollen load analysis indicated random foraging. Combined, the bee assemblage differences and pollen collection biases of the native bees likely result in more nonrandom mating in sympatry than expected from flower production alone, even though the prevalence of mixed pollen loads suggests intercytotype mating continues to occur (i.e., diploid pollen on tetraploid flowers, tetraploid pollen on diploid flowers, mixed pollen loads on either ploidy). These findings also highlight the complex and cryptic nature of polyploid plant–animal interactions and the potential widespread importance of such community‐level biotic interactions for other polyploid species.

Diploid and tetraploid flowering time differences

In previous studies investigating broader scale patterns of phenotypic and phenological differences, diploid and tetraploid L. tridentata were shown to differ subtly in flower size, pollen size, flowering phenology, whole‐plant architecture, and leaf size and to have unique environmental associations at their distributional boundaries (Laport and Minckley, 2013; Laport et al., 2013, 2016; Laport and Ramsey, 2015). Diploid and tetraploid L. tridentata diverged within the last few hundred thousand years (Laport et al., 2012, 2016), meaning the phenotypic differences between the cytotypes likely arose at the time of tetraploid formation, or over the relatively short period of time since the origin of tetraploids. For the focal plants in this study, we did not find a significant ploidy × flower number effect on bee collections, but did find that (1) tetraploid plants tend to produce more flowers than diploids (though not significantly), (2) tetraploid plants tend to have higher overall bee visitation driven primarily by A. mellifera, and (3) native bees appear to collect pollen from diploid flowers more frequently than expected from a model of random mating. Though diploid and tetraploid differences in flower size, pollen size, and flowering phenology are subtle, they may suggest a mechanism that biases bee visitation and pollen collection (Husband et al., 2016).

Our findings are concordant with recent studies of autopolyploid Heuchera grossulariifolia (Segraves and Thompson, 1999; Thompson and Merg, 2008), Chamerion angustifolium (Husband and Schemske, 2000; Kennedy et al., 2006), Libidibia ferrea (Borges et al., 2012), and Galax urceolata (Barringer and Galloway, 2017) that suggest that floral size, shape, and/or phenology differences were correlated with pollinator visitation. Moreover, these studies suggest that even the typically subtle differences among intraspecific autopolyploids may be important for mediating plant–insect interactions. We did not investigate the relationship between floral phenotype or phenology and bee visitation in this study, and it remains unclear how important the timing of flower opening or floral and pollen phenotype differences between co‐occurring diploid and tetraploid L. tridentata are for bee visitation and pollen collection. Similarly, other traits such as nectar composition, floral scent, and light reflectance in wavelengths known to be important for insect vision (e.g., UV) may also contribute to visitation biases, and remain unexplored for L. tridentata.

The bee assemblage we observed visiting sympatric diploid and tetraploid L. tridentata was taxonomically rich, comprising 61 taxa. Though broadly overlapping, the bee assemblages on diploids and tetraploids were significantly different largely due to A. mellifera occurring more frequently on tetraploids, though some native species were apparently rarer, or absent, on either diploids or tetraploids (e.g., Lasioglossum microlepoides on diploids, Perdita exlamans on tetraploids; Fig. 3; Appendix S4). Native bees also comprised a slightly larger proportion of floral visitors to diploids than tetraploids (though not significantly), suggesting native bees may prefer to visit diploids. The abundance of A. mellifera foraging on L. tridentata has previously been found not to influence the abundance or diversity of native bees on L. tridentata (Minckley et al., 2003). Yet, we cannot rule out the possibility that these visitation differences indicate A. mellifera displaces native bees from tetraploid flowers by either initiating foraging earlier in the day, by more efficiently or aggressively removing pollen and nectar from tetraploid flowers, and/or by foraging on tetraploid flowers in greater numbers. Interactions between A. mellifera and native bees can also be complex, and such interactions have been shown to increase the pollination effectiveness of A. mellifera on other species (Greenleaf and Kremen, 2006). Studies of the pollinator assemblages on sympatric cytotypes of other plant species similarly report varying degrees of pollinator overlap, but generally infer that pollinator assemblage differences contribute to intercytoype reproductive isolation in excess of that predicted by a model of random mating (Kennedy et al., 2006; Thompson and Merg, 2008; Borges et al., 2012; Roccaforte et al., 2015; Barringer and Galloway, 2017). For example, Roccaforte et al. (2015) observed that some bee species visited both diploid E. albidum and tetraploid E. mesochoreum, but at different frequencies leading to appreciable reproductive isolation where the two species co‐occur. Additional investigations involving choice experiments or experimental arrays could shed additional light on whether co‐evolutionary dynamics have shaped the resource preferences of L. tridentata pollinators (Harder and Johnson, 2009), but the differences we documented in bee assemblages are consistent with cytotype‐specific specialization, nonrandom mating between diploids and tetraploids, and potentially an easing of minority cytotype exclusion that could facilitate cytotype co‐occurrence.

At a single mixed‐cytotype site, we collected approximately half of the reported bee diversity on L. tridentata over three spring flowering seasons (~61 of 120 species; Appendix S3). The native bee fauna of arid western North America is diverse, with estimates of ~300–600 species in locations across the southwestern United States (Meiners et al., 2019). In a long‐term study of bee richness across a diversity of habitat types in the San Bernardino Valley spanning the United States–Mexico border ~160 km from the current study site, ~500 bee species have been reported (Minckley, 2008; Danforth et al., 2019, Meiners et al., 2019; R. L. Minckley and W. R. Radke, unpublished data). Nevertheless, collecting this many bee taxa on a single plant species at a single site is surprising. The richness of the generalist bee assemblage may potentially be high in this area because it occurs within the Chihuahuan–Sonoran Desert ecotone and floristic elements of both deserts co‐occur providing an abundance of pollen and nectar resource options (McLaughlin, 1986; Laport and Minckley, 2013). The co‐occurrence of diploid and tetraploid L. tridentata in this area may similarly contribute to the bee diversity by representing non‐equivalent pollen and nectar resources presented over a prolonged phenological period. It is unclear how unsampled L. tridentata bee species might contribute to the overall assemblage differences and intercytotype pollen movement in sympatry. Some of these species have ranges that do not overlap with the zone of sympatry, while others are active only at other times of the year and could have different interactions with the co‐occurring cytotypes than documented here (Hurd and Linsley, 1975; Simpson, 1977; Simpson et al., 1977). Future collecting efforts over longer periods, or targeting different phenological periods, would help resolve this uncertainty, as well as to more fully characterize the importance of the Chihuahuan–Sonoran Desert ecotone for promoting native plant and bee biodiversity.

The application of flow cytometry to investigate the pollen load compositions of individual bees revealed that pollinator assemblages alone provide incomplete information about intercytoytpe reproductive interactions in L. tridentata. Our investigations, while broadly consistent with observed bee assemblage differences, suggest individual bees have cryptic biases in pollen collection from, and movement within and between, sympatric diploid and tetraploid plants. We found that native bees had diploid pollen loads in excess of that expected under a model of random mating based upon mean flower production of the co‐occurring cytotypes (Tables 1, 2). Moreover, native bees collected from diploid plants rarely had tetraploid or mixed pollen loads, while those collected from tetraploid plants usually had mixed pollen loads (Table 2). It is not clear whether these biases evolved after the formation of tetraploids or whether they arose from an ancestral preference for diploids before the origin of tetraploids. Consistent with observations in other polyploid species (Nghiem et al., 2011; Borges et al., 2012; Barringer and Galloway, 2017), A. mellifera appears to account for most of the diploid–tetraploid reproductive interactions (comprising >50% of bees in this study). Pollen loads removed from A. mellifera often comprised pollen from both diploids and tetraploids (Table 1) suggesting the introduction of this generalist pollinator to North America within the last ~400 years could have altered intercytotype reproductive interactions. Specifically, random A. mellifera‐mediated intercytotype pollen movement may now be swamping tetraploids (the apparent minority cytotype) at this zone of sympatry with pollen from diploids (the apparent majority cytotype).

The native bee diploid pollen load bias documented here suggests diploid plants might experience a fitness advantage over sympatric tetraploid plants in the absence of A. mellifera. The apparently greater diploid abundance and spatial clustering of the cytotypes in the surveyed population may increase the encounter rate of native bees with diploid flowers, resulting in the greater observed number of diploid pollen loads on native bees. Indeed, when the apparent differences in cytotype abundance in sympatry are accounted for by multiplying the observed diploid and tetraploid flower production by the previously documented cytotype frequencies (69% diploid, 31% tetraploid; Laport and Ramsey, 2015), native bee pollen load frequencies do not differ from the random expectation ( = 1.339, P = 0.512, not shown). In contrast, the frequency of A. mellifera pollen loads are significantly different from the random expectation, with a bias toward tetraploid pollen, after accounting for apparent cytotype frequency differences in sympatry ( = 40.958, P < 0.001, not shown). Yet, the cytotypes may be more equally represented in sympatry than previously documented (~54% diploids, ~46% tetraploids; R. G. Laport, unpublished data), the foraging flight distances of solitary bees typically approximate the spatial scale of the study site (on the order of 100–300 m; Zurbuchen et al., 2010), and mixed pollen loads were removed from bees collected throughout the study site, suggesting bee visitation might not simply reflect differences in sitewide cytotype abundance or spatial clustering (Tables 1, 2; Appendix S3). Thus, regardless of how the random expectation is generated from observed flower numbers, either native bees or A. mellifera exhibit pollen load biases, conferring a potential fitness advantage to either diploids or tetraploids. These pollen load biases could contribute to a relaxation of frequency dependent selection against tetraploids (the apparent minority cytotype) despite the prevalence of bee‐derived mixed cytotype pollen loads. For example, diploid pollen appears to mostly be transferred to tetraploids in mixed pollen loads that provide some opportunity for tetraploid–tetraploid matings, and tetraploid plants produce a greater number of flowers than diploids, potentially countering the numerical advantage of diploid plants. Finally, the observed pollen load biases and mismatches are likely conservative estimates of intercytotype reproductive interactions because the dynamics of pollen collection, transfer, and deposition are more complex than our simplifying assumption that pollen collection and deposition are concordant with floral visits.

We expected native specialist species, which have co‐evolved with L. tridentata, to be most likely to exhibit cytotype specialization (Waser, 1986; Minckley et al., 1999; Lopez‐Uribe et al., 2016; Vaudo et al., 2016; Danforth et al., 2019). The greatest species richness and abundance of L. tridentata pollen specialist species has previously been shown to occur where spring flowering is least predictable, suggesting ongoing co‐evolutionary dynamics between L. tridentata and its pollinators in response to triggers for the initiation of bloom (Minckley et al., 2000). Yet, our observations indicate that native pollen generalist species exhibited similar pollen collection patterns as native pollen specialist species. However, mixed pollen loads were more often recovered from specialist bees collected on tetraploid plants than diploid plants (Table 2), consistent with diploid pollen being more frequently collected and moved onto tetraploid plants vs. tetraploid pollen being collected and moved onto diploid plants. Generalist bees had a similar number of mixed pollen loads on diploids and tetraploids, consistent with diploids and tetraploids being comparable pollen resources for these species (Table 2). The bias in diploid pollen movement onto tetraploids by specialist bees suggests a bee preference for diploid pollen, but is also consistent with classical predictions (Stebbins, 1971) and prior observations of unidirectional intercytotype gene flow in polyploid species (e.g., Sutherland and Galloway, 2017). Prior analyses of chloroplast haplotypes, paternally inherited in L. tridentata (Yang et al., 2000), indicate some of the naturally occurring tetraploid seedlings in sympatry have a diploid chloroplast haplotype, suggesting the bee visitation and pollen collection patterns documented here may be facilitating occasional introgression of diploid plastid genomes into tetraploids (Laport et al., 2016).

Prior studies focusing on the visitation dynamics of large pollinator assemblages (Roccaforte et al., 2015; Barringer and Galloway, 2017) or a few target species (Nghiem et al., 2011; Borges et al., 2012) have revealed that ecologically mediated assortative mating may often be stronger than classically assumed for polyploid species (Schluter, 2000). Yet, studies focusing only on pollinator assemblage differences may underestimate the dynamics of pollen movement between sympatric cytotypes. Recent investigations quantifying pollen deposition or pollination effectiveness of individual pollinators have shown that pollinator assemblages do not always predict realized reproductive output (Kennedy et al., 2006; Thompson and Merg, 2008). The ploidy analysis of pollen loads removed from individual bees using flow cytometry offers similar insight into an ecological component of polyploid reproductive interactions in L. tridentata that may facilitate the persistence of a mixed‐cytotype population. It is enticing to extrapolate from the results of this study to suggest the observed pollinator differences may have also played a role in the establishment of tetraploids. However, diploid and tetraploid L. tridentata likely diverged within the last few hundred thousand years and have largely evolved independently since that time (Hunter et al., 2001; Laport et al., 2012, 2016). Though this represents a relatively recent divergence, tetraploid L. tridentata is well‐established and likely offers only limited insight into understanding the population dynamics responsible for the original establishment of new polyploids (Segraves and Anneberg, 2016). Additional studies are required to better understand how the novel phenotypes of neopolyploids might “tip the scales” from extinction to persistence and spread for new cytotypes (Ramsey and Ramsey, 2014).

Confident pollen species identification from DNA content histograms is challenging without a priori knowledge of the pollen source (Kron et al., 2014). However, the bimodal or trimodal fluorescence histograms produced by pollen from L. tridentata specialist bees and the large pollen loads of A. mellifera provide some reassurance for our inference of pollen identity, pollen load composition, foraging patterns, and the exclusion of non‐L. tridentata pollen (Appendix S2). Although some of the co‐occurring plant species in the study area also likely have binucleate pollen (Brewbaker, 1967) and DNA contents similar to L. tridentata (2C2 DNA content = 1.5 pg, 2C4 DNA content = 2.4 pg; e.g., Prosopis velutina 2C DNA content = 0.86 pg, Carnegiea gigantea 2C DNA content= 2.87 pg, Fouqueria splendens 2C DNA content = 1.06 pg; Pellicer and Leitch, 2020), the bimodal fluorescence histograms produced by L. tridentata pollen and comparison to the external standards aided the exclusion of pollen loads that did not conform to expectations for L. tridentata from our analyses. Furthermore, our estimates of pollen load composition represent conservative evaluations for the presence or absence of diploid and tetraploid pollen. We did not attempt to estimate the relative proportion of diploid and tetraploid pollen in mixed pollen loads because it was difficult to obtain sufficient pollen of known cytotype to evaluate various levels of diploid–tetraploid mixing and because doublets (two adhering nuclei) may have occasionally complicated pollen load composition inference from fluorescence histograms (Kron et al., 2014; Appendix S2). It is likely that some mixed pollen loads containing low levels of pollen from one cytotype were excluded because of our approach to scoring pollen loads. A greater proportion of mixed pollen loads would indicate more intercyotype pollen movement (and gene flow) than we estimated, and the ability to detect asymmetrically mixed pollen loads should be investigated further. While the drivers of pollen load composition differences among the bee functional groups studied here remain unclear, floral and pollen morphology and chemistry should be further investigated to better understand why diploid and tetraploid pollen appear not to be equivalent resources to some bees. What does seem clear, however, is that nonrandom pollinator‐mediated pollen movement (and potentially mating) may be more common in mixed‐ploidy populations than previously thought, and should be investigated more broadly for polyploid species.

CONCLUSIONS

Over the last few decades, pollinator‐mediated assortative mating among closely related populations has been documented as an important ecological mechanism of genetic divergence and speciation (Bradshaw and Schemske, 2003; Sobel and Streisfeld, 2015). Yet, such interactions remain relatively understudied in populations exhibiting ploidal variation, despite representing ecological differentiation that simultaneously influences reproductive interactions and fitness that may facilitate overcoming minority cytotype exclusion, polyploid establishment, and cytotype divergence (e.g., Segraves and Thompson, 1999; Sobel et al., 2009; Ramsey, 2011; Glennon et al., 2012; Martin and Husband, 2013; Roccaforte et al., 2015; Husband et al., 2016). Our study adds to a growing body of research on the biodiversity implications of whole‐genome duplication by documenting pollinator visitation differences among populations differing in ploidy. In addition to differentiated bee assemblages, we revealed diploid pollen collection biases by native bees using flow cytometry analyses of collected pollen loads that favor diploid–diploid matings at a frequency above the random expectation. These bee assemblage differences and nonrandom pollen load distributions may play an important role in facilitating the continued coexistence of mixed‐cytoytpe populations and may offer at least some insight into the past establishment of tetraploid L. tridentata. Such nonrandom reproductive interactions may also be contributing to genetic divergence between diploids and tetraploids (Coyne and Orr, 2004). At the same time, mixed pollen loads and pollen–plant ploidy mismatches remain common, suggesting ongoing reproductive interactions and potential intercytotype gene flow between diploid and tetraploid L. tridentata.

Parallel investigations into whether similar patterns of bee assemblage and pollen collection biases occur between sympatric tetraploids and hexaploids would provide a more comprehensive view of the ecological aspects of polyploid reproductive interactions in L. tridentata. For example, the strength of plant–insect interactions may differ among higher ploidies (Sutherland and Galloway, 2017; O’Connor et al., 2019), and the potential for multiple origins of tetraploid and hexaploid L. tridentata (Laport et al., 2016) may set the stage for complex geographic patterns associated with both environmental and genetic/ploidal variation (Thompson, 2005). Moreover, additional investigations into relationships between flower and bee sizes, foraging flight distances, the effects of flower phenology differences, the frequency of self‐fertilization, and pollination efficiency by individual bee species would prove illuminating with respect to pollinator discrimination and the potential for intercytotype pollen movement by specialist and generalist bees. Nevertheless, the patterns of bee visitation observed here reveal the sometimes cryptic nature of important plant–insect interactions, support calls for broader recognition of polyploids as distinct units of biodiversity (Soltis et al., 2007; McIntyre and Strauss, 2017; Laport and Ng, 2017), and are consistent with assertions that unrecognized ploidal variation is important for conservation and biodiversity considerations (Severns and Liston, 2008).

AUTHOR CONTRIBUTIONS

R.G.L. designed the research and led writing with input from D.P. and R.L.M.; R.G.L. collected and analyzed the data with assistance from D.P.; R.L.M. identified bee specimens.

APPENDIX S1. R script and data for calculating Pianka’s niche overlap index.

Click here for additional data file.

APPENDIX S2. Supplemental description of flow cytometry methodology and example data.

Click here for additional data file.

APPENDIX S3. List of identified bees and accession numbers collected on sympatric diploid and tetraploid Larrea tridentata in 2014, 2015, and 2016.

Click here for additional data file.

APPENDIX S4. Abundance of bee species with >10 total occurrences on sympatric diploid and tetraploid Larrea tridentata.

Click here for additional data file.