Enhanced heterozygosity from male meiotic chromosome chains is superseded by hybrid female asexuality in termites

We compared 4,191 biallelic SNPs across 84 individuals from sexual and asexual populations of G. nakajimai (Fig. 1A and Dataset S1A). Principal coordinate analysis (PCoA) revealed three distinct clusters (Fig. 1B): 1) individuals derived from three sexual populations on small islands in southern Japan and one on the main island of Honshu (collectively referred to hereafter as sexual lineage 1 [SL1]); 2) individuals derived from two asexual populations in Shikoku and three asexual populations in Kyushu (hereafter, asexual lineage 1 [AL1]); and 3) individuals derived from three asexual populations in Shikoku (hereafter asexual lineage 2 [AL2]). Individuals from AL1 and AL2 respectively formed tight genetic clusters, indicative of a lack of genetic variation among members of each asexual lineage. In contrast, SL1 individuals were segregated into four subclusters (Fig. 1B), reflecting the four different geographically separated populations (Fig. 1A). Significant genetic differentiation between each pair of the four populations of SL1 (i.e., pairwise F_ST) was detected (range = 0.376 to 0.548, P < 0.001). On the other hand, all pairwise population F_ST values within each of the asexual lineages were lower and nonsignificant (AL1: range = 0.004 to 0.019, P = 0.245 to 0.502; AL2: range = 0.000 to 0.007, P = 0.253 to 0.535) (SI Appendix, Table S1). An analysis using STRUCTURE revealed that the genetic variability observed in individuals from field colonies of G. nakajimai was best explained using K = 10 (Ln P = −154,082) and recovered the same genetic clustering as the PCoA (Fig. 1C). The two asexual G. nakajimai lineages (i.e., AL1 and AL2) are sympatrically distributed in Shikoku (Ashizuri and Muroto populations), and we found that individuals of the two lineages can coexist within a single colony (Fig. 1B). This is explained by the fact that incipient colonies of the asexual G. nakajimai are founded by multiple queens (range = 2 to 25) (22).

All sexual populations of G. nakajimai showed negative mean inbreeding coefficient (F_IS) values (−0.228 to −0.031), with males displaying lower F_IS values than females in each of the four populations (SI Appendix, Table S2). This suggests that mechanisms exist to avoid inbreeding in males of G. nakajimai. Indeed, similar to the case for other drywood termite species, we observed meiotic chromosome chain formation in G. nakajimai (Fig. 2A), with 30 linked chromosomes found in males (i.e., kings) (n = 2) collected from a population of SL1 (Okinawa Island population). Therefore, out of a total of 17 haploid chromosomes, 15 Y + neo-Y chromosomes are expected to be inherited as a unit. A significant proportion of the genome should therefore be linked and not undergo recombination with the 15 “homologous” X + neo-X chromosomes when in the male germline (28). This is among the highest number of end-to-end linked chromosomes in male meiosis of any animal or plant (29) and expected to lead to this portion of the genome remaining heterozygous. On the other hand, X + neo-X chromosomes, when in the female germline, retain the capacity to undergo recombination with their homologous chromosomes, potentially allowing heterozygosity to become reduced in females under inbreeding. In addition, examination of meiotic chromosomes from a male (i.e., a king) (n = 1) collected from a second population of SL1 (Ogasawara Islands population) revealed a chain of 12 chromosomes plus 11 bivalents (SI Appendix, Fig. S1). Variability in neo-sex chromosome number within a species has previously been reported among different populations of drywood termite species (20, 30).

In agreement with previous hypotheses on the effects of neo-sex chromosomes (11), we found that male G. nakajimai possessed significantly higher mean levels of heterozygosity than females in sexual populations [10.5% versus 7.8% of alleles (SI Appendix, Table S3); P < 0.001, Tukey’s honestly significant difference (HSD) test] (Fig. 2B and Dataset S1B). Mean heterozygosity was up to 47.6% higher in males than females [range 18.9 to 47.6% across four sexual populations (SI Appendix, Table S3)]. However, females in both AL1 and AL2 asexual populations were found to possess significantly high heterozygosity levels (approximately fourfold) when compared with both female and male individuals from sexual populations (P < 0.001, Tukey’s HSD test) (Fig. 2B and Dataset S1B).

