Age-related variation in the oral microbiome of urban Cooper’s hawks (Accipiter cooperii)

Background Bird species worldwide are affected by trichomoniasis caused by the protist Trichomonas gallinae. In avivorous raptors such as Cooper’s hawks (Accipiter cooperii), nestlings are more susceptible than fledglings and adults. Previous research suggested a link between oral pH and susceptibility: the oral pH of fledgling and adult hawks is more than seven times more acidic than that of nestlings. We speculated that this age-specific difference in pH would correspond to age-specific differences in the oral microbiota of Cooper’s hawks. We examined the oral microbiomes of 31 healthy, wild Cooper’s hawks in Tucson, Arizona (USA). Individuals represented three age classes (nestlings, fledglings, and adults). We designed our study with multiple controls, replicated sampling, mock communities, and stringent quality-controls to address challenges that can limit the inferential quality of microbiome data sets. Results Richness of bacterial communities in oral cavities of Cooper’s hawks differed as a function of age but not as a function of sex, sampling date, or sampling location. Bacterial communities in oral cavities of nestlings differed from those of fledglings and adults, whereas communities in fledglings and adults did not differ from each other. Communities were similar in males and females and did not differ over the sampling season. Prevalence of acid-producing bacteria in fledgling and adults vs. nestlings is consistent with previous reports of age-specific variation in oral pH, but further research is needed to establish a causal link to pH levels or susceptibility to disease. Analyses of mock communities demonstrated high repeatability and showed that operon number and read abundance were highly correlated. Conclusions The oral microbiota of wild Cooper’s hawks differs between nestlings and older birds. Variation in the oral microbiome is consistent with differences in oral pH between nestlings and older individuals. Overall our study provides a first perspective on bacterial communities associated with oral cavities of a wild raptor. Electronic supplementary material The online version of this article (10.1186/s12866-019-1413-y) contains supplementary material, which is available to authorized users.


Background
Bird species worldwide are affected by diverse diseases including trichomoniasis, a disease caused by the flagellated protist Trichomonas gallinae [1,2]. Birds of prey that consume members of the family Columbidae are especially susceptible to trichomoniasis, as pigeons and doves are primary hosts for T. gallinae [3,4].
Consequently, infectious diseases such as trichomoniasis rank among the most common causes of morbidity in some birds of prey [5]. Susceptibility to trichomoniasis varies among age classes in several species of predatory birds. In general, young birds are more likely to contract the disease than older birds, often leading to high rates of nestling mortality [2,6,7]. For example, trichomoniasis was identified as the cause of mortality in ca. 40% of nestlings produced annually in a population of urban-nesting Cooper's hawks in Tucson, Arizona [2]. For susceptible Cooper's hawks, infection by T. gallinae results in lesions in the mouth and oropharynx, with subsequent spread of lesions to the esophagus and crop [8]. Reduced food consumption and intermittent regurgitation resulting from severe lesions leads to weight loss and ultimately can be fatal [9].
Previous research suggests that susceptibility of Cooper's hawks to trichomoniasis is related strongly to the age-specific pH of the oral cavity [10]. The mean pH of fluid in the oral cavity of nestling Cooper's hawks is 6.8, whereas that of fledglings and adults is 6.0-6.1 (i.e., is at least seven times more acidic) [10]. Trichomonads generally are sensitive to environmental pH (e.g., T. gallinae [11], T. gallinae and T. vaginalis [12], T. vaginalis [13], T. foetus [14]). Trichomonas gallinae thrives when pH is between 6.5 and 7.5 (optimum 7.2 [11,12]). When the pH is near optimum for trichomonads, some species (e.g., T. vaginalis) secrete a substance that facilitates invasion of host epithelial cells [15], but the ability to do so declines rapidly when pH is above or below the optimum.
Fluid in the oral cavity of Cooper's hawks becomes more acidic after birds have fledged and are nearing independence (≥50 days of age; [10]), but the reason for the change is unknown. Many animals undergo similar changes in body chemistry (i.e., a change in acidity) during maturation, and those changes often are associated with changes in their bacterial communities [16][17][18][19][20][21][22]. We speculated that the community of bacteria in the oral cavity of nestling Cooper's hawks differs from that in older hawks. Our objectives in this study were to provide a first perspective on the bacterial communities associated with the oral cavity of a wild raptor, and to compare the oral microbiomes of apparently healthy Cooper's hawks in three age classes, sampled contemporaneously from a shared geographic area.

