DOI: 10.5061/dryad.pg4f4qrzk
title: Data from: Bombus impatiens colony microsatellite genotype data for mating frequency analysis

description: The genus Bombus (bumble bees) includes approximately 265
 species, many of which are in decline in North America and Europe. To
 estimate the colony abundance of bumble bees in natural and agricultural
 habitats, sibship relationships are often reconstructed from genetic data
 with the assumption that colonies have one monandrous queen. However, some
 species such as the North American common eastern bumble bee (Bombus
 impatiens Cresson) can display low levels of polyandry, which may bias
 estimates of colony abundance based on monandrous sibship reconstructions.
 To accurately quantify rates of polyandry in wild and commercially mated
 queens of this species, we empirically estimated mating frequencies using
 a novel statistical model and genotypes from 730 bees. To genotype
 individuals, we used a highly polymorphic set of microsatellites on
 colonies established from 20 wild-caught gynes and 10 commercial colonies.
 We found multiple fathers in three of the wild colonies and three of the
 commercial colonies. This resulted in average effective mating frequencies
 of 1.075 ± 0.18 and 1.154 ± 0.25 for wild and commercial colonies,
 respectively. These findings agree with previous reports of low rates of
 polyandry for B. impatiens. Using a large empirical dataset, we
 demonstrate that assuming monandry for colony abundance estimation in
 species that violate this assumption results in an overestimation of the
 number of colonies. Our results emphasize the importance of studying
 mating frequencies in social species of conservation concern and economic
 importance for the accuracy of colony abundance estimation and for
 understanding their ecology and sociobiology.
