A total of 217,541 raw, barcoded amplicons of the V1–V2 region of the 16S rRNA gene (“pyrotags”), generated by Roche 454 FLX Titanium sequencing as described by Mattila et al. 
, were downloaded from the NCBI Sequence Read Archive (accession #DRA000526), trimmed to remove pyrosequencing adaptors, low quality base calls (<27 Phred score) and size-selected (≥350 bp) within CLC Workstation (www.clcbio.com
). When downloading the primer, barcode and adapter sequences from the Sequence Read Archive website, we noticed that the forward primer sequence deposited in the Archive had a “G” at the 12th
position substituted for an “M,” which is the IUPAC code for an “A/C” degeneracy, as described in the Lane et al. 
reference cited in Mattila et al. 
. Use of the forward primer with the “G” substitution in the raw read preprocessing, which includes a quality filtering step that excludes reads with errors in the primer sequence, resulted in loss of all reads. We changed the “G” in the primer sequence to an “M” and were able to proceed with nearly all of the reads to subsequent quality-filtering steps and sequence similarity clustering. A web-implementation of Pyrotagger 
was used to select for reads whose barcode and forward primer sequences were 100% error-free, had a Phred score of ≥27 over ≥90% of the sequence, were between 350–400 bp long and were not flagged as possible chimeras. High-quality sequence reads were clustered based on ≥97% sequence identity via uclust 
into operational taxonomic units (OTUs). Sequences representative of each OTU were taxonomically classified by blastn-based comparison to greengenes and silva databases within Pyrotagger as well as to a local NCBI non-redundant ‘nt’ nucleotide sequence database (
, downloaded June 29, 2011).
To remove noise from the data, including potential rare contaminants, we removed OTUs that did not meet the criterion of being present as at least two sequences in each of at least two samples; this resulted in elimination of only 0.14% of the reads. The resulting set of OTUs was used in diversity analyses (see below).
To assign taxonomic designations to OTUs, we used representative sequences, selected as the most frequent sequence within each OTU, in blastn searches against the NCBI non-redundant ‘nt’ nucleotide sequence database (March 29 2012). For classifying OTUs with top hits (at 99–100% identity) to a representative of one of the previously retrieved species from bees, we used the designations first applied by Babendreier et al. 
, and subsequently used by Cox-Foster et al. 
, Martinson et al. 
and Moran et al. 
. These classifications were applied if the match in blastn was 99–100% to a bee-associated bacterium and if blastn gave no other close matches. We then designated the species by examining the associated GenBank file of top matches. In the case of the previously studied bee-associated bacterial species, this designation is indicated in a “note” field for sequences from Martinson et al. 
and is given in the publication for sequences from Babendreier et al. 
In some cases, different OTUs had top matches, at 100% identity, to sequences assigned to the same species, due to the presence of several slightly different 16S rRNA sequences for a species in the database. In these cases, we grouped the OTUs into the same species, to facilitate visualization of the relative abundances of species present in each sample. We refer to these single or clustered OTU groups as species. (Previous studies have sometimes referred to these as ‘phylotypes’, in recognition that each may include distinct strains with different ecologies and gene sets). We note that sequences with <97% identity for the amplicon containing hypervariable V1–V2 regions (corresponding to these pyrotags) sometimes have >97% identity for the overall 16S rRNA, with the result that different OTUs may have top hits to sequences that would be considered the same species under a 97% criterion for species delineation. Also, we note that the 97% OTUs generated directly by clustering were used in analyses of diversity, as described below.
In order to clarify apparent discrepancies between conclusions of Mattila et al. 
and other findings for the honey bee gut community composition, we obtained representative sequences for the main taxonomic assignments as reported in of Mattila et al. 
from the corresponding author (I. Newton) of the study. We used these sequences in web-based implementations of the Ribosomal Database Classifier and SeqMatch 
(RDP database, release 10, update 28), and in blastn searches against the NCBI nucleotide database (searches performed on March 27, 2012). We recorded the SeqMatch classifications and also the top hits in the NCBI database, including the percent nucleotide identities, and compared these to the assignments and percent identities for the taxonomy assignments reported in of Mattila et al. 
. We designated species using similar criteria as used for representative sequences of our OTUs.
Our quality-filtering criteria yielded almost twice as many sequences from the Matilla et al. 
dataset, but fewer OTUs. We reanalyzed the gut-derived reads, using criteria for quality and clustering similar to those described by Mattila et al., using the Mothur platform 
with requirement of quality score >27 for 97% of the read, and length between 300 and 400 bp. OTUs were at 97% identity. Rare OTUs were excluded using the same rule as above.
We retested other conclusions of Mattila et al. 
regarding bacterial community patterns. We used our clusters, generated as described above, for these tests. We restricted our analyses to the bee gut communities, because communities in bee bread are distinct from those in guts and cannot be pooled in the same analyses. In our analyses we used the samples designated by “bgmdi” and bgsdi”, and we simplify these tags in this paper, using “M” and “S” for colonies with multiply and singly mated queens respectively, followed by a number designating the colony sampled. Colonies with multiply and singly mated queens are presumed to differ in genetic diversity among workers.
First, we tested a potential relationship between number of times the queen mated (singly versus multiply) and the gut bacterial communities, using several approaches. To determine whether colony genetic diversity affected the abundances of putative pathogens, we performed a Mann-Whitney U test on the ratio of putative pathogen numbers divided by total sequences per sample. The putative pathogen category included sequences from the Enterobacteraceae
and from M. plutonius
, similar to the pathogen category as defined in Mattila et al. 
(We agree with Mattila at al. 
that these Enterobacteriaceae
might speculatively be considered to be pathogens, mostly based on their erratic distributions among individual bees and colonies).
Second, we tested whether gut community composition or diversity differed between colonies with singly and multiply mated queens. For all community multivariate analyses, PC-ORD (version 4.25) was used 
. Analyses were conducted on the 42 OTUs that resulted from clustering the data, with elimination of extremely rare OTUs, as described below. Bacterial sequence reads per colony were standardized to the same sample size before diversity and multivariate community analyses were conducted. Standardization was carried out by randomly selecting 1230 reads (the smallest sample size) per colony using a custom Perl script. For nonparametric community structure and composition analyses, Multi-Response Permutation Procedures (MRPP) 
were used to test for differences among a priori
groups (singly and multiply mated queen colonies). Nonmetric Multidimensional Scaling (NMS) 
was used to visualize differences among honey bee gut microbiotas between bee colonies from both groups. Outlier Analysis 
identified bee colonies S5 and S7 as extreme outliers, and they were removed because such outliers can distort ordination solutions. All criteria, distance measures, and parameters chosen are similar to those in Hansen et al. 
. Briefly, a three-dimensional solution was chosen based on a combination of low stress, final instability, and P-values; 500 iterations were conducted. For ease of visualization, the two axes explaining the highest amount of variation in the dataset were presented in a 2-D figure. For alpha diversity analyses, colonies were considered subsamples of each treatment (N
10 and 12, respectively for singly and multiply mated queen colonies). Species accumulation curves were plotted with error bars (plus or minus two standard deviations) to contrast species diversity between the two treatments, and to evaluate the adequacy of dataset sample size. Subsampling was repeated 500 times for this analysis and the number of OTUs and average distance (Sorensen) for each subsample size was computed. Using EstimateS 
evenness of bacterial OTUs for each treatment was estimated using Simpson's dominance index (D) 
, and richness for each treatment was calculated using the Chao 2 richness index 
The above diversity analyses were repeated for the alternative OTU set, based on 166 retained OTUs.
Finally, we performed a Pearson's correlation analyses to test the proposed relationship between the number of Bifidobacterium cells and numbers of putative pathogens in the gut microbiota. We analyzed gut microbiota samples only for a relationship between numbers of sequences corresponding to Bifidobacterium and numbers of sequences corresponding to putative pathogens as defined above.