Termites are eusocial insects that host a range of prokaryotic and eukaryotic gut symbionts and can differentiate into a range of caste phenotypes. Soldier caste differentiation from termite workers follows two successive molts (worker–presoldier–soldier) that are driven at the endocrine level by juvenile hormone (JH). Although physiological and eusocial mechanisms tied to JH signaling have been studied, the role of gut symbionts in the caste differentiation process is poorly understood. Here, we used the JH analog methoprene in combination with the antibiotic kanamycin to manipulate caste differentiation and gut bacterial loads in Reticulitermes flavipes termites via four bioassay treatments: kanamycin, methoprene, kanamycin+methoprene, and an untreated (negative) control. Bioassay results demonstrated a significantly higher number of presoldiers in the methoprene treatment, highest mortality in kanamycin+methoprene treatment, and significantly reduced protist numbers in all treatments except the untreated control. Bacterial 16S rRNA gene sequencing provided alpha and beta diversity results that mirrored bioassay findings. From ANCOM analysis, we found that several bacterial genera were differentially abundant among treatments. Finally, follow-up experiments showed that if methoprene and kanamycin or untreated termites are placed together, zero or rescued presoldier initiation, respectively, occurs. These findings reveal that endogenous JH selects for symbiont compositions required to successfully complete presoldier differentiation. However, if the gut is voided before the influx of JH, it cannot select for the necessary symbionts that are crucial for molting. Based on these results, we are able to provide a novel example of linkages between gut microbial communities and host phenotypic plasticity.

Termites are eusocial insects living in colonies composed of three main caste phenotypes (workers, soldiers and reproductives), each with their own well-defined socio-behavioral roles (Lainé and Wright, 2003; Roisin, 2000; Wilson, 1971). Termite workers further have the unique characteristic capacity to molt into other caste phenotypes through the process known as caste differentiation, which is a form of phenotypic plasticity. Caste differentiation and development in Reticulitermes flavipes can follow either apterous/neuter or imaginal routes (Scharf et al., 2003) and can have three different developmental potentials that include status quo molts into next-instar workers, differentiation into a presoldier that later molts into a soldier (Zhou et al., 2006a), or differentiation into a supplementary reproductive individual. Research on the worker–soldier molt process is of recent interest as a non-destructive form (soldier) develops from the most economically destructive life stage (worker). Several endogenous and exogenous factors have been reported to regulate the worker–soldier molt. At the endocrine level, the sesquiterpenoid juvenile hormone (JH) acts as a hormone to differentiate workers into soldiers via an intermediate presoldier phenotype (Tarver et al., 2009; Zhou et al., 2007). In R. flavipes and most other insects, JH III is produced by a neuroendocrine gland known as the corpora allata (CA) closely associated with the brain (Yagi et al., 2005). The entire soldier differentiation process occurs in response to the elevated endogenous JH titers (Chan et al., 2011; Elliott et al., 2009; Elliott and Stay, 2008; Mao et al., 2005; Park and Raina, 2004). However, in other insects, elevated JH retains the immature phenotype through development resulting in extra-nymphal or larval instars (Nijhout et al., 2014; Smykal et al., 2014). In addition to regulating termite polyphenism, JH is also involved in vitellogenesis through inducing increased vitellogenin protein, insect diapause and other processes (Engelmann, 1984; Scharf et al., 2005; Yagi, 1976).

Besides JH, termite caste differentiation could be associated with their gut microorganisms. The termite gut is a habitat for numerous microorganisms from various taxa, including protists, Bacteria and Archaea, living in a consortium and contributing to termite physiology (Brune and Ohkuma, 2010; Ohkuma, 2008; Scharf and Tartar, 2008). Termite symbionts are crucial in cellulose/hemicellulose degradation, nitrogen fixation, acetogenesis and anti-fungal defense (Doolittle et al., 2008; Inoue et al., 2000, 1997; Peterson and Scharf, 2016a; Peterson and Scharf, 2016b). The literature is flooded with reports on the impact of gut symbionts on endocrine regulation in humans, mice and insects. The human gut bacterium Lactobacillus reuteri was reported to upregulate the neuropeptide hormone oxytocin to enhance wound-healing properties through nerve-mediated pathways (Poutahidis et al., 2013). Studies on human and mouse models have supported the conclusion that gut symbionts can affect immune and endocrine modulation in a variety of signaling pathways (e.g. Chen et al., 2015). Similarly, in some insects, gut symbionts were observed to influence the pheromones released, social behavior, host development, egg production, expression of hemolymph proteins, and hormones produced (Engl and Kaltenpoth, 2018; Lee et al., 2017). In honeybees, gut symbionts were also reported to increase body mass by mediating changes in insulin and vitellogenin signaling pathways (Zheng et al., 2017).

Several aspects of bacteria–termite symbiosis are yet to be understood. No reports are yet available on the interaction of termite gut symbionts with presoldier differentiation, which is regulated by JH produced by the endocrine system. In female bean bugs (Riptortus clavatus), the gut symbiont Burkholderia was reported to modulate the JH titer (Lee et al., 2019). Association with commensal bacteria such as Lactobacillus plantarum in Drosophila melanogaster has also been correlated with hormone signaling, providing strong evidence for the influence of gut bacteria on insect endocrine functioning (Storelli et al., 2011), while in termites, JH was reported to alter the expression of hundreds of symbiont genes in the gut (Sen et al., 2013). Another study has demonstrated that presoldier phenotypes have biased expression of genes coding for vitellogenin, a protein thought to bind JH (Scharf et al., 2003), while in other insects, gut symbionts were reported to influence the expression of vitellogenin (Zheng et al., 2017). Several studies on other organisms have further demonstrated that gut symbionts can impact hormone signaling pathways.

The goal of this study was to investigate whether termite gut bacteria can influence JH-dependent presoldier differentiation in R. flavipes, with the overarching hypothesis being that specific gut symbionts are connected to the presoldier differentiation process. Our specific objectives were to (i) study the effects of JH and antibiotic treatments on presoldier emergence and/or survivorship and protist numbers, and (ii) seek to connect specific gut bacteria with the above conditions using 16S rRNA gene sequence surveys of the gut microenvironment. Our findings provide a novel example of the influences of gut microbial symbionts on host phenotypic plasticity.

Termite bioassays

This study is based on a preliminary experiment where JH III and kanamycin were used to induce presoldier differentiation and results were similar to those obtained in the present study. We opted to use methoprene in place of JH III because it is consistently available with high purity, whereas JH III is not. Termites [Reticulitermes flavipes (Kollar 1837)] from colonies reared in the laboratory were used in bioassays. There were two phases for this experiment. In the first phase, paper towel sandwiches were pretreated with either kanamycin (50% w/v in water) or water alone, and were then fed to replicate groups of termites for 2 days. For the second phase, termites were then moved to paper towel sandwiches that were treated with either acetone or methoprene (30 µg per dish) in acetone (Chem Service Inc., West Chester, PA, USA) and dried in a fume hood prior to use. Methoprene and kanamycin were used to induce soldier caste differentiation and manipulate gut microbial load, respectively. Based on previous studies, kanamycin significantly reduced and altered the gut microbial load with minimal impacts to the host (Peterson et al., 2015), while exogenous application of methoprene has been reported to induce the genes in JH signaling pathways and similarly induce presoldier differentiation (Du et al., 2020; Tarver et al., 2012). Also, preliminary experiments with JH III demonstrated similar presoldier differentiation as seen in the present study with methoprene. Phase 1 pre-treatments lasted 2 days, which was identical for all the experiments; however, the duration of the final treatment was dependent on the type of experiment, as detailed in the following sections. Treatment abbreviations are as follows: Kan (kanamycin alone), Meth (methoprene alone), KanMeth (kanamycin in phase 1 and methoprene in phase 2) and UnCtl (untreated control), which received water in phase 1 and acetone (evaporated) in phase 2.

Presoldier initiation and survivorship assessment bioassays