Methods
We conducted the field portion of the study in Tucson, Arizona from April to July 2015. We sampled the oral microbiome of hawks from three age groups: nestlings (14-28 days old), fledglings (55-120 days old), and breeding adults (> 1 year old). We used a bal-chatri trap [23] to capture fledglings and adult males, and a dho-gaza trap [23] to capture adult females. Nestlings were taken by hand from their nests. We differentiated males from females by diameter of the tarsometatarsus [24]. All hawks were captured from nests monitored in a long-term study of population dynamics [25]. The University of Arizona's IACUC reviewed and approved all of the activities in our study that involved the capture and handling of adult, fledgling, and nestling Cooper's hawks (Protocol 12-329). Sampling locations are shown in Table 1 and Additional file 1: Figure S1.

Sampling the oral microbiome
We used a sterile foam-tipped applicator (Whatman, model WB100032) to swab the tongue and palate of each hawk. Prior to swabbing, we used scissors to trim the foam tip so it fit easily in the oral cavity. Scissors were soaked thoroughly with 95% ethanol prior to each use. We wore latex gloves and used a wooden tongue depressor to open the mouths of adult hawks. We used a new depressor for each hawk. Immediately after swabbing, we stored each foam tip in a sterile, 2 mL Eppendorf tube. We stored all tubes at − 80°C within 2 h of sampling. In total we sampled 34 hawks; of these, 31 yielded high-quality DNA extractions and were included in our analyses (Table 1).

Field controls
We used two types of controls in the field that allowed us to distinguish between bacterial associates of the oral cavity of hawks and those that might be captured incidentally during our sampling process. First, at each sampling event we retained a foam applicator that had been trimmed with scissors and exposed to the air as above, but not exposed to a hawk (field, air control). Second, we used an applicator that was trimmed as described above to swab a section of each tongue depressor that had not contacted a hawk or a person (field, wood control). These field-control swabs were processed contemporaneously with oral swabs as described below.

DNA extractions
We used half of each foam tip for DNA extraction, retaining the second half for archival purposes. We used the PowerSoil DNA Isolation Kit (Mo Bio Laboratories: Carlsbad, California, USA) to extract total genomic DNA from each half-tip. We followed the manufacturer's instructions, except that we added an initial incubation for 10 min at 65°C in lysis buffer. We used the Qubit 2.0 Fluorometer (Invitrogen: Carlsbad, California, USA) with the Qubit dsDNA HS (High Sensitivity) Assay Kit to quantify DNA concentrations following the manufacturer's instructions.

Library preparation and Illumina sequencing
We used a dual-barcoded, two-step PCR approach for high throughput amplicon sequencing on the Illumina MiSeq platform (Illumina, Inc.) [26]. In the first PCR (PCR1) we amplified the V4 region of the 16S rRNA with primer pair 515F/806R [27]. Each PCR1 primer contained a universal 22 basepair (bp) consensus sequence tag (i.e., CS1 forward and CS2 reverse), 0-5 bp for phase-shifting, a 2 bp linker, and the locus-specific primer (515F or 806R). We pooled forward or reverse PCR1 primers with different phase-shifting lengths in equimolar concentrations prior to PCR, such that each amplification consisted of a random mixture of different phase-shifting lengths. This allowed for the maximum amount of sequence diversity in the first four bases, which is critical for accurate cluster identification and color matrix estimation on the Illumina MiSeq [28]. We chose linkers and phase-shifting bases with low identity between primers and the target sequences.
In the second amplification step (PCR2) we used forward and reverse primers that each contained the complement of the CS tag, a 12 bp barcode, and the corresponding Illumina sequencing primer. The addition of barcodes in PCR2 allowed for maximum flexibility in the locus of interest or to include multiple targets in the same sequencing reaction, without needing to purchase a large number of barcoded target-specific primers.
We performed PCR1 in triplicate in 15 μl reaction volumes that each contained 0.5 μl DNA template, 7.5 μl 1X Phusion Flash High-Fidelity PCR Master Mix (Ther-moFisher Scientific, Austin, Texas, USA), 0.2 μl of 50 μM forward and reverse primers (505F and 806R Details for free-living Cooper's hawks for which the oral microbiome was sampled in Tucson, Arizona, USA. Thirty-four individuals were sampled, but high-quality DNA extractions were obtained for 31 individuals (listed here). Information for each individual includes the US Fish and Wildlife Service identification number (USFWS ID), sex, age class, capture date, capture location, capture method, diameter of the metatarsus, and the number of operational taxonomic units (OTU) of bacteria recorded from swabs of the palate and tongue. NA, not applicable: individuals were too young to have been assigned a USFWS ID, identified to sex, or measured for metatarsus diameter. Asterisk indicates one individual excluded from richness and community analyses due to low bacterial richness respectively) [27], 1 mg/mL of molecular grade bovine serum albumin (BSA; New England Biolabs, Ipswich, Massachusetts, USA), and 6.0 μl of molecular grade water. We used the following cycling protocol: an initial denaturing step at 98°C for 10 s; 25 cycles consisting of denaturation at 98°C for 1 s, annealing at 57°C for 5 s, and extension at 72°C for 20 s; and a final extension at 72°C for 1 min. We used sterile, molecular grade water instead of template for negative controls. We used SYBR Green 1 (Molecular Probes, Invitrogen) to visualize amplification on a 2% agarose gel after electrophoresis. We then pooled the three PCR1 products for each sample. From this, we diluted 5 μl of the pooled amplicons with molecular grade water to a final concentration of 1:15. We then used 1 μl of the pooled, diluted PCR1 product as the template for PCR2.
Each PCR2 reaction contained a final concentration of 1X Phusion Flash High-Fidelity PCR Master Mix, 0.075 μM barcoded primers (forward and reverse previously pooled at a concentration of 2 μM), and 0.24 mg/ mL of BSA, for a final volume of 20 μL. We used the following cycling protocol: an initial denaturing step at 98°C for 10 s; five cycles consisting of denaturation at 98°C for 1 s, annealing at 51°C for 5 s, and extension at 72°C for 15 s; and a final extension step at 72°C for 1 min.
We visualized products of PCR2 on a 2% agarose gel to verify amplification with minimal primer dimers. We quantified PCR2 products with PicoGreen and the Biotek Synergy H1 Multi-Mode Reader (Winooski, Vermont, USA). We normalized amplicons to 1 ng/μL and pooled 2 μL of each for Illumina sequencing. We purified the final amplicon pool with Agencourt AMPure XP beads (Beckman Coulter, Inc., Brea, California, USA) at a ratio of 1:1 to remove excess primers, nucleotides, salts, enzymes, and primer dimers, following the manufacturer's instructions.
We evaluated the amplicon library on a BioAnalyzer 2100 (Agilent Technologies) to determine concentration and fragment size distribution. Paired-end sequencing was performed on an Illumina MiSeq with Reagent Kit v3 (2 × 300 bp) according to protocols at the IBEST Genomics Core at the University of Idaho.

Laboratory controls
We complemented our field controls (i.e., air and wood controls) with several precautionary measures in the laboratory to limit contamination. We prepared PCR mixes in a sterile, dedicated "pre-PCR" hood, which never was exposed to amplified DNA. We decontaminated pipettes and all surfaces in the pre-PCR hood with DNA Away (Molecular Bioproducts, Inc., San Diego, California, USA) prior to each use. We used a dedicated "post-PCR" hood for all pre-sequencing steps after PCR1 (i.e., PCR1 pooling, dilutions, addition of diluted products to PCR2 master mix, and PCR2 amplicon pooling). We used DNA Away and treated the hood and all equipment with ultra-violet light for a minimum of 30 min prior to each use. We used sterile, aerosol-resistant pipette tips at all steps to minimize cross-contamination of samples, and we used separate reagents, pipettes, tips, and consumables for pre-and post-PCR setup. We pooled negative controls from PCR1 and used them as template for reactions in PCR2 to ensure that no contamination occurred during pooling or PCR2 setup. Finally, although we detected no contamination, we combined 5 μL of each PCR negative control, extraction blank, and field control in a separate pool and subjected them to the same pre-sequencing treatment as positive amplicon pools. We sequenced these negative controls in parallel with our samples.

Mock communities
To measure error rate, determine the consistency of quality among replicates, assess the correlation of operon number and read count, and filter sequence data, we amplified and sequenced bacterial mock communities in parallel with our samples. The mock communities were obtained through BEI resources (American Type Culture Collection, Manassas, Virginia, USA), the National Institute of Allergy and Infectious Diseases, and the National Institutes of Health as part of the Human Microbiome Project. We used two preparations: Genomic DNA from Microbial Mock Community B (Even, Low Concentration), v5.1 L, for 16S rRNA Gene Sequencing, HM-782D; and Genomic DNA from Microbial Mock Community B (Staggered, Low Concentration), v5.2 L, for 16S rRNA Gene Sequencing, HM-783D. Each preparation contained 20 phylogenetically diverse bacterial taxa either at even concentrations of 100,000 copies per organism per μL (even mock community) or staggered concentrations from 1000 to 1000,000 copies per organism per μL (staggered mock community). We amplified 1 μL of each mock community for PCR1 as above. To assess within-run variation, we used 1 μL of diluted PCR1 product (even and staggered) for five separate PCR2 amplifications, resulting in amplicons with five different barcodes for each mock community. Overall, we sequenced 10 replicate mock samples in parallel with oral swabs, field controls, and laboratory controls.
We used the UPARSE pipeline with USEARCH v.8.1.1861_i86linux32 [29] for all subsequent analyses. We used the fastq_eestats2 command to create a summary report showing the number of reads that would pass an expected error filter (maxEE 0.25, 0.50, 1.0) at different length thresholds (length_cutoffs 200, 250, 10). Based on this report we chose a length cutoff that would yield a sufficient number of high quality reads per sample while also retaining maximum read length. We used the fastq_filter command to trim forward (R1) reads to 200 bp (fastq_trunclen = 200) and to remove reads with greater than one error rate (max_ee = 1). Overall, 952,214 reads passed this quality filtering (94%).
We used the command derep_fulllength to dereplicate quality-trimmed reads and then removed singletons (parameters -sizeout -minuniquesize 1). We clustered dereplicated sequences (48,221 sequences) into operational taxonomic units (OTU) at 97% sequence similarity. In addition to de novo chimera checking performed during clustering, we used the RDP classifier for reference-based chimera checking of representative sequences for each OTU (16S rRNA reference database v.9). We used the command-usearch_global (id = 0.97) to map raw reads back to 478 chimera-checked OTU. From the resulting OTU table we removed 48 OTU that represented potential contaminants (i.e., OTU present in the PCR negative controls (16), extraction blanks (4), or field controls (11), or spurious OTU found in the mock community (17)). We retained 430 OTU for analysis (Additional file 2: Table  S1). Following quality control we retained high-quality data from 31 hawks (Table 1).

Data analysis
We used the RDP naïve Bayesian 16S rRNA classifier to estimate taxonomy of each OTU. We used analyses of variance (ANOVA) to compare read number among age classes (nestlings, fledglings, and adults), sampling dates (Table 1), and sampling locations (Table 1). We used a t-test to compare read number between male and female hawks. We used a Kruskal-Wallis test to compare richness among age classes. We used ANOVA to compare read number for replicate analyses of the even mock community. We used linear regression to compare read number and operon number for the staggered mock community. Normality was confirmed before all parametric tests. We used analyses of similarity (ANOSIM) based on presence-absence data (Jaccard's Index) and abundance data (Simpson's Index) to evaluate community composition, implementing two sets of analyses: one that included only the OTU represented in the data set by ≥1000 reads (55 OTU), and one that included all nonsingleton OTU (236 OTU). We used a Bonferroni correction to correct for multiple comparisons in each analysis and non-metric multidimensional scaling implemented in PAST [30] to visualize results.
We analyzed residuals from the previous analysis to determine whether richness differed as a function of sex, sampling date, or sampling location once differences due to age class were taken into account. We observed no difference in richness between male and female hawks (t 28 = − 0.3657, P = 0.7173), among sampling dates (F 9,20 = 1.3524, P = 0.2730), or among sampling locations (F 14,15 = 0.5887, P = 0.8356).
Analyses of mock communities demonstrated the consistency of results among replicates when operon counts were equal (even mock community, Additional file 3: Table S2), and a strong correlation of read number with operon number when inputs were unequal (staggered mock community, R 2 = 0.88; Additional file 1: Figure S2). The latter observation led us to use read count as a proxy for abundance in subsequent analyses, permitting analyses based on both presence/absence data and abundance.
To assess the quality of our inference, we evaluated whether differences in community structure among age classes could instead reflect differences in the sex of individuals (i.e., males vs. females) or the timing of sampling events. However, analyses restricted to fledglings and adults revealed that the oral microbiota of female and male hawks was similar (analysis of 55 OTU, as above; Jaccard's index, ANOSIM R = 0.0583, P = 0.1743; Simpson's index, R = 0.0306, P = 0.3171). Bacterial communities in these age classes did not differ as a function of sampling month (Table 1; Jaccard's Index, ANOSIM R = 0.1181, P = 0.1363; Simpson's Index, R = 0.0880, P = 0.1931). Analyses based on the full data set without singletons yielded the same results (Additional file 1: Figure  S4).
Read counts for 13 of the 55 most common OTU (i.e., 24.1%) were at least 2.5 times as common in fledglings than in nestlings ( Table 2). Twelve of these were also at least 2.5 times as common in adults than in nestlings ( Table 2). These 12 OTU included putative members of the genera Porphyromonas, Proteiniphilum, Veillonella, Peptoniphilus, Kingella, Parvimonas, and Suttonella, as well as representatives of the Corynebacteriaceae and Dermabacteraceae (Table 2). Seven additional OTU were at least 2.5 times as common in adults vs. nestlings but were not as common in the transitional fledgling stage, although sampling of the fledgling age class was relatively limited ( Table 2).

Discussion
The oral microbiota of urban Cooper's hawks differs between nestlings and older individuals. Differences in the oral microbiota could not be attributed to sampling month, location, or the sex of birds. These differences correspond to previously detected differences in the oral pH of Cooper's hawks [10]. It is unclear whether the previously documented change in oral pH between nestlings and older birds could be attributed to an increase in the abundance of acid-producing bacteria as birds mature, or whether the shift is endogenous and selects for the establishment of acidophilic and/or acid-producing bacteria in older individuals. However, we observed that communities in the oral cavities of fledgling and adult hawks have high abundances of acid-producing bacteria, providing a basis for further study.

Comparison of the oral microbiome of Cooper's hawks relative to that of other vertebrates
The oral microbiome of Cooper's hawks included bacterial phyla that are prevalent in the oral cavity of other vertebrates. The most common OTU in the present study represented Firmicutes, Actinobacteria, Proteobacteria, and Bacteroidetes. Firmicutes, Proteobacteria, and Bacteroidetes are especially common in the oral microbiome of healthy dogs [31]. These phyla also predominate in the oral microbiome of domestic cats [32,33] and humans [34], wherein Spirochaetes, Actinobacteria, and Synergisetes also are common. As in those vertebrates, Porphyromonas was particularly prevalent in the oral microbiome of Cooper's hawks, especially in older age classes ( Table 2). Previous studies on bird-associated microbiota have not focused on the oral cavity Previous work on the microbiomes of birds has focused primarily on cloacal [35], gut [36][37][38][39][40], and facial microbiomes [40]. Studies of vultures suggest an oral-gastrointestinal-fecal route for Clostridia and other microbial taxa, consistent with our observations of these bacteria in the oral cavity of Cooper's hawks [40]. Many of the dominant genera observed here are more consistently found in oral cavities of other vertebrates than in studies of other bodily microbiomes in birds [31][32][33][34][35][36][37][38][39][40].

Age-specific variation in the oral microbiota
Several studies focusing on microbial communities associated with birds have explored bacterial community composition as a function of age. For example, the gut microbiota of juvenile kakapos (Strigops habroptilus) differed from that of adults in the prevalence of Lactobacillus, but communities did not differ otherwise among individuals of different ages [41]. Cloacal microbiomes differed significantly between chicks and adults of kittiwakes (Rissa tridactyla) [35]. In poultry such as chickens (Gallus gallus domesticus), major changes in the cecal microbiome are observed several weeks after hatching relative to birds on the day of hatch [37]. Among 59 neotropical bird species, the age of individuals appears to be less important than taxonomic affiliation in predicting the composition of the gut microbiome, suggesting that age-related structure is more easily observed within bird species, rather than among them [36].
We found that nestlings differed markedly in the composition of their oral microbiomes relative to adults and fledglings. This difference is striking given the feeding behavior of this species. Early in the nestling phase, Cooper's hawks are fed by the adult female. She tears pieces of flesh from prey items delivered by the adult male and individually feeds each nestling. Late in the nestling phase, chicks are strong enough to dismantle prey delivered to the nest on their own. After fledging, prey items are delivered whole to fledglings who usually are perched near the nest. Prey species delivered to nestlings in Tucson are dominated by birds, 57% of which are members of the Columbidae [42]. Breeding adults consume the same prey items they feed their nestlings, so it is unlikely that changes in Thirteen of the 55 most common operational taxonomic units (OTU) were ≥ 2.5 times more common in fledglings (F) than in nestlings (N) (marked F > N in the fledglings vs. nestlings column, with the fold-difference in read number shown in parentheses). Of these, 12 also were more common in adults (A) than in nestlings (N) (marked A > N in the adults vs. nestlings column, with the fold-difference in read number shown in parentheses). OTU 17 was 2.3-fold more common in adults vs. nestlings and is marked with an asterisk. Seven OTU were at least 2.5 fold more common in adults than in nestlings, but were not common in fledglings (OTU 136, 168, 170, 157, 152, 160, 53, marked with ‡), potentially reflecting the relatively small number of fledglings sampled here and consistent with the lower richness observed in fledglings ( Fig. 1) diet among age classes play a role in the changes in oral microbiomes we observed.
Taxonomic composition and potential relevance to the pH of the oral cavity in Cooper's hawks Taxa such as Porphyromonas, Proteiniphilum, Parvimonas, Kingella, Suttonella, Peptoniphilus, and Veillonella were especially common in samples from fledglings and adult hawks but were less commonly observed in nestlings. Porphyromonas species are known to produce acidic products (e.g., butyric acid) [43]. Proteiniphilum strains grow at a pH similar to that of the oral microbiome of fledgling and adult hawks (ca. 6.0) and can produce acetic acid [44]. Parvimonas strains metabolize peptone and amino acids to form acetic acid [45]. Acid production is known in Kingella, Suttonella, and Peptoniphilus species [46][47][48], and Veillonella strains can produce acetic and proprionic acids [49] (also characteristic of some Corynebacteriaceae, Table 2). Although it is unclear whether these bacteria are responsible for the age-related change in pH in the oral microbiome of Cooper's hawks, the results are consistent with observed differences in pH reported previously [10]. In future work, isolating culturable members of the oral microbiota may be useful to support experiments evaluating whether it is pH per se, or microbial community composition, that correlates with the age-specific differences in susceptibility to T. gallinae observed previously in Cooper's hawks.

Methodological considerations
We designed our study to maximize quality control in the sampling, molecular analysis, and bioinformatics steps. Analyses of mock communities demonstrated the consistency of runs and a strong correlation of read number with operon number when inputs were unequal, as is anticipated for natural communities. Use of multiple field controls, replicated sampling within each age class, phylogenetically diverse mock communities prepared in two different ways, diverse measures to limit contamination in the laboratory, stringent quality-control thresholds, analyses that considered the most thoroughly sampled subset of the microbial community as well as the community as a whole, and analyses based on both presence/absence and read abundance allowed us to overcome some of the challenges that can restrict the inferential quality of microbiome data sets. Overall, our detection of robust differences in community structure among bacteria of the oral cavity in nestling hawks vs. fledglings and adults, and our testing of alternative explanations for patterns observed here (e.g., sex-related differences or differences as a function of sampling month) provide support for our conclusions.

Conclusions
Our observations of the age-structured oral microbiome of Cooper's hawks are consistent with established, age-related differences in oral pH between nestlings and older individuals. Previous work has shown that distinctive microbial communities are associated with diseases caused by trichomonads (e.g., in the human vaginal microbiome [50]). Although further study is needed to evaluate causal relationships among microbiome composition, oral pH, and trichomoniasis, our study provides a first perspective on the bacterial communities associated with the oral cavity of a wild raptor and sets the stage for further research to explore how susceptibility to parasites such as T. gallinae may be placed in a microbial community context.

Additional files
Additional file 1: Figure S1. Sampling locations in Tucson, Arizona, USA. Points are numbered by representative hawks, corresponding to hawk ID numbers shown in Table 1. Scale, 0.6 cm = 1 km. Figure S2. Strong positive correlation between read count and operon number for staggered mock community, used to confirm that read abundance could be used as a proxy for OTU abundance. Figure S3. Species-accumulation curve for the 55 most common OTU indicates thorough sampling of bacterial communities, providing a basis for statistical analyses presented in the text. Figure S4. Additional file 3: Table S2. Mock community, number of operons, and average reads observed for the even and staggered mock communities.