Termite collection
We collected 74 mature colonies and four incipient colonies, which comprise only primary reproductives and young larvae, of the dry-wood termite G. nakajimai from 15 sites across ten populations (Honshu [Kushimoto], Shikoku [Tokushima, Muroto, and Ashizuri], Kyushu [Saiki, Toi, and Sata], Amami-Oshima Island, Okinawa Island, and Ogasawara Islands, Japan) from November 2014 to March 2016 (Tables 1, 2, and 4). We very carefully dismantled nest woods and extracted all colony members using an aspirator and forceps. The reproductives (kings and queens), soldiers, workers (also called pseudergates), nymphs, alates, and young instars (see Fig. 1) from each colony were placed in a moist unwoven cloth in a 90-mm Petri dish and preserved at − 25 °C until further use. Portions of workers and nymphs from each colony were kept in the laboratory as stock colonies in 90-mm Petri dishes that contained damp chips of sliced Oregon pine wood at 25 °C under constant darkness until subsequent experiments. The sex of all collected reproductives, soldiers, and 100 workers from each colony was determined based on the configuration of the caudal sternites [23, 24] under a stereomicroscope (SZX7; Olympus, Tokyo, Japan). To compare the proportion of soldiers to other individuals of mature field colonies between asexual and sexual populations, we used general linear modeling (GLM) (JMP 9; SAS Institute, Cary, NC).
Sexing G. nakajimai reproductives, soldiers, and workers
To test the accuracy of sexing individuals using external morphology in G. nakajimai, we used reproductives, soldiers, and workers of three mature field colonies of the sexual lineage (SN150430B, NZ150526F, and HH151016F in Table 2). Under a stereomicroscope (SZX7; Olympus, Tokyo, Japan), five putative male and five putative female reproductives, five putative male and five putative female soldiers, and ten putative male and ten putative female workers were randomly chosen from each colony using the configuration of the caudal sternites [23, 24]. The true sex was then ascertained through inspection of the gonads [51].
Spermatheca analysis
To examine the insemination status of egg-laying queens, we developed a spermatheca staining method to observe inseminated sperm. Two egg-laying queens were randomly chosen from each of the six field colonies collected from asexual populations (TO150911A, TO150911B, MR150910A, MR150910B, SK150715A, and TI150728A in Table 1) and from sexual populations (IZ150430B, SN150430C, SN150430E, NZ150526C, NK150527C, and CC151015C in Table 2). Queens were dissected in phosphate-buffered saline (PBS), and their spermathecae were removed under a stereomicroscope (SZX7; Olympus). Whole spermathecae of individual queens were immersed in fixation solution (4% paraformaldehyde) for 30 min at room temperature, washed with PBST (PBS with 0.1% Triton X-100) three times (5 min each), transferred into staining solution (propidium iodide [Dojindo, Kumamoto, Japan], final concentration: 2 ng/ml in PBST), and incubated for 20 min at room temperature. The spermathecae were washed with PBST for 10 min, mounted on glass slides with mounting medium (Vectashield H-1000; Vector Laboratories, Burlingame, CA), and kept at 4 °C until analysis. To acquire fluorescent images of spermathecae, we photographed the stained spermathecae using a confocal laser scanning microscope (TCS-SP8; Leica, Wetzlar, Germany) equipped with a 100×/1.4 oil HC PL APO objective. The parameter settings of the microscope were as follows: excitation wavelength, 561 nm; laser intensity, 20–40%; gain, 20–50; and laser type, DPSS 20 mW 561 nm. To compare the rate of insemination of egg-laying queens between asexual and sexual populations, we used Fisher’s exact probability test (Statistica 10; StatSoft, Tulsa, OK).
Comparison of hatching success of unfertilized eggs between the asexual and sexual lineages
To investigate the difference in costs of parthenogenesis between the asexual and sexual lineages, we compared the hatching success of unfertilized eggs between the two lineages. Virgin alates were obtained from three stock colonies of the asexual lineage (AS141110A, AS141111I, and AS141111K in Table 1) and of the sexual lineage (IZ150430C, SN150430F, and SN150501A in Table 2), separated by sex before swarming began, and maintained in 90-mm Petri dishes containing moist unwoven clothes until they shed their wings. Then, two females, or a female and a male were randomly chosen from each colony and placed in a 52 × 76-mm glass cell that contained mixed sawdust bait blocks, as described in a previous study [41]. All treatments were replicated 20 times for each colony. The glass-cell colonies were kept at 25 °C under constant darkness for 100 days. We counted eggs and larvae by checking the glass-cell colonies every 3 days. The hatching success, calculated as percentage of eggs hatched within 100 days after colony foundation, was compared among eggs of glass-cell colonies founded by pairs of female reproductives (i.e., unfertilized eggs) of the asexual lineage, those of glass-cell colonies founded by pairs of female reproductives (i.e., unfertilized eggs) of the sexual lineage, and those of glass-cell colonies founded by pairs of female and male reproductives (i.e., fertilized eggs) of the sexual lineage using Fisher’s exact probability tests with Bonferroni correction (Statistica 10; StatSoft). Because egg protection behavior by reproductives is indispensable for egg survival, data for the glass-cell colonies in which at least one reproductive died were excluded from the analysis. Moreover, because there were no significant differences between the natal colonies and between the glass-cell colonies within egg types (unfertilized eggs of the asexual lineage, unfertilized eggs of the sexual lineage, and fertilized eggs of the sexual lineage), respectively (P > 0.05, Fisher’s exact probability test with Bonferroni correction [Statistica 10; StatSoft]), we pooled the data for both the natal colonies and the glass-cell colonies of each egg type.
Phylogenetic analyses
To infer intraspecific relationships among asexual and sexual populations of G. nakajimai, we conducted phylogenetic analyses of Glyptotermes termites. One worker randomly chosen from each of the 15 G. nakajimai collection sites (six collection sites of asexual populations and nine collection sites of sexual populations) was used for phylogenetic analyses. The sequences of other kalotermitid species obtained in this study and from GenBank were also used in the analyses. Heads of individual termites were ground in Chelex-100 resin solution (Bio-Rad, Richmond, CA) and DNA was extracted and purified in accordance with standard Chelex-based protocols [52]. Individuals were sequenced for the mitochondrial cytochrome c oxidase subunit II (COII) and the nuclear internal transcribed spacer 2 (ITS2) genes. Primer sequences, PCR conditions, and sequencing methods are detailed in previous reports (COII: [53]; ITS2: [54]). We obtained COII (624 bp) and ITS2 (444-bp) sequences. Consensus sequences were aligned using the ClustalX algorithm [55] from the BioEdit 7.0.4.1 sequence editor [56] and corrected manually. Bayesian inference was performed using MrBayes 3.1.2 [57] with the GTR + I + G model for COII and with the HKY + G model for ITS2 selected by the hierarchical likelihood ratio test (hLRT) in MrModeltest 2.3 [58]. Kalotermes flavicollis (Fabricius) (GenBank accession number KX688868) and G. fuscus Oshima (GenBank accession number KX688884) were used as outgroups for the phylogenetic analyses of COII and ITS2, respectively. Methods of phylogenetic analyses are described in a previous report [53].
We also analyzed the sequence alignments using maximum likelihood in RAxML version 7.7.1 [59]. The GTRGAMMA model was selected for the combined datasets and 100 bootstrap replicates were performed. Kalotermes flavicollis (GenBank accession number KX688868) and G. fuscus (GenBank accession number KX688884) were used as outgroups for the phylogenetic analyses of COII and ITS2, respectively.
The COII and ITS2 gene sequences obtained in this study were deposited in the DDBJ/EMBL/GenBank nucleotide sequence databases under accession numbers KX688845–KX688885.
Divergence-dating analysis
We analyzed the COII DNA sequence alignment with a relaxed molecular-clock model using the Bayesian phylogenetic software BEAST 1.8.2. Rate variation was modeled among branches using uncorrelated lognormal relaxed clocks [60], with a single model for all genes. A Yule speciation process was used for the tree prior [61], and posterior distributions of parameters, including the tree, were estimated using MCMC sampling. We performed two replicate MCMC runs, with the tree and parameter values sampled every 1000 steps over a total of 50 million generations. A maximum clade credibility tree was obtained using Tree Annotator within the BEAST software package with a burn-in of 10,000 trees. Acceptable sample sizes and convergence to the stationary distribution were checked using Tracer 1.5 [60]. The molecular clock was calibrated using four minimum age constraints based on the following fossils: Coptotermes priscus Emerson, ~ 26 million years old; Dolichorhinotermes dominicanus Schlemmermeyer and Cancello, ~ 18 million years old; Anoplotermes sp., 18 million years old; and Reticulitermes antiquus (Germar), ~ 44 million years old; (for further details, see [62]). Fossil calibrations were implemented as exponential priors on node times.
Cytological analysis
To examine the karyotypes of the asexual and sexual lineages, we used female alates of six stock colonies of the asexual lineage (TO150911B, MR150910B, AS141110A, SK150715A, TI150728A, and ST160304B in Table 1) and four stock colonies of the sexual lineage (IZ150430C, NK150527C, HD160328C, and HH151016F in Table 2), and male neotenic reproductives of two stock colonies of the sexual lineage (NK150527C and HH151016F in Table 2). Three alates or five neotenic reproductives were randomly chosen from each colony. The somatic chromosomes of alates and neotenic reproductives were observed using the lactic acid dissociation drying method described in previous studies [41, 63].
Genome size analysis
To investigate the difference in genome size between the asexual and sexual lineages, we used female workers of three mature field colonies of the asexual lineage (TO150911A, MR150910A, and ST160304B in Table 1) and of the sexual lineage (IZ150430A, NZ150526C, and HH151016D in Table 2). Using flow cytometry with propidium iodide staining, we estimated genome size (C value). Tissue for flow cytometric analysis was processed with a Cycletest PLUS DNA Reagent Kit (Becton Dickinson, San Jose, CA). All procedures were performed in accordance with the manufacturer’s supplied protocol. Five female workers were randomly chosen from each colony and their heads were ground together with the heads of female Drosophila melanogaster Meigen (strain Oregon R) (1C = 0.18 pg [64]), which served as an internal standard, with a tight-fitting pestle. Detailed descriptions of the nuclear isolation and staining are described in a previous report [65]. Stained nuclei were analyzed for DNA-PI fluorescence using an Accuri C6 Flow Cytometer (Becton Dickinson) at an excitation wavelength of 488 nm. Approximately 1000 cells were acquired for each measurement, and gating was performed using Accuri C6 software v1.0.264.21 (Becton Dickinson). Each C value (pg) was calculated based on the main peak of 2C cells (G0/G1 phase) of G. nakajimai and D. melanogaster. To compare the genome size (C value) between the asexual and sexual lineages, we used nested ANOVA. Assuming that the sexual lineage is diploid, as are all known termite species [11], a relative ploidy level was assigned to the asexual lineage using the following formula: ploidy level = mean genome size of the asexual lineage/mean genome size of the sexual lineage × 2.
Examination of sex differences in soldier morphology
To compare morphology between male and female soldiers, we used soldiers of five mature field colonies of the sexual lineage (SN150501A, CC151015B, HH151016D, HH151016F, and HH151016H in Table 2). Head width and head length were measured for all collected male and female soldiers, excluding severely damaged individuals, from each colony. Head width, head length, the head width to length ratio (head width/head length) were analyzed using two-way ANOVA (Statistica 10; StatSoft, Tulsa, OK).
Comparison of the within-colony CV of soldier head width between the asexual and sexual lineages
To investigate whether asexuality leads to stabilization of soldier head width, we compared the within-colony coefficient of variation (CV = standard deviation/mean) of soldier head width between the asexual and sexual lineages. Excluding severely damaged individuals, all collected soldiers from five mature field colonies of the asexual lineage (TO150911A, MR150217B, MR150910B, AS141110A, and AS141111I in Table 1) and of the sexual lineage (SN150501A, CC151015B, HH151016D, HH151016F, and HH151016H in Table 2) were used for comparison. Head width was measured under a stereoscope (SZX7; Olympus) using a digital imaging system (FLVFS-LS; Flovel, Tokyo, Japan). Within-colony CVs of soldier head width, adjusted for sample size differences [66], were compared between the asexual and sexual lineages using Mann–Whitney U test (Statistica 10; StatSoft).
Development of microsatellite markers
To characterize the asexual and sexual lineages of G. nakajimai genetically, we developed six polymorphic microsatellite markers for termites of the genus Glyptotermes. Genomic DNA was extracted from a pool of ten individuals from a single colony of G. fuscus collected on Okinawa Island using the DNeasy tissue Kit (Qiagen, Hilden, Germany). Microsatellite enrichment was achieved using streptavidin-coated magnetic spheres according to the protocol described in a previous study [49]. Enrichment amplicons were cloned using a pGEM T-Easy cloning kit (Promega, Madison, WI) following the manufacturer’s instructions. Inserted DNA obtained from the color-positive clones was amplified by PCR. PCR products were checked for the insert size by agarose gel electrophoresis with a 100-bp ladder, and PCR products exhibiting a single band of 400–1000 bp were sequenced. Primer sequences, PCR conditions, and sequencing methods are detailed in a previous report [49].
Primers were designed for the 18 inserts containing a microsatellite region with more than nine repeat units. Genomic DNA was extracted from each head of 15 individuals per colonies of G. fuscus collected on Okinawa Island, those of the sexual lineage of G. nakajimai (two to three colonies for each population: IZ150430A, IZ150430B, SN150430F, SN150501A, NZ150526A, NZ150526F, NK150527C, HD160328B, HD160328C, NJ150226A, YJ150226B, CC151014G, CC151015B, HH151016D, and HH151016E in Table 2) and those of the asexual lineage of G. nakajimai (three to four colonies for each population: TO150911A, TO150911B, MR150217B, MR150910B, MR150910D, AS141110A, AS141111H, AS141111K, SK150715A, SK150715B, TI150324A, TI150728A, ST160304A, ST160304B, and ST160304C in Table 1). PCR was carried out in accordance with the following protocol. Each of the 25 μL reaction mixtures contained 0.5 μM each primer (for detailed information on the primers, see Additional file 4: Table S1), 1× PCR buffer, 0.5× Q-solution, 0.5 mM MgCl2, 0.25 mM each dNTP, 0.75 U Taq DNA polymerase (Qiagen), and 1 μL of template DNA. The reactions were run according to the PCR cycle, which consisted of an initial denaturation step at 94 °C for 3 min, followed by 35 cycles of 30 s at 94 °C and 60 s at 60 °C, and one step at 72 °C for 2 min to complete the extension at the end. The PCR products were electrophoresed with the GeneScan-600 LIZ size standard (Applied Biosystems, Foster City, CA) on a 3500 Genetic Analyzer (Applied Biosystems).
Of 18 loci, polymorphisms within the Okinawa Island population of G. fuscus were found for six (Additional file 5: Table S2). Genetic diversity based on number of alleles, observed, and expected heterozygosities, and deviations from Hardy–Weinberg equilibrium for each locus in the Okinawa Island population of G. fuscus, the sexual and asexual lineages of G. nakajimai, respectively, were calculated using Genepop on the Web 4.2 [67]. When appropriate, Bonferroni correction for multiple tests (Statistica 10; StatSoft) was applied.
Mode of parthenogenesis in G. nakajimai
To assess the mode of parthenogenesis in the sexual lineage of G. nakajimai, we genotyped the primary queens and larvae in the glass-cell colonies founded by pairs of female reproductives (the asexual laboratory colonies) described above. All individuals used in this analysis were placed in vials containing 99.5% (vol/vol) ethanol and stored until DNA extraction. Heads of queens or whole larvae were ground in Chelex-100 resin solution (Bio-Rad, Richmond, CA), and DNA was extracted and purified in accordance with standard Chelex-based protocols [52]. Individuals were genotyped at two polymorphic microsatellite loci: Gly8 and Gly18 (Additional file 4: Table S1 and Additional file 5: Table S2). PCR and electrophoresis were performed as described above.