The assessment of presoldier emergence and survivorship was done to observe the differences among treatments at the whole-insect ‘holobiont’ level. There were four biological and three technical replicates each for this assessment (12 replicates total). For this bioassay, 25 worker termites were provided with kanamycin- or water-treated paper towel sandwiches in Petri dishes for 2 days (phase 1). After 2 days, 20 termites were transferred for final treatments as detailed above (phase 2) and observed for 15 days. Termites emerging from molting with transparent elongated mandibles were scored as presoldiers. Termites that were immobile and could not move their appendages when disturbed by soft forceps were considered dead. Survivorship in untreated controls, even in smaller groups, is quite high and invariable. We have tried working with larger groups in the past, and treatment-induced mortality is still a problem, coupled with an inability to score assays quickly and without causing additional, confounding assay stress. Therefore, we used small groups of workers for this experiment.

Protist counts

Each bioassay treatment described above was replicated 6 times and 35 termites were placed in each Petri dish for pre-treatment (phase 1). For final treatments (phase 2), 30 termites were transferred to Petri dishes with methoprene- or acetone-treated paper towel sandwiches for 9 days, which was the day of emergence of at least one presoldier in the methoprene treatment. The hind guts from five live termites each from assay days 0, 2, 6 and 11 were dissected. The guts were squeezed gently into 500 μl phosphate buffered saline (PBS, pH 7.4). A 10 μl suspension was used to count the protist numbers using a phase contrast microscope at 20× magnification and a Sedgwick Rafter counting cell (SPI Supplies, West Chester, PA, USA).

16S rRNA gene sequencing

DNA extraction

Bioassays with 25 worker termites were set up for DNA extraction. There were four replications each for the treatments. Similar to the phenological study, 20 termites were taken for the second phase of the experiment. Live termite workers were sampled and dissected after at least one presoldier was sighted (on day 11) in methoprene bioassay treatments. We picked day 11 for sampling DNA because the survival of termites was high and only live termite workers were dissected to acquire sufficient high-quality DNA from gut tissues. Also, a histological study of fat body (Cornette et al., 2007) revealed that the production of fat body and protein granules starts within the first 3 days after hormonal treatment. These granules disappear within a few days of presoldier development, indicating that termite physiological changes can be documented up to the point of presoldier molting. Four replications of 10 guts per treatment were sampled for DNA extraction. The 10 guts from live termite workers were dissected, placed in 200 μl PBS and processed for DNA extraction using the DNeasy® Blood and Tissue Kit (Qiagen, Valencia, CA, USA). Along with the experimental samples, water was also run through the DNA extraction protocol as an extraction control. The quality of DNA in the aliquots was assessed through gel electrophoresis and was quantified using a Nanodrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and later verified using fluorometric detection by the University of Minnesota Genomics Core. Thirty microliters (41.3–125.9 ng µl−1) of extracted DNA was sent for library preparation and sequencing, and the remaining DNA was used for 16S rRNA gene quantitative PCR (qPCR) analysis.

16S quantitative PCR
The 16S rRNA gene copy number from different treatments detailed above was assessed using two independent primer sets that included (1) primers 338F and 518R and (2) primers 515F and 806R, flanking the V3 and V4 regions of 16S rRNA gene, respectively. For primer set 1, the total reaction volume was 15 μl, containing 7.5 μl buffer (SsoFast EvaGreen® Supermix, Bio-Rad Inc., Hercules, CA, USA), 1 μl reverse and forward primers each (10 µmol l−1 final concentration), 4 μl nuclease-free water and 1.5 μl DNA template (0.31–1.48 ng). The cycling conditions used were initial denaturation at 95°C for 3 min followed by 40 cycles of denaturation at 95°C for 15 s, annealing at 55°C for 15 s, and extension at 72°C for 30 s followed by a melt curve from 65°C to 95°C with 0.5°C increments. The 16S rRNA qPCR with primer set 2 was performed by the University of Minnesota Genomics Core sequencing facility; the cycling conditions used were initial denaturation at 95°C for 5 min followed by 35 cycles of denaturation at 98°C for 20 s, annealing at 55°C for 15 s and extension at 72°C for 60 s. Final extension was at 72°C for 5 min, after which the samples were held at 4°C. DNA copy number (molecules) was determined by comparison against a DNA standard curve (see equation below), with analysis and PCR efficiency determinations by the method 10(–1/slope):
formula
(1)
where X is the amount of DNA (in ng) calculated using the equation of a standard curve, A is Avogadro's number (6.022×1023 molecules mole−1), N is the primer length, M is the average mass of one bp of dsDNA (660 g mole−1) and C is the conversion factor (109 ng g−1).
16S rRNA gene sequencing

The V3–V4 region of 16S rRNA gene was sequenced using bacterial universal primers 341F-CCTACGGGAGGCAGCAG and 806R-GGACTACHVGGGTWTCTAAT (Caporaso et al., 2011; Peterson et al., 2015). A two-step PCR was used for library preparation. The indexed PCR products were normalized to 10 ng µl−1, and 10 µl of each sample was pooled and concentrated to ∼100 µl and cleaned using solid reverse phase immobilization (SPRI) purification. The cleaned library was quantified by Qubit and fragment analyzed by Agilent Tapestation, and the library was diluted to 2 nmol l−1. The pooled samples were denatured with NaOH, diluted to 8 pmol l−1 in Illumina HT1 buffer, spiked with 15% PhiX and heat denatured at 96°C for 2 min immediately before loading. A MiSeq 600 (1/8th lane Stowaway) cycle v3 kit was used to sequence the sample by the University of Minnesota Genomics Center, MN, USA.

Sequence processing and community analysis

The paired end (301 bp) Illumina reads were processed and analyzed using the QIIME2 pipeline (Bolyen et al., 2019). Divisive Amplicon Denoising Algorithm 2 (DADA2) was used to obtain high-quality reads following the filtering, trimming, denoising, dereplicating, merging and removing chimera steps (Callahan et al., 2016). During this process, a total of 1,559,624 demultiplexed sequences were processed to obtain 1,234,554 non-chimeric high-quality sequences in QIIME2. The moving picture tutorial in QIIME2 was followed to calculate the diversity matrices. To calculate the diversity matrices (alpha and beta), the sequences were subsampled from the feature table from phylogenetic diversity analysis at the sampling depth of 2000 sequences that removed one replicate from KanMeth and four replicates of extraction control, respectively. The alpha diversity measures used were observed operational taxonomic units (OTUs), Pielou’s evenness and the Shannon index. The beta diversities, phylogenetic Unifrac distances [weighted and unweighted (Lozupone et al., 2011)] and non-phylogenetic distances [Jaccard and Bray–Curtis] were calculated and visualized using R and QIIME2. The reference database SILVA (version 132) was customized following the instructions on training the classifier for the V3/V4 region of 16S rRNA gene. Taxonomy was assigned to the data using naïve Bayes trained Silva 13,299% OTU classifier bounded by the 341F/806R primers set as the reference reads (Quast et al., 2013). Bar graphs representing the relative abundance of amplicon sequence variant (ASV) at the phylum and family levels were generated using the R package phyloseq from QIIME2-generated phylotype tables and taxonomy files (McMurdie and Holmes, 2013). Finally, analysis of composition of microbiomes (ANCOM) was used to study differential microbiota composition at phylum and genus level in different treatments via the QIIME2 plugin (Mandal et al., 2015).

Trophallaxis experiment

A final trophallaxis experiment was conducted to understand the effect of different treatments in a eusocial context on donor and recipient termites. Prior to the experiment, ‘donor’ termites were treated with Nile Blue dye 0.5% w/v (Sigma-Aldrich, St Louis, MO, USA) for 2 days. After 2 days, the donor termites were transferred into Petri dishes containing filter paper sandwiches treated with methoprene or kanamycin, with an exposure duration of 2 days. The methoprene and kanamycin treatments and the untreated controls were set on the same day. Similarly, some untreated termites were exposed to kanamycin in different Petri dishes; these termites later served as ‘recipients’ for methoprene-treated donors. The treated termite donors were then placed in new Petri dishes with untreated or kanamycin-treated recipients in the ratio of 1:1 with a total of 20 individuals. The donor–recipient combinations during the experiment were: kanamycin-treated donors+untreated recipients (Kan-D+Untreated-R), methoprene-treated donors+kanamycin-treated recipients (Meth-D+Kan-R) and methoprene-treated donors+untreated recipients (Meth-D+Untreated-R). Controls for this experiment included kanamycin (Kan control), methoprene (Meth control) and untreated (Untreated control). Owing to the total mortality of termites in some treatments, the experiment was ended after 11 days. The data on termite mortality and presoldier initiation were taken every day until 11 days, with the cumulative data from day 11 of the treatment being presented. Each treatment combination was replicated three times.