To date, all known examples of asexual reproduction in lower termites (i.e., all termites excluding the most derived family Termitidae) involve automixis with terminal fusion (31). Under this mode of reproduction (and also under automixis with gamete duplication), offspring are homozygous for a single maternal allele (if no crossing-over takes place) and are expected to contain an even number of chromosomes. The high heterozygosity displayed by asexual G. nakajimai (Fig. 2B), and the fact that they invariably possess 2n = 35 chromosomes, indicates they do not reproduce via automixis.

On the other hand, clonal reproduction via mitosis (i.e., apomixis) produces offspring that are heterozygous at loci that were also heterozygous in the mother (26, 31), and an uneven number of chromosomes can be maintained through reproduction (32, 33). To assess the mode of reproduction in G. nakajimai, we compared SNP genotypes between queens (i.e., mothers) and their larvae (i.e., offspring) in laboratory-founded colonies whose natal colony had been collected from an asexual population (Ashizuri population), where all individuals (two queens and one larva from each of three laboratory-founded colonies) were identified as AL1 by additional PCoA (SI Appendix, Fig. S2). The offspring inherited all or nearly all heterozygous SNP loci (99.5 to 100% of the loci) from the mother, whereas a small portion of homozygous SNP loci of the mothers changed to heterozygous in the offspring (0 to 0.6% of the loci) (Table 1 and Dataset S1A). Given the presence of new mutations among asexual offspring, this result is suggestive of apomixis where almost all heterozygosity is maintained. In addition, we conducted a crossbreeding experiment with the asexual and sexual G. nakajimai. As expected for apomixis where meiosis is suppressed, the hatching success of hybrid eggs of the asexual and sexual G. nakajimai was much lower than that of unfertilized eggs of the asexual G. nakajimai (P < 0.0001, Fisher’s exact probability test) (SI Appendix, Fig. S3). Only 9 of 59 hybrid eggs (15.3%) developed into larvae, possibly being triploid (infertile) individuals (Dataset S1C). Further work involving cytological observation of chromosomes during oogenesis and embryogenesis of both AL1 and AL2 is required to confirm clonal reproduction as the mode of reproduction in asexual G. nakajimai.

Under the assumption that asexual G. nakajimai individuals arise via clonal reproduction, their high levels of heterozygosity have two potential origins. One is the accumulation of mutations in allele pairs [i.e., the Meselson effect, as seen in ancient asexual animals and plants (34, 35)], another is through hybridization of genetically divergent parent taxa. In the former case, asexual individuals that show high levels of divergence from one another at nuclear loci are also expected to display divergence at mitochondrial loci. In the case of AL1 and AL2, which show clear divergence at nuclear loci (Fig. 1B), examination of a 702–base pair (bp) fragment of the mitochondrial cytochrome c oxidase subunit II (COII) gene between AL1 and AL2 revealed 100% identity (GenBank accession nos. MT387025–MT387032 and MT387033–MT387036, respectively). This suggests that high levels of heterozygosity within asexual individuals have not arisen through gradual accumulation of mutations over long periods of evolutionary time.

Instead, the origin of asexual G. nakajimai can be reasonably explained through intraspecific hybridizations between sexual lineages having different chromosome numbers, 2n = 34 and 2n = 36, respectively. Each parental lineage would have contributed n = 17 and n = 18 chromosomes respectively to their hybrid offspring, explaining the presence of 2n = 35 chromosomes in asexual individuals (22). Following initial hybridization, females would have then undergone clonal reproduction (see the subsection The Mode of Reproduction in Asexual G. nakajimai above). The complete identity of COII sequences between AL1 and AL2 individuals suggests the maternal ancestors of the two lineages are genetically similar or very closely related. Therefore, the genetic differences between AL1 and AL2 (Fig. 1B) can be primarily attributed to differences between their paternal ancestors. Interestingly, the STRUCTURE analysis at K = 2 revealed that individuals of AL2 possessed mixed genetic components from AL1 and SL1 (Fig. 1C). Given the presence of reproductive barriers between the asexual and sexual lineages (described in the subsection The Mode of Reproduction in Asexual G. nakajimai above), this implies that AL2 comprises hybrids with half ancestry from SL1 as the paternal ancestor [which has 2n = 34 chromosomes (22)] and another (unidentified) sexual lineage as the maternal ancestor (consequently having 2n = 36 chromosomes; SL2 in SI Appendix, Fig. S4). AL1 is predicted to have arisen from hybridization between this same maternal ancestor and a third unidentified sexual lineage with 2n = 34 (paternal ancestor; SL3 in SI Appendix, Fig. S4).

