Email updates

Keep up to date with the latest news and content from BMC Genomics and BioMed Central.

Open Access Research article

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

Ruchira Sen1, Rhitoban Raychoudhury1, Yunpeng Cai24, Yijun Sun25, Verena-Ulrike Lietze3, Drion G Boucias3 and Michael E Scharf1*

Author Affiliations

1 Department of Entomology, Purdue University, West Lafayette, IN, USA

2 Interdisciplinary Center for Biotechnology Research, University of Florida, Gainesville, FL, USA

3 Entomology and Nematology Department, University of Florida, Gainesville, FL, USA

4 Current Address: Research Center for Biomedical Information Technology, Shenzhen Institutes of Advance Technology, Chinese Academy of Sciences, Shenzhen, China

5 Current Address: Department of Microbiology and Immunology & New York State Center of Excellence in Bioinformatics and Life Sciences, The State University of New York at Buffalo, Buffalo, NY 14203, USA

For all author emails, please log on.

BMC Genomics 2013, 14:491  doi:10.1186/1471-2164-14-491


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2164/14/491


Received:11 March 2013
Accepted:9 July 2013
Published:19 July 2013

© 2013 Sen et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Termites are highly eusocial insects and show a division of labor whereby morphologically distinct individuals specialize in distinct tasks. In the lower termite Reticulitermes flavipes (Rhinotermitidae), non-reproducing individuals form the worker and soldier castes, which specialize in helping (e.g., brood care, cleaning, foraging) and defense behaviors, respectively. Workers are totipotent juveniles that can either undergo status quo molts or develop into soldiers or neotenic reproductives. This caste differentiation can be regulated by juvenile hormone (JH) and primer pheromones contained in soldier head extracts (SHE). Here we offered worker termites a cellulose diet treated with JH or SHE for 24-hr, or held them with live soldiers (LS) or live neotenic reproductives (LR). We then determined gene expression profiles of the host termite gut and protozoan symbionts concurrently using custom cDNA oligo-microarrays containing 10,990 individual ESTs.

Results

JH was the most influential treatment (501 total ESTs affected), followed by LS (24 ESTs), LR (12 ESTs) and SHE treatments (6 ESTs). The majority of JH up- and downregulated ESTs were of host and symbiont origin, respectively; in contrast, SHE, LR and LS treatments had more uniform impacts on host and symbiont gene expression. Repeat “follow-up” bioassays investigating combined JH + SHE impacts in relation to individual JH and SHE treatments on a subset of array-positive genes revealed (i) JH and SHE treatments had opposite impacts on gene expression and (ii) JH + SHE impacts on gene expression were generally intermediate between JH and SHE.

Conclusions

Our results show that JH impacts hundreds of termite and symbiont genes within 24-hr, strongly suggesting a role for the termite gut in JH-dependent caste determination. Additionally, differential impacts of SHE and LS treatments were observed that are in strong agreement with previous studies that specifically investigated soldier caste regulation. However, it is likely that gene expression outside the gut may be of equal or greater importance than gut gene expression.

Keywords:
Metagenomics; Soldier head extract; Microarray; Caste differentiation; Live soldier; Live reproductive

Background

The successful maintenance of social insect colonies relies on the efficiency of workers. Functions like foraging, cleaning, brood care, colony defense, etc. are carried out by individuals, which are often differentiated into worker sub-castes that each performs specific duties. Such differentiation is either morphological, where physical features determine the tasks (polyphenism), or age-based, where workers carry out different functions at different ages (polyethism) [1]. In both types of caste differentiation, juvenile hormone (JH) is known to play a significant role in most social insects [2]. Experimental application of JH to workers of different social insects has been shown to induce two major types of changes; such as (i) stimulating young workers to carry out functions that are usually performed by older individuals [3] or (ii) inducing workers to go through physical changes and differentiate from one phenotype to another [4].

Termite colonies contain one or more pairs of reproductives and a large number of non-reproductives that are morphologically differentiated into workers, nymphs and soldiers. Workers and nymphs carry out foraging, brood care and cleaning while soldiers, with bigger and stronger mandibles, are dedicated to colony defense. As hemimetabolous insects, termites go through several juvenile instars before reaching adulthood; in each of these juvenile instars JH is presumed to perform its stereotypical “status quo” function [2]. In lower termites, including Reticulitermes flavipes, helper castes are juvenile stages composed of both workers, which are eyeless and wingless, and nymphs, which are immature imagoes [5]. Nymphs and workers diverge from undifferentiated “larvae” after the second instar. Workers have three alternate developmental trajectories that include (i) status quo molts into workers, (ii) molts into presoldiers (followed by terminal soldier differentiation), or (iii) molts into apterous neotenic reproductives. Nymphs, conversely, have two alternate developmental trajectories that include molts into (i) brachypterous neotenic reproductives or (ii) adult imagoes that eventually become primary reproductives that start incipient colonies. Aside from two caste-regulatory genes and two soldier-derived primer pheromones linked to presoldier caste regulation [5-7], relatively little is known about the molecular mechanisms of caste differentiation in termites and the factors that initiate this process.

Different juvenile hormone analogues (JHAs) have been tested for their effect on caste differentiation in many species of termites through various bioassay methods [8]. Most of these experiments demonstrated increased production of presoldiers, intercastes and pseudoimagoes and increased or decreased production of neotenic reproductives. The application of JH also resulted in defaunation of protozoan symbionts and bacterial endosymbionts, inhibition of ecdysis and feeding, atrophication of the prothoracic gland, and apparent toxicity [8]. It can thus be hypothesized that JH and host-symbiont interactions within the gut may impact caste differentiation processes.

In R. flavipes and other species of the genus Reticulitermes, JH and JHAs predominantly induce presoldier differentiation [7,9-12]. Moreover, Elliott and Stay [13] have suggested that, in R. flavipes, workers that are destined to become presoldiers have a 2.5-fold higher JH titer compared to workers destined to become neotenic reproductives. Park and Raina [14,15] also reported an increased titer of JH in presoldiers and new soldiers of a closely related species, Coptotermes formosanus. Soldiers are known to have inhibitory effects on presoldier differentiation mediated by soldier-derived primer pheromones in lower termites [16-18]. Two primer pheromone candidates, ɣ-cadinene (CAD) and its aldehyde ɣ-cadinenal (ALD) have been identified in R. flavipes soldier head extract (SHE), which, when applied in combination with JH, increases soldier caste differentiation [7,12]. However, SHE alone does not impact caste differentiation, survivorship, or any other aspect of worker biology [7,12,18]. Also, two fat body-expressed hexamerin-encoding genes (Hex-1 and Hex-2) play a key role in maintaining a developmental status quo in workers by being JH-inducible, sequestering JH and thereby promoting high worker caste proportions [11,19-22].

In addition to endocrine effects, social effects on gene expression have been investigated to some extent in R. flavipes. For example, being held with soldiers increases levels of the primer pheromone ALD by 10× in workers [7] and such workers are less likely to undergo presoldier formation [7,18]. However, with respect to reproductive effects, while some phenotypic impacts have been noted in R. speratus[37], in R. flavipes the impact on workers of being held with live reproductives is not known. Termite biology, however, is also influenced by the numerous symbionts that are harbored in the termite gut. In R. flavipes, these symbionts consist of both pro- and eukaryotes and include over 5,000 ribotypes of prokaryotes [23] and 11–12 different protists [24]. The protists, especially, have been shown to be in a nutritional symbiosis with R. flavipes[25]. Despite this emerging understanding, no studies have yet focused on the worker termite gut or its resident symbionts as potential molecular determinants of caste differentiation. Since the lower-termite gut environment is centrally important to nutrition, physiology and symbiotic relationships with protists and bacteria, we hypothesized that gene expression in this environment is substantially altered in response to caste-regulatory factors and caste composition. Raychoudhury et al.[26] recently tested the effects of different diets on gene expression of both the termite gut and protist symbionts using an oligonucleotide cDNA microarray. Here we use the same protocol to test the effects of hormones (i.e., JH), primer pheromones (SHE) and the presence of live reproductives (LR) and live soldiers (LS) on host and symbiont gene expression in the gut of R. flavipes workers. We show that JH has the most pronounced effect while SHE treatment and the presence of live reproductives or soldiers have much lower impacts on gene expression in the gut environment. We further investigated the expression level of a few selected upregulated and downregulated genes (as found by microarray) in termites independently treated with JH, SHE and the JH + SHE combination. Lastly, we report details on a previously un-annotated 50-kDa midgut protein gene that was the most significantly JH-upregulated EST.