Statistical analyses

Statistical analyses were performed in the QIIME2 interface and R (Bolyen et al., 2019; https://www.r-project.org/). The Kruskal–Wallis test was used to compare presoldier emergence, total mortality, average number of protists and alpha diversities among the treatment groups and trophallaxis experiment. Pairwise comparisons were made using the Mann–Whitney U-test and were corrected for Benjamini–Hochberg false discovery rate. Beta diversity matrices were subjected to permutational multivariant analysis of variance (PERMANOVA) and the PERMDISP test of homogeneity of dispersions to study statistical significance among groups. The graphs for 16S rRNA gene sequencing data were generated and visualized using R and Microsoft Excel. Canonical correspondence analysis (CCA) (ter Braak, 1986) was conducted to study the relationship between the bacterial genera and other variables – protist number, presoldier initiation and total mortality. CCA was implemented with the R package ‘vegan’. For the CCA analysis, four replicates from the presoldier initiation, total mortality and protist numbers were used to study their association with the taxa abundance.

Presoldier differentiation, survivorship and protist numbers significantly differed among treatments

The effect of treatments on cumulative percentage of presoldiers initiated (H=37.625, d.f.=3, P<0.001), cumulative percentage mortality (H=27.027, d.f.=3, P<0.001) on day 17 and average protist numbers per gut (H=19.048, d.f.=3, P<0.001) on day 11 were significant based on Kruskal–Wallis tests. Based on pairwise Mann–Whitney comparisons (P<0.05), the total percentage of presoldiers emerged was significantly higher in the Meth than in the KanMeth treatment, while no presoldiers emerged in the UnCtl and Kan treatments (Fig. 1A). From the same bioassays, survivorship of workers+presoldiers was also recorded. Total mortality in the KanMeth treatment was significantly higher than in any other treatments, and as expected, the UnCtl treatment had the lowest mortality. Mortality levels in the individual Kan and Meth treatments were not significantly different (Fig. 1B). Some variability in the mortality of termites from different replicates in the Kan, KanMeth and Meth treatments was observed, but this variability was not sufficient to overcome the statistically significant trends as reported. The variability could be due to the variation in group response to the treatments. Similar group variation in the mortality of termites against JH analog was observed by Hrdý et al. (2006) while treating different termite colonies with the same compound. The average number of protists per termite gut on day 11 was significantly higher in the UnCtl treatment in comparison with other treatments, while the Kan and KanMeth treatments had similar protist numbers (Fig. 1C). The average number of protists in the Meth treatment was significantly higher than in the Kan and KanMeth treatments, but still was greatly reduced in comparison with the UnCtl treatment. The average mortality of termites in all treatments on the day of protist and DNA sampling (day 11) was less than 50% (Fig. S1).

Fig. 1.

Bioassay and protist count results. (A) Cumulative percentage presoldiers initiated (12 replications, sample size=20), (B) cumulative percentage mortality (12 replications, sample size=20) and (C) average number of protists per termite gut (6 replications, sample size=5). Boxplots show the median (horizontal line in the box), interquartile range (IQR, the box) and samples within 1.5 times IQR (vertical lines); samples outside 1.5 times IQR are shown as small circles. P-values for global Kruskal–Wallis comparisons across treatments are shown at the top of each graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in graphs with the same letters are not significantly different (Mann–Whitney U-test; P<0.05).

Fig. 1.

Bioassay and protist count results. (A) Cumulative percentage presoldiers initiated (12 replications, sample size=20), (B) cumulative percentage mortality (12 replications, sample size=20) and (C) average number of protists per termite gut (6 replications, sample size=5). Boxplots show the median (horizontal line in the box), interquartile range (IQR, the box) and samples within 1.5 times IQR (vertical lines); samples outside 1.5 times IQR are shown as small circles. P-values for global Kruskal–Wallis comparisons across treatments are shown at the top of each graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in graphs with the same letters are not significantly different (Mann–Whitney U-test; P<0.05).

Confirmation of antibiotic impacts on gut bacterial populations

The two independent primer sets, 338F–518R and 515F–806R, used for qRT-PCR showed highly similar results, indicating that the Kan treatment had an expected result of reducing bacterial loads (∼2.5× decrease compared with UnCtl), but also unexpectedly that the Meth (∼1.2× decrease compared with UnCtl) and KanMeth treatments (∼2× decrease compared with UnCtl) had reduced bacterial loads as well.

Bacterial alpha diversity significantly differed among treatments

Rarified 2000 sequence reads obtained from 16S rRNA gene sequencing were used to calculate alpha and beta diversity measures. Alpha diversity was estimated using observed OTUs for species richness (Fig. 2A), Pielou’s evenness (Fig. 2B) and the Shannon diversity index (Fig. 2C). The effect of different treatments on observed OTUs (H=12.829, d.f.=3, P=0.005), Pielou’s evenness (H=12.567, d.f.=3, P=0.005) and the Shannon index (H=12.567, d.f.=3, P=0.005) were significant based on Kruskal–Wallis tests. Based on pairwise Mann–Whitney comparisons (P<0.05), the UnCtl treatment demonstrated significantly higher species richness and evenness, followed by the Meth treatment, whereas the Kan and KanMeth treatments were not significantly different from each other. Similar results were obtained for average numbers of protists estimated per gut (Fig. 1C), suggesting a close connection between some bacteria and protist symbiont communities.

Fig. 2.

Alpha diversity matrices from 16S surveys. (A) Observed operational taxonomic unit (OTU) numbers, (B) Pielou’s evenness and (C) the Shannon diversity index. Boxplots show the median (horizontal line in the box), IQR (the box) and samples within 1.5 times IQR (vertical lines). Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=10). P-values for global Kruskal–Wallis comparisons across treatments are shown at the top of each graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in the graphs with the same letters are not significantly different (P<0.05).

Fig. 2.

Alpha diversity matrices from 16S surveys. (A) Observed operational taxonomic unit (OTU) numbers, (B) Pielou’s evenness and (C) the Shannon diversity index. Boxplots show the median (horizontal line in the box), IQR (the box) and samples within 1.5 times IQR (vertical lines). Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=10). P-values for global Kruskal–Wallis comparisons across treatments are shown at the top of each graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in the graphs with the same letters are not significantly different (P<0.05).

Bacterial beta diversity significantly differed among treatments

PERMANOVA statistics show that the weighted UniFrac analysis was statistically significant (pseudo-F=36.101, P=0.001), which means the effect of treatments was significant in creating the dissimilarities among the different treatment groups (Fig. S2). Also based on the PERMDISP test, the bacterial community structures among the replicates were not determined to be significantly different (PERMDISP test, F=0.7168, P=0.224). Pairwise PERMANOVA results further indicate all treatments differed in terms of community structure for weighted Unifrac. Similar results were observed for beta diversity matrices (Bray–Curtis and Jaccard). However, for the unweighted Unifrac analysis, the Kan and KanMeth treatments were not significantly different, whereas other pairwise comparisons were significant.

Bacterial taxonomic composition differed among treatments

The taxonomic composition and relative abundance of major bacterial phyla at abundances greater than 0.1% varied greatly among treatments but was consistent among the replicates within treatments (Fig. 3). At the phylum level, in the UnCtl treatment, Spirochaetes, Firmicutes and Bacteroidetes were the most abundant (∼70%). In the Meth treatment, the most abundant phyla were Proteobacteria and Bacteroidetes (∼55%). In the KanMeth treatment, the most abundant phyla were Proteobacteria and Actinobacteria (∼85%). Finally, the Kan treatment was dominated by the phylum Actinobacteria (∼75%). Similar to phylum level results, family level relative abundances also varied widely among the treatments.

