Introduction

Reproductive division of labor is a hallmark of social insects, where individuals follow different developmental trajectories resulting in distinct morphological castes (Wilson 1971). Increasing evidence in the last decade for heritable influences on division of labor has put an end to the assumption that social insect broods are fully totipotent and environmental factors alone determine caste (Smith et al. 2008; Schwander et al. 2010). The relative influence of environmental and heritable effects (genetic or epigenetic effects) on caste can vary among species along a continuum from strict environmental caste determination to almost pure genetic caste determination.

Termite colonies are typically founded by a monogamous pair of primary reproductives (one king and one queen) derived from alates (winged adults). In Reticulitermes termites, larvae differentiate into either the nymph-alate pathway (sexual pathway) or the worker pathway (neuter pathway; fig. 1). Most nymphs develop into alates and fly out to establish new colonies. Neotenics (also called secondary reproductives or supplementary reproductives) are produced from within the colony on the death of the primary king and/or queen and take over reproduction in the colony (Vargo and Husseneder 2009). Neotenics can develop either from nymphs to become “nymphoid” reproductives or from workers to become “ergatoid” reproductives (fig. 1). In the termite Reticulitermes speratus and other related species, the parental phenotypes (worker- or nymph-derived reproductives) influence the caste fate of the offspring (Hayashi et al. 2007; Kitade et al. 2010). These laboratory-based studies demonstrated that termite caste determination is not purely environmental but involves some heritable components. As a candidate mechanism by which parental phenotypes influence caste determination of offspring, Hayashi et al. (2007) proposed a genetically influenced caste determination (GCD) model (fig. S1). The GCD model assumes the existence of an x-linked caste-determining locus wk; females have two x chromosomes and therefore three possible genotypes at wk (wk^AA, wk^AB, or wk^BB), and males have one x and one y chromosome (wk^AY or wk^BY; fig. S1). Offspring with wk^AB and wk^AY genotypes develop into workers, and offspring with genotypes wk^AA and wk^BY develop into nymphs. The genotype wk^BB is lethal. The important prerequisite of this GCD model is that ergatoid (worker-derived) reproductives are very common and that breeding by ergatoid reproductives (wk^AB and/or wk^AY) is required for nymph production and therefore alate production (fig. 2a).

In contrast to the GCD model, by directly genotyping the royals (kings and queens) and other members of mature R. speratus colonies collected in the field, Matsuura et al. (2009) demonstrated that queens use parthenogenesis to produce nymphoid (nymph-derived) neotenic queens, which supplement and then replace the primary queens. In this reproductive system, named asexual queen succession (AQS; fig. 2b), parthenogenetically produced individuals are destined to develop into neotenics and both the primary queens and the neotenic queens produce other colony members through sexual reproduction (Matsuura et al. 2009; Matsuura 2017).

It is apparent that AQS is contradictory to any genetic models including the GCD model proposed by Hayashi et al. (2007; fig. 2). Genetic models assume a genotype-phenotype linkage; that is, nymph-derived king and queen pairs (carrying nymph genotype[s]) produce worker genotype(s) (and thus worker phenotype). To maintain this genotype-phenotype linkage, the nymph-derived king and queen pairs must be replaced by worker-derived reproductives before producing nymphs. In the AQS system, however, both workers and nymphs are produced by the same combination of parental phenotypes—that is, nymph-derived kings and queens. Alate production by AQS colonies breaks the genotype-phenotype linkage, which is required for the maintenance of the GCD system (fig. 2). The discrepancy between genetic models and the AQS system can be explained by one of two possibilities: (1) there is geographic variation in AQS, where AQS can be found in a limited area and colonies of non-AQS populations have a GCD system; (2) the heritable influence of parental phenotypes is not due to genetic effects but results from other heritable factors.

In the latter case, a new epigenetic model can explain the heritable influence of parental phenotypes in termites that is consistent with AQS. Studies investigating heritable effects on caste determination might have largely overlooked epigenetic inheritance, although transmission to offspring of factors other than DNA sequences including epigenetic states can also affect offspring phenotype (Chong and Whitelaw 2004; Bonduriansky and Day 2008; Jablonka and Raz 2009). Recent studies of genomic imprinting show that genes can be differentially marked in egg and sperm and that inheritance of these epigenetic marks cause genes to be expressed in a parental origin–specific manner in the offspring (Sha 2008; Ferguson-Smith 2011). We here develop a genomic imprinting model, in which queen- and king-specific epigenetic marks influence caste determination of offspring, to understand the heritable influence of parental phenotypes and AQS in a single framework.

In this study, we first conducted wide geographic sampling of R. speratus colonies throughout its range in Japan and performed genotyping of royals and other members of each colony to identify kin structure. Genetic models predict that nymph-producing colonies have worker-derived reproductives, while our genomic imprinting model allows colonies with nymph-derived reproductives to produce both workers and nymphs. We found that R. speratus colonies are AQS throughout its range, where almost all nymph-producing colonies have only nymph-derived kings and queens. Therefore, instead of genetic models, we developed a genome imprinting model that accounts for all known laboratory experimental data on caste differentiation in Reticulitermes species as well as their reproductive systems in nature. We also examined caste development of offspring under the presence of founder males and females in incipient colonies to determine the relative importance of genomic imprinting and environmental factors in caste determination in R. speratus.

Material and Methods

Genomic Imprinting Model

We developed a novel genomic imprinting model for caste predisposition in R. speratus and other Reticulitermes termites that closely matches empirical data from both lab and field studies. Sex-specific epigenetic marks (hereafter called “epimarks”) of the parents can influence the phenotype of the offspring through incomplete erasure of the epimarks in the germ line (Rice et al. 2016). Our genomic imprinting model assumes that heritable epimarks (i.e., imprinting) increase in strength with the sexual maturation of the parents (fig. 3a). The sex- and phenotype-specific imprinting, which can influence the expression of both sex and growth regulatory genes, can be inherited, resulting in parent-specific gene expression in offspring (fig. 3b). The sex-specific imprinting effect of a parent of a certain phenotype on the development of offspring is given by

$$\begin{array}{}\text{(1)}& {\mathit{I}}_{\mathrm{sex},\mathrm{phenotype}}={\delta }_{\mathrm{phenotype}}\left(s_{\mathrm{sex}},g_{\mathrm{sex}}\right),\end{array}$$

where sex is either male (m) or female (f) and phenotype is either ergatoid (E), nymphoid (N), or alate-derived primary reproductive (A). The values s_sex and g_sex are imprinting effects on the expression of sex regulatory genes and growth regulatory genes, respectively. The value δ_phenotype is the relative strength of imprinting effects of each parental phenotype. In the model, caste phenotype of the parents influences the imprinting level; nymphoid and primary reproductives, which develop from the reproductive pathway, obtain higher imprinting levels than ergatoids, which develop from the neuter pathway (figs. 1, 3a). Therefore, the values of δ_phenotype can be predicted to be $0<{\delta }_{\mathrm{E}}<{\delta }_{\mathrm{N}}={\delta }_{\mathrm{A}}=1$, although the range of δ_E is not theoretically constrained to be less than 1. The effects of maternal imprinting by a nymphoid queen I_f,N and a primary queen I_f,A on the development of offspring are given by

$$\begin{array}{}\text{(2)}& {\mathit{I}}_{\mathrm{f},\mathrm{N}}={\mathit{I}}_{\mathrm{f},\mathrm{A}}=\left(s_{\mathrm{f}},g_{\mathrm{f}}\right).\end{array}$$

Because paternal imprinting acts antagonistically toward maternal imprinting with respect to the expression of sex regulatory genes, the effects of paternal imprinting by a nymphoid king I_m,N and a primary king I_m,A are

$$\begin{array}{}\text{(3)}& {\mathit{I}}_{\mathrm{m},\mathrm{N}}={\mathit{I}}_{\mathrm{m},\mathrm{A}}=\left(-s_{\mathrm{m}},g_{\mathrm{m}}\right).\end{array}$$

The effects of imprinting by an ergatoid queen I_f,E and an ergatoid king I_m,E are

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{I}}_{\mathrm{f,E}}=& {\delta }_{\mathrm{E}}\left(s_{\mathrm{f}},g_{\mathrm{f}}\right),\\ {\mathit{I}}_{\mathrm{m,E}}=& {\delta }_{\mathrm{E}}\left(-s_{\mathrm{m}},g_{\mathrm{m}}\right).\end{array}\end{array}$$

Let $s_{\mathrm{m}}={\mathit{ks}}_{\mathrm{f}}$ and $g_{\mathrm{m}}={\mathit{kg}}_{\mathrm{f}}$, where k is the relative strength of paternal imprinting to maternal imprinting, and let $g_{\mathrm{f}}={\mathit{rs}}_{\mathrm{f}}$, where r is the relative strength of parental imprinting of growth factors to sexual factors. The effect of parental imprinting on the development of offspring is expressed as the sum of maternal and paternal imprinting. In the case of mating between a female nymphoid (fN) and a male nymphoid (mN), the effect of genomic imprinting I_fNmN is

$$\begin{array}{}\text{(5)}& {\mathit{I}}_{\mathrm{fNmN}}={\mathit{I}}_{\mathrm{f},\mathrm{N}}+{\mathit{I}}_{\mathrm{m},\mathrm{N}}=\left(s_{\mathrm{f}},g_{\mathrm{f}}\right)+\left(-s_{\mathrm{m}},g_{\mathrm{m}}\right)=\left\{\left(1-k\right)s_{\mathrm{f}},\left(1+k\right){\mathit{rs}}_{\mathrm{f}}\right\}\end{array}$$

and that of mating between a female nymphoid (fN) and a male ergatoid (mE) is

$$\begin{array}{}\text{(6)}& {\mathit{I}}_{\mathrm{fNmE}}={\mathit{I}}_{\mathrm{f},\mathrm{N}}+{\mathit{I}}_{\mathrm{m},\mathrm{E}}=\left(s_{\mathrm{f}},g_{\mathrm{f}}\right)+{\delta }_{\mathrm{E}}\left(-s_{\mathrm{m}},g_{\mathrm{m}}\right)=\left\{\left(1-{\delta }_{\mathrm{E}}k\right)s_{\mathrm{f}},\left(1+{\delta }_{\mathrm{E}}k\right){\mathit{rs}}_{\mathrm{f}}\right\}.\end{array}$$

In AQS species, in which females are capable of diploid thelytoky (parthenogenetic production of diploid female offspring), the effects of genomic imprinting on the parthenogenetic daughters produced by a female nymphoid I_fNfN and a female ergatoid I_fEfE are

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{I}}_{\mathrm{fNfN}}=& 2{\mathit{I}}_{\mathrm{f,N}}=2\left(s_{\mathrm{f}},g_{\mathrm{f}}\right)=2\left(s_{\mathrm{f}},{\mathit{rs}}_{\mathrm{f}}\right),\\ {\mathit{I}}_{\mathrm{fEfE}}=& 2{\mathit{I}}_{\mathrm{f,E}}=2{\delta }_{\mathrm{E}}\left(s_{\mathrm{f}},g_{\mathrm{f}}\right)=2{\delta }_{\mathrm{E}}\left(s_{\mathrm{f}},{\mathit{rs}}_{\mathrm{f}}\right).\end{array}\end{array}$$

Caste fate of offspring is determined by the relative expression of growth regulatory genes (ε_g: growth factors) to sex regulatory genes (ε_s: sexual factors); growth without sexual development channels individuals to the neuter caste (i.e., workers), and growth with sexual development directs individuals to the reproductive caste (i.e., nymphs; fig 3b). This can be modeled such that offspring with ε_g and ε_s values below the threshold line (male: ${\epsilon }_{\mathrm{g}}=-{\epsilon }_{\mathrm{s}}$; female: ${\epsilon }_{\mathrm{g}}={\epsilon }_{\mathrm{s}}$) develop into nymphs (fig. 3b).

Here, the effect of genomic imprinting on the sex and growth regulatory genes is denoted by $\mathit{I}=\left(i_{\mathrm{s}},i_{\mathrm{g}}\right)$. In the model, the developmental pathway of individuals is primarily determined by imprinting but also subjected to random variation due to developmental noise, that is, phenotypic variation among individuals with identical genotypes living in an identical environment (Scheiner et al. 1991). We incorporated the developmental noise into the model as follows: gene expression levels have uniform distribution within a circular region ($\mathrm{radius}=e$) with center I on the ${\epsilon }_{\mathrm{s}}-{\epsilon }_{\mathrm{g}}$ plane (fig. 3b). Thus, the proportion of female offspring with imprinting $\mathit{I}=\left(i_{\mathrm{s}},i_{\mathrm{g}}\right)$ destined to develop into nymphs is

$$\begin{array}{}\text{(8)}& {\phi }_{\mathrm{f}}\left[\mathit{I},e\right]=P\left\{{\epsilon }_g<{\epsilon }_s|{\left({\epsilon }_s-i_s\right)}^2+{\left({\epsilon }_g-i_{\mathrm{g}}\right)}^2<e^2\right\}.\end{array}$$

This can be calculated by

$$\begin{array}{}\text{(9)}& {\phi }_{\mathrm{f}}\left[\mathit{I},e\right]=\{\begin{array}{cl}0& (d\ge e,i_g\ge i_s)\\ A_s/A_c& (d<e,i_g\ge i_s)\\ 1-A_s/A_c& (d<e,i_g<i_s)\\ 1& \left(\text{otherwise}\right)\end{array},\end{array}$$

where d is the distance from the center I to the threshold line (${\epsilon }_{\mathrm{g}}={\epsilon }_{\mathrm{s}}$), A_c is the area of the circle with center I and radius e, and A_s is the area of the circular segment bounded by the circle and the threshold line.

Similarly, the proportion of male offspring with imprinting $\mathit{I}=(i_{\mathrm{s}},i_{\mathrm{g}})$ destined to develop into nymphs is

$$\begin{array}{}\text{(10)}& {\phi }_{\mathrm{m}}\left[\mathit{I},e\right]=P\{{\epsilon }_g<-{\epsilon }_s|{({\epsilon }_s-i_s)}^2+{({\epsilon }_g-i_g)}^2<e^2\},\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\phi }_{\mathrm{m}}[\mathit{I},e]=\{\begin{array}{cl}0& (d\ge e,i_g\ge -i_s)\\ A_s/A_c& (d<e,i_g\ge -i_s)\\ 1-A_s/A_c& (d<e,i_g<-i_s)\\ 1& \left(\mathrm{otherwise}\right)\end{array}.\end{array}$$

We can set e constant (=1) since the magnitude of the noise is determined in relation to the strength of imprinting.

We tested whether our model can explain the observed caste differentiation patterns of offspring with various combinations of parental phenotypes (nymphoid or ergatoid) in Reticulitermes species (R. speratus [Hayashi et al. 2007], Reticulitermes kanmonensis, Reticulitermes okinawanus, and Reticulitermes yaeyamanus [Kitade et al. 2010]). We explored the parameter sets with the greatest explanatory power by calculating root mean square errors for the rate of nymphs produced by each combination of parental phenotypes. The parameter ranges are shown in table 1. We present a sensitivity analysis of the model, testing the effects of the s_f, δ, k, and r values on the caste fate of offspring (figs. S2–S5).

			Species
Symbol	Definition	Range	R. spe	R. kan	R. yae	R. oki
sf	Maternal imprinting effects on sex regulatory genes	1.0–10.0	1.9	10.0	8.8	5.6
r	Relative strength of parental imprinting of growth factors to sexual factors	.01–1.00	.27	.88	.87	.98
k	Relative strength of paternal imprinting to maternal imprinting	.10–10.0	.85	.13	.19	.70
δE	Relative strength of imprinting effects of ergatoid to nymphoid	.01–1.00	.32	.02	.02	.09
RMSE	…	…	.151	.045	.049	.025

Note. We calculated the probabilities of individuals developing into nymphs across the range of parameters to explore parameter sets that can explain the empirical data of each species. The most explainable parameter sets are determined by using root mean square errors (RMSE) for the proportion of nymphs produced by each combination of parental phenotypes. Ranges and most fitted parameters are shown. R. kan = Reticulitermes kanmonensis; R. oki = Reticulitermes okinawanus; R. spe = Reticulitermes speratus; R. yae = Reticulitermes yaeyamanus.

View Table Image

During the colony-founding stage, the founder males and females strongly inhibit their first brood from differentiation into nymphs and encourage worker differentiation to obtain an initial labor force (Oster and Wilson 1978). We incorporated the effects of the founder males and females on offspring into our genomic imprinting model. The effects of queen presence φ_f,A and king presence φ_m,A, likely through primer pheromones, are

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathbf{\phi }}_{\mathrm{f,A}}=& (-s_{\mathrm{P}},g_{\mathrm{P}})=\eta (-s_{\mathrm{f}},{\mathit{rs}}_{\mathrm{f}}),\\ {\mathbf{\phi }}_{\mathrm{m,A}}=& (s_{\mathrm{P}},g_{\mathrm{P}})=\eta (s_{\mathrm{f}},{\mathit{rs}}_{\mathrm{f}}).\end{array}\end{array}$$

Here, we set the relative strength of s_P and g_P to s_f and g_f as $\eta (\eta =s_{\mathrm{P}}/s_{\mathrm{f}}=g_{\mathrm{P}}/g_{\mathrm{f}},0<\eta <1)$. These effects on offspring are sex specific; that is, φ_f,A acts only on female offspring and φ_m,A acts only on male offspring. In a female-female (FF) colony, the effect of queen presence on female offspring is 2φ_f,A. These effects are normally moderated after the emergence of first-brood workers. We predicted caste differentiation patterns with the effect of parent presence using the parameters determined for R. speratus in the above calculations. All calculations were performed using R, version 3.1.3 (https://cran.r-project.org).

Results

Rejection of Genetic Influence Hypothesis

Our wide-area assessment showed that Reticulitermes speratus exhibits no geographic variation in its reproductive system (fig. 4a). We performed microsatellite analysis on the 20 representative colonies collected from throughout its range in Japan to determine parentage of the neotenic queens and other colony members. Our genetic data showed that all R. speratus colonies collected in all areas exhibited AQS; that is, neotenic queens are produced by parthenogenesis of primary queens (fig. 4a; S2 data set). Among these 20 colonies, all colonies had parthenogenetically produced secondary queens, where nearly all of the secondary queens (158 of 165) were produced by parthenogenesis (see the S2 data set for raw genotyping data). Therefore, it is clear that geographic variation cannot explain the inconsistency of the genetic model (Hayashi et al. 2007) and AQS (Matsuura et al. 2009).

We found that among the 162 kings collected from 114 mature colonies in the field, 104 were primary and 58 were nymphoid, while no ergatoid king was found in the field. Among the 6,824 queens, 6,812 (99.82%) were nymphoid and six were primary (alate-derived), while only six (0.088%) were ergatoid (fig. 4b; S1 data set). Moreover, all of these mature colonies had nymphs. These results are clearly inconsistent with the prediction of the genetic influence hypothesis that nymph-producing colonies have ergatoid reproductives.

For example, the GCD model (Hayashi et al. 2007) assumes a link between genotype and phenotype so that offspring with wk^AB and wk^AY genotypes develop into workers and offspring with genotypes wk^AA and wk^BY develop into nymphs (fig. S1). The important prerequisite of this GCD model is that primary reproductives must be replaced by ergatoid (worker-derived) reproductives (wk^AB and/or wk^AY) before nymph and alate production (fig. 2a). Our field data refuted this assumption as described above. Hence, the GCD model cannot explain the caste system of R. speratus. This requires an alternative model to explain why parental phenotypes (nymphoid or ergatoid) strongly influence the caste fate of offspring (Hayashi et al. 2007; Kitade et al. 2010).

Genomic Imprinting Model Results

Sensitivity analysis of our genomic imprinting model showed that each of the parameters (s_f, r, k, and δ_E) influences the caste fate of the offspring depending on the combination of parental phenotypes (figs. S2–S5) and found the parameter combinations with the greatest explanatory power for the empirical data (table 1). The prediction of our genomic imprinting model matched the empirical data for caste differentiation patterns in the AQS species R. speratus (Hayashi et al. 2007) and the non-AQS species Reticulitermes kanmonensis, Reticulitermes yaeyamanus, and Reticulitermes okinawanus (Kitade et al. 2010; fig. 5a) when the model parameters are allowed to vary among species (table 1). In addition, the GCD model does not explain the differentiation pattern of R. okinawanus (Kitade et al. 2010), whereas our genomic imprinting model successfully accounts for the observed pattern.

Discussion

We conclude that our genomic imprinting model is a viable model to explain the influence of parental phenotypes on the caste predisposition of the offspring, while the assumption of genetic models was clearly rejected by our field data. In Reticulitermes species, both parental phenotype and environment affect caste differentiation patterns. To harmoniously incorporate these factors, our model assumes that the expression of sexual regulatory genes and growth regulatory genes in offspring are affected by heritable factors (i.e., parent sex– and phenotype-specific imprinting) and by socio-environmental factors (e.g., pheromones and nutritional condition). In our model, the relative expression of sexual and growth regulatory genes determines the developmental pathway of individuals. Males and females differentiated into either the worker pathway or the nymph pathway and establish their de novo phenotype-specific epimarks, which can be transgenerationally inherited by their offspring through eggs and sperm (fig. 3a). Therefore, our model can reasonably incorporate both heritable factors and environmental factors into the life cycle of these termites and matches the reproductive system in nature as well as the results of laboratory experiments.

The presumed evolutionary pathway leading to our genomic imprinting system is plausible and parsimonious. With the development of reproductive division of labor, sexual dimorphism of the reproductive castes should increase; that is, kings and queens should have greater differences in their sexual traits so as to be specialized in sperm and egg production, respectively. Development of sexual dimorphism requires sex-specific gene expression, and sexual dimorphism is often linked to the establishment of sex-specific epigenetic marks, which are involved in canalizing sexual development (D’Ávila et al. 2010; Ferguson-Smith 2011; Glastad et al. 2016). Sex-specific epimarks would contribute to king- and queen-specific gene expression and thus increase reproductive output. Incomplete erasure of these epimarks in the germ line would lead to transgenerational carryover of the sex-specific epimarks, which could act antagonistically in the offspring influencing their sexual development (Rice et al. 2016). It is unnecessary to invoke novel mechanisms for the evolution of our genomic imprinting system. Hence, our model matches well with the reproductive systems of a broad range of termite species.

Note that recent studies identified three AQS species in higher termites (Termitidae), which are phylogenetically distant from the lower termite genus Reticulitermes (Fougeyrollas et al. 2015, 2017; Fournier et al. 2016). The termitid species with AQS exhibit different mechanisms of parthenogenesis: Cavitermes tuberosus uses automixis with terminal fusion (Fournier et al. 2016), resulting in homozygous female parthenogens as in Reticulitermes, whereas Embiratermes neotenicus and Silvestritermes minutus use automixis with central fusion, resulting in heterozygous parthenogens (Fougeyrollas et al. 2015, 2017). Although there is no genetic difference between the founder queen, which is derived from an alate, and her parthenogenetic daughters in the latter species, the parthenogenetic daughters (carrying only maternal chromosomes) also have a higher propensity to develop into neotenic queens than sexually produced daughters (carrying paternal and maternal chromosomes). Thus, our genomic imprinting model also matches the situation in the higher termite species. In addition, only our genomic imprinting model can reasonably explain the evolution of AQS in termites. The origin of AQS requires the simultaneous evolution of two components: the capacity for parthenogenesis and the higher propensity of parthenogens to develop into neotenic queens than sexually produced daughters (Matsuura 2017). Our model provides a mechanism for the simultaneous evolution of parthenogenesis and the developmental caste bias of parthenogens.

The basic idea of our genomic imprinting model, that maternal and paternal imprinting act on the offspring’s sexual development in different directions, is in accordance with the sexual antagonism theory of genomic imprinting (Day and Bonduriansky 2004; Rice et al. 2016). An example of the sexual antagonism of parental imprinting can be seen in the genomic imprinting sex determination mechanism in some haplodiploid insects (Dobson and Tanouye 1998; Gadagkar 2000). One novel prediction of genomic imprinting is that the neuter termite worker is generated by the counterbalanced actions of maternal and paternal imprinting on the offspring’s sexual development. Another leading evolutionary explanation for genomic imprinting is kinship theory, which predicts that differences in inclusive fitness between mothers and fathers can lead to conflicts between maternal and paternal genes in their offspring (Haig 2000; Queller 2003; Burt and Trivers 2006; Kronauer 2008; Drewell et al. 2012; Kocher et al. 2015; Galbraith et al. 2016). Our genomic imprinting model and the kinship theory of genomic imprinting are not mutually exclusive since the former deals with the mechanism of caste determination and the latter concerns evolutionary dynamics. However, it is important to note that our genomic imprinting model can explain the observed caste predisposition in termites without assuming any kin conflict between kings and queens.

Our genomic imprinting model provides qualitative predictions for future studies of the molecular mechanisms by which parental phenotypes influence the caste fate of offspring. Our prediction of the genomic imprinting parameters indicates that in all of the study species the strength of maternal imprinting should be greater than paternal imprinting ($k<1$; table 1), suggesting that females have higher levels of epimarks (e.g., higher DNA methylation or histone modification). Similarly, nymphoids have stronger imprinting effects on the offspring than ergatoids (${\delta }_{\mathrm{E}}<1$; table 1). In this article, we proposed a simplified model to explicitly demonstrate the influence of genome imprinting on caste determination. To obtain quantitative predictions for interspecific comparisons of parameter values, however, we would need to explore other factors that may affect caste differentiation. For instance, the effect of phenotypes on imprinting strength might differ between sexes; that is, δ can be different between males and females. The parameter r, the relative strength of parental imprinting of growth factors to sexual factors, can also be different between sexes or among phenotypes. In AQS species the original queen is replaced by multiple parthenogenetically produced daughter queens; that is, queens are genetically immortal. It is known that sex-asymmetric inbreeding leads to female-biased sex allocation in AQS species (Kobayashi et al. 2013), although the mechanism regulating alate sex ratios remains to be explored. Our model implies that genomic imprinting influences nymph differentiation of male and female offspring, which can act as an important factor regulating sex allocation. Further studies enabling us to assign more precise quantitative values to the parameters would contribute to our understanding of the mechanisms underlying the evolution of life-history traits including sex ratios.

Our study, together with the growing evidence of nongenetic inheritance or parent-of-origin effect in many contexts in various organisms including the eusocial Hymenoptera (Oldroyd et al. 2014; Kocher et al. 2015; Galbraith et al. 2016), calls for a reevaluation of existing data on genetic caste determination from the viewpoint of genomic imprinting effects in social insects. Our findings lay the groundwork for further investigation of the role of genomic imprinting on the developmental fate of offspring across the social insects.

Associate Editor: Madeleine Beekman

Editor: Yannis Michalakis