Other termite species that exhibit facultative asexual reproduction are known to produce eggs without openings for sperm entry (micropyles) (36). We counted the number of micropyles of eggs collected from a field colony of each of two populations of the asexual G. nakajimai, as well as those collected from a field colony of each of two populations of the sexual G. nakajimai. We found that all examined eggs of the asexual G. nakajimai possessed a substantial number of micropyles (Tokushima population [which probably contains only individuals of AL2, as mentioned in the subsection The Evolution of Asexuality in G. nakajimai below]: 48.29 ± 4.26 SEM, range = 32 to 66, n = 7; Saiki population [which probably contains only individuals of AL1, as mentioned in the subsection The Evolution of Asexuality in G. nakajimai below]: 46.50 ± 4.10 SEM, range = 35 to 61, n = 8) (SI Appendix, Fig. S5 A and B and Dataset S1D). In addition, no significant difference in the number of micropyles was observed between eggs of the asexual and sexual G. nakajimai (colony: F₂, ₂₆ = 0.13, P = 0.88; reproductive type: F₁, ₂₆ = 2.80, P = 0.11; nested ANOVA with colonies nested within reproductive types) (SI Appendix, Fig. S5B). Thus, the evolution of asexuality in G. nakajimai cannot be explained by the production of eggs without micropyles.

To estimate when AL1 and AL2 originated, we performed pairwise comparisons of SNP genotypes between individuals within each lineage. Of the 4,191 loci, 103 and 145 were the largest number of SNP differences between individuals within AL1 and AL2, respectively (Dataset S1A). Given that the mode of reproduction in the asexual lineages of G. nakajimai appears to be apomixis (see the subsection The Mode of Reproduction in Asexual G. nakajimai above), the number of generations of the two lineages can be roughly calculated as the largest number of SNP differences between individuals within a lineage divided by the number of new SNP mutations per generation, divided by two. The above-mentioned comparison of SNP genotypes between the mothers and their offspring in laboratory-founded colonies showed that the number of new SNP mutations in one generation was about 11.33 (calculated from the data in Table 1). Thus, the calculated generation numbers of AL1 and AL2 are 9 and 13, respectively. In drywood termites, new queens are produced after colony maturation, which requires about 4 y (37), and the reported maximum of queen lifespan is 14 y (38). As a result, the estimated ages of AL1 and AL2 were 18 to 63 [i.e., (9 × 4 − 14)/2] and 26 to 91 [i.e., (13 × 4 − 14)/2] years, respectively. These results suggest that the two asexual lineages originated recently (within the last few hundred years). To further examine the maternal origin of AL1 and AL2, we sequenced a 144-bp fragment of the mitochondrial A+T-rich region. We detected three changes across this region between AL1 and AL2 (GenBank accession nos. MT387011–MT387017 and MT387018–MT387024, respectively). Based on intraspecific rates of insect mitochondrial evolution, we estimated the divergence time of the maternal lineage of AL1 and AL2 as 104,000 to 333,000 y ago (39). This suggests that the maternal ancestors of AL1 and AL2 might have had different sequences in their A+T-rich regions.

An analysis of contemporary gene flow between populations revealed evidence for migration between asexual, but not sexual, populations (Fig. 3 and SI Appendix, Table S4). These results indicate that the Tokushima and Sata populations may be the primary source for other populations in AL1 and AL2, respectively. Notably, the Ashizuri and Muroto populations were predicted to have received migrants from both asexual lineages, in contrast to the presence of only one type of asexual lineage in each of other populations (Fig. 3). These results, in combination with the predicted origins of both AL1 and AL2 within the last few hundred years (described in the above paragraph), suggest that human movement of one or more sexual lineages of G. nakajimai may have led to hybridization events and the appearance of asexual lineages. Based on our widespread sampling across the breadth of the distribution of G. nakajimai, it appears that these asexual lineages have replaced two of their predicted sexual ancestors (i.e., SL2 and SL3; SI Appendix, Fig. S4) on the mainland of Japan, since the only sexual lineage we detected was SL1 (which was found only at the southernmost part of Honshu [Kushimoto]). We hypothesize that such replacement of sexual lineages by asexual lineages would have been facilitated by the high levels of heterozygosity in asexual lineages compared with sexual lineages (inferred from our comparison of heterozygosity levels in AL1 and AL2 with SL1; Fig. 2B), despite the presence of neo-sex chromosomes in the sexual lineages. The twofold rate of production of females by asexuals compared with sexuals is another advantage that would promote the spread of the former.

Hybridization between closely related social insect lineages has been shown to have unusual outcomes with regard to the production of different castes within colonies. In Pogonomyrmex harvester ants (Hymenoptera: Formicidae), it has led to the genetic determination of the queen caste and worker offspring with high heterozygosity (40, 41). In G. nakajimai, all colony members possess relatively high heterozygosity in relation to their sexual relatives, and caste determination appears unaffected as a result of hybridization.

We collected 17 mature colonies of G. nakajimai from four sexual populations (Honshu [Kushimoto], Amami-Oshima Island, Okinawa Island, and Ogasawara Islands, Japan) and six asexual (all-female) populations (Shikoku [Ashizuri, Muroto, and Tokushima] and Kyushu [Sata, Toi, and Saiki], Japan) from November 2014 to May 2021. The colonies were transported back to the laboratory with colonized wood. The nest woods were dismantled, and all colony members (reproductives [queens and kings], soldiers, workers, nymphs, alates, and young instars) were extracted using an aspirator and forceps. The eggs were also collected if they were present. Individuals from each colony were placed in a moist unwoven cloth in a 90-mm Petri dish and preserved at −25 °C until sexing was carried out. The sex of individuals was determined based on the configuration of the caudal sternites (22) under a stereomicroscope (SZX7; Olympus). 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.

We conducted high-throughput genome-wide SNP genotyping of individuals from sexual and asexual populations of G. nakajimai. In total, 5 female and 5 male workers were randomly chosen from each of the four field colonies of sexual populations collected in Kushimoto (colony code: IZ150430A), Amami-Oshima Island (colony code: NK150527C), Okinawa Island (colony code: HD160328C), and Ogasawara Islands (colony code: CC151014G); 10 female workers were randomly chosen from each of the four field colonies of asexual populations collected in Ashizuri (colony code: AS141111K), Muroto (colony code: MR150217B), Sata (colony code: ST160304C), and Toi (colony code: TI150728A); two female workers randomly chosen from each of the two field colonies of asexual populations collected in Tokushima (colony code: TO150911B) and Saiki (colony code: SK150715A); and two queens and one larva from each of the three laboratory-founded colonies whose natal colony had been collected in Ashizuri (colony code: AS141111C) (details of the laboratory-founded colonies are described in the subsection Investigation of the Mode of Asexual Reproduction below) were used for genotyping. The termite individuals were preserved in 99.5% (volume/volume) ethanol for genotyping. DNA was extracted from the whole body (excluding gut) of each individual using a High Pure PCR Template Preparation Kit (Roche). Genotyping was performed by Diversity Arrays Technology Pty. Ltd. using DArTsEq (50–53). Four methods of complexity reduction were tested in the Glyptotermes termites and double digestions with PstI-SphI method were selected. Further genotyping methodology details are published elsewhere (54). Approximately 152,000 sequences per barcode per sample were identified and used in marker calling. After quality filtering using the R package “dartR v0.93” (55), our data yielded 4,191 SNPs (average call rate 100%, average reproducibility rate 100%) (Dataset S1A).

To visualize genetic similarities and differences among individuals and populations, we generated a PCoA for individuals from the field colonies using the R package “dartR v0.93” (55).

To investigate patterns of population structure and admixture among populations, we performed a Bayesian clustering analysis of the SNP data using STRUCTURE version 2.3.4 (56) implemented in parallel through StrAuto 1.0 (57). Markov chain Monte Carlo simulations were performed under the assumption of 1 to 10 genetic clusters (K), with 10 replicates of 500,000 iterations for each value of K and with 10% burn-in. All analyses allowed admixture and independent allele frequencies. The Markov chains reached convergence, and alpha values were stable after 200,000 iterations. Owing to the known problem of inferring population clustering from ΔK (58, 59), the optimal K value was inferred using a hierarchical approach by sequential STRUCTURE analyses of clusters identified at each step (60). The results of each replicate of K were summarized using CLUMPAK version 1.1.2 (61) and STRUCTURE HARVESTER (web) version 0.6.94 (62) to obtain marginal likelihoods. Bar plots were generated using DISTRUCT version 1.1 (63).

To estimate the direction and magnitude of contemporary gene flow among populations, we analyzed the SNP data using a Bayesian approach (64) in BayesAss version 3 (65). Each run was 8 × 10 steps, with a burn-in of 2 × 10⁷ steps and sampling every 8,000 steps. The mixing parameter of ΔA (allele frequencies) was optimized at 0.6 to ensure appropriate acceptance rates.

Based on the evidences of the SNP analyses that the asexual G. nakajimai contains two lineages (i.e., AL1 and AL2), we further compared the percentage of heterozygous loci within individuals between males of SL1, females of SL1, females of AL1, and females of AL2 (Dataset S1B) using nested ANOVA followed by Tukey’s HSD test (Statistica 10; StatSoft). Percent data were arcsine transformed prior to analysis. In addition, we performed the following analyses. We measured pairwise population F_ST for SL1, AL1, and AL2 by analysis of molecular variance with 9,999 permutations in GENALEX 6.5 (66). We calculated mean F_IS values for males and females in each sexual population, for males in each sexual population, for females in each sexual population, and for females in each asexual population using GENALEX 6.5 (66).

To examine the male mitotic and meiotic karyotypes of the sexual G. nakajimai, we used two primary kings from two field colonies collected from one of the sexual populations, Okinawa Island (colony code: NJ210511A and NJ210511B), and a neotenic king from one of the field colonies from Ogasawara Islands. The mitotic and meiotic chromosomes of these individuals were successfully observed using the lactic acid dissociation drying method (modified from refs. 67 and 68). Demecolcine (colcemid) was used to block cells in metaphase. The chromosomes of kings from Okinawa Island were stained with DAPI and observed with a confocal microscope (FV1000; Olympus). The chromosomes of a king from Ogasawara Islands were stained with 3% Giemsa and observed with an optical microscope (TBR-1; Yashima Optical).

Based on the evidence of the SNP analyses that the asexual G. nakajimai contains two lineages, we compared mitochondrial A+T-rich and COII sequences between them. The extracted DNA of 14 and 12 individuals (including at least one individual of each lineage from each population when present) that genotyped as described in the subsection Genome-Wide SNP Analyses above were used for A+T-rich and COII sequencing, respectively. A fragment of A+T-rich was amplified by PCR using the following custom primer set, modified from ref. 69: Forward primer (5′-TATTTTGGTGGTGGTTGGTGCAC-3′), reverse primer (5′-CCTACAAACACAATAACAFC-3′). PCR for A+T-rich was performed on a MyCycler thermal cycler system (Bio-Rad) with initial denaturation at 95 °C for 2 min, followed by 35 cycles of denaturing at 94 °C for 30 s, annealing at 52 °C for 30 s, extension at 72 °C for 90 s, and a final extension at 72 °C for 5 min. A fragment of COII was amplified by PCR using the primer set TL2-J-3037 (5′-ATGGCAGATTAGTGCAATGG-3′) and TK-N-3785 (5′-GTTTAAGAGACCAGTACTTG-3′) (70). The PCR for COII consisted of initial denaturation at 94 °C for 1 min, followed by 35 cycles of denaturing at 94 °C for 30 s, annealing at 50 °C for 30 s, extension at 72 °C for 40 s, and a final extension at 72 °C for 3 min. PCR products for A+T-rich and COII were sequenced in both directions in a commercial sequencing facility (Macrogen Inc.), and forward and reverse chromatograms were edited using BioEdit 7.0.4.1 (71) and resulted in a 144 nucleotide sequence and a 702 nucleotide sequence, respectively. The A+T-rich and COII sequences obtained in this study were deposited in the DNA Data Bank of Japan/European Molecular Biology Laboratory/GenBank nucleotide sequence databases under accession nos. MT387011–MT387036.

Divergence time was estimated based on the A+T-rich sequences and intraspecific rates of mitochondrial evolution (39).

To count the number of micropyles of eggs, we used all collected eggs of two field colonies of asexual populations (7 eggs of a colony from Tokushima [colony code: TO150911B] and 8 eggs of a colony from Saiki [colony code: SK150715A]) and those of sexual populations (5 eggs of a colony from Kushimoto [colony code: SN150430C] and 10 eggs of a colony from Amami-Oshima Island [colony code: NZ150526A]) of G. nakajimai. The number of micropyles of eggs was counted under scanning electron microscope (VE-8800; Keyence). To compare the numbers of micropyles of eggs between the asexual and sexual G. nakajimai (Dataset S1D), we used nested ANOVA followed by Tukey’s HSD test.