Results

JH-presoldier induction bioassays

Presoldier induction bioassays were conducted on one Florida colony used for microarray studies (B1) and another from Indiana used in “follow up” qPCR bioassays (WI-9). After a limited 24-hr exposure period identical to that employed in microarray and follow up bioassays, the B1 and WI-9 colonies had average ± std. deviation presoldier formation of 90.0 ± 7.1% and 91.1 ± 11.5% after 25 days (n = 4-5 per colony). No presoldier formation (0%) occurred in acetone controls.

Numbers of differentially expressed transcripts and their annotation

A total of 10,990 distinct transcripts, represented in 15,208 array positions, were studied in this microarray experiment (see Methods for details). To visualize the array data, we calculated a fold change ratio for every array position based on its normalized average intensity in each treatment (i.e., JH, SHE, LS, LR) divided by its normalized average intensity in the control (acetone) treatment. We selected array positions that had significant fold ratios (P <0.05) and mean log2 fold change ratios ≤ −0.25 and ≥0.25 (emphasized in the Volcano plots in Figure 1) which corresponded to actual fold changes of <0.84 and >1.19, respectively. Array positions with ≥0.25 mean log2 fold change ratio represented ESTs which were upregulated by the various treatments, and array positions with ≤ −0.25 mean log2 fold change ratio represented ESTs which were downregulated by the various treatments.

thumbnailFigure 1. Volcano plots illustrating microarray results. Shown are numbers of significantly (P <0.05) upregulated (blue; ≥1.19-fold) and downregulated (red; ≤0.84-fold) microarray positions from four treatments that included juvenile hormone (A), soldier head extract (B), live reproductives (C) and live soldiers (D). The fold ratios on the x-axis represent the ratio of EST abundance for each of the four treatments relative to untreated controls. The numbers of contigs formed from ESTs of host or symbiont origin are given in each panel.

Effects of juvenile hormone III (JH)

A total of 179 and 981 JH up- and downregulated microarray positions were identified. After contiging the upregulated and downregulated ESTs separately (90% similarity level), we found 96 and 12 unique sequences of host and symbiont origin, respectively, in the upregulated set; in the downregulated set, we found 44 and 349 unique sequences of host and symbiont origin, respectively, along with 6 sequences of mixed origin (Table 1, Figure 1). The term “EST” is used hereafter to refer to both contigs and non-contiguous orphan EST “singletons”. Significantly more upregulated ESTs were from the host termite; whereas, significantly more downregulated ESTs were from protist symbionts (G test for independence, G = 242.7, P <0.001). Out of 108 JH upregulated ESTs, 74 could be annotated by BLASTx against the nr database through BLAST2GO, and 332 out of 399 downregulated ESTs could be annotated (Tables 2 and 3, Additional file 1: Table S1A, S1B).

Table 1. Number of upregulated and downregulated ESTs in different treatments compared to acetone controls

Table 2. KEGG pathways upregulated by JH treatment

Table 3. KEGG pathways downregulated by JH treatment

Additional file 1. Identity, fold change and gene ontology terms for passing host and symbiont genes from JH microarrays (Additional file 1: Table S1-S4).

Format: DOCX Size: 155KB Download fileOpen Data

Effects of SHE (soldier head extract)

Three and 5 SHE up- and downregulated array positions, respectively, were identified. Each of the 3 upregulated ESTs was unique, i.e., they did not form any contig at 90% similarity, whereas the 5 downregulated ESTs contiged into 3 unique sequences. In the upregulated set, 2 out of 3 ESTs were of host origin while all 3 downregulated contig sequences were of symbiont origin (Table 1, Figure 1). Only 1 host EST out of the 3 upregulated ESTs and 3 symbiont downregulated ESTs could be annotated with BLASTx searches of the nr database (Additional file 1: Table S2A, S2B).

Effects of live reproductives (LR)

Ten and 3 LR up- and down regulated array positions were identified. Only 2 upregulated host ESTs formed contigs while all other ESTs were unique. In the upregulated set, 6 (out of 9) ESTs were of host origin while all 3 sequences in the downregulated set were of symbiont origin (Table 1, Figure 1). Among the upregulated ESTs, 4 host ESTs and 1 symbiont EST could be annotated by BLASTx against the nr database (Additional file 1: Table S3A); whereas all 3 downregulated symbiont ESTs had nr database matches (Additional file 1: Table S3B).

Effects of live soldiers (LS)

The presence of live soldiers had a more pronounced impact on worker gut gene expression. A total of 37 and 6 LS up- and downregulated array positions, respectively, were identified (Figure 1). These 43 array positions contiged into 12 host and 6 symbiont-derived up- and downregulated sequences, respectively. The ESTs in the downregulated set did not form any contigs and only 1 out of these 6 was of host origin. Among the upregulated ESTs, 6 and 4 host and symbiont ESTs, respectively, could be annotated by BLASTx against the nr database through BLAST2GO. In the downregulated set, only 2 symbiont ESTs set could be annotated (Additional file 1: Table S4A and B).

Gene ontology and KEGG pathway analyses with BLAST2GO and DAVID

Gene Ontology (GO) terms were obtained through BLAST2GO and were performed in the three categories Molecular Function, Cellular Location and Biological Process for the JH, SHE, LR and LS datasets, with results for each passing EST shown in Additional file 1: Table S1-S4. In addition, Molecular Function GO results were tallied across all 4 treatment categories and are summarized in Table 4. As expected, the highest numbers of GO-Molecular Function terms were in the JH up- and downregulated categories (411 and 91 terms). In most other treatments categories, numbers of passing ESTs had proportional numbers of associated GO-MolecularFunction terms; however, the SHE downregulated category, which had only 2 passing ESTs (dna replication licensing factor mcm7 and serine/threonine-protein kinase mph1), had 19 GO-Molecular Function terms.

Table 4. Summary of molecular function gene ontology (GO) terms obtained for annotated expressed sequence tags (ESTs) with BLASTx for treatments with juvenile hormone (JH), soldier head extract (SHE), exposure to live soldiers (LS) and exposure to live reproductive (LR)

Next, KEGG pathway analysis was carried out on the JH dataset to understand the general effects of JH (Tables 2 and 3) on cellular metabolism in the gut transcriptome. Only the JH dataset was examined because of the relatively large numbers of genes available as compared with the SHE, LR and LS datasets. All annotated upregulated ESTs were of host origin, while all but one of the downregulated ESTs were of symbiont origin. Aldehyde dehydrogenase and peroxidase enzymes were abundant in the upregulated dataset and found to be involved in various pathways that included intermediary and amino acid metabolism and cuticle biosynthesis (Table 2). The downregulated KEGG dataset was much more diverse but most notably included enzymes and pathways linked to terpenoid biosynthesis (e.g., HMG Co-A reductase; Table 3).

Finally, the program DAVID [27] was used for further bioinformatic analyses of ESTs with identifiable homologs in Drosophila melanogaster. DAVID utilizes the background annotation information from various sequenced genomes to assign enrichment scores. Out of the 74 and 332 annotated up- and downregulated ESTs, respectively, 50 and 156 D. melanogaster homologs were obtained. (Additional file 1: Table S1A, and S1B). This result suggests JH-related pathways are well conserved between D. melanogaster and R. flavipes. We present three sets of analyses with JH treatment which show enrichment of GO-Molecular Function terms in the sequences that were downregulated from both host and symbiont libraries and upregulated from the host library. From the symbiont library there were only three annotated ESTs that were upregulated, which proved insufficient for enrichment scores. Additional file 1: Table S1C shows the GO-Molecular Function terms that were enriched from the upregulated host fraction. In general, molecular functions related to iron binding, peptidase activity and cuticle formation were enriched. However, molecular functions related to peptidase activity were also enriched among the downregulated host ESTs (Additional file 1: Table S1D). Clearly, ESTs with similar molecular functions were both up- and downregulated by JH treatment; however, further details cannot be elucidated by the present analyses. The downregulated sequences from the symbiont fraction, conversely, show a clear trend of shutting down vital cellular functions like translation, transcription and other enzymatic activities (Additional file 1: Table S1E). Whether this reflects the capability of the protists to detect the presence of JH and the ensuing molt of the workers towards non-feeding pre-soldier phenotypes, or other host functions which suppress their cellular activities, remains to be investigated.

