Transcriptomics on Social Interactions in Termites: Effects of Soldier Presence

The organization of social insect colonies requires sophisticated mechanisms to regulate caste composition according to colony demands. In termites, the soldier caste is responsible for the inhibition of soldier differentiation, but the mechanism underlying the regulation of soldier differentiation is still unclear. In this study, we performed transcriptome analyses to identify genes expressed in workers that fluctuated in the presence of soldiers in the subterranean termite Reticulitermes speratus. First, soldier differentiation was artificially induced via juvenile hormone (JH) application, and the inhibitory effects of soldier differentiation on soldier presence were evaluated. Second, transcriptomes were prepared from workers with or without soldiers under JH treatment, and expression analyses were performed to identify differentially expressed genes (DEGs) for each treatment. The expression levels of several DEGs were verified by quantitative real-time PCR. The results indicated that only a small number of DEGs were upregulated by the presence of soldiers. A homology search of DEGs and gene ontology (GO) analysis of the DEGs showed that some genes were responsible for the regulation of hormone levels, social interaction, and response to xenobiotic substances, suggesting that they could be involved in developmental arrest and pheromonal regulation in workers. Moreover, GO analysis indicated that the expression of many genes, including those involved in hormone metabolic processes, fluctuated with JH application. Suppression of soldier differentiation in the presence of soldiers could be accomplished by the expression of a large number of genes required for soldier differentiation.


INTRODUCTION
Social insects have multiple phenotypes (castes) in the same colony. The optimal composition of each caste and cooperation among individuals are needed for regular colony growth and maintenance (Wilson, 1971). Termites are major social insect groups, and their castes are normally discriminated as workers, soldiers, and reproductives (Roisin, 2000;Korb and Thorne, 2017). Termite caste differentiation is thought to be regulated by caste-specific gene expression through environmental stimuli that may trigger various hormonal conditions (Noirot, 1991;Watanabe et al., 2014). Because a distinctive sterile soldier caste is crucial for termite sociality and evolution, soldier differentiation has been extensively studied during the last two decades (Miura and Maekawa, 2020). Genes and cascades involved in specific weapon formations have been analyzed, and several important factors, including hormone-related and toolkit genes, have been clarified (e.g., Zhou et al., 2006;Toga et al., 2012;Masuoka et al., 2018;Sugime et al., 2019). The expression patterns of these factors can be affected by environmental stimuli via inter-individual (i.e., social) interactions (Watanabe et al., 2014). However, the effects of social interactions on gene expression and endocrine status have not yet been elucidated (Miura and Maekawa, 2020).
Termite soldiers are differentiated from workers via an intermediate stage called a presoldier. Presoldier and soldier formation is known to be inhibited by soldiers themselves and promoted by the reproductive caste in the same colony (Watanabe et al., 2014). Indeed, for the inhibitory mechanism by soldier existence, candidate primer pheromones [soldier inhibiting factors, SIFs (Park and Raina, 2005)], transferred from soldiers to workers, were identified in some species of Reticulitermes [γ-cadinenal and (-)-β-elemene; Tarver et al., 2009;Tarver et al., 2011;Mitaka et al., 2017]. These chemicals are terpenoids and are probably defense substances released from the frontal gland, which is developed in soldiers of phylogenetically apical termite lineages (Miura and Maekawa, 2020). Previous studies have also shown that live soldiers or whole extracts of soldiers had strong inhibitory effects on soldier differentiation (Tarver et al., 2011;Mitaka et al., 2017). Moreover, similar inhibitory effects were observed in the extracts of soldiers from more basal termite species (without frontal glands), which possess physical weapons such as enlarged mandibles and head capsule (Korb et al., 2003). Consequently, SIFs may consist of complex and lineage-specific compounds. However, because information about the internal changes of workers are lacking, it is unclear how the degree of inhibitory effects of SIFs is determined.
The percentage of soldiers is usually very small in natural colonies (Howard and Haverty, 1981), and thus soldier differentiation is not frequently observed compared to those of workers. However, soldier differentiation can be artificially induced by juvenile hormone (JH) or JH analog treatments with workers (Watanabe et al., 2014;Miura and Maekawa, 2020). In R. speratus, artificial presoldier induction rates were significantly decreased by the presence of soldiers compared to those without soldiers (Watanabe et al., 2011). Importantly, JH titers of workers with soldiers were significantly lower than those of workers without soldiers just 5 days after JH treatment. If workers are reared without soldiers, only small numbers of workers differentiate to presoldiers and soldiers (e.g., about 7% for 16 days in Coptotermes formosanus; Park and Raina, 2005). However, we never identify the soldier-destined workers before the presoldier molt normally in the mature colony. Consequently, to determine effectively whether gene expression changes are affected by the existence of soldiers (probably using SIFs), an artificial presoldier induction method is considered useful, and R. speratus workers should be investigated within 5 days after JH treatment.
Here, internal transcriptomic changes in JH-treated workers caused by coexisting soldiers were analyzed in R. speratus. The genome sequence and transcriptome of each caste/caste differentiation have been clarified for this species (Shigenobu et al., 2022;Saiki et al., submitted manuscript). Transcriptomes were prepared from worker individuals with or without soldiers within 5 days after JH treatment. Based on the differentially expressed genes (DEGs) observed and specific gene ontology (GO) terms detected, the molecular basis of physiological changes in workers with coexisting soldiers is discussed.

Termites
Mature termite colonies (total: four) were collected in Furudo, Toyama Prefecture, Japan, in 2014. The nest logs were transported to the laboratory and kept in a plastic case under constant darkness. Sixth and seventh instar workers were selected for the following experiments based on the number of antennal segments and body size (Tsunoda et al., 1986;Takematsu, 1992). All experimental insects were maintained in an incubator at 25 • C under constant darkness for at least 3 days before use.

Dish Assays for Induction of Presoldier Differentiation
Dish assays were performed in accordance with the procedure described previously (Tsuchiya et al., 2008;Watanabe et al., 2011;Masuoka et al., 2013). Briefly, filter papers (55 mm diameter; Advantec, Japan) were treated with 20 and 40 µg JH III (Santa Cruz Biotechnology, Dallas, TX, United States) dissolved in 200 µL acetone. Filter papers treated with acetone alone were used as controls. After the acetone evaporated, each filter paper was moistened with approximately 450 µL of distilled water and placed in a 65 mm Petri dish. Then, 20 workers were exposed to each filter paper along with 0 or 10 soldiers. Each category was replicated four times using individuals sampled from four different colonies. All dishes were kept in an incubator at 25 • C in constant darkness. The number of dead individuals and differentiated presoldiers was checked every 24 h. If a dead worker was detected, it was immediately removed from the dish. If a soldier died, a live soldier was added from a separate dish kept in the incubator. On day 16, the induction rates of newly molted presoldiers were compared between dishes with 0 or 10 soldiers. Presoldier differentiation rates (mean ± S.D. values) were calculated from dishes replicated four times (20 workers per dish) and evaluated by the generalized Wilcoxon test (80 individuals in each JH III concentration) using the statistical software Mac Statistical Analysis, version 1.5 (Esumi, Japan). Statistical significance was set at P < 0.05.

Dish Assays for RNA Extraction
Dish assays for RNA extraction were performed using the method described in Section "Dish Assays for Induction of Presoldier Differentiation." Three days after treatment with 40 µg JH III (or acetone treatment as control), all workers in each dish were fixed with liquid nitrogen and stored at −80 • C until RNA extraction.

RNA Extraction, Library Preparation, and Sequencing
Total RNA was extracted from four categories [(1) workers with 0 soldiers 3 days after acetone treatment, (2) workers with 10 soldiers 3 days after acetone treatment, (3) workers with 0 soldiers 3 days after JH treatment, (4) workers with 10 soldiers 3 days after JH treatment] using an SV Total RNA extraction kit (Promega, Madison, WI, United States). RNA extracted from 20 individuals without guts was used for each library. Four replicates derived from four different colonies (biological quadruplicates) were prepared (four libraries × four categories = 16 libraries). The amounts of total extracted RNA and DNA were quantified using a Qubit fluorometer (Life Technology, Eugene, OR, United States) and the quality was confirmed using an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, United States). Total RNA (500 ng) was used for cDNA synthesis and purification using a low-throughput protocol with a TruSeq Stranded RNA LT Kit (Illumina, San Diego, CA, United States). A halfscale reaction of the standard protocol was used for library preparation. The quality and quantity of cDNA were validated using an Agilent 2100 bioanalyzer and a KAPA qPCR SYBR Green PCR Kit (GeneWorks, Thebarton, SA, Australia). Sixteen libraries were sequenced by single-end sequencing (101 bp) using Hiseq1500 (Illumina, San Diego, CA, United States). All reads were deposited in the DDBJ Sequence Read Archive (DRA) database under accession number DRA013774.

Identification of Differentially Expressed Genes
For each read, nucleotides with a low-quality score at the sequence ends and adapter sequences were removed using SolexaQA ver. 2.5 (Cox et al., 2010), and the cutadapt program version 1.2.1 (Martin, 2011), respectively. These reads were mapped against the R. speratus genome (Rspe OGS1.0; Shigenobu et al., 2022) using the TopHat program version 2.0.13 (Kim et al., 2013) with default settings.
The counting of reads and detection of DEGs were performed using Cufflinks pipeline version 2.2.1 (Trapnell et al., 2010). In this program, the count data were normalized and analyzed using the same algorithm implemented in DESeq (Anders and Huber, 2010). These counts were scaled using the median of the geometric means of fragment counts across all libraries. After normalization, pairwise comparisons were performed using cuffdiff command with "-frag-bias-correct" and "-multi-readcorrect" options (Roberts et al., 2011;Trapnell et al., 2013). We used the following four schemes for comparison: (i) without JH and 0 soldiers vs. without JH and 10 soldiers, (ii) without JH and 0 soldiers vs. 40 µg JH and 10 soldiers, (iii) without JH and 10 soldiers vs. 40 µg JH and 10 soldiers, and (iv) 40 µg JH and 0 soldiers vs. 40 µg JH and 10 soldiers. Homology searches of DEGs obtained in (i) and (iv) (the same JH concentration and relatively small numbers of DEGs; see results) were performed by BLASTX using the NCBI non-redundant protein database (nr) (run on March 2022), and the top hit sequence was obtained for each DEG.
We also conducted a generalized linear model (GLM) analysis. Because Cufflinks pilelined did not generate raw count data, we re-mapped and counted our RNA-seq data. Briefly, we used the Bowtie program (Langmead et al., 2009) for mapping and RSEM v1.3.3 software (Li and Dewey, 2011) to estimate the relative abundance and expected read counts of all genes. The count data were normalized by the TMM (Trimmed Mean of M values) method in the edgeR software package (Robinson and Oshlack, 2010;McCarthy et al., 2012). The adjusted count values were then used for the DEG analysis. Two model designs were used for DEG detection. In the first design, we considered three explanatory factors: soldier presence, JH treatment, and these interactions (Soldier presence × JH treatment). In the second design, in addition to these three factors, we used colony information as another explanatory factor. The threshold of statistical analysis is false discovery rate (FDR) cutoff of 0.05.

Annotation and Gene Ontology Enrichment Analysis
Since GO terms were only assigned to model organisms, we identified the fruit fly ortholog of each termite gene. All predicted termite amino acid sequences were searched against fruit fly (Drosophila melanogaster) amino acid sequences using BLASTP. The complete amino acid sequence datasets for the fruit flies were downloaded from Flybase version FB2012_04. 1 BLASTP searches were performed using termite genes as queries, and a 10 −5 e-value cutoff was used against the fruit fly dataset. The top hit proteins were defined as orthologs of the focal termite gene.
Gene ontology enrichment analysis was performed using clusterProfiler software package (version 2.4.3, R version 3.3.3; Yu et al., 2012). We performed an enrichment test for GO terms by assuming a hypergeometric distribution. To prevent high FDR due to multiple tests, we also estimated the q-values for FDR control (Storey, 2002). For these analyses, we used a list of DEGs identified by pairwise comparisons. To identify enriched GO biological processes, we conducted an enrichment analysis of these DEGs, in which p < 0.1 and q < 0.05 were used as strict cutoff values.

Quantitative Real-Time PCR
To validate the transcriptome analysis, qPCR was performed for the identified DEGs (a total of three genes, see results). Individuals were sampled from five different colonies collected in Toyama Prefecture, Japan from 2019 to 2020. The dish assay was performed following the method described above, and the presoldier induction rates were confirmed 16 days after the treatment. Total RNA was extracted from workers with 0 or 10 soldiers 3 days after treatment with 40 µg JH III using ISOGEN II (Nippon Gene, Tokyo, Japan). A total of 17-20 live workers (whole bodies) were used in each sample and homogenized using a Bead Mill 4 (Thermo Fisher, Waltham, MA, United States). Five replicates derived from five different colonies (biological quintuplicates) were prepared for each category. The amounts of RNA and DNA were quantified using a Qubit fluorometer, and the purity and quantity of the extracted RNA were determined by spectroscopic measurements at 230, 260, and 280 nm using a NanoVue spectrophotometer (GE Healthcare Bio-Sciences, Tokyo, Japan). DNase-treated RNA (2 µg per sample) was transcribed using High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher). Quantitative PCR (200 nM of each primer) was performed in biological quintuplicates using PowerUP SYBR Green Master Mix (Thermo Fisher) and QuantStudio 3 Real-Time PCR System (Thermo Fisher). According to the previous study (Miyazaki et al., 2021), to determine an internal control gene, the suitability of six reference genes, EF1-alfa (accession no. AB602838), NADH-dh (AB602837), beta-actin (AB520714), GstD1 (gene ID: RS001168), EIF-1 (RS005199), and RPS18 (RS015150), were evaluated using GeNorm (Vandesompele et al., 2002) and NormFinder (Andersen et al., 2004). Specific primers were designed against each gene sequence using Primer3Plus (Untergasser et al., 2007;Supplementary Table 1). Statistical analysis, Welch's t-test, was performed using the statistical software Mac Statistical Analysis version 3.0 (Esumi, Tokyo, Japan).

Presoldier Induction Rates
As shown in a previous study (Watanabe et al., 2011), presoldier induction rates after both 20 and 40 µg JH III treatments were significantly reduced by soldier presence (Figure 1). The JH titer (endogenous + applied JH III) levels in workers with 0 soldiers were shown to be higher than those with 10 soldiers 5 days after JH III treatment (Watanabe et al., 2011(Watanabe et al., , 2014. Because we tried to effectively detect gene expression changes in workers before the reduction of JH III titer levels, total RNA was extracted from workers 3 days after 40 µg or without JH III treatments.

Numbers of Differentially Expressed Genes
We obtained different numbers of DEGs for the four pairwise comparisons (Figure 2). The number of DEGs between 40 µg and without JH III treatments in workers with 0 soldiers was 3,089. Of these genes, 1,516 and 1,573 were upregulated in 40 µg and without JH III-treated workers, respectively (Figure 2;  Supplementary Table 2). Similarly, the number of DEGs between 40 µg and without JH III treatments in workers with 10 soldiers was 3,014. Of these genes, 1,567 and 1,447 were upregulated in 40 µg and without JH III-treated workers, respectively (Figure 2;  Supplementary Table 3). In contrast, there were only 44 DEGs between soldier absence (26 upregulated genes) and soldier presence (18 upregulated genes) in acetone-treated workers (Figure 2; Supplementary Table 4). The expression levels of 171 genes were significantly different between soldier absence (142 upregulated genes) and soldier presence (29 upregulated genes) in JH III-treated workers (Figure 2; Supplementary  Table 5). In particular, the expression of a large number of genes (approximately 3,000) in workers was fluctuated by JH/acetone treatments. These numbers of DEGs were similar to, but slightly larger than, those between JH analog (methoprene) and acetone-treated workers in Coptotermes formosanus (2,547 unigenes) (Du et al., 2020). Reticulitermes and Coptotermes are phylogenetically closely related to each other (Bucek et al., 2019), and inferred genome sizes (800-900 Mb) and total gene numbers (approximately 13,000-15,000) of both species are very similar (Maekawa et al., 2022). Consequently, these discrepancies may be due to the differences in treated chemicals [JH homolog (this study), synthetic juvenoid (Du et al., 2020)] and/or RNA-sequencing (RNA-seq) methods [genome-based assay (this study), de novo transcriptome assembly-based assay (Du et al., 2020)]. Most of these DEGs were also identified by GLM analysis (Supplementary Tables 6, 7), suggesting that our analysis effectively listed genes which are related to soldier differentiation.

Quantitative Real-Time PCR
We selected three genes [cytochrome P450 (Cyp4c3; RS013835), alpha-amylase (Amy-p; RS006137), and multidrug resistance protein (Mdr49; RS002816)], which were listed as DEGs in JH III-treated workers with 0 and 10 soldiers (Supplementary Table 5). We focused on this category because it was expected that there would be large variations in gene expression levels among samples. RNA-seq analysis indicated that all these three genes were upregulated in workers with 10 soldiers. A suitable reference gene (NADH-dh) was selected for real-time qPCR analysis using GeNorm and NormFinder (Supplementary Table 8). Real-time qPCR analysis showed that high expression levels of the three genes examined were observed in JH III-treated workers with 10 soldiers (soldier presence), compared to those with 0 soldiers (soldier absence) (Figure 3). Although the statistical support of RS002816 expression levels between soldier presence and absence was weak (p = 0.106; Welch's t-test), it should be noted that the q value of RS002816 (0.012) was higher than that of the remaining two genes (0.007; Supplementary Table 5). Overall, the qPCR analysis supported the reliability of the RNA-seq results.

Differentially Expressed Genes Between Soldier Presence and Absence
We focused on DEGs specifically observed in workers with 0 and 10 soldiers under the same JH concentrations. Note that all DEGs discussed below were detected by both designes of GLM analysis (Supplementary Tables 6, 7), except that Cyp4c3 (RS013835) was not observed only in the first design. In the acetone treatment, three takeout protein genes (RS013762, RS014812, and RS010477) were highly expressed in workers with 0 soldiers, whereas a JH-inducible protein gene (RS007835) was Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 1 | Presoldier induction rates (mean ± S.D., biological quadruplicates) 16 days after the 20 or 40 µg juvenile hormone (JH) III application in Reticulitermes speratus. In both applications, presoldiers were more frequently differentiated from workers with 0 soldiers (soldier absence, "Sol = 0") than those with 10 soldiers (soldier presence, "Sol = 10") (generalized Wilcoxon test, *P < 0.05).
FIGURE 2 | Numbers of differentially expressed genes (DEGs, bold italic) between each category. Small italic numerals indicate the numbers of upregulated genes in each category. Gene numbers between juvenile hormone (JH) III and acetone treatments (3,014 in soldier presence, 3,089 in soldier absence) were much larger than those between soldier presence and absence (171 in JH treatments, 44 in acetone treatments). Numbers of significant gene ontology (GO) terms of the DEGs are indicated in squares. Table 4). Takeout proteins are normally able to bind JH because they possess the JH binding protein domain conserved in insects (Noriega et al., 2006;Chamseddin et al., 2012). Although the functions of JH-inducible and takeout proteins are unanalyzed, there is a possibility that these work on the role of JH sequestration from the hemolymph of workers, such as JH binding proteins (e.g., hexamerin, Zhou et al., 2006). Alternatively, these proteins may be crusial for the change of the JH sensitivity in workers, which is similar to the case in Pheidole ants (Wheeler and Nijhout, 1984;Nijhout, 2003). In any case, these results support that JH levels in workers are affected by the three-day interaction with soldiers. Note that other JH binding proteins (including hexamerin) and JH biosynthetic/degradation genes, all of which were annotated in previous literature (Shigenobu et al., 2022), were not detected in our transcriptome analysis.

highly expressed in workers with 10 soldiers (Supplementary
Two chemosensory protein genes (RS000584 and RS010442) were highly expressed in workers with 0 soldiers (Supplementary Table 4) and were recognized as RspeCSP1 and RspeCSP7, respectively (Shigenobu et al., 2022). Previous RNA-seq analysis in R. speratus indicated that both genes were highly expressed in the head compared to the remaining part of the body of workers (RspeCSP7) or workers and soldiers (RspeCSP1) (Shigenobu et al., 2022). Although their precise expression sites (antennae or other head parts) should be clarified, they may be involved in SIFrelated social communication between soldiers and workers.
In the JH III treatment, the functions of some DEGs (total: 52) could not be identified based on the BLASTX search (hypothetical or uncharacterized proteins, no hit; Supplementary Table 5). Most DEGs shown in the acetone treatment described above (Supplementary Table 4) were not detected, but many genes involved in antimicrobial and xenobiotic responses and digestive enzymes were observed in workers with 10 soldiers [e.g., Mdr49 (RS002816), prolixicin antimicrobial protein (AttD, RS000201), laccase2 (RS004166), and Amy-p (RS006137)] and those with 0 soldiers [e.g., C-type lysozyme-3 (RS003406), toll-like receptor 6 (RS012895), and venom protease-like (RS012253)]. Interestingly, fatty acyl-CoA reductase (RS002448), a well-known soldier-specific gene (Maekawa et al., 2022), was highly expressed in workers with 10 soldiers. As this gene may be involved in the production of soldier-specific cuticular hydrocarbon (CHC) profiles (Wu et al., 2018;Maekawa et al., 2022), the presence of soldiers can induce changes in worker CHC profiles. Cyp4c3 (RS013835) was also highly expressed in workers with 10 soldiers. Because CYP4 is involved in the last step of CHC biosynthesis in D. melanogaster (Qiu et al., 2012), it may also support the above hypothesis. The soldier ratio applied in this study was quite high, and further analysis should be performed in colonies with an appropriate soldier ratio under natural conditions (about 2% in R. flavipes; FIGURE 3 | Quantitative real-time PCR expression levels of the three genes (mean ± S.D., biological quintuplicates) in the 40 µg juvenile hormone (JH) III-treated workers with 10 soldiers (Soldier presence, "Sol = 10") and 0 soldiers (Soldier absence, "Sol = 0"). Each value was normalized to the expression levels of NADH-dh (Supplementary Table 8). Asterisks above the bars indicate significant differences (*P < 0.05, **P < 0.01; Welch's t-test). Howard and Haverty, 1981) to clarify the general tendency of this hypothesis.
Finally, in JH III treatment, high expression levels of some members of the multigene family [beta-glucosidase (RS004143) and lipocalin (RS013912); Shigenobu et al., 2022] were observed in workers with 10 soldiers (Supplementary Table 5). Many paralogs of both genes have been identified in the R. speratus genome, and these paralogs showed different expression patterns among castes (Shigenobu et al., 2022). Beta-glucosidase is essential for cellulose digestion in termite workers (Watanabe and Tokuda, 2010), but RS004143 was highly expressed in the thorax and abdomen of soldiers and reproductives (Shigenobu et al., 2022). Thus, RS004143 may have a different role other than wood digestion in R. speratus. It is interesting to note that betaglucosidase (called Neofem2) is involved in queen-recognition pheromones that probably function in the suppression of reproductive emergence in Cryptotermes secundus (Korb et al., 2009;Korb, 2016). Similarly, lipocalin is a member of the protein transporter family, but molecular phylogeny showed that RS013912 was closely related to soldier-specific protein 1 (SOL1) identified in Hodotermopsis sjostedti (Shigenobu et al., 2022). SOL1 may function as a signaling molecule for defensive social interactions among colony members (Miura et al., 1999;Miura, 2005). We suggest that both RS004143 and RS013912 are involved in chemical communication, and their expression changes in workers are affected by different social circumstances, with or without soldiers.

Gene Ontology Enrichment Analysis of Differentially Expressed Genes
We performed GO enrichment analysis of the DEGs observed in each category (Figure 2). The number of GO terms detected between 40 µg and without JH III treatments in workers with 0 soldiers was 65 (Supplementary Table 9). Similarly, the number between 40 µg and without JH III treatments in workers with 10 soldiers was 60 (Supplementary Table 10). More than half of these terms (total: 38) were common (bold italic terms in Supplementary Tables 9, 10), and four out of 38 terms (GO: 0016053, 0046394, 0072330, and 0032787; metabolic and biosynthetic processes of some molecules) were specifically observed during the workerpresoldier molt (Saiki et al., submitted manuscript). However, specific GO terms involved in hormone levels (e.g., GO: 0010817, 0042446, and 0045455) were detected only in workers with 10 soldiers (Supplementary Table 10). Moreover, the number of GO terms significantly detected between soldier presence and absence in the acetone-treated workers was only three (Supplementary Table 11); all of which were related to the regulation of hormone levels. These results clearly indicate that the JH levels in workers are affected by the presence of soldiers.
Finally, a total of 10 significant GO terms were observed between the presence and absence of soldiers in JH-treated workers (Supplementary Table 12). These terms included metabolic processes of some molecules (GO: 0006022, 1901071, and 0006040), chitin and cuticle development (GO: 0006030 and 0040003), and response to xenobiotic substances (GO: 0009617, 0042742, and 0050830). These results may be explained by the effect of the presence of soldiers on developmental arrest during JH-induced presoldier differentiation, which is generally accompanied by specific cuticle development via the tyrosine metabolic pathway (Masuoka et al., 2013;Masuoka and Maekawa, 2016).

CONCLUSION
Most workers treated with commercial JH III are differentiated into presoldiers, and living soldiers strongly inhibit the presoldier differentiation by rapid JH decrease soon after interaction with workers (Watanabe et al., 2011(Watanabe et al., , 2014. This study aimed to understand the gene expression profiles of workers affected by coexisting soldiers. RNA-seq analysis supported that worker JH levels are affected by the presence of soldiers, probably by the functions of JH binding and inducible proteins, regulatory factors of social interaction, and response to xenobiotic substances. Further gene function analyses of candidate targets are needed to clarify this possibility.