To investigate the mode of reproduction in the asexual G. nakajimai, we genotyped the primary queens and larvae in the laboratory-founded colonies. Virgin female alates were obtained from a colony collected from one of the asexual populations, Ashizuri (colony code: AS141111C). Colonies of asexual populations are founded by more than two female alates (young queens) (22), probably due to the necessity of grooming partners that would be essential to survive in a pathogen-rich environment because termites cannot clean the whole of their body through self-grooming (67). Therefore, two virgin female alates were randomly chosen from the colony and placed in 35-mm Petri dish that contained layers of a filter paper and two damp chips of Oregon pine wood (22.5 × 22.5 × 4 mm), as described in a previous study (72). This procedure was replicated 20 times. The laboratory-founded colonies were kept at 25 °C under constant darkness for 500 d. Although 17 of 20 laboratory-founded colonies could not survive for 500 d, two queens and one larva (all surviving individuals) were obtained from each of the remaining three laboratory colonies (I to III). The six queens and three larvae were genotyped as detailed in the subsection Genome-Wide SNP Analyses above. Using the SNP data, we calculated the percentage of SNP identity between individuals. The percentage of SNP identity between an offspring (larva) and its two possible mothers (queens) was compared, and the queen with the highest genetic similarity to an individual larva was the inferred mother. The percentage of heterozygous SNPs in an individual larva (for the SNPs where the inferred mother was heterozygous) was calculated and then the observed proportion of heterozygosity in the offspring was compared with the expected proportion of heterozygosity in candidate modes of asexual reproduction.

Based on the evidence of the SNP analyses that the asexual G. nakajimai contains two lineages (i.e., AL1 and AL2), we further conducted an additional PCoA for individuals both from the field colonies and the laboratory-founded colonies using the R package “dartR v0.93” (55) to determine whether individuals of the laboratory-founded colonies belong to AL1 or AL2.

To investigate the possibility of hybridization between the asexual and sexual G. nakajimai, we performed a crossbreeding experiment. Virgin alates were obtained from two stock colonies of asexual populations collected in Muroto (colony code: MR150910D) and Sata (colony code: ST160304C) and of sexual populations collected in Kushimoto (colony code: IZ150430A) and Ogasawara Islands (colony code: HH151016D), separated by sex before swarming began, and maintained in 90-mm Petri dishes containing moist unwoven clothes until they shed their wings (i.e., dealates). Then, individual dealates were randomly chosen from each colony and assigned to either pairs of a female from an asexual population and a male from a sexual population (FM pairs) or pairs of females from an asexual population (FF pairs), where FM pairs consisted of four different combinations (F_MR150910DM_IZ150430A, F_MR150910DM_HH151016D, F_ST160304CM_IZ150430A, and F_ST160304CM_HH151016D) and FF pairs consisted of two different combinations (F_MR150910DF_MR150910D and F_ST160304CF_ST160304C). Each combination was replicated 10 times. Pairs were placed in a 52 × 76-mm glass cell that contained mixed sawdust bait blocks, as described in a previous study (67). The glass-cell colonies were kept at 25 °C under constant darkness for 100 d. We counted eggs and larvae by checking the glass-cell colonies every 3 d. The hatching success, calculated as percentage of eggs hatched within 100 d after colony foundation, was compared among eggs of glass-cell colonies founded by FM pairs and those of glass-cell colonies founded by FF pairs using Fisher’s exact probability tests with sequential 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. In addition, we genotyped the reproductives (primary queens and kings) and newborn larvae in glass-cell colonies founded by FM pairs at two polymorphic microsatellite loci (Gly8 and Gly18) as described previously (22), and data for asexual offspring in the colonies of FM pairs were excluded from the analysis (SI Appendix, Table S5). Moreover, because there were no significant differences between the combinations and between the glass-cell colonies within pair types (i.e., FM pairs and FF pairs), respectively (P > 0.05, Fisher’s exact probability test with sequential Bonferroni correction [Statistica 10; StatSoft]), we pooled the data for both the combinations and the glass-cell colonies of each pair type and compared the hatching success between eggs of FM pairs (i.e., hybrid eggs of the asexual and sexual G. nakajimai) and those of FF pairs (i.e., unfertilized eggs of the asexual G. nakajimai) (Dataset S1C) using Fisher’s exact probability tests.