Candidate genes of interest

The JH microarray dataset contained 501 significantly differentially expressed EST contigs (Table 1, Additional file 1: Table S1A, B). The most highly JH up- and downregulated of these ESTs encoded homologs of a host 50 kDa midgut protein (50MGP), which was upregulated 2.9-fold with JH treatment, and a protist symbiont cysteine synthase a, which was downregulated 3.1-fold with JH treatment. Other highly upregulated genes in the JH dataset included nli interacting factor-likephosphatase, Arylsulfatase B, and several apolipoproteins, chymotrypsins and serine proteases. Two cytochrome P450 homologs from families 6 and 4 were upregulated in the JH dataset and five others from families 6 and 4 were downregulated. Genes related to cuticle formation were upregulated by JH, including larval and pupal cuticle proteins, resilin, Tyramine beta hydroxylase and Dopamine N-acetyl transferase. The JH dataset also contained a number of upregulated host genes with links to phosphate related post-translational modification, e.g., nli interacting factor-likephosphatase, seven de-phosphorylating phosphatases, seventeen kinases (predicted to add phosphate groups), and an insulin receptor homolog with predicted kinase activity.

SHE arrays revealed only 6 ESTs with significant differential expression. The most highly SHE up- and downregulated transcripts were an un-annotated host gene upregulated 2.0-fold and a symbiont serine/threonine-protein kinase mph1 homolog that was 1.4-fold downregulated. Also, a DNA replication factor, dna replication licensing factor mcm7, was downregulated 1.2-fold by SHE treatment.

LR arrays revealed 12 significantly differentially expressed ESTs. The two most highly LR upregulated genes were both from the host and had ~2.2-fold increased expression; the first was a serine protease 13 homolog, and the second had no significant database matches. The most LR downregulated gene (1.6-fold) encoded a symbiont linker histone h1 and h5 family protein.

LS arrays revealed 25 differentially expressed ESTs. With 2.7-fold up-regulation, a host-derived venom allergen 3-like homolog with predicted protease functions was the most highly LS upregulated EST. Five additional upregulated ESTs had predicted carbohydrate-active and immune functions; these ESTs included two alpha amylase homologs (2.0- to 2.4-fold), two lysozyme homologs (1.4 and 1.5-fold), and a C-type lectin homolog (1.5-fold). The most downregulated ESTs in the LS dataset included 1 host and 1 symbiont EST with no significant database matches.

50 kDa midgut protein: cDNA sequence and amino acid translations

As noted above, the 50MGP gene was the most highly JH upregulated gene identified in this study. The full length 50MGP cDNA was assembled from six overlapping EST contigs [28] and verified by database searches and independent resequencing as described under Methods. The 50MGP cDNA and translated amino acid sequences (Genbank Accession No. KC751537) are provided in Additional file 2: Figure S1. The cDNA sequence contains at least 60 base pairs (bp) of 5′ untranslated region (UTR) before the ATG start codon, a 1419-bp open reading frame, and a 3′ UTR of 133 bp between the “taa” termination codon and the first base of the poly-A tail. The 3′ UTR also contains an “aataa” polyadenylation signal 20 bp upstream of the poly-A tail.

Additional file 2. 50 kDa midgut protein translation (Additional file 2: Figure S1) and HPLC analysis of soldier head extract (Additional file 2: Figure S2).

Format: DOCX Size: 4.5MB Download fileOpen Data

The translated 50MGP protein sequence had 473 amino acids (AA) and begins with a 19-AA secretory signal peptide “MKTQAILIAAVALLLGTEG”, indicating the mature protein is soluble and secreted. The mature protein without signal peptide contains 454 AA with a predicted mass of 49.9 kDa and isoelectric point (pI) of 7.09. The translated AA sequence has 17 predicted phosphorylation sites on 10 Ser, 4 Thr and 3 Tyr residues (see gray shading in Additional file 2: Figure S1). There are no predicted glycosylation sites. The 50MGP protein has predicted function as a cell-envelope enzyme with lyase and/or ligase activity. Gene ontology categories predicted for the full length AA sequence, from most to least probable, are “immune response”, “stress response” and “signal transduction.”

A ClustalW alignment of homologous 50MGP proteins from various insects is shown in Figure 2. The alignment contained all homologs that could be identified by BLASTx searches of the GenBank nr database as of Nov 2012. Included in the alignment were homologs sharing 20-35% AA identity from R. flavipes, the sandfly Phlebotomous papatasi, the bark beetle Dendroctonus ponderosae, and the red flour beetle Tribolium castaneum (two homologs). There were 38 invariant AA residues in the alignment, with P. papatasi sharing the most identity with R. flavipes (24.5%), followed by D. ponderosae (22.8%), and T. castaneum (22.2 and 21.6%). All homologs are equally rich in phosphorylation sites. The sandfly 50MGP is most similar to the termite 50MGP protein; it shares a secretion signal peptide, predicted cell envelope interaction, and predicted ligase/lyase, immune, stress response and signal transduction functions with the termite protein.

thumbnailFigure 2. A ClustalW multiple alignment of homologous 50 kDa Midgut Protein sequences from various insects, including R. flavipes (current study), Phlebotomous papatasi (GenBank accession No. ABV44742), Dendroctonus ponderosae (AEE61518) and Tribolium castaneum (EFA 07930, XP976444). Gray shaded amino acids are identical to the R. flavipes sequence. Signal peptides for the R. flavipes and P. papatasi sequences are enclosed in boxes.

Validation of microarray results by qRT-PCR

Correlation analysis

A subset of 52 ESTs was used to validate relative expression levels determined from microarray hybridization data (Figure 3). This subset included 18 and 22 JH up- and downregulated ESTs, 1 and 2 SHE up- and downregulated ESTs, 3 and 6 LR up- and downregulated ESTs, and 6 LS upregulated ESTs. Template cDNA used in these qRT-PCR reactions was reverse-transcribed from the original RNA samples used for microarray hybridizations. A test of correlation was carried out between the 2-ΔΔCT values (obtained from qPCR CT values) and microarray fold change values of the JH up- and downregulated ESTs. The ΔΔCT values were positively correlated (Spearman Rank Correlation, Rs = 0.874, P <0.001) with array fold change values as expected since higher fold change indicates presence of more transcripts, which provides smaller CT values. The sample sizes for other treatments were too small for conducting statistical tests; however, ΔΔCT values for all but one tested ESTs showed similar negative correlation trends (i.e., there were negative ΔΔCT values for upregulated ESTs and positive ΔΔCT values for downregulated ESTs, Additional file 3: Table S5).

thumbnailFigure 3. Significant positive correlation between fold change and threshold cycle (CT) values (Spearman rank correlation, Rs= 0.874, N = 40, P <0.001). Each data point represents 2-ΔΔCT of quantitative real-time PCR performed on a subset of sequences to verify the robustness of microarray results. Red and blue spots represent down- and upregulated transcripts, respectively. The black spots represent transcripts that were neutral to the treatment.

Additional file 3. Summary of ΔΔCT values for repeat bioassay qPCRs (Additional file3: Table S5) and the list of primers used for qPCR validations (Additional file3: Table S6).

Format: DOCX Size: 30KB Download fileOpen Data

“Follow-up” bioassays investigating select candidate genes in JH, SHE and JH + SHE treatments

We also reassessed the microarray results using qRT-PCR with cDNA samples from “follow-up” bioassays that exposed new worker termites to JH and SHE for 1 day. Additionally, a combination treatment of JH + SHE that was not included in microarray studies was evaluated. SHE extracts were prepared as described previously and verified by HPLC [7]; Additional file 2: Figure S2). ESTs tested in this experiment included 50MGP (two ESTs; FL636982 and FL636656), Apolipoprotein d (FL640421), Radial Spoke Protein (FL643521), Unknown Ribosomal Protein (FL644772), DNA Replication Licensing Factor (FL644436), and SoldierSpecific ProteinNtSp1-like” (FL637031). First, a correlation analysis comparing follow-up bioassay results showed a high degree of correlation with microarray results (Rs = 0.874; Figure 4A). Second, expression levels for all ESTs examined showed generally opposite effects between JH and SHE treatments, i.e., when transcript levels were upregulated by JH they were downregulated by SHE, and vice-versa (Figure 4B-4G). With the exception of the DNA Replication Licensing Factor (FL644436) and Soldier Specific ProteinNtSp1-like” ESTs, expression levels in JH + SHE treatments were always intermediate between JH and SHE treatments (Figure 4B-4F). Also, the two 50MGP ESTs showed highly similar expression profiles (Figure 4B & C). However, only the two 50MGP ESTs showed significant variation among treatments (Kruskal-Wallis test, P < 0.05); nevertheless, we avoid over-interpreting this result due to our small sample sizes (3 biological replicates). The latter 50MGP results also provide additional support that these two ESTs represent portions of the same cDNA.