Fig. 3.

Relative abundance (from abundance >0.001), measured as the percentage of the total bacterial abundance for each treatment. Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=10). The index indicates the phylum level identification of bacterial species.

Fig. 3.

Relative abundance (from abundance >0.001), measured as the percentage of the total bacterial abundance for each treatment. Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=10). The index indicates the phylum level identification of bacterial species.

Based on the ANCOM test at the phylum level, the bacterial phylum Elusimicrobia was the most significantly abundant among treatments. In contrast, at the genus level, the following top 13 bacteria were significantly differentially abundant among treatments (centered log ratio, clr>200): uncultured delta proteobacterium (Proteobacteria), uncultured Synergistetes (Synergistetes), Lactococcus (Firmicutes), Propionivibrio (Proteobacteria), uncultured Epsilonproteobacteria bacterium (Epsilonbacteraeota), Tyzzerella 3 (Firmicutes), Treponema (Spirochaetes), Treponema 2 (Spirochaetes), Candidatus Endomicrobium (Elusimicrobia), unidentified Firmicutes (Firmicutes), Stenoxybacter (Proteobacteria), uncultured Bacteroidetes bacterium (Bacteroidetes) and Cerasicoccus (Verrucomicrobia). The observed log abundances tied to different genera are shown in Fig. S3, and the list of significant genera based on ANCOM hypothesis test with their respective clr and W scores is shown in Table S1.

CCA confirms significant factors contributing to the differences among the treatments

The CCA model indicated that the variation in bacterial species abundance across treatments was significantly correlated with presoldier initiation, total mortality and protist numbers (chi-square, F3,11=1.7603, P=0.006). The CCA plot shows that the number of protists and presoldiers initiated was strongly associated with the UnCtl and Meth treatments, respectively, whereas total mortality was strongly associated with Kan and KanMeth treatments (Fig. 4). Axes 1, 2 and 3 explain 38.6%, 12.03%, and 7.8% of the total constrained variation, respectively.

Fig. 4.

Canonical correspondence analysis model performed with relative abundance of bacterial taxa and environmental variables – protist numbers, presoldier initiated and total mortality. Overall significance of the model: chi-square, F3,11=1.7603, P(>F)=0.006, with first, second and third axes representing 38.6%, 12.03% and 7.8% of the total constrained variation, respectively. Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=N/A).

Fig. 4.

Canonical correspondence analysis model performed with relative abundance of bacterial taxa and environmental variables – protist numbers, presoldier initiated and total mortality. Overall significance of the model: chi-square, F3,11=1.7603, P(>F)=0.006, with first, second and third axes representing 38.6%, 12.03% and 7.8% of the total constrained variation, respectively. Total replications for Kan, Meth and UnCtl were 4 while KanMeth had 3 replications (sample size=N/A).

Significant effects of trophallaxis on presoldier initiation and termite survivorship

The effect of trophallaxis ‘rescue’ treatments on cumulative percentage of presoldiers initiated (H=10.692, d.f.=3, P=0.013) and cumulative percentage mortality (H=16.111, d.f.=5, P=0.006) were significant based on Kruskal–Wallis tests. No presoldiers emerged in the Meth-D+Kan-R and untreated control treatments. All the other treatment combinations that did not result in presoldier emergence were removed from the presoldier initiation data analysis. Based on pairwise Mann–Whitney comparisons (P<0.05), the total percentage of presoldiers emerged in meth controls was significantly higher than the treatment combination Meth-D+Untreated-R termites (Fig. 5). From the same bioassays, survivorship of termites was also recorded. Total mortality in Kan control and Meth-D+Kan-R were significantly higher than in other treatment combinations. The mortalities in Kan-D+Untreated-R and Meth-D+Untreated-R were not statistically different from one another (Fig. S4). Collectively, these results show that presoldier emergence is higher in methoprene-treated individuals held with control individuals containing intact microbiomes than in methoprene-treated individuals held with kanamycin-treated individuals.

Fig. 5.

Trophallaxis experiment on presoldier initiation study showing cumulative percentage presoldiers initiated over the course of 11-day assays (3 replications, sample size=20). Boxplots show the median (horizontal line in the box), IQR (the box) and samples within 1.5 times IQR (vertical lines). P-values for global Kruskal–Wallis comparison across treatments are shown at the top of the graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in the graph with the same letters are not significantly different (Mann–Whitney U-test; P<0.05). Treatment abbreviations are defined in the Materials and Methods.

Fig. 5.

Trophallaxis experiment on presoldier initiation study showing cumulative percentage presoldiers initiated over the course of 11-day assays (3 replications, sample size=20). Boxplots show the median (horizontal line in the box), IQR (the box) and samples within 1.5 times IQR (vertical lines). P-values for global Kruskal–Wallis comparison across treatments are shown at the top of the graph. Significance levels for pairwise comparisons are indicated by lowercase letters. The individual treatment groups in the graph with the same letters are not significantly different (Mann–Whitney U-test; P<0.05). Treatment abbreviations are defined in the Materials and Methods.

This study experimentally quantified the effect of gut symbiont communities on the phenomenon of JH-mediated presoldier differentiation in termites. The study was conducted using methoprene, whose chemical structure closely approximates that of JH (Charles et al., 2011). Methoprene is also consistently available in a highly pure form, while JH III is not, which makes independent verification of findings for JH III difficult. These findings provide a novel example of gut microbial community influence on host phenotype. Four different treatments were used to investigate possible roles of bacterial and protist gut symbionts in the caste differentiation process. These treatments included kanamycin alone (Kan), kanamycin+the JH analog methoprene (KanMeth), methoprene alone (Meth) and an untreated control (UnCtl). The effect of these treatments on presoldier formation, survivorship and protist abundance were observed through phenotypic bioassays and protist quantification. Treatment effects on bacterial abundance and composition were studied through the use of quantitative PCR and16S rRNA gene sequence analysis, which together were used to infer alpha and beta diversities and differential abundance through data analysis procedures.

The JH analog methoprene demonstrates a similar effect as endogenous JH would on a termite worker, resulting in significant presoldier initiation (Zhou et al., 2006b). Kanamycin was previously shown to reduce both protist and bacterial symbiont load in the termite gut (Peterson et al., 2015), and a similar effect was observed in this experiment as well. Therefore, the higher percentage of presoldiers emerged in the Meth compared with the KanMeth treatment suggests that gut symbionts along with JH play a crucial role in presoldier initiation and their survival. In the survivorship assessment among different treatments, the higher percent mortality of termites in the KanMeth treatment could be linked to the double stress from antibiotic and methoprene/JH. That means exposure to kanamycin prior to methoprene reduced the gut microbial load (possibly the gut microbes needed to survive stress during the molting process), resulting in the significantly high mortality observed in phenological bioassays. Additionally, the mortality observed with the Meth-D+Kan-R treatment in our rescue experiment further supports that presoldier initiation within colonies is lethal to termites lacking intact gut microbiomes. The trophallaxis experiment further support that the microbiome is relevant for surviving the presoldier molt in termites.

The mortality observed in the Meth treatment was also higher than the mortality observed with JH III by Hrdý et al. (2006). The possible explanation of this higher mortality could be due to the strong affinity of methoprene at the JH receptor or the stress caused during and after presoldier molting that is normally alleviated by the presence of specific gut symbionts and trophallaxis exchange with other colony members. Another cause of bioassay mortality could be the excessive proportion of presoldiers initiated in our study, which are unable to eat or be fed by workers exposed to methoprene in Petri dishes. Lastly, similar mortality percentages were observed by Howard and Haverty (1978) while using a similar concentration of methoprene as ours. Protist numbers were significantly reduced in termites prior to molting in Meth treatments (Nalepa, 2017; Sen et al., 2013), while some protists were still intact. Termite gut undergoes defaunation during the molting process causing a reduction in the numbers of gut symbionts. This defaunation or reduced protist number could be another cause for the significant mortality we observed in Meth treated individuals. Similar to previous findings (Peterson et al., 2015), protist numbers were reduced in the Kan treatment, which corresponds well with the observed reduction of several bacterial taxa (Spirochaetes, Elusimicrobia, Firmicutes, Bacteroidetes) that are known to maintain symbiotic associations with protists.

