Loss of males from mixed-sex societies in termites

Glyptotermes nakajimai contains multiple asexual populations

To confirm the presence of asexual populations of G. nakajimai, we investigated the colony composition of 74 mature colonies from ten populations. Colonies contained multiple reproductives and up to 2300 non-reproductive members (soldiers, workers, nymphs, alates, and young instars; see Fig. 1). In lower termites, including the family Kalotermitidae, individuals (except young instars) usually exhibit sexual dimorphism in their external morphology, such as the female-specific character of having an elongated seventh sternite [23, 24] (but see [25]). Dissections revealed that our sexing of G. nakajimai individuals on the basis of the configuration of the caudal sternites was 100% accurate in reproductives (n = 30), soldiers (n = 30), and workers (n = 60) (Additional file 1: Figure S1), confirming that the sex of these castes of G. nakajimai can be easily determined from sternite morphology. In the Shikoku and Kyushu populations, all reproductives (n = 194), soldiers (n = 316), and workers (n = 3700) collected from 37 colonies were female (Fig. 2a and Table 1). In contrast, all 37 colonies collected from the other populations (Honshu, Amami-Oshima Island, Okinawa Island, and Ogasawara Islands) contained both male and female reproductives (females/[females + males] = 0.53 ± 0.02 SEM [n = 548]), soldiers (females/[females + males] = 0.42 ± 0.03 SEM [n = 713]), and workers (females/[females + males] = 0.49 ± 0.01 SEM [n = 3700]) (Fig. 2a and Table 2).

		No. reproductives (f/[f + m])			No. soldiers (f/[f + m])	No. workers (f/[f + m])
Colony codea	Location	A-derived	Neo.	Total
	Shikoku, Japan
TO150911A	Tokushima, Tokushima	6/6	1/1	7/7	17/17	100/100
TO150911B	Tokushima, Tokushima	2/2	1/1	3/3	4/4	100/100
MR150216B	Muroto, Kochi	1/1	1/1	2/2	9/9	100/100
MR150216C	Muroto, Kochi	1/1	1/1	2/2	3/3	100/100
MR150216E	Muroto, Kochi	2/2	0/0	2/2	8/8	100/100
MR150216G	Muroto, Kochi	5/5	0/0	5/5	14/14	100/100
MR150216H	Muroto, Kochi	3/3	1/1	4/4	11/11	100/100
MR150216I	Muroto, Kochi	6/6	0/0	6/6	15/15	100/100
MR150216J	Muroto, Kochi	5/5	1/1	6/6	12/12	100/100
MR150216K	Muroto, Kochi	3/3	2/2	5/5	12/12	100/100
MR150216M	Muroto, Kochi	2/2	0/0	2/2	8/8	100/100
MR150217B	Muroto, Kochi	5/5	5/5	10/10	17/17	100/100
MR150217C	Muroto, Kochi	2/2	1/1	3/3	4/4	100/100
MR150217D	Muroto, Kochi	2/2	2/2	4/4	13/13	100/100
MR150910A	Muroto, Kochi	12/12	1/1	13/13	14/14	100/100
MR150910B	Muroto, Kochi	7/7	1/1	8/8	22/22	100/100
MR150910D	Muroto, Kochi	16/16	0/0	16/16	13/13	100/100
AS141110A	Ashizuri, Kochi	10/10	1/1	11/11	17/17	100/100
AS141110B	Ashizuri, Kochi	5/5	0/0	5/5	1/1	100/100
AS141110D	Ashizuri, Kochi	6/6	3/3	9/9	6/6	100/100
AS141111A	Ashizuri, Kochi	3/3	2/2	5/5	1/1	100/100
AS141111H	Ashizuri, Kochi	3/3	0/0	3/3	6/6	100/100
AS141111I	Ashizuri, Kochi	2/2	0/0	2/2	17/17	100/100
AS141111J	Ashizuri, Kochi	2/2	0/0	2/2	1/1	100/100
AS141111K	Ashizuri, Kochi	8/8	2/2	10/10	16/16	100/100
	Kyushu, Japan
