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, respectively20,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 previously23 (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 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 research23, 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 days23, 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.
Figure 1
(A) Artificial induction method of each molt and time schedule for sampling of individuals examined. Each molt required a similar period of time. Ten individuals were collected from 4 different periods (worker or nymph, pre-GP, GP and molt stages) in each molt. To collect individuals before the gut purge (pre-GP), 2 or 3 individuals were sampled daily and mixed during the same period. Gut-purged individuals (GP-0, 1, 2, 3 and 4) were equally mixed for the GP individuals. Molt individuals were sampled within 24 h after each molt. (B) Examples of gut-purged workers (GP workers; arrowheads) with yellowish-white abdomens. Non-GP workers were used for the artificial induction treatment described in (A).
Induction of presoldier differentiation by JH III application
Induction of presoldier differentiation was performed as previously described23,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, nymph-nymphoid) × 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 FastQC29. The adaptor sequences, low-quality sequences (< 20 quality scores), and too short reads (< 50 bp) were removed from the sequence data using Trimmomatic v0.3530. Trimmed sequence data were mapped to the R. speratus genome data (gene prediction model RspeOGS1.019) using Tophat v2.1.0 with Bowtie2 v2.2.331. Counting reads were performed using featureCounts32.
DEG analysis
Gene expression levels were compared using the generalized linear model (GLM) approach implemented in the edgeR Bioconductor package33. 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) method34. 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 soldier-specific morphogenesis, mainly occurred in the head parts during the worker-presoldier molt11. Independent real-time qPCR analysis was used to compare gene expression levels in heads (pre-GP period) between worker-presoldier 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: RS00116819), eukaryotic initiation factor 1A (EIF-1, RS00519919), and ribosomal protein S18 (RPS18, RS01515019) was evaluated using GeNorm37 and NormFinder38 software. Specific primers were designed against each gene sequence using Primer3Plus39 (Table S2). We performed the Levene’s test on the variance equality using R ver. 3.3.3 (available at https://cran.r-project.org/). Statistical analysis (two-way analysis of variance (two-way ANOVA)) was performed using the Mac Statistical Analysis ver. 2.0 (Esumi, Tokyo, Japan).
Gene ontology (GO) enrichment analysis
Translated protein sequences of the R. speratus genome data (gene model RspeOGS1.019) were subjected to BlastP searches against the ‘nr’ database of NCBI with an e-value threshold of 1e-4. In addition, protein domain searches using InterProScan version 5.1740 with the default settings were performed for those proteins. The results of BlastP and InterProScan searches were analyzed using the Blast2GO v2.8 software41 to provide functional annotations including GO terms. GO terms with at least 30 genes were retained for the enrichment analysis. Enrichment of DEGs in each of the GO terms were examined using Fisher’s exact test implemented with the ‘fisher.exact’ function of the ‘Exact2 × 2’ package42 in R.
KEGG enrichment analysis
To assign the RspeOGS1.0 genes to KEGG pathways, we performed the BlastKOALA analysis43. 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).