thumbnailFigure 4. qPCR reassessment for seven ESTs after treatments with JH, JH + SHE and SHE from the follow-up bioassay. A) Shows a positive correlation of qPCR results from follow-up bioassays against microarray fold change results for the six JH influenced ESTs (Spearman rank correlation, Rs = 0.89, df = 5, P = 0.02). The empty circle represents an EST that was downregulated by both JH and SHE. B-H) Show expression profiles for three JH upregulated ESTs (B, C and D), two JH downregulated ESTs (E, F and H), one SHE downregulated EST (G) and one SHE upregulated EST (H). The expression level is presented as 2-ΔΔCT, calculated from qPCR cycle threshold value. The microarray fold change values are provided in the parentheses for each EST. FL644436 (G) was downregulated by both JH (fold change: 0.543) and SHE (fold change: 0.811).

Discussion

Overview

This study used a microarray-based approach to compare the effects of a morphogenetic hormone (JH), soldier-derived primer pheromones (SHE), live reproductives (LR) and live soldiers (LS) on worker gut gene expression. We focused on worker termites because in lower termites like R. flavipes (i) workers compose >90% of colonies, (ii) their guts house both eukaryotic and prokaryotic symbionts, (iii) they are responsible for the majority of lignocellulose digestion, and (iv) workers are totipotent juveniles that retain the capacity to differentiate into both soldier and reproductive caste phenotypes [1,20,29,30]. Our central hypothesis is that the worker gut and its eukaryotic symbionts are potential molecular determinants of termite caste differentiation, eusocial polyphenism and, ultimately, social structure.

We found that variable numbers of expressed genes were affected among the four treatment categories (Table 1). JH treatment had the largest impact on gut gene expression (501 total ESTs affected), followed by LS (24 ESTs), LR (12 ESTs) and SHE (6 ESTs). The JH induced change is consistent with JH’s well-established ability to induce worker-to-soldier caste differentiation in lower termites [4], as well as presoldier induction assays in which 90-91% JH-induced presoldier differentiation was observed on a sub-sample of two colonies (see Results). Interestingly, the majority of JH upregulated genes were from the host termite, whereas the majority of JH downregulated genes were from protist symbionts. While the complement of upregulated host genes seems to provide insights into caste-regulatory physiology (discussed below), the down-regulation of symbiont genes by JH seems more likely an indicator of protist susceptibility to host hormones or a purging of symbionts in response to rising hormone titers [8,31]. The comparatively weak impact of SHE on gut gene expression is consistent with its lack of impacts on caste differentiation when applied alone [12,18]. Similarly, the greater impact on gene expression in LS treatments is consistent with the more pronounced impacts of live soldiers on limiting JH-dependent caste differentiation [7,18]. Finally, the unexpectedly small number of genes impacted in LR treatments is consistent with worker-derived apterous neotenic reproductives being relatively uncommon in Reticulitermes colonies [32], as well as mounting evidence of a genetic (rather than environmental) basis for reproductive differentiation in lower termites [33-36]. However, it may also be possible the 12 genes identified in LR treatments are mediating volatile primer pheromone signals coming from Reticulitermes neotenic females [37].

As shown by volcano plots (Figure 1), the SHE, LR and LS treatments had a number of upregulated array positions with large magnitudes of change but a lack of statistical significance. These profiles imply physiological heterogeneity among the test worker populations for individuals that are responsive to SHE, reproductive and soldier-based cues, and suggest a need for targeted research to better understand physiological variability among individuals in termite colonies. Also, as a contrast to directly comparing numbers of differentially expressed genes among the four treatment categories, GO analyses were performed to gain insights based on predicted cellular location, biological process and molecular function of differentially expressed genes. In particular, GO-Molecular Function analyses revealed some notable trends (Table 4). First, despite having only two passing genes, the SHE downregulated category had 19 GO-Molecular Function terms suggesting broad pleiotropic impacts by these two passing genes. Also, the numbers of GO-Molecular Function terms in the SHE downregulated and LS upregulated categories are consistent with previous bioassay results showing contrasting impacts by these two treatments (i.e., SHE + JH induces higher levels of caste differentiation, and live soldiers limit JH-dependent caste differentiation [7,12,18]).

Regarding the small expression fold changes obtained in the current study, we conclude this to be a result of the microarray-based platform that was used. This conclusion is based on the larger fold-change expression magnitudes obtained using qPCR relative to the microarray platform in the current study (Tables 2 and 3, Figure 4), as well as two previous studies that compared feeding impacts on gut gene expression using microarrays and quantitative pyrosequencing (showing higher expression by qPCR and pyrosequencing [26,39]). In the sections that follow, we further discuss candidate gene trends, related genes passing in multiple treatment categories, insights into termite social regulation and symbiosis, and overall conclusions.

Candidate gene trends

The most highly JH up- and downregulated ESTs encoded homologs of a host 50 kDa midgut protein50MGP” and a protist symbiont cysteine synthase a. The cysteine synthase a transcript was also the most highly abundant transcript identified previously with paper (cellulose) feeding [26], which was the control condition (with acetone) in the current study. Cysteine synthases catalyze production of acetate, which is an important metabolic intermediate in the termite gut [38]. The 50MGP gene is novel to the current study and discussed under Conclusions.

A majority of JH downregulated ESTs were of symbiont origin and are considered indicators of general symbiont decline in response to JH treatment. Most notably, the list of downregulated symbiont genes includes five GHF7 cellulases. Because GHF7 cellulases play important roles in R. flavipes digestion [26,39], this trend suggests compromised digestive capabilities in association with JH-induced caste morphogenesis. The remaining JH downregulated symbiont ESTs are not discussed hereafter but can be found in Table 3.

The second and third most highly upregulated ESTs in the JH dataset were of host origin and included a nli interacting factor-likephosphatase with predicted protein dephosphorylation function and an Arylsulfatase B homolog with predicted functions in cleaving sulfate groups from glycosaminoglycans, including N-acetyl monosaccharides that compose gut integument [40]. Other noteworthy genes in the JH upregulated dataset include two Apolipoproteins and three chymotrypsin/serine proteases with potential roles in lipid- and protease-based developmental signaling similar to JH responsive proteases identified previously [41,42]. Additionally, two cytochromeP450 homologs from families 6 and 4 were upregulated in the JH dataset and five others from the same families were downregulated. Interestingly, the cytochrome P450 Cyp15F1, which has been implicated in JH-dependent presoldier induction in R. flavipes and is inducible in the gut by wood feeding [6,26], was not impacted by JH treatment in the current study.

ESTs related to cuticle formation were also upregulated with JH treatment, including larval and pupal cuticle proteins[21,43], resilin[44], and enzymes potentially regulating melanin biosynthesis such as Tyramine beta hydroxylase and Dopamine N-acetyl transferase[45,46]. Possibly related to the high number of phosphorylation sites predicted on the translated 50MGP protein, the JH dataset also contained a number of upregulated ESTs with links to phosphate-related post-translational modification; for example, (i) the nli interacting factor-likephosphatase noted above, (ii) other phosphatases potentially involved in de-phosphorylation, (iii) kinases potentially involved in protein phosphorylation, and (iv) an insulin receptor homolog [47] with predicted kinase activity.

The most highly SHE upregulated EST was a host gene with no significant GenBank nr database match; however, translated searches of the GenBank EST database revealed numerous un-annotated homologs of this EST in other termite and cockroach cDNA libraries. The most SHE downregulated EST was a symbiont serine/threonine-protein kinase mph1 homolog related to mitotic function [48]. Also, an important DNA replication factor, DNA replication licensing factor mcm7[49], was downregulated by SHE treatment. These two SHE downregulated genes had a large number GO-Molecular Function terms associated with them, suggesting broad pleiotropic impacts.