Alpha diversity matrices from the 16S rRNA gene sequencing demonstrated that the UnCtl treatment had the highest species evenness and richness, which was expected. Other treatment groups, however, had reduced alpha diversity in ways that were not entirely expected. The reductions could be due to ‘SPRY’ genes, which have been documented in association with reduced gut symbiont numbers in termites when challenged with hormones and antibiotics (Scharf et al., 2017). The significantly higher species evenness and richness in the Meth compared with the Kan and KanMeth treatments strengthens our above conclusion, i.e. that certain gut bacteria are required for successfully completing the presoldier transition, which is well supported by both the presoldier differentiation and trophallaxis experiments. Even though the KanMeth treatment resulted in some presoldier emergence, the total presoldier emergence was reduced significantly in comparison with the Meth treatment. A likely explanation for this difference is that gut symbionts were selected at two levels by kanamycin and JH in the KanMeth treatment. Bacterial diversity was further investigated among treatments using beta diversity matrices. Weighted UniFrac analysis demonstrated that different treatment groups had microbial communities different from each other. Similar results were observed for Bray–Curtis and Jaccard analyses, indicating that each treatment resulted in significantly different bacterial community structures (Lozupone et al., 2007). Overall, from the alpha and beta diversity matrices, and supporting statistical analyses, it is clear that different hormonal and antimicrobial treatments selected for different gut symbiotic load and content. This determination provided rationale to look further into the types of bacteria selected by the different treatments as detailed below.

Several protists, Archaea and Bacteria inhabit the small micro-environment inside a termite gut and provide synergistic physiological benefits to the host (Ohkuma, 2008; Scharf and Peterson, 2021). In this experiment, the major bacterial phyla observed were Spirochaetes, Firmicutes and Bacteroidetes, comprising more than 70% of abundance similar to the previous findings from 16S V5–V6 sampling of the R. flavipes gut lumen (Boucias et al., 2013). Shifts in both bacterial composition and abundance were observed among different treatments, following a hologenome theory of evolution, i.e. the artificially introduced stress was able to change the gut microbial community (Zilber-Rosenberg and Rosenberg, 2008). Most notably, Kan and KanMeth treatments lacked Spirochaetes and Elusimicrobia and had a reduction in Firmicutes. In wood-feeding termites, Spirochaetes are associated with acetogenesis and cellulose-degrading protists, whose absence could be detrimental to the survival of termite. Reduction in xylose released from lignocellulosic diet materials was also observed when the abundances of Elusimicrobia and Firmicutes were reduced (Leadbetter et al., 1999; Lucey and Leadbetter, 2014; Otani et al., 2014; Peterson et al., 2015). This evidence from prior studies highlights the involvement of these bacterial phyla in lignocellulose digestion and ultimately in the survival of termites, which is in agreement with the current findings. Spirochaete bacteria are ectosymbionts to flagellate protists (Iida et al., 2000), whereas Bacteroidetes, Elusimicrobia and Firmicutes are endosymbionts to protists (Fröhlich and König, 1999; Noda et al., 2007, 2005; Ohkuma et al., 2007; Stingl et al., 2005). Because the Kanamycin treatment eliminated or reduced the endosymbionts and ectosymbionts, the disruption of the bacteria–protist symbiosis likely resulted in significant reductions in protist counts and, ultimately, decreased termite survivorship (Peterson et al., 2015; Stephens and Gage, 2020; Treitli et al., 2019).

Based on ANCOM analysis, the genus Cerasicoccus (Verrucomicrobia) differed significantly by being completely absent from the Meth, KanMeth and Kan treatments. The taxa Propionivibrio (Proteobacteria), uncultured Epsilonproteobacteria bacterium (Epsilonbacteraeota), Tyzzerella 3 (Firmicutes), Treponema (Spirochaetes), Treponema 2 (Spirochaetes), Candidatus Endomicrobium (Elusimicrobia) and unidentified Firmicutes (Firmicutes) demonstrated the highest abundance in the UnCtl treatment, while these taxa were completely absent or highly reduced in the Kan/KanMeth treatments and reduced in the Meth treatment, suggesting their importance in the survival and successful completion of the presoldier differentiation process. In contrast, the taxa uncultured delta proteobacterium (Proteobacteria), uncultured Synergistetes (Synergistetes) and Lactococcus (Firmicutes) were equally abundant in the Meth and UnCtl treatments but were completely absent in the Kan and KanMeth treatments. And although the taxa Stenoxybacter (Proteobacteria) and uncultured Bacteroidetes bacterium (Bacteroidetes) were slightly reduced in the UnCtl compared with the Meth treatment, they were completely absent in the Kan and KanMeth treatments, suggesting their importance in the successful completion of the presoldier differentiation process. The significantly high mortality observed in the Kan and KanMeth treatments can be associated with gut stress owing to the disappearance or reduction of these taxa. Another possible explanation based on the trophallaxis experiment could be that when termites were treated with kanamycin and methoprene, they were not able to provide the symbionts to each other needed for survival, resulting in high mortality; this possibility also explains the lower or lack of presoldier differentiation in these treatments and provides further eusocial context to our findings.

The CCA analysis also confirms the strong association of total mortality with the Kan and KanMeth treatments. The reduction of specific protist flagellate-associated bacterial groups (Elusimicrobia, Firmicutes, Spirochaetes) with the Meth treatment is in agreement with the fact that termite guts are voided prior to molting (Nalepa, 2017; Sen et al., 2013) and explains the mortality that might be present in the colony during the differentiation, which could be one of the reasons why there are typically under 3% soldiers in R. flavipes colonies. Based on the presoldier emergence results from the trophallaxis experiment, it is also possible that normally faunated termites have inhibitory effects on presoldier differentiation, which may also contribute to the low numbers of soldiers present in R. flavipes colonies. In humans, bacterial symbionts are documented to impact host lipid metabolism, the immune system and neuroendocrine pathways, which in turn can affect the nervous system and brain function through microbiota–gut–brain axis pathways (Forsythe and Kunze, 2013; Krishnan et al., 2015). Another study in early life stages of mice demonstrated that the acquisition of gut bacteria enhanced stress responses (Sudo et al., 2004). In the present study, the combined results from survivorship, protist numbers, 16S rRNA gene sequence data and the trophallaxis experiment indicated the reduced stress response capabilities of termites resulting from reduced gut symbiont abundance and diversity. In particular, the stress owing to kanamycin and methoprene was reduced by holding treated individuals with faunated termites in the trophallaxis experiment, as would occur under natural colony conditions. Based on the tightly linked associations among different experimental outcomes, it is clear that gut symbionts and their relative proportions play important roles during the soldier caste differentiation process. Prior experimental evidence reveals how JH is able to select for a differential abundance of fat-body-related CYP4 and hexamerin genes and influence their expression during presoldier initiation (Zhou et al., 2006a,b, 2007). Subsequent studies have shown how JH dramatically reduces protist flagellate gene expression and populations specifically in the gut (Scharf et al., 2017; Sen et al., 2013). The current study builds on these prior findings by showing that JH selects for gut symbionts at certain abundance levels, which in turn contributes to successful completion of the presoldier differentiation process.

Lastly, the differential bacterial compositions we identified here between the Kan and KanMeth treatments are one of the advancements provided by this study. Specifically, we found that JH decreases protists and associated bacteria, while the remaining bacteria help emerging presoldiers survive the transition. However, if the gut microbial balance is disrupted before JH application, termites cannot withstand the physiological stress caused by the molting process. These impacts of JH on the gut microbial system are now disentangled by our findings. As the current study demonstrates, gut symbionts under the influence of JH, besides having importance in termite physiological functions such as digestion, reproduction, immunity, nitrogen fixation and acetogenesis (Cleveland, 1923; Rosengaus et al., 2011; Scharf et al., 2011), are directly linked to presoldier differentiation and emergence and are important in the survival of termites under colony conditions. Moreover, the composition and abundance of these gut microbes may predict the phenotype of the termite host. However, the specific mechanisms involved await determination. Unraveling these mechanisms and the roles of specific bacteria is likely to provide important new insights into termite caste differentiation, and more broadly into the processes related to gut microbial control over host phenotypic plasticity.

We thank the University of Minnesota Genomics Center for library preparation and sequencing of the 16S rRNA gene; Hannah Stewart, Ethan Hillman and Brittany Peterson for assay development in earlier phases of this research; Dr Timothy Johnson and Manoj Pandey for assistance with the data analysis; and Zachery Wolfe for laboratory assistance.

Author contributions

Conceptualization: M.E.S.; Methodology: R.S.; Validation: R.S.; Formal analysis: R.S., C.H.N.; Investigation: R.S.; Data curation: R.S., C.H.N.; Writing - original draft: R.S.; Writing - review & editing: C.H.N., M.E.S.; Visualization: R.S.; Funding acquisition: M.E.S.

Funding

This research was supported by a Ross Graduate Fellowship to R.S., USDA-NIFA-Hatch support to M.E.S. (project no. 1010572), and the O.W. Rollins/Orkin Endowment in the Department of Entomology at Purdue University.

Data availability

Sequences from this study are available through the NCBI Sequence Read Archive database under the BioProject accession number PRJNA693690 and the individual sequences are SRR13495705SRR13495724.

Bolyen
,
E.
,
Rideout
,
J. R.
,
Dillon
,
M. R.
,
Bokulich
,
N. A.
,
Abnet
,
C. C.
,
Al-Ghalith
,
G. A.
,
Alexander
,
H.
,
Alm
,
E. J.
,
Arumugam
,
M.
,
Asnicar
,
F.
et al. 
(
2019
).
Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2
.
Nat. Biotechnol.
37
,
852
-
857
.
Boucias
,
D. G.
,
Cai
,
Y.
,
Sun
,
Y.
,
Lietze
,
V. U.
,
Sen
,
R.
,
Raychoudhury
,
R.
and
Scharf
,
M. E.
(
2013
).
The hindgut lumen prokaryotic microbiota of the termite Reticulitermes flavipes and its responses to dietary lignocellulose composition
.
Mol. Ecol.
22
,
1836
-
1853
.
Brune
,
A.
and
Ohkuma
,
M.
(
2010
).
Role of the termite gut microbiota in symbiotic digestion
. In
Biology of Termites: a Modern Synthesis
(ed.
D.
Bignell
,
Y.
Roisin
and
N.
Lo
), pp.
439
-
475
.
Dordrecht
:
Springer
.
Callahan
,
B. J.
,
McMurdie
,
P. J.
,
Rosen
,
M. J.
,
Han
,
A. W.
,
Johnson
,
A. J. A.
and
Holmes
,
S. P.
(
2016
).
DADA2: high-resolution sample inference from Illumina amplicon data
.
Nat. Method.
13
,
581
-
583
.
Caporaso
,
J. G.
,
Lauber
,
C. L.
,
Walters
,
W. A.
,
Berg-Lyons
,
D.
,
Lozupone
,
C. A.
,
Turnbaugh
,
P. J.
,
Fierer
,
N.
and
Knight
,
R.
(
2011
).
Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample
.
Proc. Natl. Acad. Sci. USA
108
,
4516
-
4522
.
Chan
,
K. K.
,
Abel
,
D. S.
and
Stay
,
B.
(
2011
).
Fine structure of corpora allata of castes with different rates of juvenile hormone production in the termite Reticulitermes flavipes
.
Arth. Struct. Dev.
40
,
26
-
38
.
Charles
,
J.-P.
,
Iwema
,
T.
,
Epa
,
V. C.
,
Takaki
,
K.
,
Rynes
,
J.
and
Jindra
,
M.
(
2011
).
Ligand-binding properties of a juvenile hormone receptor, methoprene-tolerant
.
Proc. Natl. Acad. Sci. USA
108
,
21128
-
21133
.
Chen
,
J.
,
Li
,
Y.
,
Tian
,
Y.
,
Huang
,
C.
,
Li
,
D.
,
Zhong
,
Q.
and
Ma
,
X.
(
2015
).
Interaction between microbes and host intestinal health: modulation by dietary nutrients and gut–brain–endocrine–immune axis
.
Nat. Method.
16
,
592
-
603
.
Cleveland
,
L. R.
(
1923
).
Symbiosis between termites and their intestinal protozoa
.
Proc. Natl. Acad. Sci. USA
9
,
424
.
Cornette
,
R.
,
Matsumoto
,
T.
and
Miura
,
T.
(
2007
).
Histological analysis of fat body development and molting events during soldier differentiation in the damp-wood termite, Hodotermopsis sjostedti (Isoptera, Termopsidae)
.
Zoolog. Sci.
24
,
1066
-
1074
.
Doolittle
,
M.
,
Raina
,
A.
,
Lax
,
A.
and
Boopathy
,
R.
(
2008
).
Presence of nitrogen fixing Klebsiella pneumoniae in the gut of the Formosan subterranean termite (Coptotermes formosanus)
.
Bioresour. Technol.
99
,
3297
-
3300
.
Du
,
H.
,
Tong
,
R. L.
,
Huang
,
X.
,
Liu
,
B.
,
Huang
,
R.
and
Li
,
Z.
(
2020
).
Methoprene-induced genes in workers of Formosan subterranean termites (Coptotermes formosanus Shiraki)
.
Insects
11
,
71
.
Elliott
,
K. L.
and
Stay
,
B.
(
2008
).
Changes in juvenile hormone synthesis in the termite Reticulitermes flavipes during development of soldiers and neotenic reproductives from groups of isolated workers
.
J. Insect Physiol.
54
,
492
-
500
.
Elliott
,
K. L.
,
Chan
,
K. K.
,
Teesch
,
L.
,
Clor
,
O.
and
Stay
,
B.
(
2009
).
Identification of Phe-Gly-Leu-amide type allatostatin-7 in Reticulitermes flavipes: its localization in tissues and relation to juvenile hormone synthesis
.
Peptides
30
,
495
-
506
.
Engelmann
,
F.
(
1984
).
Regulation of vitellogenesis in insects: the pleiotropic role of juvenile hormones
. In
Biosynthesis, Metabolism and Mode of Action of Invertebrate Hormones
(ed.
J.
Hoffmann
and
M.
Porchet
), pp.
444
-
453
.
Berlin, Heidelberg
:
Springer
.
Engl
,
T.
and
Kaltenpoth
,
M.
(
2018
).
Influence of microbial symbionts on insect pheromones
.
Nat. Prod. Rep.
35
,
386
-
397
.
Forsythe
,
P.
and
Kunze
,
W. A.
(
2013
).
Voices from within: gut microbes and the CNS
.
Cell. Mol. Life Sci.
70
,
55
-
69
.
Fröhlich
,
J.
and
König
,
H.
(
1999
).
Rapid isolation of single microbial cells from mixed natural and laboratory populations with the aid of a micromanipulator
.
Syst. Appl. Microbiol.
22
,
249
-
257
.
Howard
,
R. W.
and
Haverty
,
M. I.
(
1978
).
Defaunation, mortality and soldier differentiation: concentration effects of methoprene in a termite
.
Sociobiology
3
,
73
-
77
.
Hrdý
,
I.
,
Kuldová
,
J.
,
Hanus
,
R.
and
Wimmer
,
Z.
(
2006
).
Juvenile hormone III, hydroprene and a juvenogen as soldier caste differentiation regulators in three Reticulitermes species: potential of juvenile hormone analogues in termite control
.
Pest Manag. Sci.
62
,
848
-
854
.
Iida
,
T.
,
Ohkuma
,
M.
,
Ohtoko
,
K.
and
Kudo
,
T.
(
2000
).
Symbiotic spirochetes in the termite hindgut: phylogenetic identification of ectosymbiotic spirochetes of oxymonad protists
.
FEMS Microbiol. Ecol.
34
,
17
-
26
.
Inoue
,
T.
,
Murashima
,
K.
,
Azuma
,
J.-I.
,
Sugimoto
,
A.
and
Slaytor
,
M.
(
1997
).
Cellulose and xylan utilisation in the lower termite Reticulitermes speratus
.
J. Insect Physiol.
43
,
235
-
242
.
Inoue
,
T.
,
Kitade
,
O.
,
Yoshimura
,
T.
and
Yamaoka
,
I
. (
2000
).
Symbiotic associations with protists
. In
Termites: Evolution, Sociality, Symbioses, Ecology
(ed.
T.
Abe
,
D.E.
Bignell
and
M.
Higashi
). pp.
275
-
288
.
Dordrecht
:
Springer
.
Krishnan
,
S.
,
Alden
,
N.
and
Lee
,
K.
(
2015
).
Pathways and functions of gut microbiota metabolism impacting host physiology
.
Curr. Opin. Biotechnol.
36
,
137
-
145
.
Lainé
,
L. V.
and
Wright
,
D. J.
(
2003
).
The life cycle of Reticulitermes spp. (Isoptera: Rhinotermitidae): what do we know?
. Bull. Entomol. Res.
93
,
267
.
Leadbetter
,
J. R.
,
Schmidt
,
T. M.
,
Graber
,
J. R.
and
Breznak
,
J. A.
(
1999
).
Acetogenesis from H2 plus CO2 by spirochetes from termite guts
.
Science
283
,
686
-
689
.
Lee
,
J. B.
,
Park
,
K. E.
,
Lee
,
S. A.
,
Jang
,
S. H.
,
Eo
,
H. J.
,
Am Jang
,
H.
,
Kim
,
C. H.
,
Ohbayashi
,
T.
,
Matsuura
,
Y.
,
Kikuchi
,
Y.
et al. 
(
2017
).
Gut symbiotic bacteria stimulate insect growth and egg production by modulating hexamerin and vitellogenin gene expression
.
Dev. Comp. Immunol.
69
,
12
-
22
.
Lee
,
J.
,
Kim
,
C.-H.
,
Am Jang
,
H.
,
Kim
,
J. K.
,
Kotaki
,
T.
,
Shinoda
,
T.
,
Shinada
,
T.
,
Yoo
,
J.-W.
and
Lee
,
B. L.
(
2019
).
Burkholderia gut symbiont modulates titer of specific juvenile hormone in the bean bug Riptortus pedestris
.
Dev. Comp. Immunol.
99
,
103399
.
Lozupone
,
C. A.
,
Hamady
,
M.
,
Kelley
,
S. T.
and
Knight
,
R.
(
2007
).
Quantitative and qualitative β diversity measures lead to different insights into factors that structure microbial communities
.
Appl. Environ. Microbiol.
73
,
1576
-
1585
.
Lozupone
,
C.
,
Lladser
,
M. E.
,
Knights
,
D.
,
Stombaugh
,
J.
and
Knight
,
R.
(
2011
).
UniFrac: an effective distance metric for microbial community comparison
.
ISME J.
5
,
169
-
172
.
Lucey
,
K. S.
and
Leadbetter
,
J. R.
(
2014
).
Catechol 2, 3–dioxygenase and other meta–cleavage catabolic pathway genes in the ‘anaerobic’ termite gut spirochete Treponema primitia
.
Mol. Ecol.
23
,
1531
-
1543
.
Mandal
,
S.
,
Van Treuren
,
W.
,
White
,
R. A.
,
Eggesbø
,
M.
,
Knight
,
R.
and
Peddada
,
S. D.
(
2015
).
Analysis of composition of microbiomes: a novel method for studying microbial composition
.
Microb. Ecol. Health Dis.
26
,
27663
.
Mao
,
L.
,
Henderson
,
G.
,
Liu
,
Y.
and
Laine
,
R. A.
(
2005
).
Formosan subterranean termite (Isoptera: Rhinotermitidae) soldiers regulate juvenile hormone levels and caste differentiation in workers
.
Ann. Entomol. Soc. Am.
98
,
340
-
345
.
McMurdie
,
P. J.
and
Holmes
,
S.
(
2013
).
phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data
.
PLoS ONE
8
,
e61217
.
Nalepa
,
C. A.
(
2017
).
What kills the hindgut flagellates of lower termites during the host molting cycle?
Microorganisms
5
,
82
.
Nijhout
,
H. F.
,
Riddiford
,
L. M.
,
Mirth
,
C.
,
Shingleton
,
A. W.
,
Suzuki
,
Y.
and
Callier
,
V.
(
2014
).
The developmental control of size in insects
.
Wiley Interdisciplinary Reviews Dev. Biol.
3
,
113
-
134
.
Noda
,
S.
,
Iida
,
T.
,
Kitade
,
O.
,
Nakajima
,
H.
,
Kudo
,
T.
and
Ohkuma
,
M.
(
2005
).
Endosymbiotic Bacteroidales bacteria of the flagellated protist Pseudotrichonympha grassii in the gut of the termite Coptotermes formosanus
.
Appl. Environ. Microbiol.
71
,
8811
-
8817
.
Noda
,
S.
,
Kitade
,
O.
,
Inoue
,
T.
,
Kawai
,
M.
,
Kanuka
,
M.
,
Hiroshima
,
K.
,
Hongoh
,
Y.
,
Constantino
,
R.
,
Uys
,
V.
,
Zhong
,
J.
et al. 
(
2007
).
Cospeciation in the triplex symbiosis of termite gut protists (Pseudotrichonympha spp.), their hosts, and their bacterial endosymbionts
.
Mol. Ecol.
16
,
1257
-
1266
.
Ohkuma
,
M.
(
2008
).
Symbioses of flagellates and prokaryotes in the gut of lower termites
.
Trends Microbiol.
16
,
345
-
352
.
Ohkuma
,
M.
,
Sato
,
T.
,
Noda
,
S.
,
Ui
,
S.
,
Kudo
,
T.
and
Hongoh
,
Y.
(
2007
).
The candidate phylum ‘Termite Group 1’ of bacteria: phylogenetic diversity, distribution, and endosymbiont members of various gut flagellated protists
.
FEMS Microbiol.Ecol.
60
,
467
-
476
.
Otani
,
S.
,
Mikaelyan
,
A.
,
Nobre
,
T.
,
Hansen
,
L. H.
,
Koné
,
N. G. A.
,
Sørensen
,
S. J.
,
Aanen
,
D. K.
,
Boomsma
,
J. J.
,
Brune
,
A.
and
Poulsen
,
M.
(
2014
).
Identifying the core microbial community in the gut of fungus-growing termites
.
Mol. Ecol.
23
,
4631
-
4644
.
Park
,
Y. I.
and
Raina
,
A. K.
(
2004
).
Juvenile hormone III titers and regulation of soldier caste in Coptotermes formosanus (Isoptera: Rhinotermitidae)
.
J. Insect Physiol.
50
,
561
-
566
.
Peterson
,
B. F.
,
Stewart
,
H. L.
and
Scharf
,
M. E.
(
2015
).
Quantification of symbiotic contributions to lower termite lignocellulose digestion using antimicrobial treatments
.
Insect Biochem. Mol. Biol
59
,
80
-
88
.
Peterson
,
B. F.
and
Scharf
,
M. E.
(
2016a
).
Lower termite associations with microbes: synergy, protection, and interplay
.
Front. Microbiol.
7
,
422
.
Peterson
,
B. F.
and
Scharf
,
M. E.
(
2016b
).
Metatranscriptome analysis reveals bacterial symbiont contributions to lower termite physiology and potential immune functions
.
BMC Genom
17
,
772
.
Poutahidis
,
T.
,
Kearney
,
S. M.
,
Levkovich
,
T.
,
Qi
,
P.
,
Varian
,
B. J.
,
Lakritz
,
J. R.
,
Ibrahim
,
Y. M.
,
Chatzigiagkos
,
A.
,
Alm
,
E. J.
and
Erdman
,
S. E.
(
2013
).
Microbial symbionts accelerate wound healing via the neuropeptide hormone oxytocin
.
PLoS ONE
8
,
e78898
.
Quast
,
C.
,
Pruesse
,
E.
,
Yilmaz
,
P.
,
Gerken
,
J.
,
Schweer
,
T.
,
Yarza
,
P.
,
Peplies
,
J.
and
Glöckner
,
F. O.
(
201
3).
The SILVA ribosomal RNA gene database project: improved data processing and web-based tools
.
Nucleic Acids Res.
41
,
D590
-
D596
.
Roisin
,
Y
. (
2000
).
Diversity and evolution of caste patterns
. In
Termites: Evolution, Sociality, Symbioses, Ecology
(ed.
T.
Abe
,
D. E.
Bignell
and
M.
Higashi
), pp.
95
-
119
.
Dordrecht
:
Springer
.
Rosengaus
,
R. B.
,
Zecher
,
C. N.
,
Schultheis
,
K. F.
,
Brucker
,
R. M.
and
Bordenstein
,
S. R.
(
2011
).
Disruption of the termite gut microbiota and its prolonged consequences for fitness
.
Appl. Environ. Microbiol.
77
,
4303
-
4312
.
Scharf
,
M. E.
and
Peterson
,
B. F.
(
2021
).
A century of synergy in termite symbiosis research: linking the past with new genomic insights
.
Ann. Rev. Entomol.
66
,
23
-
43
.
Scharf
,
M. E.
and
Tartar
,
A.
(
2008
).
Termite digestomes as sources for novel lignocellulases
.
Biofuel Bioprod. Biorefin
2
,
540
-
552
.
Scharf
,
M. E.
,
Wu-Scharf
,
D.
,
Pittendrigh
,
B. R.
and
Bennett
,
G. W.
(
2003
).
Caste-and development-associated gene expression in a lower termite
.
Genome Biol.
4
,
R62
.
Scharf
,
M. E.
,
Ratliff
,
C. R.
,
Wu-Scharf
,
D.
,
Zhou
,
X.
,
Pittendrigh
,
B. R.
and
Bennett
,
G. W.
(
2005
).
Effects of juvenile hormone III on Reticulitermes flavipes: changes in hemolymph protein composition and gene expression
.
Insect Biochem. Mol. Biol.
35
,
207
-
215
.
Scharf
,
M. E.
,
Karl
,
Z. J.
,
Sethi
,
A.
and
Boucias
,
D. G.
(
2011
).
Multiple levels of synergistic collaboration in termite lignocellulose digestion
.
PLoS ONE
6
,
e21709
.
Scharf
,
M. E.
,
Cai
,
Y.
,
Sun
,
Y.
,
Sen
,
R.
,
Raychoudhury
,
R.
and
Boucias
,
D. G.
(
2017
).
A meta-analysis testing eusocial co-option theories in termite gut physiology and symbiosis
.
Commun. Integr. Biol.
10
,
e1295187
.
Sen
,
R.
,
Raychoudhury
,
R.
,
Cai
,
Y.
,
Sun
,
Y.
,
Lietze
,
V.-U.
,
Boucias
,
D. G.
and
Scharf
,
M. E.
(
2013
).
Differential impacts of juvenile hormone, soldier head extract and alternate caste phenotypes on host and symbiont transcriptome composition in the gut of the termite Reticulitermes flavipes
.
BMC Genom.
14
,
491
.
Smykal
,
V.
,
Daimon
,
T.
,
Kayukawa
,
T.
,
Takaki
,
K.
,
Shinoda
,
T.
and
Jindra
,
M.
(
2014
).
Importance of juvenile hormone signaling arises with competence of insect larvae to metamorphose
.
Dev. Biol.
390
,
221
-
230
.
Stephens
,
M. E.
and
Gage
,
D. J.
(
2020
).
Single-cell amplicon sequencing reveals community structures and transmission trends of protist-associated bacteria in a termite host
.
PLoS ONE
15
,
e0233065
.
Stingl
,
U.
,
Radek
,
R.
,
Yang
,
H.
and
Brune
,
A.
(
2005
).
“Endomicrobia”: cytoplasmic symbionts of termite gut protozoa form a separate phylum of prokaryotes
.
Appl. Environ. Microbiol
71
,
1473
-
1479
.
Storelli
,
G.
,
Defaye
,
A.
,
Erkosar
,
B.
,
Hols
,
P.
,
Royet
,
J.
and
Leulier
,
F.
(
2011
).
Lactobacillus plantarum promotes Drosophila systemic growth by modulating hormonal signals through TOR-dependent nutrient sensing
.
Cell Metab.
14
,
403
-
414
.
Sudo
,
N.
,
Chida
,
Y.
,
Aiba
,
Y.
,
Sonoda
,
J.
,
Oyama
,
N.
,
Yu
,
X. N.
,
Kubo
,
C.
and
Koga
,
Y.
(
2004
).
Postnatal microbial colonization programs the hypothalamic–pituitary–adrenal system for stress response in mice
.
J. Physiol.
558
,
263
-
275
.
Tarver
,
M. R.
,
Schmelz
,
E. A.
,
Rocca
,
J. R.
and
Scharf
,
M. E.
(
2009
).
Effects of soldier-derived terpenes on soldier caste differentiation in the termite Reticulitermes flavipes
.
J. Chem. Ecol
35
,
256
-
264
.
Tarver
,
M. R.
,
Florane
,
C. B.
,
Zhang
,
D.
,
Grimm
,
C.
,
Lax
,
A. R.
and
Traniello
,
J.
(
2012
).
Methoprene and temperature effects on caste differentiation and protein composition in the Formosan subterranean termite, Coptotermes formosanus
.
J. Insect. Sci.
12
,
18
.
ter Braak
,
C. J.
(
1986
).
Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis
.
Ecology
67
,
1167
-
1179
.
Treitli
,
S. C.
,
Kolisko
,
M.
,
Husník
,
F.
,
Keeling
,
P. J.
and
Hampl
,
V.
(
2019
).
Revealing the metabolic capacity of Streblomastix strix and its bacterial symbionts using single-cell metagenomics
.
Proc. Natl. Acad. Sci. USA
116
,
19675
-
19684
.
Wilson
,
E. O.
(
1971
).
The Insect Societies
.
Cambridge, MA
:
Harvard University Press
.
Yagi
,
S
. (
1976
).
The role of juvenile hormone in diapause and phase variation in some lepidopterous insects
. In
The Juvenile Hormones
(ed.
L. I.
Gilbert
), pp.
288
-
300
.
Boston, MA
:
Springer
.
Yagi
,
K. J.
,
Kwok
,
R.
,
Chan
,
K. K.
,
Setter
,
R. R.
,
Myles
,
T. G.
,
Tobe
,
S. S.
and
Stay
,
B.
(
2005
).
Phe-Gly-Leu-amide allatostatin in the termite Reticulitermes flavipes: content in brain and corpus allatum and effect on juvenile hormone synthesis
.
J.Insect Physiol.
51
,
357
-
365
.
Zheng
,
H.
,
Powell
,
J. E.
,
Steele
,
M. I.
,
Dietrich
,
C.
and
Moran
,
N. A.
(
2017
).
Honeybee gut microbiota promotes host weight gain via bacterial metabolism and hormonal signaling
.
Proc. Natl. Acad. Sci. USA
114
,
4775
-
4780
.
Zhou
,
X.
,
Oi
,
F. M.
and
Scharf
,
M. E.
(
2006a
).
Social exploitation of hexamerin: RNAi reveals a major caste-regulatory factor in termites
.
Proc. Natl. Acad. Sci. USA
103
,
4499
-
4504
.
Zhou
,
X.
,
Song
,
C.
,
Grzymala
,
T. L.
,
Oi
,
F. M.
and
Scharf
,
M. E.
(
2006b
).
Juvenile hormone and colony conditions differentially influence cytochrome P450 gene expression in the termite Reticulitermes flavipes
.
Insect Mol. Biol.
15
,
749
-
761
.
Zhou
,
X.
,
Tarver
,
M. R.
and
Scharf
,
M. E.
(
2007
).
Hexamerin-based regulation of juvenile hormone-dependent gene expression underlies phenotypic plasticity in a social insect
.
Development
134
,
601
-
610
.
Zilber-Rosenberg
,
I.
and
Rosenberg
,
E.
(
2008
).
Role of microorganisms in the evolution of animals and plants: the hologenome theory of evolution
.
FEMS Microbiol. Rev.
32
,
723
-
735
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information