Comparison of gene expression profiles among caste differentiations in the termite Reticulitermes speratus

Termite castes express specialized phenotypes for their own tasks and are a good example of insect polyphenism. To understand the comprehensive gene expression profiles during caste differentiation, RNA-seq analysis based on the genome data was performed during the worker, presoldier, and nymphoid molts in Reticulitermes speratus. In this species, artificial induction methods for each molt have already been established, and the time scale has been clarified. Three different periods (before the gut purge (GP), during the GP, and after the molt) were discriminated in each molt, and two body parts (head and other body regions) were separately sampled. The results revealed that many differentially expressed genes (head: 2884, body: 2579) were identified in each molt. Based on the independent real-time quantitative PCR analysis, we confirmed the different expression patterns of seven out of eight genes in the presoldier molt. Based on the GO and KEGG enrichment analyses, the expressions of genes related to juvenile hormone titer changes (e.g., JH acid methyltransferase), nutrition status (e.g., Acyl-CoA Delta desaturase), and cell proliferation (e.g., insulin receptor), were shown to specifically fluctuate in each molt. These differences may have a crucial impact on caste differentiation. These data are important resources for future termite sociogenomics.

Termites, ants and bees constitute the major social insect groups. The division of labor, including reproduction, maintains the social insect colony and increases overall colony productivity 1 . Castes are morphologically and behaviorally specialized for specific tasks, and they are a good example of insect polyphenism in response to various environmental factors 2,3 . Termite castes are classified as workers, soldiers, or reproductives [4][5][6] . Workers feed and care for the brood and have highly developed digestive organs whereas soldiers are differentiated from workers by an intermediate stage (presoldier) and are specialized for colony defense with sclerotized heads and weapons (e.g., enlarged mandibles). Primary reproductives, who establish new colonies, and neotenic reproductives, who take over reproduction in their natal nests, make up the reproductive caste 6 . Both have functional gonads for reproduction. Caste differentiation to form different phenotypes is accomplished through the expression of different gene sets in response to environmental stimuli 7 .
The candidate gene approach has primarily been used in molecular analyses of termite caste differentiation 8,9 . These analyses were particularly focused on soldier differentiation because, in many species, juvenile hormone (or JH analog) treatments could be used to induce soldiers from workers 10,11 . Some important genes, such as hexamerin 8 , insulin receptor 12 , distal-less 9 and deformed 13,14 have been found to be specifically expressed at a certain

Materials and methods
Termites used for RNA-seq analysis. Termites were collected from five colonies (#A, B, C, D, and E) for transcriptome sequencing. All colonies were collected in Furudo, Toyama Prefecture, Japan, between June and September 2013. All colonies were maintained in plastic cases at 25 °C in constant darkness until the induction of the worker-worker, worker-presoldier, and nymph-nymphoid molts (see below). The 5th stage nymphs (N5 nymphs) and the 4th-5th stage workers (old-age workers) were collected from the colonies based on the shape of wing buds, body size, and antennal segments. Nymphs were distinct from workers by the possession of wing puds, and N5 nymphs and old-age workers had 16 and 15-17 antennal segments, respectively 20,28 .

Induction of worker-worker molt by 20E application.
To obtain individuals at multiple times during the worker-worker molt, we used the 20-hydroxyecdysone (20E) application method, as described previously 23 (Fig. 1A). Old-age workers were collected from each of three different colonies (#A, B, and C) and kept overnight amidst moist red-colored thin papers (Goshikizuru No.14, Goukaseishi Co., Ltd., Aichi, Japan). Gut-purged workers (GP workers) were identified as having yellowish-white abdomens (Fig. 1B). Non-GP workers were used for the following analyses. Red colored thin paper was treated with 40 µg 20E dissolved in 400 µL acetone and placed in a 65 mm Petri dish with 20 non-GP workers. Ten Petri dishes (containing 200 workers) from each colony were used to collect individuals for RNA extraction. 20E and control-acetone-treated Petri dishes were also used to measure the induction rate (n = 1-3 dishes, 20 individuals per dish). In addition, to obtain natural  www.nature.com/scientificreports/ GP workers, 200 workers were kept in a 90 mm Petri dish with no 20E or acetone treatments. Because of low frequency of natural molting events and limitation of the space in the incubator, we used a 90 mm Petri dish for this analysis. All Petri dishes were kept in an incubator at 25 °C. Over a 2-week period, all dishes were checked every 24 h for dead individuals. According to previous research 23 , GP workers appeared approximately 5 days after 20E treatments. During the worker-worker molt, four main periods were identified: (1) before the 20E treatment (worker), (2) just before the gut purge (pre-GP worker), (3) during the gut purge (GP worker), and (4) just after the molt (molt worker). For the pre-GP workers, we randomly collected 10 individuals from multiple dishes each day following the 20E treatment (total 4 days). GP workers, including natural GP workers, that emerged on the same day were transferred to a new 40 mm dish (up to 10 individuals). These workers (referred to as GP-0 workers) were kept in an incubator at 25 °C. Because the period for gut purging was approximately 5 days 23 , GP-0, -1, -2, -3, and -4 workers (i.e., GP workers 0, 1, 2, 3 and 4 days after the emergence, respectively) were collected, respectively (n = 10 per day). For the molt workers, individuals were collected within 24 h of the worker molt. For RNA extraction, 10 individuals were used per period. For the pre-GP workers, 2, 3, 3, and 2 individuals were randomly selected at 1-4 days after 20E treatment, respectively, and mixed for the same period. For the GP workers, GP-0, -1, -2, -3, and -4 workers were equally mixed (Fig. 1A). Each individual was dissected on ice, with the head and other parts of the body (including thorax and abdomen with guts; hereafter body) separated, immediately frozen in liquid nitrogen, and stored at -80 °C until RNA extraction.
Induction of presoldier differentiation by JH III application. Induction of presoldier differentiation was performed as previously described [23][24][25] . Old-age workers collected from each of the three different colonies (#A, B, C) were kept overnight with wetted red-colored thin paper. Non-GP workers were used for the following analyses ( Fig. 1). Red colored thin paper was treated with 80 µg JH III (Santa Cruz Biotechnology Inc., California, USA) dissolved in 400 µL acetone and placed in a 65 mm Petri dish with 20 non-GP workers. Ten Petri dishes (containing 200 workers) from each colony were used to collect individuals for RNA extraction. We also prepared JH and control-acetone-treated dishes to measure the induction rate (n = 1-3 dishes, 20 individuals per dish). All Petri dishes were kept in an incubator at 25 °C. Over a 2-week period, all dishes were checked every 24 h for dead individuals. As the method described for worker-worker molt ("Induction of worker-worker molt by 20E application" section), pre-GP workers (individuals collected at 1-4 days after JH treatments), GP workers (GP-0, -1, -2, -3, and -4 workers), and molt presoldiers (individuals collected within 24 h of the presoldier molt) were sampled. For RNA extraction, 10 individuals were used of each period. Each sample was dissected to separate the head and body on ice, immediately frozen in liquid nitrogen, and stored at -80 °C until RNA extraction.

Induction of nymphoid differentiation by isolation from the natal nest.
Nymphoid differentiation was induced by isolation from the natal nest, according to the method of Saiki and Maekawa (2011) 26 . The 5th stage nymphs (N5 nymphs) were collected from two different colonies (#D, E) and kept overnight with wetted red-colored thin paper. Non-GP nymphs were used for subsequent analyses (Fig. 1). Red thin paper was placed in a 65 mm Petri dish with 20 non-GP nymphs prepared from colony #D. We prepared > 11 dishes (> 220 nymphs) to collect individuals for RNA extraction. Since a large number of nymphs (> 600) were collected from colony #E, red-colored paper (90 mm) was placed in a 90 mm Petri dish with approximately 200 non-GP nymphs from this colony. We prepared three dishes (> 600 nymphs) to collect individuals for RNA extraction. Because of the sampling difficulties and large number of nymphs in colony #E, we doubled the number of Petri dishes for biological replications using nymphs from colony #E. All Petri dishes were kept in an incubator at 25 °C. Over a 2-week period, all dishes were checked every 24 h for dead individuals. As the method described for worker-worker molt ("Induction of worker-worker molt by 20E application" section), the N5 nymphs, the pre-GP nymphs (individuals collected at 1-4 days after isolation), GP nymphs (GP-0, -1, -2, -3, and -4 nymphs), and molt nymphoids (individuals collected within 24 h of the nymphoid molt) were used. For RNA extraction, 10 individuals were used in each period. Each sample was dissected on ice to separate the head and body, immediately frozen in liquid nitrogen, and stored at 80 °C until RNA extraction.
Total RNA extraction. We prepared 72 categories: 3 molts (worker-worker, worker-presoldier, nymphnymphoid) × 4 periods (workers or nymphs, pre-GP, GP, and molt individuals), × two body parts (head and body), × three colonies (in the case of nymphoid differentiation, triplicates were prepared from two colonies; see "Induction of nymphoid differentiation by isolation from the natal nest" section) (10 individuals were used per category). Total RNA was isolated from each category using the SV Total RNA Isolation System (Promega, Madison, WI, USA). DNA was digested with RNase-free DNase I for 20 min at 37 °C. The quantity and quality of the extracted RNA were checked using a NanoVue spectrophotometer (GE Healthcare Bio-Sciences, Tokyo, Japan), Qubit 2.0, fluorometer (Life Technologies, Darmstadt, Germany), and Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA).
cDNA synthesis for HiSeq 2500 and RNA-seq analysis. cDNA libraries for RNA-seq were prepared using the Truseq Stranded mRNA LT kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. First-and second-strand cDNA synthesis, adaptor ligation, and amplification were carried out. After cDNA preparation, the quantity and quality of cDNA were checked using the KAPA qPCR SYBR green PCR kit (Geneworks, Thebarton, Australia) and Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA). Twelve cDNA libraries were pooled in equal quantities, and six tubes from 72 libraries were prepared. All libraries were subjected to single-end sequencing of 150 bp fragments on Hiseq2500 (Illumina, San Diego, CA). Sequence read quality was determined using FastQC 29 . The adaptor sequences, low-quality sequences (< 20 quality scores), and too short reads (< 50 bp) were removed from the sequence data using Trimmomatic v0.35 30 . DEG analysis. Gene expression levels were compared using the generalized linear model (GLM) approach implemented in the edgeR Bioconductor package 33 . For the DEG analysis, we kept only genes with at least one count per million (CPM) in at least three samples. Then, the DEG analysis was carried out using the read count data normalized with the trimmed mean of M-values (TMM) method 34 . We compared gene expression levels (head or body) during each molting process at four different periods (worker/nymph, pre-GP, GP, and molt individuals). We selected DEGs with an FDR cut-off of 0.01 and log2 fold change cut-off of 1. To test the effect of artificial 20E treatment, the number of DEGs was compared between 20E-induced GP workers and natural GP workers using edgeR.
Real-time quantitative PCR (qPCR) analysis. To validate the DEGs list obtained in "DEG analysis" section, we focused on the genes that were highly expressed during the worker-presoldier molt (total eight; Fig. S1 and Table S1). We especially focused on those genes; some of which could be involved in the soldierspecific morphogenesis, mainly occurred in the head parts during the worker-presoldier molt 11 . Independent real-time qPCR analysis was used to compare gene expression levels in heads (pre-GP period) between workerpresoldier and worker-worker molts. In October 2020, three new colonies were found in Himi and Yatsuo, Toyama Prefecture. Following the methods described above, each colony received both 20E and JH treatments, and workers were collected at three time points: early (day 1-2), middle (day 2-3) and late (day 3-4). Four heads were placed into one tube, and four tubes were prepared for each category. RNA extraction was carried out using the ReliaPrep RNA Tissue Miniprep System (Promega, Madison, USA). RNA and DNA were quantified using a Qubit fluorometer, and RNA purity and quantity were measured using a NanoVue spectrophotometer. DNase-treated RNA was used for cDNA synthesis using the High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher). Quantitative real-time PCR was performed using a Thunderbird SYBR qPCR Mix (TOYOBO, Japan) and a QuantStudio 3 Real-Time PCR System (Thermo Fisher). To determine an internal control gene, the suitability of six reference genes, EF1-alfa (accession no. AB602838) 35 , NADH-dh (AB602837) 35 , beta-actin (AB520714) 36 , glutathione-S-transferase 1 (GstD1, gene ID: RS001168 19 ), eukaryotic initiation factor 1A (EIF-1, RS005199 19 ), and ribosomal protein S18 (RPS18, RS015150 19 ) was evaluated using GeNorm 37 and NormFinder 38 software. Specific primers were designed against each gene sequence using Primer3Plus 39 (Table S2). We performed the Levene's test on the variance equality using R ver. 3.3.3 (available at https:// cran.r-proje ct. org/). Statistical analysis (two-way analysis of variance (two-way ANOVA)) was performed using the Mac Statistical Analysis ver. 2.0 (Esumi, Tokyo, Japan).

KEGG enrichment analysis.
To assign the RspeOGS1.0 genes to KEGG pathways, we performed the BlastKOALA analysis 43 . Enrichment of DEGs in each of the GO terms was examined using Fisher's exact test, as described above (see "Gene ontology (GO) enrichment analysis" section).

Results and discussion
Induction of each molt. In all colonies, neither worker nor presoldier molts occurred during the controlacetone treatments within 2 weeks (Table S3). The induction rates of the worker and presoldier molts were 73.3-85% and 66.7-90.0% (average), respectively. Because sufficient numbers of the N5 nymphs could not be obtained from multiple colonies (see "Induction of nymphoid differentiation by isolation from the natal nest" section), induction rates of nymphoids were not measured in this study. A previous study 26 reported that the average nymphoid induction rate was approximately 35%.
Sequencing and mapping. The minimum and maximum read counts in all 72 RNA-seq libraries were 12,592,720 and 19,447,093, respectively. Low quality and short sequences (quality score < 20, < 50 bp sequences) were removed, and the remaining reads were counted. The minimum and maximum counts of trimmed sequences were 11,886,368 and 18,346,386, respectively. These RNA-seq reads were used to map the R. speratus genome data (gene model OGS1.0) 19 . The minimum and maximum mapping rates were 57.50 ± 1.57% and 86.23 ± 1.46% (average ± S.D. of triplicates), respectively (Table S4). The rates observed in this study were slightly lower than but comparable with those published of the damp-wood termite Z. nevadensis (65.0-98.58%) 18 . Low mapping rates (< 70%) were observed of the bodies of workers (57.5%) and pre-GP workers (60.6-62.9%) (Table S4), most likely due to the presence of symbiotic protists in their hind guts [44][45][46] . Even in the worker bodies, however, high mapping rates were observed after the gut purge (82.2-85.5%), owing to the loss of symbionts from their hind guts. We excluded the data probably derived from symbiotic protists and other microorganisms in the following analysis, because these data might not map to the genome sequence, as previously shown in Z. nevadensis 18 . GP workers and the natural GP workers (FDR < 0.01), when all genes were analyzed (11,884 and 12,481 genes for head and body, respectively). Because presoldier differentiation is usually quite rare in mature colonies under natural conditions, most natural GP worker castes are composed of individuals that are at the stage before the worker molt. These results suggest that artificial 20E treatment could be used to induce more developmentally synchronized individuals during and after the worker molt. In the mature colony, the timing of worker molts is usually different from each individual, and we can never collect the pre-GP workers without the 20E treatment.

Numbers of DEGs during each molt.
The numbers of differentially expressed and co-expressed genes in each molt are shown in the Venn diagram (Fig. 2). This is the first evidence to show the DEGs among all termite caste differentiations, and many co-expressed genes (614 and 1025 for heads and bodies, respectively) responsible for three completely different molts in R. speratus. The number of overlapping genes between the worker and presoldier molts (head: 614 + 718, body: 1025 + 663) were larger than those between the worker and nymphoid molts (head: 614 + 269, body: 1025 + 67) and between the presoldier and nymphoid molts (head: 614 + 20, body: 1025 + 75). This may be due to the different initial stages (N5 nymphs) of the nymphoid molt, compared to those of the other molts (old-age workers).

Verification of DEGs in worker-presoldier molt.
We concentrated on some genes that were highly expressed in the heads of JH-treated workers during the pre-GP period (Table S1). Based on the stability values obtained by GeNorm and NormFinder software (Table S5) Relatively higher gene expression levels were observed in JH-treated than in 20E-treated workers, and significant differences were detected in seven out of eight genes examined ( Fig. 3; two-way ANOVA, P < 0.05). These results validated the DEGs list derived from the RNA-seq analysis. Specific morphogenesis of (pre)soldiers (e.g., frontal gland formation in the head part 11 ) is probably initiated during the pre-GP period 23 , and thus the function of the identified genes should be studied in-depth to facilitate a full understanding of the proximate mechanism of soldier differentiation.
Gene ontology (GO) terms specifically observed during each molt. The number of GO terms specifically observed during the worker-worker molt (odds ratio > 1) was much larger than that of the other two molts (Tables 1, 2, 3). In particular, the expression levels of genes involved in chemical transport and transporter activity specifically fluctuated during the worker-worker molt (16 out of 32 specific GO terms; Table 1 and S6). In contrast, the expression of genes involved in the chemical compound biosynthetic or metabolic processes specifically fluctuated during the worker-presoldier molt (12 out of 17 specific GO terms; Table 2). For example, in the isoprenoid metabolic processes (GO: 0006720), expression levels of juvenile hormone acid methyltransferase (JHAMT) (gene ID RS007861), and predicted genes of cytochrome P450 305a1 (RS000980, RS000985, and RS013788) fluctuated significantly (Table S7). Although this GO term was not found in the head (including the corpora allata), expression changes in these genes may be related to JH titer changes during the presoldier molt. In the same GO term, expression levels of predicted genes of geranylgeranyl pyrophosphate synthase (GGPPS) (RS007482 and 100014) were also shown to be significantly changed. A previous study found that several GGPPS paralogs were specifically expressed in the frontal gland cells of R. speratus soldier heads 35 .   47 , these findings may be related to the production of chemicals involved in both defensive (a toxic secretion to the enemies) and non-defensive (regulation of soldier differentiation and antimicrobial activity) roles [48][49][50] . Moreover, in the amine metabolic processes (GO: 0009308, 0044106 and 0006576), expression levels of N-b-alanyldopamine (NBAD) synthase (= ebony, RS003967) and predicted genes of Aromatic-L-amino-acid decarboxylase (= DDC, RS006642) and henna (RS010135) were significantly fluctuating (Table S7). Henna is an enzyme that catalyzes the reaction from phenylalanine to tyrosine in Locusta migratoria 51 while DDC and ebony are important enzymes in the tyrosine metabolic (TM) pathway required for the cuticular tanning of insects 52,53 . Changes in the expression of these genes were consistent with the results obtained during soldier differentiation in Z. nevadensis 54 . However, the expression levels of genes related to chitin-based cuticle development (GO ID 0040003), including laccase2 (= Lac2, RS002050) and predicted genes of some cuticle proteins (RS004549, RS004555, RS004558, RS004559, and RS013966), fluctuated significantly during the nymph-nymphoid molt (Table 3 and S8). In Z. nevadensis, high Lac2 expression was observed not only during the presoldier-soldier molt, but also during the nymph-alate molt 54 . Although alates have well-developed black cuticles, nymphoids also have dark pigmentation when compared to nymphs, particularly in the head region 55,56 . Thus, some genes involved in cuticular formation may be highly expressed during nymph-nymphoid molt.
Kyoto encyclopedia of genes and genomes (KEGG) terms specifically observed during each molt. Nine specific KEGG terms (odds ratio > 1) were observed during the worker-worker molt (Table 4).
From the head, more than half of the genes in the insect hormone (JH and 20E) biosynthesis pathway (KEGG:  (Table S1 and Fig. S1). Gene expression levels (mean ± S.E., n = 4) were compared between worker-presoldier molt (gray line) and worker-worker molt (black line). JH-and 20E-treated individuals collected at three time points [early (day 1-2), middle (day 2-3) and late (day 3-4)] during pre-GP period, and workers in the colony (day 0), were used for the analysis. All data are consistent with the use of parametric statistics by the Levene's test (p values are shown in the "Verification of DEGs in worker-presoldier molt" section). Statistical results of twoway ANOVA are described in each box (*p < 0.05, **p < 0.01). Significant differences among time points were only observed in RS002081, and no interaction was detected between treatments and time points in all genes. www.nature.com/scientificreports/ ko00981) were identified as DEGs. These genes were predicted to be juvenile hormone epoxide hydrolase (JHEH) 1-like (RS011542), JHAMT (RS007861), and predicted genes of cytochrome P450 (two mitochndrial (P450 302a1 (RS012246) and 315a1 (RS010451)) and three nuclear genes (P450 18a1 (RS002863), 306a1 (RS002862) and 307a1 (RS010514)) ( Table S9). JHEH is an enzyme that inactivates JH. In a previous study, JHEH hydrolyzed JH III to JH III-diol in the cat flea Ctenocephalides felis 57 . Thus, changes in the expression of these JH-related genes in the head may be related to changes in JH titer during the worker-worker molt. Moreover, in the body region, the expression levels of many genes involved in the fatty acid biosynthesis pathway (ko00061) fluctuated (Table 4). In the biosynthesis of unsaturated fatty acids and PPAR signaling pathways (ko01040 and ko03320), the expression levels of 14 predicted Acyl-CoA Delta (11) desaturase, which is an important gene for lepidopteran sex pheromone synthesis (derived from fatty acids) 58 , significantly fluctuated in heads (Table S9). Thus, Acyl-CoA Delta (11) desaturase could be involved in the synthesis of some releaser pheromones or related substances derived from worker heads. Alternatively, expression patterns of genes related to biosynthesis and/or metabolism of fatty acids may reflect the changes of nutrition status during a molt. A biochemical analysis of workers prior to the molt should be performed to clarify these possibilities. During the worker-presoldier molt, a total of nine specific KEGG terms (odds ratio > 1) were observed ( Table 4), most of which were metabolic pathways of amino acids (histidine, phenylalanine, and tryptophan) www.nature.com/scientificreports/ (ko00340, ko00360, and ko00380). These results suggest that nutrition status of individuals is remarkably changed during the molting event. The other remarkable KEGG terms were the terpenoid backbone biosynthesis (ko00900) and longevity regulating pathway (ko04212), both of which were specifically observed in the head.
In the former pathway, expression levels of 10 genes fluctuated significantly, including a total of eight predicted geranylgeranyl pyrophosphate synthase (GGPPS) genes (Table S10). These results were not contradictory to the GO analysis ("Verification of DEGs in worker-presoldier molt" section). These findings also suggest that there are numerous GGPPS paralogs in the R. speratus genome that are specifically expressed in soldier frontal gland cells, as shown in the previous study 35 . Recent genome analysis of R. speratus clarified this possibility. Total 13 GGPPS paralogs were observed in the same scaffold and some of these were highly expressed in soldier heads 19 . Further research is required to determine whether they play a role in the biosynthesis of soldier defensive chemicals during the presoldier molt. In the longevity regulating pathway, expression levels of the predicted insulin-like receptor (RS007018) were found to fluctuate in heads during the presoldier molt (Table S10). A previous study found that insulin receptor was strongly expressed in mandibular epithelial tissues, and that RNAi treatment caused soldier-specific morphogenesis, including mandibular elongation, in Hodotermopsis sjostedti 12 .
The present results also suggest that insulin/insulin-like signaling activity and related specific cell proliferation may activate drastic morphological reorganization during presoldier differentiation. Only one specific KEGG term was observed in the heads during nymph-nymphoid molt (Table 4). Five out of seven DEGs were predicted to be angiotensin-converting enzyme (ACE) (RS009620, RS009621, RS009623, RS009070, and RS009071) (Table S11). Several insects' brains showed immune responses to ACE, suggesting that insect ACE is involved in the activation and/or inactivation of peptide hormones in the nervous system 59 . Furthermore, several ACE paralogs have been identified in the genome of Bombyx mori, and highly expressed

Conclusion
This study aimed to develop a comprehensive list of the molecular mechanisms underlying termite caste differentiation. We performed both GO and KEGG enrichment analyses using the differential gene expression profiles of each molt (worker, presoldier, and nymphoid molts) of R. speratus. Based on the timing of initiation of the gut purge (GP), we identified three developmental stages during each molt (pre-GP, GP, and after the molt). The results revealed that a large number of genes were differentially expressed during each molt and/or were differentially expressed in each body part. In each molt, the expression of genes, related to JH titer changes, nutrition status, and cell proliferation, specifically fluctuated. These differences may have a crucial impact on caste differentiation. To understand the regulatory mechanisms of caste differentiation in termites, more functional and detailed expression analyses of the genes listed here are required. It should be noted that some gene expressions, especially observed in the worker-presoldier molt, may be affected by the artificial hormone treatment. Further comparative transcriptome analysis should be performed between artificial conditions described in this paper and natural conditions, the latter of which has been reported in other species 17 www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.