description: <strong>Specimen Collection and
 Rearing</strong>   
 <em>Wild B. impatiens</em> gynes were collected between April
 and May of 2018 from three locations in Pennsylvania, USA: University
 Park, Newport, and Landisville. Twenty colonies established from these
 wild-caught gynes were then reared following standard protocols (Treanore
 et al. 2021). Colonies were held in a walk-in incubator under constant
 darkness, at  28–30 ºC, 60% RH, and supplied ad libitum with a 60% sucrose
 solution and honey bee-collected pollen. Colony maintenance and offspring
 collection were performed under red light. To quantify mating frequencies
 in wild colonies, we collected the first 10 emerging workers, and if
 possible, 10 additional workers approximately 3-4 weeks after to capture
 any differences in paternity in the later emerging workers. The queen was
 removed and frozen along with the rest of the workers after collections of
 the second cohort. Additional bees with unknown emergence dates were also
 collected for a total of 13-32 bees per wild colony. For the assessment of
 mating frequencies in managed colonies, ten commercial colonies were
 purchased from Koppert Biological Systems (Howell Michigan, USA).
 Commercial colonies were purchased and left in the field until late summer
 when the entire colony was removed from the field and frozen. The original
 queen and 14-21 workers per colony were removed for microsatellite
 genotyping analysis.   
 <strong>Quality Assessment of Microsatellite Loci
 Population Genetic Analysis</strong>
   We used 11 microsatellite markers
 previously developed for <em>B. impatiens </em>and other
 Bombus spp (Supp 2). Because the family relationships within our dataset
 violate assumptions of Hardy-Weinberg equilibrium, we used the genetic
 data published by McGrady et al (2021) to estimate parameters of allele
 frequencies of the source population for the 11 loci included in this
 study. We calculated the number of alleles per locus and evaluated the
 probability of linkage disequilibrium using the package poppr (v2.9.3,
 Kamvar, Tabima and Grünwald, 2014) in R Statistical Software (v4.1.2; R
 Core Team 2021). We performed null allele analysis according to
 Chakraborty et al. (1994) using the package PopGenReport (v3.0.7, Adamack
 and Gruber 2014). To evaluate deviations from Hardy-Weinberg Equilibrium,
 we performed a test with chi-squared comparisons using the R package pegas
 (v1.1, Paradis, 2010). Given that we included 11 loci for all tests, we
 used the corresponding Bonferroni adjusted α value of 0.0045 for Hardy
 Weinberg analysis and 0.0009 for linkage disequilibrium to correct for
 multiple comparisons.   
 <strong>Microsatellite
 Genotyping</strong>  
 After pinning, the right mesothoracic leg from each bee was
 removed and placed in a 96-well PCR plate for DNA extraction. A total of
 150 μl Chelex® 100 (5%, in milli-q H20) and 5 μl of 10 mg/mL Proteinase K
 were added to each well of the plate containing a leg. Each sample was
 then heated in a Mastercycler® pro thermocycler for 60 minutes at 55 °C,
 15 minutes at 99 °C, 1 minute at 37 °C, and 15 minutes at 99 °C. Extracted
 DNA was stored for up to 3 weeks in a 4°C refrigerator. We amplified 11
 microsatellite loci with different concentrations optimized for multiplex
 amplification (Supp. 2). Polymerase chain reactions were performed using a
 Mastercycler® pro thermocycler with the following cycle settings: 3 min 30
 s initial denaturation at 95 °C, followed by 30 cycles of 30 s at 95 ºC, 1
 m 15 s at an annealing temperature of 55 ºC and 45 s at 72 ºC, then a
 final 15-minute extension at 72 ºC with a hold temperature at 15 ̊C.
 Following amplification, PCR products were diluted 1:10 with molecular
 water, then sized on an ABI3730XL® DNA Analyzer (Applied Biosystems) with
 GeneScan LIZ 500 internal size standard (Applied Biosystems) at the Penn
 State Genomics Core Facility.  
 <strong>Paternity Reconstruction and Power
 Analysis</strong>  
 Fragment results were imported into Geneious Prime v.2022.0.1
 (Biomatters Ltd) for genotyping and alleles were scored manually. Dropout
 and mistyping error rates were estimated jointly with the posterior
 probability of the number of fathers on a colony-by-colony basis. To
 jointly estimate paternity and error rates while including known queen
 genotypes, we developed a Bayesian approach to paternity reconstruction
 based on the genotyping error model in Wang (2004). Briefly, the model
 assumes that observed genotypes are perturbed from true genotypes by two
 classes of errors: dropouts (ɛ1 where a heterozygous genotype appears
 homozygous) and mistyping (ɛ2 where a given allele appears to be another
 allele). The likelihood is calculated by summing (integrating) over
 possible maternal, paternal, and offspring genotypes conditional on a
 paternity assignment. We used a nonparametric prior (a Dirichlet process)
 on full sibling groups and developed a Gibbs sampler that targets the
 joint posterior distribution of error rates and paternity assignments. The
 model and Gibbs sampler are detailed in Supplementary Methods and are
 implemented in an R package (https://github.com/nspope/paternity_DP). To
 assess the accuracy of the method under realistic conditions, we applied
 it to 100 simulations for each combination of number of fathers (1 to 6)
 and error rates (0.01, 0.05, and 0.1 for both classes of errors but
 consistent across loci). Each simulation used parental genotypes that were
 sampled from the observed maternal allele frequencies from wild colonies;
 and produced observed genotypes for 20 offspring and a mother. In each
 simulation, one father was assigned the majority of offspring and each
 additional father was allocated only one offspring, as this represents the
 most challenging conditions for distinguishing genotyping errors from
 distinct paternal genotypes. To determine the impacts of different numbers
 of loci on the accuracy in jointly estimating paternity and error rates,
 100 simulations were repeated with 2, 4, 8, and 16 loci. To simplify the
 analysis, each locus was assigned a conservative number of 3 alleles at
 equal frequencies and a fixed error rate of 0.05 for both dropout and
 mistyping.   To reconstruct
 the paternity of our focal wild and commercial colonies, we applied our
 paternity reconstruction method to the empirical genotypes from each
 colony, after removing bees/loci with &gt;50% missing genotypes,
 monomorphic loci, and “drifter” bees. These ‘drifters’ are individuals who
 originated from a different mother colony and likely escaped their
 original colony and returned to a different colony during colony
 maintenance procedures. Because the genotypes of the queens were known,
 the identification of drifters was accomplished by including a second
 unobserved queen in the model. After paternity reconstruction, we
 calculated effective mating frequency– defined as where pi is the
 proportion of each offspring that belongs to father i (Starr 1984) – by
 averaging over 1,000 Markov chain Monte Carlo (MCMC) samples for each
 colony. All following statistical tests were computed in R v.4.1.1 (R Core
 Team 2020).   
 <strong>Colony Abundance Estimation From Sampled Worker
 Genotypes</strong>  
 To determine the influence of polyandrous or monandrous
 assumptions when recreating sibships from genetic data, we reanalyzed a
 previously published dataset of 6,306<em> B. impatiens</em>
 workers collected from 30 agricultural field sites in Pennsylvania (USA)
 (McGrady et al. 2021) that used the same loci and approach for genotyping
 used in our study. For each site, the genotypes of the workers were
 imported into COLONY v.2.0.6.5 (Jones and Wang 2010; Wang 2004) to
 reconstruct sibship relationships among individuals. We analyzed the data
 twice in COLONY: once with the assumption of female monandry, and once
 with the assumption of female polyandry, as the software only provides a
 binary parameter for mating frequency. For monandry, the number of
 full-sibling groups was recorded as the number of detected colonies, and
 the number of individuals comprising each colony was recorded. In order to
 simplify the analysis, we did not reject any full-sibling groups based on
 COLONY’s statistical probability of accuracy. For polyandry, the estimated
 number of mothers was recorded as the number of detected colonies to
 include half-siblings. Additionally, the number of workers comprising each
 colony was also recorded. Mating frequencies were calculated from the
 output of COLONY with polyandry to compare it with our empirical
 estimation of mating frequencies in <em>B. impatiens</em>. A
 Wilcoxon rank-sum test was used to determine whether mating frequency
 differed between monandrous and polyandrous colonies. Total colony
 abundance was estimated in R with the package CAPWIRE (Miller et al. 2005)
 for both monandrous and polyandrous colonies from each site. A Welch’s
 paired t-test was performed in R v.4.1.1 (R Core Team 2020) to determine
 if the number of detected and estimated colonies differed when calculated
 under assumptions of monandry and polyandry.
description: # Data from: *Bombus impatiens* colony microsatellite genotype data for
 mating frequency analysis
 [https://doi.org/10.5061/dryad.pg4f4qrzk](https://doi.org/10.5061/dryad.pg4f4qrzk) To accurately quantify rates of polyandry in wild and commercially mated queens of *Bombus impatiens*, the common eastern bumblebee, we empirically estimated mating frequencies using a novel statistical model and genotypes from 730 bees. To genotype individuals, we used a highly polymorphic set of microsatellites on colonies established from 20 wild-caught gynes and 10 commercial colonies. We found multiple fathers in three of the wild colonies and three of the commercial colonies. This resulted in average effective mating frequencies of 1.075 ± 0.18 and 1.154 ± 0.25 for wild and commercial colonies, respectively. These findings agree with previous reports of low rates of polyandry for *B. impatiens*. Using a large empirical dataset, we demonstrate that assuming monandry for colony abundance estimation in species that violate this assumption results in an overestimation of the number of colonies. Our results emphasize the importance of studying mating frequencies in social species of conservation concern and economic importance for the accuracy of colony abundance estimation and for understanding their ecology and sociobiology. ## Genotype data and description of the data and file structure The first sheet in the Excel file describes the microsatellite genotypes for each individual used in this study. Columns A and B describe the unique ID for each insect as well as the sample ID. Columns C-L are colony and individual information used to create the sample IDs. In column E, "Bee_type" WR stands for "worker random" describing workers of an unknown emergence date, W1 describes the early emerging workers and W2 describes the later emerging workers. Column I "Type_num" describes the number of individuals per worker/queen type. Columns M and N describe the plate number and the location on the well plate of each sample. Column O contains any important notes regarding each sample. Columns Q-T describe the dates that each step was performed as well as the PlateID for the sequencing facility. The final 11 columns U-AE each represent a microsatellite locus, and the genotypes for each individual at that locus. The numbers represent the length of the microsatellite repeat and these represent genotypes (2 for diploid females). In the Genotype data sheet, columns D, F, H, J, and L containing “_” are only there to separate the data and to be used in an Excel formula to form the SampleID in column B. The colors of these columns do not have any significance except to make reading the columns easier for the researcher. For any data cells that are left blank, this is because the data is not available. The 2nd tab, "AlleleFreqs" describes the frequency of each allele for each microsatellite locus.    This genotype data might be useful for future paternity analyses. These genotypes may be uploaded into any analysis software, including R, to evaluate relationships between individuals. However, this data may not be useful for population analyses since many of the individuals are related.