SK150715A	Saiki, Oita	2/2	3/3	5/5	3/3	100/100
SK150715B	Saiki, Oita	2/2	0/0	2/2	1/1	100/100
TI150324A	Toi, Miyazaki	3/3	0/0	3/3	7/7	100/100
TI150324B	Toi, Miyazaki	2/2	0/0	2/2	1/1	100/100
TI150324C	Toi, Miyazaki	0/0	3/3	3/3	4/4	100/100
TI150324D	Toi, Miyazaki	4/4	0/0	4/4	3/3	100/100
TI150728A	Toi, Miyazaki	2/2	1/1	3/3	7/7	100/100
ST150322C	Sata, Kagoshima	2/2	0/0	2/2	7/7	100/100
ST150323L	Sata, Kagoshima	2/2	0/0	2/2	3/3	100/100
ST160304A	Sata, Kagoshima	4/4	0/0	4/4	4/4	100/100
ST160304B	Sata, Kagoshima	15/15	0/0	15/15	13/13	100/100
ST160304C	Sata, Kagoshima	4/4	0/0	4/4	2/2	100/100
Total		161/161	33/33	194/194	316/316	3700/3700
Mean (SEM)				100 (0)	100 (0)	100 (0)

All collected reproductives, soldiers, and 100 workers from each colony were sexed. A-derived: alate-derived; Neo: neotenic; f/[f + m], females/[females + males]

^aNumbers in colony codes indicate the dates when colonies were collected. For example, colony TO150911A was collected on 11 September 2015

		No. reproductives (f/[f + m])			No. soldiers (f/[f + m])	No. workers (f/[f + m])
Colony codea	Location	A-derived	Neo.	Total
	Honshu, Japan
IZ150430A	Izumo, Kushimoto, Wakayama	13/22	0/1	13/23	2/11	53/100
IZ150430B	Izumo, Kushimoto, Wakayama	3/5	0/1	3/6	2/7	49/100
IZ150430C	Izumo, Kushimoto, Wakayama	1/2	0/1	1/3	6/17	62/100
SN150430B	Misaki, Kushimoto, Wakayama	31/53	2/3	33/56	19/30	58/100
SN150430C	Misaki, Kushimoto, Wakayama	36/68	11/25	47/93	17/29	54/100
SN150430D	Misaki, Kushimoto, Wakayama	5/10	1/1	6/11	7/9	52/100
SN150430E	Misaki, Kushimoto, Wakayama	3/6	3/4	6/10	2/7	68/100
SN150430F	Misaki, Kushimoto, Wakayama	3/6	0/1	3/7	4/8	47/100
SN150430G	Misaki, Kushimoto, Wakayama	13/25	2/2	15/27	5/12	50/100
SN150501A	Misaki, Kushimoto, Wakayama	5/8	10/14	15/22	14/33	49/100
	Amami-Oshima Is., Japan
NZ150526A	Naze, Kagoshima	3/9	3/3	6/12	2/8	44/100
NZ150526C	Naze, Kagoshima	1/2	3/4	4/6	4/12	60/100
NZ150526F	Naze, Kagoshima	8/19	0/0	8/19	6/15	55/100
NK150527A	Setouchi, Kagoshima	0/1	3/4	3/5	3/16	49/100
NK150527C	Setouchi, Kagoshima	3/6	0/0	3/6	13/22	53/100
	Okinawa Is., Japan
HD150225A	Hedo, Okinawa	3/8	0/0	3/8	6/17	41/100
HD150225B	Hedo, Okinawa	4/7	7/14	11/21	5/37	40/100
HD150225C	Hedo, Okinawa	1/2	3/4	4/6	4/6	43/100
HD150225D	Hedo, Okinawa	2/3	1/1	3/4	3/8	42/100
HD160328B	Hedo, Okinawa	2/4	0/0	2/4	0/9	48/100
HD160328C	Hedo, Okinawa	2/4	0/0	2/4	5/15	59/100
NJ150226A	Nakijin, Okinawa	0/1	2/6	2/7	1/4	27/100
YJ150226A	Yagajijima, Okinawa	2/4	0/0	2/4	6/11	45/100
YJ150226B	Yagajijima, Okinawa	1/3	0/0	1/3	2/7	48/100
	Ogasawara Isis., Japan
CC151014A	Chichijima, Tokyo	22/45	1/2	23/47	8/19	52/100
CC151014B	Chichijima, Tokyo	1/2	0/0	1/2	1/10	51/100
CC151014G	Chichijima, Tokyo	0/1	5/5	5/6	7/13	42/100
CC151014H	Chichijima, Tokyo	3/5	0/0	3/5	10/19	43/100
CC151015A	Chichijima, Tokyo	1/3	2/3	3/6	14/21	43/100
CC151015B	Chichijima, Tokyo	7/14	6/10	13/24	17/36	38/100
CC151015C	Chichijima, Tokyo	2/4	0/0	2/4	12/20	53/100
CC151015H	Chichijima, Tokyo	1/3	0/0	1/3	10/14	53/100
HH151016A	Hahajima, Tokyo	0/0	2/3	2/3	8/13	47/100
HH151016D	Hahajima, Tokyo	6/10	0/0	6/10	38/83	56/100
HH151016E	Hahajima, Tokyo	12/23	0/0	12/23	4/37	24/100
HH151016F	Hahajima, Tokyo	11/25	0/0	11/25	29/44	46/100
HH151016H	Hahajima, Tokyo	11/21	1/2	12/23	12/34	60/100
Total		222/434	68/114	290/548	308/713	1804/3700
Mean (SEM)				0.53 (0.02)	0.42 (0.03)	0.49 (0.01)

All collected reproductives, soldiers, and 100 workers from each colony were sexed. A-derived, alate-derived; Neo., neotenic; f/[f + m], females/[females + males]

^aNumbers in colony codes indicate the dates when colonies were collected. For example, colony IZ150430A was collected on 30 April 2015

To rule out the possibility that cryptic males exist in the all-female populations (Shikoku and Kyushu), we examined sperm storage among egg-laying queens by observing the spermathecae (sperm storage organs) stained by propidium iodide. Spermathecae of all examined queens (n = 12) in the all-female populations of G. nakajimai were empty, whereas those of all examined queens (n = 12) in the other mixed-sex populations were filled with sperm (P < 0.0001, Fisher’s exact probability test) (Fig. 2a). This is, to our knowledge, the first demonstration of queens in natural colonies that lack sperm in their spermathecae. In all the cases, we could find in the literature every termite queen in natural colonies that has been checked has had sperm in her spermatheca, regardless of the presence or absence of the kings (e.g., [26, 27]).

We also compared the hatching success of unfertilized eggs between all-female and mixed-sex populations. As expected, the all-female populations exhibited high rates of hatching success of unfertilized eggs (112 of 134 eggs [83.6%] hatched) (Fig. 3). Notably, some unfertilized eggs from the mixed-sex population also hatched (7 of 193 eggs [3.6%] hatched), indicating the presence of tychoparthenogenesis (occasional development of eggs without fertilization) [1, 28, 29] in sexual colonies of this species (Fig. 3). No significant difference in hatching success was observed between unfertilized eggs of the all-female population (112 of 134 eggs [83.6%] hatched) and fertilized eggs of the mixed-sex population (117 of 127 eggs [92.1%] hatched) (P = 0.12, Fisher’s exact probability test with Bonferroni correction) (Fig. 3). Overall, our results clearly indicate complete asexuality in the Shikoku and Kyushu populations, in contrast to the other normal sexual populations (Fig. 2b).

A single origin of asexuality in G. nakajimai

To elucidate the evolutionary relationships among asexual and sexual populations of G. nakajimai, we conducted phylogenetic analyses of Glyptotermes termites based on two independently evolving markers: mitochondrial cytochrome c oxidase subunit II (COII) and nuclear internal transcribed spacer 2 (ITS2) sequences. The monophyly of G. nakajimai was unequivocally supported in each analysis (COII: Bayesian posterior probability [BPP] = 1.00; RAxML bootstrap support [RBS] = 100%, ITS2: BPP = 1.00; RBS = 100%) (Fig. 4, Additional file 2: Figure S2). All sequences of COII and ITS2 from six collection sites across the range of asexual populations were identical (COII: 100% identity, ITS2: 100% identity). Likewise, COII and ITS2 sequences from nine collection sites across the range of sexual populations were highly similar to each other (COII: 99%–100% identity, ITS2: 96%–100% identity), and the monophyly of all sexual populations was well-supported (COII: BPP = 1.00; RBS = 93%, ITS2: BPP = 0.97; RBS = 100%) (Fig. 4, Additional file 2: Figure S2). Given that all other known species of termites, including members of the genus Glyptotermes, are sexual, these results demonstrate a single origin of asexuality within G. nakajimai. The estimated divergence time between the two lineages was 14.1 million years ago (95% confidence interval = 8.1–22.4 million years ago [Mya]) (Fig. 5). Whether asexuality has persisted in G. nakajimai for this period of time is not known. It is possible that the most closely related sexual relatives of the asexual G. nakajimai were not represented among our samples (perhaps because they are extinct, or found in areas we did not collect from), and that asexuality evolved much more recently than 14.1 Mya. Further work is required to investigate this issue.

Previous work has shown that the external morphology and cuticular hydrocarbon profiles of G. nakajimai representatives across its range are indistinguishable [21], including the asexual populations reported here. To further investigate potential fine-scale differences between the asexual and sexual lineages, we performed karyotyping. The sexual lineage consistently displayed 2n = 34 chromosomes (female: n = 12, male: n = 12), while the asexual lineage consistently harbored a complement of 2n = 35 (n = 18) and contained a trisomy, most probably of chromosome 16 (Fig. 6). To confirm our observations that the asexual lineage has an extra chromosome, we compared genome size (based on the C value or nuclear DNA mass) between the asexual and sexual lineages. As expected, the genome size of the asexual lineage was significantly higher than that of the sexual lineage (colony: F_{4, 24} = 0.80, P = 0.54; lineage: F_{1, 24} = 55.76, P < 0.0001; nested ANOVA with colonies nested within lineages [the asexual lineage and the sexual lineage]), although the ploidy levels were similar between the two lineages (Table 3). These results are consistent with the two lineages containing different chromosome profiles (Fig. 6), in addition to their molecular sequence profiles (Fig. 4, Additional file 2: Figure S2).

Lineage	C-value (pg)a	Relative ploidy level
Asexual lineage	1.27 ± 0.00****	2.07
Sexual lineage	1.23 ± 0.00	2

Data for the asexual lineage was compared with the sexual lineage using nested ANOVA. ****, P < 0.0001

^aValues are mean ± SEM (n = 15)

Greater uniformity of head size in all-female soldiers of asexual colonies

Phenotypic variations among individuals within a colony often provide opportunities for efficient task partitioning, such as sexual specialization in tasks in animal societies [13, 15, 30, 31]. However, such variations among individuals performing the same task may lower task-efficiency when individual variation reduces group performance. The most common defense mechanism in lower termites, including the genus Glyptotermes, is phragmosis, where soldiers with plug-like heads block tunnels connecting chambers, thus preventing enemies from invading the nest [32, 33]. Tunnel width within a nest is relatively uniform, and stabilizing selection acts on soldier head width for efficient phragmotic defense in these termites, because soldiers with narrower head widths would tend to allow enemies to pass through more easily, while those with greater head width would tend to clog tunnels within the nest [34, 35]. When sampling G. nakajimai in the field, we found that individuals of several ant species carried off struggling workers when logs infested with G. nakajimai were broken open. However, soldiers typically retreated to the small tunnels between chambers and plugged the tunnels with their heads to prevent the ants from invading intact parts of the nest (i.e., phragmotic defense) (Fig. 7a).

In the sexual lineage of G. nakajimai, male soldiers were significantly smaller than female soldiers in head width (colony: F_4,191 = 5.16, P < 0.001, sex: F_1,191 = 6.86, P < 0.01, two-way ANOVA), as well as in head length (colony: F_4,191 = 10.67, P < 0.0001, sex: F_1,191 = 13.21, P < 0.001, two-way ANOVA) (Additional file 3: Figure S3). Stabilizing selection appears to act on soldier head width because the head width to length ratio differed significantly between the sexes (colony: F_4,191 = 12.46, P < 0.0001, sex: F_1,191 = 12.15, P < 0.001, two-way ANOVA), such that differences in head width between the sexes were reduced (Additional file 3: Figure S3).

To investigate whether asexuality leads to increased stabilization of soldier head width, we compared the within-colony coefficient of variation (CV) of soldier head width between the asexual and sexual lineages. The asexual lineage had significantly smaller within-colony CVs of soldier head width than those of the sexual lineage (P = 0.012, Mann–Whitney U test) (Fig. 7b), indicating reduced intracolonial variation in soldier head width in the asexual lineage. This is presumably due to a loss in sexual dimorphism in the asexual lineage, although the loss of genetic diversity resulting from parthenogenesis may also contribute. We also compared the proportion of soldiers to other individuals in mature field colonies between the asexual and sexual lineages. Mature field colonies of the asexual lineage had significantly lower proportions of soldiers to other individuals than those of the sexual lineage (P < 0.0001, GLM) (Fig. 7c), suggesting the number of soldiers required for colony defense in the asexual lineage appears to be lower than that in the sexual lineage. A reduction in the number of soldiers, which require sibling care because they are unable to feed themselves, is likely to allow investment elsewhere in the colony. Thus, our results hint at the possibility of increased defensive efficiencies arising from the greater morphological uniformity of soldiers in all-female asexual colonies. It is possible that such efficiencies contributed to the persistence and spread of the asexual lineage.

If the sexual lineage produced only female soldiers (or only male soldiers), intracolonial variation in soldier head width could also be reduced. Indeed, sexual specialization of soldiers is common in many termite taxa with sterile workers [13, 15]. However, such sexual specialization is extremely rare in dry-wood termites (Kalotermitidae) [13, 15]. This can be partly explained by their caste-developmental pathway, through which all colony members finally develop into alates, neotenic reproductives (only a small proportion of individuals), or sterile soldiers [13, 14] (Fig. 1). Sexual specialization of soldiers would therefore lead to skewed sex ratios of other colony members within colonies, most importantly alates. This would potentially lead to a failure of some alates from the colony to successfully breed.

Pre-adaptive conditions for the evolution of maleless societies

When the costs of asexuality are not high, a mutation causing asexual reproduction is expected to rapidly spread within a sexual population due to the twofold reproductive advantage of asexual reproduction. Under this scenario, asexual females produce only female offspring, while sexual females continue to produce both male and female offspring [36]. Assuming that approximately equal numbers of offspring are produced by both asexual and sexual colonies, asexual females will increase in frequency in the population and eventually replace sexual females entirely unless mutations that suppress asexuality arise.

We hypothesize the existence of a number of key traits in the ancestors of asexual G. nakajimai, which permitted them to overcome the barriers required for complete loss of males from their mixed-sex societies. Firstly, we hypothesize that the sexual ancestors of asexual G. nakajimai were likely to have been pre-adapted to overcome at least some developmental constraints associated with parthenogenesis. This is because parthenogenetic eggs have to be activated without sperm, and centrioles need to be inherited only from mothers during reproduction [5]. This hypothesis is supported by our finding of tychoparthenogenesis (occasional development of eggs without fertilization) in the sexual lineage, whereby 7 of 193 unfertilized eggs (3.6%) successfully developed into larvae (Fig. 3). Tychoparthenogenesis is thought to provide an important pathway to the evolution of parthenogenesis [29], but is also relatively common in many animal groups, regardless of the presence or absence of asexual lineages [1, 28]. Therefore, tychoparthenogenesis is unlikely to be the only characteristic of the sexual ancestors of G. nakajimai that facilitated the evolution of all-female colonies.

The presumed presence of multi-queen colony foundation in the ancestor of asexual G. nakajimai is also likely to have removed one potential barrier to the spread of the asexual phenotype: the expected failure of many colonies founded by single asexual queens [41]. Termites cannot clean their own bodies by self-grooming, with the exception of their antennae. Mutual grooming by reproductives during incipient colony foundation is therefore thought to play a key role in disease avoidance [41]. During the evolution of termite all-female asexual societies, females that founded colonies with other females (rather than alone) are likely to have been at a significant advantage, both from the benefits of grooming, as well as faster colony growth due to enhanced reproductive output.