The two most highly LR upregulated ESTs were from the host termite; one with no database matches and the other was a serine protease 13 homolog. The most downregulated EST in the LR dataset encoded a symbiont Linker histone h1 and h5 family protein involved in DNA binding and regulation of chromatin structure [50].

Finally, the most highly LS upregulated EST was a venom allergen 3-like homolog, similar to a venom protease from the ant Solenopsis invicta[51]. Other highly upregulated host ESTs from LS arrays had links to carbohydrate hydrolysis/binding and immune function. These ESTs included alpha amylase homologs, lysozyme homologs from GHF 13 and 22, and a C-type lectin homolog from CBM Family 13 [52]. The most downregulated ESTs in the LS dataset included one host and one symbiont EST, neither of which had significant database matches.

Similar ESTs passing in multiple microarray treatment categories

Kinases were a highly represented gene family in the JH dataset and other kinases also appeared in the LS, LR and SHE datasets. Of these kinases, only a serine/threonine-protein kinase mph1 homolog appeared in more than one dataset (JH and SHE), and was downregulated in both cases. An Arylsulfatase B homolog was highly induced in the JH dataset, and was also upregulated in the LS dataset. A c-type lectin precursor presumably involved in carbohydrate binding was upregulated in the JH and LS datasets. A number of protease, peptidase and/or chymotrypsin homologs were differentially expressed among the JH, LR and LS datasets; however, while each of these datasets had highly differentially expressed protease-related ESTs, the composition among each dataset was unique. Nonetheless, it is noteworthy that JH-responsive proteases have previously been identified in insect guts [41]. Lastly, host alpha amylase homologs were upregulated in the LS and LR datasets, but all were distinct. Thus, in summary, major functional protein categories sampled across treatment categories include kinases, proteases, and various carbohydrate active proteins (amylases, arylsulfatases, and lectins).

New insights into termite social regulation and symbiosis

Although JH and SHE may not be ingested by termites in nature, the simple technique of placing them on JH or SHE-treated diet ensured intake of JH and SHE by the termites. We expect that JH and SHE are ingested and absorbed by the cuticles of the termites exposed to treated filter paper. Also, because trophallaxis plays such a large role in termite sociality [1], this study in part, investigates what could be happening if termites were to acquire JH by trophallaxis [53]. We found significant JH-dependent expression changes in gut-associated genes from both the host termite and protist symbionts, as well as significant long-term presoldier induction in JH bioassays conducted on a sub-sample of colonies used in the current study (see Results). These correlated gene expression and phenotypic changes were elicited after just one day of treatment, illustrating the strong impacts of JH even after a short exposure. JH is a known caste-regulatory hormone, and termites are known to go through morphological changes after experimental JH exposure [4]; however, the mechanisms by which JH influences gene expression in termites are not well understood. This study takes us one step further in this interpretation. We have shown that JH generally up-regulates host gene expression and generally down-regulates that of protist symbionts. We do not expect that JH directly influences protist gene expression but it can influence the purging of gut contents [8,31], and our results might be indicative of a general decline of protist populations (dissections, performed post-hoc, revealed that some JH treated termites indeed lacked any gut content or had purged guts; results not shown). However, the changes brought upon by the absence of protists in the termite gut also have to be considered. The most significant effect in this respect would likely be on digestion, implying that an absence of regular nutrients could play a role in caste regulation [4,11,54]; for example, several protist GHF7 cellulases were downregulated by JH exposure (Table 3). It should also be noted that our cDNA libraries were made with only gut-associated ESTs; therefore the differentially expressed genes shown here must be only a fraction of the changes triggered by JH treatment in the whole body of termites.

We also found that exposure to SHE and LS triggered expression changes for comparatively small numbers of genes relative to JH. Previous studies showed that SHE and LS by themselves do not elicit any significant impacts on caste differentiation [7,12,18], and the current results are consistent with those findings. Also, relatively small numbers of symbiont genes were downregulated by SHE and LS treatments, supporting the idea that protist populations are not dramatically affected by these factors and thus do not mediate their signals.

Conclusions

This research took a microarray-based approach to test the hypothesis that the worker termite gut and its eukaryotic protist symbionts are potential molecular determinants of caste differentiation. This study builds on a prior study by Tarver et al. [18] that investigated whole-body expression of 49 host genes in response to JH, SHE and LS treatments. Four treatments were investigated in the current study that included JH, SHE, LR and LS. With respect to JH, our hypothesis is strongly supported. However, for SHE, LR and LS treatments our hypothesis is not well supported. These results suggest that, if they have impacts at all, the SHE, LR and LS treatments might (i) be acting more substantially beyond 1-day exposure periods [18], (ii) be acting in other body regions or tissues than gut, and/or (iii) have minimal impacts on expression of genes that influence caste differentiation and caste homeostasis. Alternatively, as suggested by GO-Molecular Function analyses, some of the small numbers of gut genes impacted by SHE, LR and LS treatments may actually have broad pleiotropic impacts. Despite these vastly different impacts on gut gene expression among treatments, major functional categories of responsive genes sampled across treatment categories encoded integumental proteins, kinases, proteases, and various carbohydrate-active proteins (amylases, arylsulfatases, and lectins).

This study also revealed a novel 50MGP cDNA as the most highly JH-inducible transcript in the R. flavipes gut. Multiple ESTs for this gene assembled into an apparent full-length cDNA that shares many common features with other host termite cDNAs [28,39]. While the function of the translated 50MGP is unknown, it is predicted to have a large number of phosphorylation sites, a well-defined secretory signal sequence, and share a predicted JH-binding region with a homologous protein from the bark beetle Dendroctonus ponderosae[55]. Other known insect homologs of this protein are from the sand fly Phlebotomous papatasi[56] and the red flour beetle Tribolium castaneum[57]. The translated 50 kDa protein has predicted immune, stress response and signal transduction functions; future functional studies will investigate links to these processes and others.

Finally, previous studies testing the JH + SHE combination revealed unexpected increases in soldier caste differentiation relative to treatments with JH alone, whereas treatments with SHE alone had no such impacts [7,12,18]. In the current study, we conducted “follow-up” bioassays and qRT-PCR on 6 significant genes from JH and SHE microarrays to investigate possible interactions. Interestingly, in these treatments JH and SHE had opposite impacts, and in most cases the JH + SHE combination resulted in gene expression intermediate between JH and SHE alone. These results and those of Tarver et al. [18] suggest next-generation transcriptome sequencing as a potentially informative approach for investigating whole-body gene expression impacts by JH + SHE treatments, as well as individual SHE components. This study in general provides important new insights into molecular determinants underlying termite caste polyphenism and homeostasis, including symbiont population decline in association with the caste differentiation process.

Methods

Colony maintenance and bioassays

All colonies were verified as R. flavipes by mitochondrial 16S rRNA sequencing as described by Szalanski et al. [58] (data not shown). Colonies were maintained in darkness in sealed plastic boxes containing wet pine wood shims and brown paper towels, within an environmental chamber kept at 22°C and 60-100% relative humidity (RH). Bioassays were conducted in darkness at 27°C and 60-100% RH. The termites used for microarray analysis originated from five established laboratory colonies at the University of Florida, Entomology and Nematology Department, in Gainesville, FL: (1) B1#1 (established 05/20/2009); (2) B2 (06/03/2010); (3) K2 (07/11/2007); (4) K5 (08/02/2008); and (5) K9 (06/29/2010). Bioassays for microarray studies were performed in August-September, 2010. Additionally, JH-presoldier induction bioassays were conducted on workers from one Florida colony used for microarray studies (B1#1) and another from Indiana (WI-9). These bioassays used 100 μg JH III (Sigma-Aldrich; Milwaukee, WI) per assay dish or acetone alone for controls and followed the methods of Tarver et al. [7,12], except only a 24-hr exposure was used. Six and four independent replicates per treatment were done for the WI-9 and B1#1 colonies. Finally, we conducted a second “follow-up” bioassay experiment to see the effect of JH, SHE and JH along with SHE on selected genes. Termites used in this bioassay originated from a single colony collected from Purdue University, West Lafayette, IN, USA (“WI-9” established in 2012).

Hormonal microarrays

