Symbiont densities vary across Cephalotes varians development – qPCR assays
In our first assessment of microbiome variation across Cephalotes ants, we utilized 16S rRNA-based qPCR assays to measure bacterial titer in four colonies of C. varians. Our two-way ANOVA detected no significant effect of colony on normalized bacterial density (F-value = 0.742; p=0.529; Additional file 2: Table S13). In contrast, this same ANOVA model revealed a significant effect of developmental stage/caste designations (F-value = 95.08; p < 2 x 10−16), suggesting variation in symbiont density between some combination of eggs; larvae at different stages of development; prepupae; pupae; soldiers; queens; or adult workers of varying age, as approximated by pigmentation measures.
Tukey’s post-hoc HSD tests allowed us to ascertain which of these stages or castes harbored varying bacterial loads. Eggs and young larvae (<2) had the lowest numbers of bacteria across the four studied C. varians colonies (Fig. 2a), with average 16S rRNA gene copy numbers (per ng total DNA) respectively estimated at 4.34 × 103 and 5.08 x 103. Unnormalized qPCR measures were comparable between eggs and our blank, negative control samples (8.63 x 102 vs. 5.00 x 102, respectively) suggesting that C. varians ants are nearly sterile upon hatching. The first detectable change in bacterial titer unfolded as larvae transitioned from small (<2 mm) to intermediate size (2-3 mm), with 16S rRNA copy number/ng DNA increasing to 2.82 x 104 (p-value for Tukey’s HSD<2 vs. 2−3 = 0.040; Additional file 2: Table S13). Titers increased again as larvae aged further, reaching 2.50 x 105 for individuals >3mm in length (p-value for Tukey’s HSD2−3 vs. >3 = 0), before falling drastically at the pre-pupal stage to 3.17 x 103 copies/ng DNA (p-value for Tukey’s HSD>3 vs. pre−pupae = 0). There was no detectable change in titer into the pupal stage, with normalized 16S rRNA copy number estimates hovering at 2.49 x 103 (p-value for Tukey’s HSDpre−pupae vs. pupae = 0.999). But callow workers (pixel value ≥ 60) showed a spike in symbiont density, with normalized copy number estimates rising to 1.31 x 105 (p-value for Tukey’s HSDpupae vs. adult−pixel>60 = 0). In addition, workers at this young adult stage showed some of the highest variability detected for symbiont titer (Fig. 2a). Adult workers at the next age class (pixel value = 40-60) showed an additional, but smaller increase in bacterial titer, with 5.45 x 105 16S rRNA copies/ng DNA (p-value for Tukey’s HSDadult−pixel≥60 vs. adult−pixel−40−60 = 2 x 10−7). Further aging was not detectably associated with titer changes, and workers at these older ages appeared to show greater consistency in 16 rRNA copy number when compared to callows (Fig. 2a; Additional file 2: Table S13).
The most mature workers (pixel value = 0-20) showed pigmentation overlap with one alate queen, one wingless queen, and ten soldiers (Additional file 2: Table S4), suggesting similar age. Comparisons of bacterial titers among these adult categories suggested little difference in symbiont quantities, with 16S rRNA copy numbers/ng DNA equaling 3.57 x 105 (alate) and 1.14 x 105 (wingless) for the single queens, and averages of 4.03 x 105 and 4.18 x 105 for soldiers and workers. Accordingly, p-values from pairwise Tukey’s HSD comparisons were non-significant, ranging from 0.997 – 1.0 (Additional file 2: Table S13).
Symbiont presence/absence varies across Cephalotes development – standard PCR assays
After discovering the above-described developmental time course for symbiont proliferation, and ascertaining that all female castes harbor symbionts during adulthood, we used standard PCR amplification with universal 16S rRNA gene primers as a crude approximation of symbiont presence/absence across 522 ants spanning various castes and stages of C. varians and twelve additional Cephalotes species (Fig. 2b; Additional file 2: Table S3). Eggs from three species (n=13 samples) failed to amplify with universal bacterial PCR primers. But amplification did occur for 20.6% (n=34) of young larvae (i.e. normalized size <2) from across four Cephalotes species. In contrast, 74.6% (n=47) of middle-aged larvae from 11 species (i.e. standardized size class: 2-3) yielded 16S rRNA PCR amplification, as did 94.5% (n=73) of older larvae (i.e. standardized size class: >3) from 12 cephalotines. All 37 pre-pupae, sampled across seven Cephalotes, failed to amplify with universal bacterial primers, while amplification was seen for only 1/57 pupae obtained from three host species. At the adult stage, all wingless queens (n=18, from 6 species), alate queens/gynes (n=7, sampled from 4 species), and soldiers (n=57, from 11 species) yielded PCR amplification in these 16S rRNA assays, as did 177/178 workers (from 13 species).
Whole ant community composition across castes and stages of Cephalotes varians
With a sequence of bacterial acquisition, proliferation, loss, and reacquisition being suggested by the aforementioned data, we next set out to understand the composition of symbiotic bacterial communities across castes and development. We began with a study of C. varians, subjecting 145 samples with PCR-detectable bacterial communities to amplicon sequencing of the bacterial 16S rRNA gene using Illumina MiSeq technology. We removed one library where the number of reads was less than 1,000, and focused on the remaining 144. After removing contaminants, chimeras, and non-bacterial sequences, 1,597,177 quality-controlled sequences remained for further analysis, an average of 11,092 sequences per library (see unique sequence, OTU, and oligotype tables in Additional file 2: Tables S6-S14).
Permutational multivariate ANOVA statistics indicated that some combination of species composition and/or relative abundances of bacterial community members differed significantly across developmental stages of C. varians (Fig. 3a; Additional file 2: Table S14; Bray-Curtis: F = 129.617, R2 = 0.415, P = 0.001; Jaccard: F = 89.087, R2 = 0.342, P = 0.001), and also among colonies (Bray-Curtis: F = 7.813, R2 = 0.075, P = 0.001; Jaccard: F = 6.329, R2 = 0.073, P = 0.001). Across all larval libraries from four C. varians colonies, six dominant 97% OTUs, including Rhizobiales (OTU001 and OTU004), Lactobacillales (OTU007, OTU022 and OTU025), and Enterobacteriales (OTU006), accounted for an average of 87.1% of reads. In very young larvae (i.e. normalized body length < 2), Rhizobiales OTU004 was the only dominant bacterium. Another type of Rhizobiales, from OTU001, became common subsequently, beginning in middle-aged larvae (i.e. normalized body length 2-3) and continuing into the oldest larval stages (i.e. normalized body length >3). Lactobacillales OTUs 007, 022, and 025, as well as Enterobacteriales OTU006, were the most abundant taxa in larvae assigned to the oldest age class, but were common also in some middle-aged C. varians larvae.
With the exception of one mature queen, discussed below, microbiomes of C. varians adults were quite different from those of larvae (e.g. Fig. 3a). To begin, Lactobacillales and Enterobacteriales were extremely rare in workers, soldiers, and one alate queen/gyne, while Rhizobiales were ubiquitous, but only modestly abundant. Instead, adults harbored OTUs found rarely in larvae, including highly abundant OTUs from the orders Opitutales (OTU003) and Xanthomonadale (OTU002 and OTU005), seven frequently present OTUs from the Burkholderiales, and individual OTUs from the orders Pseudomonadales, Flavobacteriales, Campylobacterales, and Sphingobacteriales. The dominant OTUs from these orders accounted for an average of 93.95% of the sequence reads in each library.
We re-ran our permutational multivariate ANOVA analysis using a 97% OTU table containing samples from the most darkly pigmented and, hence, most mature workers and soldiers (pixel value <20). Our assessment of caste, however, revealed no differences among these two groups, whether beta diversity was measured through the Bray-Curtis or Jaccard index (Additional file 2: Table S14; Bray-Curtis: F = 2.376, R2 = 0.043, P = 0.074; Jaccard: F = 2.023, R2 = 0.029, P = 0.107). Our remaining caste sampling lacked sufficient replication for statistics. But anecdotally, a virgin winged queen/gyne (colony YH100) harbored a similar gut bacterial community in comparison to its sibling workers and soldiers, while the previously mentioned mature, wingless queen (colony YH091) was dominated by Rhizobiales OTU004, similar to some young larval samples from the same colony (Fig. 3a).
Principal Coordinates Analysis (PCoA) on the full 97% OTU table reinforced the above results (Fig. 3b). For example, gut communities of larvae were separated from those of workers and soldiers, along the first PCoA axis (46.43%). In fact, only two adults harbored microbiomes clustering with microbiomes of larvae, including a worker dominated by Lactobacillales OTU007, and the previously described wingless queen. Matching the compositional patterns from Fig. 3a, the symbiotic community of the alate queen/gyne clustered along the first axis with all worker samples and all but one of the above-described solider samples.
Bar graphs suggested high microbiome variability and occasionally unusual composition among the youngest, most lightly pigmented workers (Fig. 3a). To pursue this pattern, we ran multivariate homogeneity of variance tests separately in the three C. varians colonies with adequate sampling across stages (Fig. 3c). We found that the average distance to the centroid for the youngest worker group (pixel value = 40-60), and occasionally members of the next youngest cohort (pixel value = 40-60), was significantly higher than for the older worker age groups, in two of the three colonies with sufficient sample sizes (Additional file 2: Table S15). While still somewhat preliminary, these patterns suggest higher variability in microbiome composition during early adulthood.
Gut localization of larvae-associated symbionts in Cephalotes varians
Microbes dominating amplicon sequencing libraries from gut tissue and solid gut contents of three older C. varians larvae were highly similar to the communities seen in the above-described amplicon sequencing of whole-larva DNA extractions (Additional file 2: Figure S2 vs, Fig. 3a). Highly abundant were Rhizobiales OTUs 001 & 004, Lactobacillales OTU007, and – in one case – Enterobacteriales OTU006. Present at lower titers were OTUs from adult-enriched symbionts in the Xanthomonadales and Burkholderiales. Similarly present at lower relative abundance was OTU025 from the Lactobacillales. Sequences from additional OTUs made up <10% from each library. OTU004 (Rhizobiales) was consistently more abundant in all three gut tissue vs. gut content comparisons. The remaining, aforementioned OTUs showed similar relative abundances across the two sample types. Given the small number of dissected samples, these trends require further study.
Whole ant microbiome composition across castes and stages of 13 Cephalotes species – 97% OTU distributions
We further examined whether patterns of community structure found in C. varians were consistent across the Cephalotes genus, expanding our study to include 12 additional species. A total of 209 further 16S rRNA sequence libraries were generated toward this end, with an average of 8714 sequences per library after quality and contaminant filtering. Our remaining analyses included 192 of these QC’d libraries (i.e. all those except for libraries from 16 callow adults and 1 pupa). Focusing on 97% OTU clusters, a permutational multivariate Analysis of Variance revealed significant differences in the bacterial communities between larvae and adults (Additional file 2: Table S16; Adonis permANOVA tests: p < 0.05) for all Cephalotes species, except C. targionii. This was true for both Bray-Curtis and Jaccard measures of beta-diversity.
To characterize the taxa driving these differences across Cephalotes development, we constructed a heatmap showing 40 abundant 97% OTUs in all 336 whole ant sequence libraries (Fig. 4), in addition to a heatmap illustrating average 97% OTU relative abundance for each stage and caste per colony (Additional file 1: Fig. S3). Similar to patterns for C. varians, larval samples of different Cephalotes species harbored overlapping, yet distinct symbiotic communities compared to adults, with most adult-enriched symbionts exhibiting rarity, or low-abundance, in larvae (Additional file 2: Table S17). Using a relative abundance threshold of 0.00294 to ascertain presence/absence (explained below), and focusing on OTUs found in least 30% of the 154 sampled adult workers (Additional file 1: Fig. S4), we found that the most depleted OTUs came from the Opitutales (OTU003 – 98.1% of workers vs. 8.7% of larvae), Comamonadaceae (OTU020 – 74.0% of workers vs. 8.7% of larvae), Xanthomonadales (OTU005 – 70.8% of workers vs. 1.9% of larvae), Pseudomonadales (OTU013 – 69.5% of workers vs. 4.9% of larvae), Flavobacteriales (OTU028 – 46.8% of workers vs. 1.0% of larvae), Campylobacterales (OTU021 – 40.3% of workers vs. 3.9% of larvae), and Sphingobacteriales (OTU030 – found in 33.1% of workers vs. 0% of 103 larvae). Showing more intermediate patterns were Alcaligenaceae OTUs 010 (93.5% of workers vs. 35.0% of larvae), 012 (97.4% of workers vs. 40.8% of larvae), and 016 (83.8% of workers vs. 19.4% of larvae), along with Xanthomonadales OTU002 (74.0% of workers vs. 17.5% of larvae).
Using the same abundance threshold for presence/absence assignment, and focusing on OTUs found in at least 10% of larval samples (Additional file 1: Fig. S4; Additional file 2: Table S17), OTUs enriched in larvae included Lactobacillales OTUs 007 (3.2% of workers vs. 36.9% of larvae), 022 (2.6% of workers vs. 28.2% of larvae), 025 (0% of workers vs. 35.9% of larvae), 082 (0% of workers vs. 11.7% of larvae), and 086 (0% of workers vs. 10.7% of larvae). Enterobacteriales OTU006 (1.3% of workers vs. 45.8% of larvae) and OTU033 (0% of adults vs. 13.6% of larvae) showed a similar trend, as did Actinobacteria OTUs 062 (Actinomycetales; 0.6% of workers vs. 10.7% of larvae) and 080 (Bifidobacteriales; 0% of workers vs. 16.5% of larvae).
Distinct from these patterns were trends for Rhizobiales symbionts, which were common in both adults and larvae. OTU001 was found, for instance, in 100% of workers and 79.6% of larvae. OTU004 was in 76.6% of workers and 83.5% of larvae. Though less prevalent, OTUs 017 (15.6% of workers vs. 19.4% of larvae) and 053 (11.0% of workers vs. 17.5% of larvae) showed similar patterns.
Maximum and average relative abundance of OTUs across development revealed qualitatively similar trends, indicating that adult-enriched symbionts with low- to intermediate-prevalence across larvae are usually found at low relative abundance within these immatures (Additional file 1: Fig. S4; Additional file 2: Table S17). Two Alcaligenaceae OTUs, 010 and 012, provided interesting exceptions – with maximum relative abundances in larvae (0.309 and 0.157) either exceeding or approaching those in workers (0.121 and 0.188). Rhizobiales, with their dual-stage prevalence, showed similar patterns, with greater average and maximum relative abundance in larvae vs. adults, as exemplified by OTUs 001, 004, and 053, but not OTU017. Larvae-enriched symbionts were, conversely, most often rare in adults when detected, with the exceptions of OTU007 (Lactobacillales) and OTU006 (Enterobacteriales) – both reaching sporadically high maximum relative abundance in a small number of workers (0.506 and 0.578).
Many of the above general trends of larval vs. adult differences – gleaned from adult workers – were recapitulated when comparing alate queens/gynes and soldiers against larval immatures (Figs. 3; Additional file 1: Fig. S3; Additional file 2: Table S17). But while these adult castes harbored fairly similar microbiomes there were some key exceptions. In C. depressus and C. pusillus, for instance, Lactobacillales were seemingly enriched in soldiers relative to workers. Furthermore, in the former species a rare bacterium from the Neisseriales was also common in soldiers but rare in workers. In contrast, Campylobacterales appeared more abundant in workers than in soldiers in 5 out of 7 colonies harboring this microbe, suggesting interesting avenues for future study.
When expanding our adult caste comparisons to include five mature queens, sampled from five Cephalotes species, we saw considerably greater microbiome divergence (Additional file 1: Fig. S3). With microbiomes bearing little resemblance to those of other adult counterparts, symbiotic communities of wingless queens instead, resembled those of some larvae, with dominance by a subset of Rhizobiales (i.e. OTUs 004 and 017). Overall, mature queen microbiomes had lower diversity, with an apparent absence of many of the adult-enriched OTUs present in their mature offspring. Representing a slight exception to the Rhizobiales-enrichment trend was a queen from C. depressus, possessing such bacteria at lower titers, but showing enrichment, instead, for the same Neisseriales and Lactobacillales found in abundance within their soldiers (see above). The mature queens from C. pusillus and C. maculatus colonies uniquely harbored Wolbachia, in addition to Rhizobiales. This was a conspicuous pattern for the latter given this symbiont’s apparent absence from the queen’s presumed larva, soldier, and worker offspring.
Whole ant microbiome composition across castes and stages of 13 Cephalotes species – oligotype distributions
To understand strain-level diversity in Cephalotes-associated bacterial communities across developmental stages and castes, we looked at 16S rRNA oligotype distribution within the 40 focal symbiont OTUs from Fig. 4 (Additional file 1: Fig. S5). Although slowly-evolving genes, like 16S rRNA, obscure higher strain-level resolution, we detected several distinct oligotypes within most dominant OTUs. Similar to our 97% OTU result, permutational multivariate Analysis of Variance revealed that larvae and adults from 11 out of 12 species harbor significantly different bacterial communities at this narrower clustering width (Additional file 2: Table S16; Adonis permANOVA tests: p < 0.05).
Despite these differences, we found identical oligotypes across siblings from these stages for a handful of trans-developmentally shared OTUs, including those from the Rhizobiales (OTUs 001, 004, 017, 023, 042 and 056). Sibling workers, soldiers and winged queens from the same colony similarly shared oligotypes, suggesting a homogenizing effect of colony on microbiome composition. Microbiomes of wingless queens remained non-diverse at the oligotype-level. And those sharing 97% OTUs with colony-mates shared oligotypes with them, as well. The sole exception came for a single oligotype belonging to OTU004, found in the mature queen, and n=1 larva of C. varians colony YH091, but in no other members of that social unit (Additional file 1: Fig. S5).
Whole ant microbiome composition across castes and stages of 13 Cephalotes species – enterotyping analysis inferred from 97% OTU composition
With our results suggesting impacts of development and caste on microbiome composition, we aimed to further characterize symbiotic bacteria of Cephalotes ants through enterotyping analysis on our 336 whole ant Illumina sequencing libraries. Results suggested the presence of three community types, or “enterotypes” (Fig. 5a), with Calinski-Harabasz values reaching a peak of 0.829 at n=3 clusters, compared to CH values of 0.662 and 0.781 under models with n=2 and n=4 (Additional file 2: Table S11). Enterotype 1 was mainly comprised of microbiomes from workers, soldiers, and 14 out of 16 alate queens/gynes (Fig. 5c). Indicator OTU analyses revealed that this cluster was dominated by 97% OTUs from the Opitutales (n=1 OTU), Burkholderiales (n=6 OTUs), Xanthomonadales (n=2 OTUs), Flavobacteriales (n=1 OTU), Campylobacterales (n=1 OTU), and Sphingobacteriales (n=1 OTU), with the latter taxon being concentrated in C. varians (Fig. 5b; Additional file 2: Table S12).
Enterotype 2 included microbiomes of larvae, four out of five mature/wingless queens, and two out of sixteen alate queens/gynes, in addition to very small numbers of workers and soldiers (Fig. 5c). Characteristic of this cluster were Rhizobiales OTUs 004 and 017 (Fig. 5b). Microbiomes assigned to the third enterotype were mostly from larvae. However, one hailed from a wingless queen/gyne, while a few belonged to soldiers and workers. Rhizobiales OTU001 was identified as an indicator OTU for this enterotype, as was Enterobacteriales OTU006, along with three Lactobacillales OTUs (Fig. 5b; Additional file 2: Table S12). Estimating larval developmental stage through body length measures, we noticed that enterotype 3 was dominated by microbiomes of older (i.e. larger) larval samples, while the mature/wingless queen-including enterotype 2 was primarily comprised of young larvae (Fig. 5c).
Correlations between host phylogeny and cephalotine microbiomes differ across enterotypes
Principal coordinates analysis plots (Fig. 6) suggested that host phylogeny plays a significant role in predicting community structure for ants hosting enterotype 1-clustering microbiomes, but not clearly those clustering with enterotypes 2 or 3. With an interest in understanding whether cephalotine microbiome composition is predicted by evolutionary history, we investigated the impact of host phylogeny and geographic location on community beta-diversity using partial Mantel tests. Bray-Curtis and Jaccard Index dissimilarity matrices computed for enterotype 1-classifying samples, were significantly correlated with the host phylogenetic distance matrix, at 97%, 98%, and 99% OTU clustering widths, and at the oligotype level, after controlling for the effect of geographic distance (P ≤ 0.001; Table 1).
Table 1
The r value and p value of Mantel's test showing the influence of host phylogeny in controling of geography or the influence of geography in controling of host phylogeny for each enterotype.
Enterotype ID | Effect | Controlling for | 97% OTU | 98% OTU | 99% OTU | Oligotype |
Partial R2 | P value | Partial R2 | P value | Partial R2 | P value | Partial R2 | P value |
Based on Bray-Curtis distance matrix |
Enterotype1 | Host phylogeny | Geography | 0.457 | 0.001 | 0.515 | 0.001 | 0.496 | 0.001 | 0.519 | 0.001 |
Geography | Host phylogeny | 0.297 | 0.015 | 0.447 | 0.003 | 0.284 | 0.01 | 0.157 | 0.141 |
Enterotype2 | Host phylogeny | Geography | 0.079 | 0.537 | 0.196 | 0.143 | 0.111 | 0.353 | 0.130 | 0.081 |
Geography | Host phylogeny | 0.132 | 0.220 | 0.429 | 0.005 | 0.170 | 0.151 | 0.176 | 0.15 |
Enterotype3 | Host phylogeny | Geography | 0.017 | 0.563 | 0.026 | 0.858 | 0.068 | 0.66 | 0.045 | 0.776 |
Geography | Host phylogeny | -0.129 | 0.200 | -0.034 | 0.816 | -0.066 | 0.688 | 0.338 | 0.025 |
Based on Jaccard distance matrix |
Enterotype1 | Host phylogeny | Geography | 0.646 | 0.001 | 0.676 | 0.001 | 0.656 | 0.001 | 0.636 | 0.001 |
Geography | Host phylogeny | 0.477 | 0.005 | 0.547 | 0.001 | 0.452 | 0.003 | 0.363 | 0.008 |
Enterotype2 | Host phylogeny | Geography | 0.363 | 0.004 | 0.434 | 0.004 | 0.562 | 0.001 | 0.262 | 0.039 |
Geography | Host phylogeny | 0.033 | 0.801 | 0.252 | 0.035 | 0.048 | 0.69 | 0.131 | 0.328 |
Enterotype3 | Host phylogeny | Geography | -0.006 | 0.959 | 0.193 | 0.225 | 0.375 | 0.021 | 0.372 | 0.02 |
Geography | Host phylogeny | 0.204 | 0.190 | 0.222 | 0.176 | 0.301 | 0.051 | 0.257 | 0.106 |
Geographic distance explained a smaller proportion of community dissimilarity for enterotype 1. Specifically, when using Bray-Curtis distances the average r-squared value across clustering widths was 0.296 for geography compared versus 0.497 for phylogeny. For Jaccard index-based analyses, the average r-squared value for geography was 0.460 compared to 0.654 for phylogeny. But like phylogeny, geography did show a significant correlation with community dissimilarity at every OTU clustering width, along with the oligotype level, for Jaccard index analyses. Significant correlations between geographic distance and community dissimilarity were evident for 97%, 98%, and 99% OTUs, but not the oligotype level, when analyzing Bray-Curtis distances (P < 0.05; Table 1).
For enterotype 2, host phylogenetic distance was significantly correlated with Jaccard dissimilarity after controlling for geographic distance, at all sequence clustering widths (P < 0.05; Table 1). These correlations were not, however, significant at any level of sequence clustering when using the Bray-Curtis dissimilarity matrix. Average r-squared values for the effect of phylogeny were slightly lower, for Jaccard-based analyses on enterotype 2 (0.405) compared to those on enterotype 1 (0.654). For this second enterotype, correlations between community dissimilarity and phylogeny-controlled geographic distances were only significant at the 98% clustering level.
For enterotype 3, host phylogenetic distance and community dissimilarity were not significantly correlated after accounting for geographic distance, except for a significantly positive correlation for the Jaccard distance computed at the 99% OTU and oligotype levels. Average r-squared values, however, were again comparatively low versus those from enterotype 1 (0.374 vs. 0.646), when averaged across these significant clustering widths. In addition, the correlations between geographic distance and community dissimilarities for enterotype 3 were only significant at the oligtype level for Bray-Curtis measures (Table 1).
We experimented with partial Mantel tests on enterotype 3, repeating analyses with only specialized, adult-enriched bacteria (defined in subsequent sections) or with non-specialized bacteria likely acquired from the environment. Focusing on 99% and oligotype levels, and beta-diversity metrics computed with the Jaccard index, we found that phylogenetic distance remained a significant predictor of microbiome dissimilarity when only specialized bacteria were included. For these results r-squared values of significant findings ranging from 0.335 – 0.502. When only non-specialists were included prior to Jaccard distance matrix calculation, phylogeny was not significant (Additional file 2: Table S18). Combined with results discussed further below, these results suggest a strong effect of host phylogeny on microbiomes made up primarily by specialized Cephalotes symbionts (enterotypes 1 & 2), but not on those comprised largely of free-living bacteria (enterotype 3).
Phylogenetic affinities of abundant bacterial species from Cephalotes larvae
To initiate a study on the evolutionary origins of larva-enriched, and in some cases, newly identified Cephalotes-associated symbionts, we used maximum likelihood analysis for phylogenetic inference on abundant sequences from the Rhizobiales (Analysis #1), Enterobacteriales (Analysis #2), and Lactobacillales (Analysis #3). Our Rhizobiales phylogeny indicated that 16S rRNA gene sequences from such microbes grouped into two monophyletic, Cephalotes ant-specific lineages (Fig. 7a). Symbiont strains from Rhizobiales OTU001 formed a distinct clade, which was sister to a ‘crown group’ clade comprised primarily of symbionts from other ants. Nested within this crown group was a second lineage comprised exclusively of Cephalotes associates, including OTUs 004, 017, 023, 042, and 056.
Our maximum likelihood phylogeny of Lactobacillales bacteria (Fig. 7b) revealed that 16S rRNA gene sequences from the more common OTU007 (sister to Lactobacillus homohiochii and L. fructivorans) and OTU022 (sister to two uncultured bacteria) formed two distinct and well-supported Cephalotes-specific clades, separated from free-living bacteria by modestly long branches. While dominant in Cephalotes with enterotype 3-like microbiomes, bacteria from these larvae-dominant lineages were found in 23 of the 24 adult ants harboring Lactobacillales. Symbiotic strains from OTU025 Lactobacillales, also common in larvae, but seemingly absent from adults, did not form a host-specific clade, instead showing close relatedness to free-living Enterococcus bacteria.
In contrast to the patterns for Rhizobiales and Lactobacillales, 16S rRNA gene sequences from three Enterobacteriales OTUs, dominant in Cephalotes larvae with enterotype 3-like microbiomes, did not form host-specific clades, grouping instead with a range of free-living bacteria from this order (Fig. 7c). In particular, sequences from OTU006 showed relatedness to Erwinia, Enterobacter, Salmonella, Pantoea, and Citrobacter. OTU033 showed phylogenetic affinity to Serratia. The more distantly related OTU063 appeared related to Morganella.
Identifying host-specialized vs. environmentally acquired bacteria – BLAST and phylogenetic analyses
Frequent relatedness to free-living, non-cephalotine associated bacteria suggested that several of the bacteria common within larval microbiome are environmentally acquired, unlike the vast majority of microbes from adults. To identify and compare the fractions of stage-enriched microbiomes that are likely acquired from the environment (i.e. highly related to free-living bacteria) vs. seemingly specialized and transmitted through evolved mechanisms for ancient timespans (e.g. Lanan et al. [41]), we began by performing BLASTn searches on the most common (“representative”) sequence found in each of the 154 97% OTUs exceeding pre-defined thresholds (i.e. 0.00294 in at least 2 libraries, or 0.05 in one; Additional file 2: Table S19). Fifty-one of these yielded a top hit to a bacterium previously determined to fall within a cephalotine-specific clade [7, 62, 66]. Across these instances the average sequence identity was 98.0%.
Of the remaining 103 OTUs we performed phylogenetic analysis on 67 of the most abundant (Additional file 1: Figs. S6-10), using inclusion criteria, and BLASTn strategies for identifying close relatives that are described in the supplementary figure legends. Analyses from the Actinobacteria (Additional file 1: Fig. S6) – found commonly in larvae (Additional file 1: Fig. S3) – revealed that OTU062 was related to previously identified gut- and head-associated sequences from Cephalotes. Showing some relatedness to an uncultured bacterium from rainwater, and more distant relatedness to Agromyces and Curtobacterium, we do not yet see strong evidence that this lineage is a long-term specialist given that its only identified cephalotine host is, thus far, C. varians from the Florida Keys. Furthermore, its 98.4% sequence identity to a non-cephalotine bacterium (accession #: KY874657) puts it at odds with most cephalotine-specific bacteria, previously shown to exhibit an average of 93.3% relatedness to their closest non-cephalotine match [65]. Only two other monophyletic clusters of cephalotine ant-derived sequences were observed on this tree. While both were distributed across multiple Cephalotes species, they were also comprised of only two unique sequences, showing close relatedness to bacteria from other habitats. These trends argue against high likelihood of ancient and specialized relationships between cephalotines and Actinobacteria, once again contrasting with results for most adult-enriched bacteria [7, 62, 66].
Using similar logic, we found minimal evidence for additional cephalotine-specialization among the Lactobacillales (Additional file 1: Fig. S7). Beyond the potentially specialized OTUs 007 and 022, only one other cephalotine-specific clade was detected. Consisting of two sequences from OTU085 with 98.81% relatedness to non-cephalotine bacteria (accession #: AB934560.1), this is also unlikely to be an ancient or widespread specialist.
In contrast, the three unique sequences from the ambiguously classifying Pseudomonadales OTU0046, were nested within a cephalotine-specific lineage, albeit with low bootstrap support (62%; Additional file 1: Fig. S8). Hailing from C. wheeleri and C. goniodontus, they showed relatedness to symbionts from three other Cephalotes species, including C. setulifer, C. rohweri, and C. atratus. OTU116, from the Bacteroidetes had only 85.4% identity to its top cephalotine-derived BLAST hit. Initially implicated as a plausible specialist, our phylogenetic analyses did place it near cephalotine-specific clades of bacteria in the Sphingobacteriales. But it was not part of a host-specific monophyletic unit, grouping instead with rumen-associated bacteria (Additional file 1: Fig. S9).
The phylogeny inferred for Analysis #8 grouped sequences from three Flavobacteriales OTUs (60, 72, and 119) into a monophyletic clade with 85% bootstrap support (Additional file 1: Fig. S10). Collectively found in four Cephalotes species (Additional file 1: Fig. S3), and with members showing only ~92.0% identity to their nearest non-cephalotine relatives (Additional file 2: Table S18), this lineage can be plausibly retained as a cephalotine-specialized clade. Sequences from the occasionally abundant Neisseriales in OTUs 019 and 069 were included in this same phylogeny. While only 94.5% identical to other, non-cephalotine derived bacteria (Additional file 2: Table S18), we did not see strong evidence for specialization due to non-monophyly of these two OTUs and confinement of each to a single host species (Additional file 1: Fig. S3). Also of interest were sequences from the Oxalobacteriaceae, including OUT041, found in three monophyletic Mexican Cephalotes, and OTU073, present in the early-branching C. unimaculatus (Dominican Republic). While related, monophyly of these sequences was interrupted by Olivibacterium flavum, a bacterium from a non-cephalotine source. We saw no additional evidence for host-specific lineages on this phylogeny.
Proportions of adult vs. larval microbiomes comprised by host-specialist bacteria
Using our phylogenetic results we expanded our list of cephalotine-specialized symbionts to include Lactobacillales (OTUs 007 and 022), Pseudomonadales (OTU046), and Flavobacteriales (OTUs 060, 072, 119) OTUs. With this finalized classification scheme in hand (Additional file 2: Table S19), we computed the proportion of each individual ants’ microbiome comprised by specialized bacteria, assessing these values across stages and castes from each of our 18 studied colonies (Additional file 1: Fig. S11). In adults, microbiomes of pre- and non-reproductive castes contained similarly high fractions of bacteria from cephalotine-specialized clades, i.e. those found exclusively in Cephalotes and/or Procryptocerus ants. Specifically, the median summed relative abundances for such specialized bacteria were, 0.984, 0.983 and 0.987 across across alate queens/gynes, soldiers, and workers, respectively. Averages for alate queens trended slightly downward compared to the other two castes (0.869 AQ, 0.940 S, and 0.956 W) due to the high prevalence of putatively non-specialized Wolbachia in the former caste for C. goniodontus.
In contrast to these patterns, symbionts from cephalotine-specialist clades comprised a median of 0.469 of mature/wingless queens’ microbiomes, and an average of 0.661. These lower fractions were driven by an abundance of Wolbachia or Neisseriales in a subset of these n=5 samples, and by the paucity or absence of specialized symbiont OTUs in the Xanthomonadales, Burkholderiales, Pseudomonadales (i.e. assigning to the genus Ventosimonas), Campylobacterales, Flavobacteriales, Gracillibacteria, and Sphingobacteriales. Along similar lines, Cephaloticoccus symbionts from the Opitutales – typically the most dominant microbes across adult ants – were abundant in only one mature queen. Indeed, in focusing on the non-Rhizobiales portion of wingless queens’ specialized core microbiomes, we saw that these, and the few other, remaining specialists made up a median of 0 and an average of 0.109 of their microbiomes. Proportions comprised by these non-Rhizobiales specialized bacteria were, contrastingly, high for alate queens/gynes, soldiers, and workers. In particular, these respective averages equaled 0.756, 0.717, and 0.719, illustrating their less Rhizobiales-dominated composition. It is interesting to note that most Rhizobiales found in mature/wingless queens hailed from only one of the two deep-branching cephalotine-associated lineages in this order (Fig. 7), with an absence of OTU001, found in all other castes and stages of every colony in our study (Additional file 1: Fig. S3).
Bearing some resemblance to wingless queens, the fractions of the larval microbiome made up by bacteria from cephalotine-specialist lineages was considerably lower than that for alate queens/gynes, soldiers, and workers, with a median and average of 0.700 and 0.675. Removal of Rhizobiales from consideration pushed these values down to 0.162 and 0.206. Additional removal of Lactobacillales OTUs 007 and 022 yielded median and average values of 0.023 and 0.060 for the specialized fraction of the larval microbiome, suggesting the rarity of other core adult-enriched specialist symbionts at this immature stage. Of notable exception were five colonies where the remaining specialized, adult-enriched symbionts comprised over 0.070 of the symbiotic community. In these cases we found a tendency toward modest abundance of specialized Burkholderiales OTUs from the Alcaligenaceae, which comprised between 0.061 to 0.279 of the total microbiome. Specialized Xanthomonales were generally rare, though they did constitute an average relative abundance of 0.065 in one colony from C. setulifer. Non-specialized portions of the larval microbiome were most often made up of Lactobacillales from beyond OTUs 007 and 022, Enterobacteriales, and sporadically abundant bacteria, including those from the orders Actinomycetales and Bifidobacteriales in the Actinobacteria (Additional file 1: Fig. S3). Bacteria from this phylum made up an average of 0.019 of the larval microbiome, reaching a maximum average of 0.121 in larvae of C. maculatus.
Assessing microbiome variation across larval development
With: 1) established definitions for specialist bacterial lineages (Additional file 2: Table S19); 2) identified patterns of high larval abundance of non-specialized symbionts from the Actinobacteria, Enterobacteriales, and Lactobacillales (Additional file 1: Fig. S3); 3) discovery of two divergent specialized Rhizobiales linages (Fig. 7a); and 4) suggestions that microbiomes of young vs. old larvae often classify into distinct enterotypes (Fig. 5c), we assessed trends of symbiont shifts across larval development by binning individuals into 11 normalized size categories, used to infer the relative age of 103 larvae. ANOVA statistics, focused on arcsin transformed relative abundances from each larval sequence library (Additional file 2: Tables S20-21) revealed significant differences in the titers of Rhizobiales OTU001 across larval development (df=10, sum of squares=0.42865, F-value=2.0602, p=0.03573), with a trend of increasing relative abundance with increasing larval size – used here as our proxy for age (Additional file 1: Fig. S12). Results were highly significant for the second cephalotine-specific Rhizobiales lineage containing OTUs 004, 017, 023, 042, and 056 (df=10, sum of squares=7.9636, F-value=5.5194, p-value=2.15 e-06); but in this case, relative abundance decreased with size. A combined pool of specialized and non-specialized Lactobacillales showed a trend of greater prevalence in larger and, hence, presumably older larvae, with similarly strong statistical support for a relationship between size and titer (df=10, sum of squares = 2.2033, F-value=3.6587, p-value=3.87 e-04). While non-specialized Enterobacteriales showed a similar trend of increasing relative abundance with size, larval size categories did not have a significant effect in our ANOVA model (df=10, sum of squares=0.7507, F-value=1.63, p-value=0.1103), suggesting a need for future study on this taxon. No other bacterial groups considered here showed significant differences in arcsin transformed relative abundance across these ranked size classes (Additional file 2: Table S21).
Cophylogeny analysis – Cephalotes species vs. Cephaloticoccus symbionts
Finally, to assess whether the earlier described phylosymbiosis signal might extend from cospeciation between turtle ants and their symbionts, we initiated what will be a deeper, more systematic investigation into host-symbiont phylogenetic congruence (Fig. 8; Additional file 1: Fig. S13). For this purpose, we used sequence data from a prior metagenomic study [7], inferring phylogenies of adult-enriched Cephaloticoccus (Opitutales) symbionts from separate analyses of two non-laterally transferred protein-coding genes. Bootstrap support vales for the nodes in each phylogeny were generally high, exceeding 80 for 10/16 internal nodes of the uvrB tree and for 9/16 internal nodes of the rpoB tree.
We compared these two phylogenies against the host Cephalotes phylogeny inferred by Price et al. [93], using Jane4. Both symbiont gene trees exhibited a greater degree of phylogenetic congruence in relation to the host phylogeny than expected by chance (Fig. 8; Additional file 1: Fig. S13). In particular, the event-based mapping of the Jane4 program estimated 10 cospeciation events compared to 7 host switches for the rpoB vs. host tree comparison. Additionally, the overall ‘cost’ of our inferred rpoB tree – 18 – was substantially lower than the all costs estimated for 1000 randomized datasets (average randomized cost = 27.82). Results were similar for the uvrB tree – with 9 cospeciation events compared to 7 host switching events, and an actual cost of 19 compared to an average of 27.74 for randomized datasets. Among the several nodes suggestive of cospeciation events, two pairs of sister symbiont taxa recapitulated sister taxa relationships on the Cephalotes phylogeny in ways that were not confounded by geography. This was true for both gene trees. Specifically, symbionts from C. pusillus (Brazil) and C. spinosus (Peru) grouped together with high bootstrap support (uvrB = 95; rpoB = 100), as did symbionts from C. minutus (Brazil) and C. simillimus (Peru) (uvrB = 80; rpoB = 100). In a distinct instance, the sister pairing of Peru-collected C. angustus and C. umbraculatus was nearly recapitulated on the rpoB gene tree, with the symbiont clade’s monophyly being broken up only by an apparent horizontal transfer from this symbiont lineage to the Brazil-derived C. clypeatus. Collectively, these patterns suggest that while the symbionts used for this analysis came from ants in a small number of geographic locales, host evolutionary history appears to be a geographically robust predictor of symbiont gene tree topology for at least one cephalotine-specific, adult-enriched symbiont.