The final traits we hypothesize that may have permitted the loss of males in G. nakajimai are those generally associated with the life histories of dry-wood termites. Dry-wood termites are single-site nesters, living in a single piece of dead wood that serves both as nest and as food [42]. This means that individual members do not need to forage outside the nest, reducing the risk of exposure to pathogens and parasites. Importantly, dry-wood nesters can easily disperse over water through wood-rafting. This is evidenced by the presence of many species, including G. nakajimai, in coastal areas and on remote islands [20, 21, 43, 44]. According to the Red Queen hypothesis, sexual reproduction is favored because it helps the population to co-evolve with specialist parasites and pathogens [45]. Long-distance dispersal events might have released sexually reproducing ancestors of asexual G. nakajimai from selection pressures favoring sexual reproduction. This is because expansion into new areas may allow founders to escape from parasites and pathogens [46, 47]. Moreover, following dispersal to new habitats, asexuality could be adaptive because it is an effective way of circumventing the challenges associated with low population densities, such as inbreeding depression and the inability to find mates [48]. Indeed, many asexual lineages of ants are unusually widespread geographically [9].

The mode of parthenogenesis in G. nakajimai

The transition from sexual to parthenogenetic reproduction may have significant negative consequences as a result of increased homozygosity and inbreeding depression [5]. Such increases in homozygosity may occur under two forms of automictic parthenogenesis: gamete duplication or terminal fusion of gametes produced during meiosis. On the other hand, heterozygosity may be maintained in offspring under alternative forms of automictic parthenogenesis: central fusion of gametes produced during meiosis, or apomixis (in which eggs cells are essentially clones produced via mitosis) [5]. Indeed, the mode of thelytokous parthenogenesis in asexual lineages of hymenopteran social insects is either automixis with central fusion or apomixis [8, 9]. The presence of an extra chromosome in all examined members of the G. nakajimai asexual lineage (2n = 35, vs 2n = 34 in members of the sexual lineage; Fig. 6) is suggestive of apomixis, as opposed to automixis. This is because an uneven number of chromosomes is likely to lead to pairing problems during the first stage of meiosis. Apomixis, on the other hand, would lead to a consistent number of chromosomes being passed to offspring [5]. We performed a preliminary investigation of the mode of parthenogenesis in G. nakajimai by genotyping individuals from both the asexual and sexual lineages (including offspring produced through tychoparthenogenesis from the sexual lineage; Fig. 3) at six polymorphic microsatellite markers developed for Glyptotermes termites (Additional file 4: Table S1). In the sexual lineage of G. nakajimai, offspring produced by tychoparthenogenesis were homozygous for a single maternal allele at two polymorphic loci (Gly08 and Gly18; Additional file 5: Table S2) (Additional file 6: Table S3). This pattern is suggestive of either (a) automixis with terminal fusion, where offspring are homozygous for a single maternal allele at all loci if no crossing-over takes place, as reported in other lower termites [16, 49], or (b) automixis with gamete duplication where offspring are completely homozygous for a single maternal allele at all loci, as suggested in a study of higher termites [16]. Unfortunately, the six markers we analyzed showed no polymorphism in the asexual lineage of G. nakajimai (Additional file 5: Table S2). Therefore, further work involving deeper sampling of the genomes of G. nakajimai representatives is required to test the hypothesis that the mode of parthenogenesis in the asexual lineage is apomixis, as suggested by its uneven number of chromosomes (Fig. 6). If our hypothesis is correct, such sampling of the asexual and sexual lineages could also determine whether heterozygous genotypes in the asexual lineage (if present) are derived from the conservation of initial allelic diversity present in the sexual ancestor. Under this scenario, the asexual lineage would have evolved through a transition from normal sexual reproduction directly to apomixis, as is known in Timema stick insects [6]. Alternatively, any heterozygous genotypes in the asexual lineage could be derived from new mutations, indicating that the asexual lineage evolved through an automictic step (as appears to be the case in tychoparthenogenetic G. nakajimai (Additional file 6: Table S3)).

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 (G₀/G₁ 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 MgCl₂, 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.