The hormonal bioassays were conducted as previously described by Tarver et al.[12,18]. Twenty workers (pseudergates) were used per assay and given moist filter paper discs (Whatman #3 filter discs, GE Healthcare Bio-Sciences Corp., Piscataway, NJ) as food. The filter paper discs received either JH III or SHE or acetone (for control). JH III (93% purity; Sigma; St. Louis, MO) was applied on one filter paper disc per assay and at 150 μg per disc in 150 μl acetone. SHE was prepared by homogenizing soldier heads in acetone with a Tenbroeck glass homogenizer and applied at two soldier head equivalents per disc in 150 μl acetone. Acetone (150 μl) was applied on discs for control assays. After applying the solution, acetone was evaporated from the discs for 20 minutes at room temperature prior to each assay. The filter papers were then moistened with 150 μl of distilled water and placed in 3.5-cm diameter Petri dishes with the termites.

Social treatment microarrays

In social treatments, groups of twenty workers were maintained with either two soldiers or two neotenic reproductives originating from the same colony. Moist filter paper discs were offered as food. Termites were maintained for one day in 3.5-cm diameter Petri dishes in darkness at 27°C and 60-100% RH.

Gut extraction and RNA isolation

Both social and hormonal treatments were conducted for 24 hr, and the workers were then cold-immobilized, surface-sterilized by a serial rinse in 0.3% sodium hypochlorite (once) and sterilized water (twice), and dissected on Parafilm® to collect digestive tracts including salivary glands. Digestive tracts were transferred into RLA Lysis Buffer (Promega, Madison, WI) and stored at −70°C until RNA isolation. RNA Extraction and cDNA Synthesis was done according to Raychoudhury et al. [26].

Microarray hybridizations

Experiments were designed after MIAME guidelines, and microarray data obtained in these studies were deposited at ArrayExpress [http://www.ebi.ac.uk/arrayexpress webcite (accession number E-MTAB-1417)]. A type II microarray [59] design used with a common-reference strategy [60]. The common reference consisted of a normalized blend of all RNA samples included in the experiment. This common reference was co-hybridized against each replicate sample on single microarrays. Dye swaps [59] were performed between replicate samples and references to check for potential dye impacts on spot intensity. Twenty-five total microarray hybridizations were performed and consisted of each of five colonies treated with JH, SHE, acetone and exposed to live reproductives (LR) and live soldiers (LS).

Statistical analyses

The Matlab statistics toolbox was used for statistical analysis of the intensity data of 25 hybridizations from five different treatments [JH, SHE, exposure to live reproductives (LR), live soldiers (LS) and acetone (A) for control]. Before comparative analysis, the individual signal intensity values obtained from the microarray probes were log-transformed (using 2 as the base) and normalized among all individual samples included in the study. Normalization was accomplished by scaling the individual log-transformed signal intensities so that each dataset had comparable lower, median and upper quartile values [61]. After the data were normalized, Student’s t-tests were used to make probe-by-probe comparisons among each treatment and control (JH vs. A, SHE vs. A, LR vs. A and LS vs. A). In each comparison, a p-value and fold change was computed for all microarray loci. In addition to p-values, q-values were computed [62]. While the p-value measures the minimum statistical false positive rate incurred when setting a threshold for test significance, the q-value measures the minimum false discovery rate incurred when calling that test significant [62]. A volcano plot for each comparison was generated that displays the negative log10-transformed p-value versus log2-transformed fold change for each array locus (Figure 1).

Bioinformatic analyses

Contig generation: All significantly differentially expressed array positions that met the fold change criteria in each bioassay were selected and processed through Sequencher (Gene Codes Corporation, Ann Arbor, MI) with a minimum match percentage of 95 to generate contigs. The generated contigs and the remaining orphan sequences were used for further analyses.

Gene Ontology analysis using BLAST2GO: The selected contigs and the orphan sequences were analyzed using the program BLAST2GO [63] for identification and annotation. By using the inbuilt BLASTx algorithm, these sequences were used as queries in BLASTx searches against the GenBank non-redundant (nr) database with an e-value cut-off of ≤ 1e-03 (last performed 06, 2012). The putative identification, annotation, and Gene Ontology (GO) terms [64] for the sequences were also obtained through BLAST2GO. JH, SHE, LR and LS influenced sequences are listed in supporting information Additional file 1: Tables S1-S4.

Validation of microarray fold change data by quantitative real-time PCR

The fold change data from the microarray results were validated by performing sets of quantitative real-time PCRs (qRT-PCR) with a CFX-96 Real-time System (Bio-Rad, Hercules, CA) using the SYBR green detection method (SensiMix SYBR & Fluorescein one-step PCR reagent; Bioline, Taunton, MA). Fifty-two different sequences (Additional file 3: Table S6) with varying degrees of fold change were used to design primer sequences using the web-based tool Real-time Design (http://www.biosearchtech.com/realtimedesign webcite). The housekeeping gene lim-1 was used as a reference gene [18,26]. Two μl of total RNA (from aliquots of 10 ng/μl) were taken from the original mRNA pools used for microarray hybridizations from all five colonies (5 treatments each) to synthesize cDNA using the iScript cDNA kit (Bio-Rad, Hercules, CA). Triplicate qRT-PCR reactions were performed for each of the 25 sets of cDNA samples, along with a no-cDNA negative control, across the 52 primer sets (Additional file 3: Table S6). Cycling conditions were an initial step of 95°C for 3 minutes followed by 39 cycles of 95°C for 20 seconds, 56°C for 45 seconds and 68°C for 50 seconds. Quantification was performed by first generating a standard curve of primer amplification efficiency, using whole-gut cDNA from colony B1#1 with a five-fold dilution series, and then extrapolating the experimental samples onto the curve. Each triplicate sample was averaged to one data point for ease of graphical representation. The mean delta threshold cycle (ΔCT) was calculated for each data point by subtracting it from the average CT values of lim-1. Then a ΔΔCT value was calculated by subtracting average acetone data point from JH, SHE, LR or LS (see formula below for JH). These ΔΔCT values were plotted against the corresponding fold change levels from the microarray studies, and a correlation test was conducted.

<a onClick="popup('http://www.biomedcentral.com/1471-2164/14/491/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2164/14/491/mathml/M1">View MathML</a>

j = number of biological replicates,

i = number of technical replicates,

P = given primer, lim1 = lim1 primer;

JHi = CT value of the ith technical replicate from the JH treated termite gut cDNA,

Ai = CT value of the ith technical replicate from the acetone treated termite gut cDNA

“Follow-up” bioassays

Follow-up bioassays were conducted with JH III and freshly made SHE. The SHE was analyzed on HPLC to confirm whether the chemical peaks, and thus overall composition, matched with previous studies of Tarver et al.[7]. As Additional file 2: Figure S2 shows, the composition of SHE matched with that used in prior work demonstrating SHE transfer and efficacy by Tarver et al.[7]. Filter paper discs were soaked with one of the following: 200 μl acetone (control), JH III (65% purity; Santa Cruz Biotechnology Inc, Santa Cruz, CA) at 100 μg/μl of acetone, SHE (2 head equivalents in 100 μl acetone) and a combination of JH III (100 μg/100 μl of acetone) and SHE (2 head equivalent in 100 μl acetone). Then all filter papers were air dried and moistened with 150 μl of distilled water before offering to the termites in 6.5-cm diameter Petri dishes. Twenty workers were used per assay with three replicates for each treatment. Individual guts were extracted into RNA lysis buffer, and RNA was extracted as described above. cDNAs were produced using the iScript cDNA synthesis kit (Bio-Rad, Hercules, CA). These cDNAs were used for qPCR with lim1 and seven other primers that were selected due to very high or very low fold change ratios from JH and SHE microarray data (Figure 4). The ΔCT values (for control and treatment) acquired from the qPCR were used to calculate the 2-ΔΔCT values to compare the results for different treatments.

50 kDa Midgut protein sequencing and sequence analysis

The full-length cDNA sequence for the 50MGP gene was assembled from five overlapping clones/ESTs from a previously described gut cDNA library (GenBank accession nos. FL639806, FL636982, FL638011, FL638525 and FL637656 from Tartar et al.[28]). The contig sequence (GenBank no. KC751537) was verified by alignments having >99% identity with seven ESTs obtained from different R. flavipes phenotype cDNA libraries [65]. Portions of the ORF sequence were further verified by PCR amplification of cDNA fragments spanning nucleotides 41–537, 366–1104, 510–807 and 656–1126, using the forward and reverse primers described in Additional file 3: Table S6. Database searches were performed at NCBI (http://www.ncbi.nlm.nih.gov/ webcite) using BLASTn and BLASTx under default settings. Peptide translations were made using SDSC Biology Workbench (http://workbench.sdsc.edu/ webcite). Signal peptide, N-glycosylation, phosphorylation, protein functional category and enzyme class analyses were performed using the SignalP, NetNGlyc 1.0, NetPhos and ProtFun servers available at http://www.cbs.dtu.dk/services/ webcite. Protein alignments were generated using ClustalW in the Lasergene software package (DNAstar; Madison, WI).

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

RS and RR did the bioinformatic analyses, ran qRT-PCRs, surveyed the library fidelity, did statistical analyses and drafted the paper. MES conceived, designed and conducted the experiments and participated in drafting the manuscript. VUL and DGB participated in the original design of the experiment and conducted bioassays and RNA isolation for the microarray. YC and YS analysed the microarray data. All authors read and approved the final manuscript.

Acknowledgements

We thank: Jesse Hoteling for assistance with termite rearing; Ameya Gondhalekar for HPLC analysis of SHE; and Amit Sethi, Andres Sandoval, Brittany Peterson and Zachary Karl for helpful discussions. This research was supported by USDA-NIFA-AFRI grants 2009–05245 and 2010-65106-30727 to MES and DGB, and by the O.W. Rollins/Orkin Endowment at Purdue University.

References

  1. Wilson EO: The Insect Societies. Cambridge: Belknap Press of Harvard University Press; 1971. OpenURL

  2. Hartfelder K: Insect juvenile hormone: from "status quo" to high society.

    Braz J Med Biol Res 2000, 33:157-177. PubMed Abstract | Publisher Full Text OpenURL

  3. Smith CR, Toth AL, Suarez AV, Robinson GE: Genetic and genomic analyses of the division of labour in insect societies.

    Nat Rev Genet 2008, 9:735-748. PubMed Abstract | Publisher Full Text OpenURL

  4. Miura T, Scharf ME: Molecular basis underlying caste differentiation in termites. In Biology of Termites: A Modern Synthesis. Edited by Bignell DE, Roisin Y, Lo N. New York: Springer; 2011:211-253. OpenURL

  5. Zhou X, Song C, Grzymala TL, Oi FM, Scharf ME: Juvenile hormone and colony conditions differentially influence cytochrome P450 gene expression in the termite Reticulitermes flavipes.

    Insect Mol Bio 2006, 15:749-761. Publisher Full Text OpenURL

  6. Tarver MR, Coy MR, Scharf ME: Cyp15F1 - a novel cytochrome P450 gene linked to juvenile hormone-dependent caste differention in the termite Reticulitermes flavipes.

    Arch Insect Biochem Physiol 2012, 80:92-108. PubMed Abstract | Publisher Full Text OpenURL

  7. Tarver MR, Schmelz EA, Scharf ME: Soldier caste influences on candidate primer pheromone levels and juvenile hormone-dependent caste differentiation in workers of the termite Reticulitermes flavipes.

    J Insect Physiol 2011, 57:771-777. PubMed Abstract | Publisher Full Text OpenURL

  8. Howard RW, Haverty MI: Comparison of feeding substrates for evaluating effects of insect growth-regulators on subterranean termites.

    J Georgia Entomol Soc 1979, 14:3-7. OpenURL

  9. Hrdy I, Krecek J: Development of superfluous soldier induced by juvenile-hormone analogs in termite, Reticulitermes lucifugus santonensis.

    Insectes Soc 1972, 19:105-109. Publisher Full Text OpenURL

  10. Scharf ME, Ratliff CR, Hoteling JT, Pittendrigh BR, Bennett GW: Caste differentiation responses of two sympatric Reticulitermes termite species to juvenile hormone homologs and synthetic juvenoids in two laboratory assays.

    Insectes Soc 2003, 50:346-354. Publisher Full Text OpenURL

  11. Scharf ME, Buckspan CE, Grzymala TL, Zhou X: Regulation of polyphenic caste differentiation in the termite Reticulitermes flavipes by interaction of intrinsic and extrinsic factors.

    J Exp Biol 2007, 210:4390-4398. PubMed Abstract | Publisher Full Text OpenURL

  12. Tarver MR, Schmelz EA, Rocca JR, Scharf ME: Effects of soldier-derived terpenes on soldier caste differentiation in the termite Reticulitermes flavipes.

    J Chem Ecol 2009, 35:256-264. PubMed Abstract | Publisher Full Text OpenURL

  13. Elliott KL, Stay B: 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 2008, 54:492-500. PubMed Abstract | Publisher Full Text OpenURL

  14. Park YI, Raina AK: Juvenile hormone III titers and regulation of soldier caste in Coptotermes formosanus.

    J Insect Physiol 2004, 50:561-566. PubMed Abstract | Publisher Full Text OpenURL

  15. Park YI, Raina AK: Regulation of juvenile hormone titers by soldiers in the Formosan subterranean termite, Coptotermes formosanus.

    J Insect Physiol 2005, 51:385-391. PubMed Abstract | Publisher Full Text OpenURL

  16. Mao LX, Henderson G, Liu YX, Laine RA: Formosan subterranean termite soldiers regulate juvenile hormone levels and caste differentiation in workers.

    ESA Annals 2005, 98:340-345. OpenURL

  17. Wiltz BA, Henderson G, Chen J: Effect of naphthalene, butylated hydroxytoluene, dioctyl phthalate, and adipic dioctyl ester, chemicals found in the nests of the Formosan subterranean termite on a saprophytic Mucor sp. , Zygomycetes: Mucorales.

    Environ Entomol 1998, 27:936-940. OpenURL

  18. Tarver MR, Zhou X, Scharf ME: Socio-environmental and endocrine influences on developmental and caste-regulatory gene expression in the eusocial termite Reticulitermes flavipes.

    BMC Mol Biol 2010, 11:28. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  19. Zhou X, Tarver MR, Bennett GW, Oi FM, Scharf ME: Two hexamerin genes from the termite Reticulitermes flavipes: Sequence, expression, and proposed functions in caste regulation.

    Gene 2006, 376:47-58. PubMed Abstract | Publisher Full Text OpenURL

  20. Zhou XG, Oi FM, Scharf ME: Social exploitation of hexamerin: RNAi reveals a major caste-regulatory factor in termites.

    Proc Nat Acad Sci USA 2006, 103:4499-4504. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Zhou X, Tarver MR, Scharf ME: Hexamerin-based regulation of juvenile hormone-dependent gene expression underlies phenotypic plasticity in a social insect.

    Development 2007, 134:601-610. PubMed Abstract | Publisher Full Text OpenURL

  22. Scharf ME, Ratliff CR, Wu-Scharf D, Zhou X, Pittendrigh BR, Bennett GW: Effects of juvenile hormone III on Reticulitermes flavipes: changes in hemolymph protein composition and gene expression.

    Insect Biochem Mol Bio 2005, 35:207-215. Publisher Full Text OpenURL

  23. Boucias DG, Cai Y, Sun Y, Lietze VU, Sen R, Raychoudhury R, Scharf ME: The hindgut lumen prokaryotic microbiota of the termite Reticulitermes flavipes and its responses to dietary lignocellulose composition.

    Mol Ecol 2013, 22:1836-1853. PubMed Abstract | Publisher Full Text OpenURL

  24. Lewis JL, Forschler BT: Protist communities from four castes and three species of Reticulitermes.

    ESA Annals 2004, 97:1242-1251. OpenURL

  25. Scharf ME, Karl ZJ, Sethi A, Sen R, Raychoudhury R, Boucias DG: Defining host-symbiont collaboration in termite lignocellulose digestion: "The view from the tip of the iceberg".

    Comm Int Biol 2011, 4:761-763. OpenURL

  26. Raychoudhury R, Sen R, Cai Y, Sun Y, Lietze VU, Boucias DG, Scharf ME: Comparative metatranscriptomic signatures of wood and paper feeding in the gut of the termite Reticulitermes flavipes.

    Insect Mol Bio 2013, 22:155-171. Publisher Full Text OpenURL

  27. Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.

    Nature Proto 2009, 4:44-57. OpenURL

  28. Tartar A, Wheeler MM, Zhou X, Coy MR, Boucias DG, Scharf ME: Parallel metatranscriptome analyses of host and symbiont gene expression in the gut of the termite Reticulitermes flavipes.

    Biotech Biofuels 2009, 2:25. BioMed Central Full Text OpenURL

  29. Cleveland LR: Symbiosis between termites and their intestinal protozoa.

    Proc Natl Acad Sci USA 1923, 9:424-428. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  30. Scharf ME, Tartar A: Termite digestomes as sources for novel lignocellulases.

    BioFPR 2008, 2:540-552. OpenURL

  31. Brugerolle G, Radek R: Symbiotic protozoa of termites. In Soil Biology vol 6. Edited by Varma A, Konig H. Berlin: Springer-Verlag; 2006:243-269. OpenURL

  32. Thorne BL: Termite terminology.

    Sociobiology 1996, 28:253-263. OpenURL

  33. Hayashi Y, Lo N, Miyata H, Kitade O: Sex-linked genetic influence on caste determination in a termite.

    Science 2007, 318:985-987. PubMed Abstract | Publisher Full Text OpenURL

  34. Lo N, Hayashi Y, Kitade O: Should environmental caste determination be assumed for termites?

    Am Nat 2009, 173:848-853. PubMed Abstract | Publisher Full Text OpenURL

  35. Matsuura K, Vargo EL, Kawatsu K, Labadie PE, Nakano H, Yashiro T, Tsuji K: Queen succession through asexual reproduction in termites.

    Science 2009, 323:1687. PubMed Abstract | Publisher Full Text OpenURL

  36. Vargo EL, Labadie PE, Matsuura K: Asexual queen succession in the subterranean termite Reticulitermes virginicus.

    Proc Royal Soc B 2012, 279:813-819. Publisher Full Text OpenURL

  37. Matsuura K, Himuro C, Yoko T, Yamamoto Y, Vargo EL, Keller L: Identification of a pheromone regulating caste differentiation in termites.

    Proc Nat Acad Sci USA 2010, 107:12963-12968. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  38. O'Brien RW, Breznak JA: Enzymes of acetate and glucose-metabolism in termites.

    Insect Biochem 1984, 14:639-643. Publisher Full Text OpenURL

  39. Sethi A, Slack JM, Kovaleva ES, Buchman GW, Scharf ME: Lignin-associated metagene expression in a lignocellulose-digesting termite.

    Insect Biochem Mol Bio 2013, 43:91-101. Publisher Full Text OpenURL

  40. Hook M, Kjellen L, Johansson S, Robinson J: Cell-surface glycosaminoglycans.

    Ann Rev Biochem 1984, 53:847-869. PubMed Abstract | Publisher Full Text OpenURL

  41. Bian G, Ralkhel AS, Zhu J: Characterization of a juvenile hormone-regulated chymotrypsin-like serine protease gene in Aedes aegypti mosquito.

    Insect Biochem Mol Bio 2008, 38:190-200. Publisher Full Text OpenURL

  42. Malik ZA, Amir S, Venekei I: Serine proteinase-like activity in Apolipophorin III from the hemolymph of desert locust, Schistocerca gregaria.

    Arch Insect Biochem Physiol 2012, 80:26-41. PubMed Abstract | Publisher Full Text OpenURL

  43. Willis JH: Structural cuticular proteins from arthropods: Annotation, nomenclature, and sequence characteristics in the genomics era.

    Insect Biochem Mol Bio 2010, 40:189-204. Publisher Full Text OpenURL

  44. Lyons RE, Wong DC, Kim M, Lekieffre N, Huson MG, Vuocolo T, Merritt DJ, Nairn KM, Dudek DM, Colgrave ML, Elvin CM: Molecular and functional characterisation of resilin across three insect orders.

    Insect Biochem Mol Bio 2011, 41:881-890. Publisher Full Text OpenURL

  45. Sugumaran M, Semensi V, Dali H, Nellaiappan K: Oxidation of 3,4-Dihydroxybenzyl alcohol-A sclerotizing precursor of cockroach ootheca.

    Arch Insect Biochem Physiol 1991, 16:31-44. PubMed Abstract | Publisher Full Text OpenURL

  46. Andersen SO: Insect cuticular sclerotization: A review.

    Insect Biochem Mol Bio 2010, 40:166-178. Publisher Full Text OpenURL

  47. Emlen DJ, Warren IA, Johns A, Dworkin I, Lavine LC: A mechanism of extreme growth and reliable signaling in sexually selected ornaments and weapons.

    Science 2012, 337:860-864. PubMed Abstract | Publisher Full Text OpenURL

  48. Yamagishi Y, Yang CH, Tanno Y, Watanabe Y: MPS1/Mph1 phosphorylates the kinetochore protein KNL1/Spc7 to recruit SAC components.

    Nat Cell Biol 2012, 14:746-752. PubMed Abstract | Publisher Full Text OpenURL

  49. Li C, Jin J: DNA replication licensing control and rereplication prevention.

    Prot Cell 2010, 1:227-236. Publisher Full Text OpenURL

  50. Fan L, Roberts VA: Complex of linker histone H5 with the nucleosome and its implications for chromatin packing.

    Proc Nat Acad Sci USA 2006, 103:8384-8389. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  51. Hoffman DR: Fire ant venom allergy.

    Allergy 1995, 50:535-544. PubMed Abstract | Publisher Full Text OpenURL

  52. Cantarel BL, Coutinho PM, Rancurel C, Bernard T, Lombard V, Henrissat B: The carbohydrate-active enZymes database, CAZy,: an expert resource for glycogenomics.

    Nucl Acids Res 2009, 37:D233-D238. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  53. Henderson G, et al.: Primer pheromones and possible soldier caste influence on the evolution of sociality in lower termites. In Pheromone Communication in Social Insects. Edited by Vander Meer RK, Breed MD, Espelie KE, Winston ML. Boulder: Westview Press; 1998:314-329. OpenURL

  54. Nalepa CA: Nourishment and the evolution of termite sociality. In Nourishment and Evolution in Insect Societies. Edited by Hunt JH, Nalepa CA. Boulder: Westview Press; 1994:507-514. OpenURL

  55. Keeling CI, Henderson H, Li M, Yuen M, Clark EL, Fraser JD, Huber DP, Liao NY, Docking TR, Birol I, Chan SK, Taylor GA, Palmquist D, Jones SJ, Bohlmann J: Transcriptome and full-length cDNA resources for the mountain pine beetle, Dendroctonus ponderosae Hopkins, a major insect pest of pine forests.

    Insect Biochem Mol Bio 2012, 42:525-536. Publisher Full Text OpenURL

  56. Ramalho-Ortigão M, Jochim RC, Anderson JM, Lawyer PG, Pham VM, Kamhawi S, Valenzuela JG: Exploring the midgut transcriptome of Phlebotomus papatasi: comparative analysis of expression profiles of sugar-fed, blood-fed and Leishmania major-infected sandflies.

    BMC Genomics 2007, 8:300. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  57. Richards S, The Tribolium Genome Sequencing Consortium, et al.: The genome of the model beetle and pest Tribolium castaneum.

    Nature 2008, 452:949-955. PubMed Abstract | Publisher Full Text OpenURL

  58. Szalanski AL, Austin JW, Ovens CB: Identification of Reticulitermes spp. from south central United States by PCR-RFLP.

    J Econ Entomol 2003, 96:1514-1519. PubMed Abstract | Publisher Full Text OpenURL

  59. Gibson G, Muse SV: A primer of genome science. Sunderland: Sinauer; 2002. OpenURL

  60. Eisen MB, Brown PO: DNA arrays for analysis of gene expression.

    Meth Enzymol 1999, 303:179-205. PubMed Abstract | Publisher Full Text OpenURL

  61. Quackenbush J: Microarray data normalization and transformation.

    Nat Genet 2002, 32:496-501. PubMed Abstract | Publisher Full Text OpenURL

  62. Storey JD: A direct approach to false discovery rates.

    J Royal Stat Soc B 2002, 64:479-498. Publisher Full Text OpenURL

  63. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa A: High-throughput functional annotation and data mining with the Blast2GO suite.

    Nucl Acids Res 2008, 36:3420-3435. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  64. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene Ontology: tool for the unification of biology.

    Nat Genet 2000, 25:25-29. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  65. Steller MM, Kambhampati S, Caragea D: Comparative analysis of expressed sequence tags from three castes and two life stages of the termite Reticulitermes flavipes.

    BMC Genomics 2010, 11:463. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL