Canalized gene expression during development mediates caste differentiation in ants

Abstract

Ant colonies are higher-level organisms consisting of specialized reproductive and non-reproductive individuals that differentiate early in development, similar to germ-soma segregation in bilateral Metazoa. Analogous to diverging cell lines, developmental differentiation of individual ants has often been considered in epigenetic terms but the sets of genes that determine caste phenotypes throughout larval and pupal development remain unknown. Here, we reconstruct the individual developmental trajectories of two ant species, Monomorium pharaonis and Acromyrmex echinatior, after obtaining >1,400 whole-genome transcriptomes. Using a new backward prediction algorithm, we show that caste phenotypes can be accurately predicted by genome-wide transcriptome profiling. We find that caste differentiation is increasingly canalized from early development onwards, particularly in germline individuals (gynes/queens) and that the juvenile hormone signalling pathway plays a key role in this process by regulating body mass divergence between castes. We quantified gene-specific canalization levels and found that canalized genes with gyne/queen-biased expression were enriched for ovary and wing functions while canalized genes with worker-biased expression were enriched in brain and behavioural functions. Suppression in gyne larvae of Freja, a highly canalized gyne-biased ovary gene, disturbed pupal development by inducing non-adaptive intermediate phenotypes between gynes and workers. Our results are consistent with natural selection actively maintaining canalized caste phenotypes while securing robustness in the life cycle ontogeny of ant colonies.

Fig. 1. Homologous caste differentiation trajectories in two ant species.

a, M. pharaonis developmental trajectories reconstructed from 568 individual transcriptomes covering all developmental stages and visualized with the Fruchterman–Reingold algorithm in an undirected network (Methods), so that individuals with similar transcriptomic profiles cluster together. Shading of connecting lines reflects the strength of mutual transcriptome correlations, ranging from 0.8 (light grey) to 1.0 (black). Symbol colours differentiate between embryos within eggs (white and blue), larvae (yellow), pupae (red) and adults (brown). Images courtesy of L.P.. b, The first three PC axes for individual transcriptomes of A. echinatior and M. pharaonis, from first instar larvae to adults, constructed from 979 individual transcriptomes across the two ant species, using 7,838 one-to-one orthologous genes. Upper panel: probability density function (PDF) of PC1 values separating the two ant species. Lower panel: PC2 and PC3 jointly distinguish between developmental stages and caste phenotypes across individual transcriptomes, plotted separately for the two ant species. Lines connect the median PC values for each caste across development stages, showing that individuals follow very similar trajectories regardless of species identity. Symbol colours and shapes as in a. c, Between-species transcriptome similarities comparing M. pharaonis (n_gynes = 171; n_workers = 154) and A. echinatior (n_gynes = 139; n_workers = 319) (left) and M. pharaonis and D. melanogaster females (n = 123) (right), based on stage-specific Spearman correlation coefficients and plotted separately for gynes (red) and workers (blue). Between-species similarities peaked in third instar larvae where gynes and workers were similar but similarities were always lower in earlier and later stages, particularly for workers. The box plot shows the median (centre line), 25% and 75% quartiles (boxes), outermost values (whiskers) and data points (overlapping with boxes and whiskers).

Fig. 2. BPA predicts individual caste phenotypes independent of external morphological traits.

a, BPA predictions of early caste identity in M. pharaonis, showing that the first two PC projections from first instar larvae (n = 54) match reproductives (an unknown mix of gynes and males) and workers among second instar larvae (n = 66). Colours of first instar sample symbols reflect the predicted probability to be a reproductive individual (left panel) and identify individuals with morphologically validated caste and sex in second instar larvae while also visualizing body length (right panel). b, HCR-FISH staining for two of the best predictor genes of caste in first instar M. pharaonis larvae based on presence (+) or absence (−) of gonad tissue, that is LOC105839887 (red; expressed in fat bodies) and Smyd3 (purple; expressed in gonads, white arrow), indicating that transcriptome-wide BPA assignments as gynes were correct because these two genes can only be detected in individuals with a germline (left Gonad+ panel). Germline presence was independently checked by vasa expression (green), showing that these transcripts were always undetectable in individuals without a germline (right Gonad− panel). The housekeeping gene Tubulin (yellow) was stained as a positive control (right images in both panels). c, BPA predictions of early caste identity in A. echinatior larvae, showing that the first two PC projections from second instar larvae (n = 84) matched the third instar individuals (n = 67) with known caste morphology. Individual samples are coloured according to their predicted probability of being gyne or worker in the second instar (left panel) while visualizing individual body lengths via the symbol diameters. Known gynes and workers in the third instar are represented by triangles and squares (right panel) confirming BPA segregation. The inserted epifluorescence microscope hair images refer to the ventral thorax region of a typical predicted second instar gyne (left) or worker (right).

Fig. 3. Individual transcriptomes quantify caste canalization and its regulation by the JH signalling pathway.

a, Developmental potential scores (−1 < Δ<1) reflecting caste commitment for gyne (positive score) and worker (negative score) phenotypes across developmental stages of M. pharaonis as body length increases, corresponding to the representative images towards the right (courtesy L.P.). Δ values are based on normalized difference of transcriptomic distances between target individuals and an average gyne (calculated as the mean transcriptome across gynes within the same stage) and an average worker in the subsequent stage. A positive Δ value means that transcriptomic distance between a target individual and the average gyne in the next stage is less than between the target individual and the average worker in the next stage, so the target individual has a gyne-biased developmental potential. Caste identities were known from morphological characters, except for first instar larvae where they were inferred via BPA. b, Caste- and body length-specific expression of JH and ecdysone (20E) pathways, coloured according to gyne-biased (red) or worker-biased (blue) expression in third instar larvae (upper panel). The expression of jhamt and E93 showed both caste-specific (qualitative) and body size-specific (quantitative) thresholds (vertical dashed lines; estimated with a threshold regression model (Methods)) (lower panel). Note that the downregulation of jhamt (terminating the JH biosynthesis) preceded the upregulation of E93 to actively initiate metamorphosis (see Extended Data Fig. 5a for expression patterns of other genes in this pathway). c, In third instar larvae of intermediate body length, feeding with JHA increased worker pupal body length (two-sided t-test; P = 0.0004; d.f. = 16.38; t = 4.44; n = 16 for both control and treatment groups), while precocene I decreased gyne pupal body length (two-sided t-test; P = 0.0001; d.f. = 50.70; t = −4.13; n = 30 and 28 for control and treatment groups, respectively). The box plot shows the median (centre line), 25% and 75% quartiles (boxes), outermost values (whiskers) and data points (jitter plot overlapping with boxes and whiskers).

Fig. 4. Freja and other major caste canalizing genes in M. pharaonis.

a, Overall gene-specific canalization scores (left) with the ten most strongly gyne-biased (top centre) and worker-biased canalized genes (bottom centre) visualized in heatmaps and colour coded by their tissue expression in D. melanogaster (right). b, HCR-FISH staining of vasa (middle, green) and the top canalized gene, Freja (right, yellow) in an adult gyne ovariole (left), showing Freja expression in ovarian follicle cells. c, Compared to GFP-RNAi controls (black, n = 51), Freja-RNAi-treated third instar gyne larvae (red, n = 58) produced pupae with reduced antennal scape length and head width (two-sided t-tests; d.f. = 102.47 and 93.63 and t = −5.43 and −5.61, respectively; both P < 10⁻⁶; left) and more frequent abnormal wing development or complete lack of wings (16/66, 24.2%) compared to the control group (4/55, 7.3%) (P = 0.01; two-sided Fisher exact test; right). d, Compared to GFP-RNAi controls (n = 34), Freja-RNAi-treated adult gynes (n = 33) induced lower mass and number of oocytes (two-sided t-tests; d.f. = 55.31 and 57.51, t = −2.40 and −3.03 and P = 0.02 and 0.004, respectively; Extended Data Fig. 7e), impairing overall reproductive functionality. The box plot shows the median (centre line), 25% and 75% quartiles (boxes), outermost values (whiskers) and data points (overlapping with boxes and whiskers).

Extended Data Fig. 1. Developmental transcriptomes in ants.

a, Transcriptomic developmental trajectories in A. echinatior, based on Spearman rank correlation similarity in gene expression across individual transcriptomes. Trajectories were constructed and visualized with a similar procedure as in Fig. 1a, except for A. echinatior having three worker subcastes, which exhibited increasing transcriptomic divergence across the pupal and adult stages. b, Between-stage transcriptomic similarity matrix in M. pharaonis (upper panel) and A. echinatior (lower panel), based on the mean values of within-group and between-group correlation coefficients. In M. pharaonis, correlation coefficients within the same stage were always higher than between stages, and transcriptomic similarities consistently clustered by adjacent developmental stages. In A. echinatior, correlation coefficients within the same caste were higher than within the same stage, consistent with morphological caste differences being more substantial in A. echinatior.

Extended Data Fig. 2. Prediction of caste identities with BPA.

a, The Backward Progressives Algorithm (BPA) predicts caste identities in previous developmental stages, using non-validated transcriptomes at the target stage (S_t) to construct PCAs and then projects these onto the subsequent stage (S_t+1) where caste identities are known. The BPA then uses the known caste labels at S_t+1 to identify PC axes that are associated with confirmed caste identity and uses linear discriminant analysis to train a predictive model, assuming that the PC axes at S_t+1 are also associated with caste identities at S_t as expected under developmental continuity. The trained model then predicts caste identities at S_t, after which it assumes these predicted caste identities to be real and initiates a next round to predict caste identities at stage S_t-1. This process continues until the prediction likelihoods at S_t-n become too low to be informative. b, BPA predicted the sex of sampled individuals among 1^st instar (44 hour) larvae in D. melanogaster (left panel). While sex can be distinguished in 2^nd instar (50 hour) larvae (right panel) via genotyping after simultaneous DNA and RNA extraction, the biomass of 1^st instar larvae was too small to perform such simultaneous extractions. By examining the proportion of reads that mapped to the Y chromosome of Drosophila, we found that predicted males in 1^st instar larvae had a significantly higher proportion of reads mapped to the Y chromosome, confirming our prediction. Individual samples are coloured according to their predicted probability to be female or male in the 1^st-instar and symbols were sized according to the number of reads mapped to Y chromosome per million reads (YPM). c, Validating BPA on samples of developmental stages with known caste identities used individual M. pharaonis transcriptomes of individuals with distinct morphology. The table presents prediction accuracies for each target stage, calculated as the ratio of the number of correctly assigned individuals and the total number of individuals sampled at each stage, using two alternative approaches: 1. Independent: Predicting caste at each targeted stage (S_n) using the observed (true) caste labels as training stage (S_n+1) to examine the prediction accuracy when the training caste identities are in fact known from morphological information. Here, the ratio in each stage reflects the accuracy of BPA in each stage. Progressive: Predicting caste starts from the late pupal stage using adults of known caste identity as training data. Here, BPA is then performed progressively using the predicted caste labels in late pupae to predict the caste identity in early pupae. This process was repeated recurrently until the 2^nd larval instar. As the first step of BPA constructs a PCA from target stage data, we also compared the accuracies between PCAs obtained from whole transcriptomes and PCAs obtained from caste DEGs at the subsequent stage (training data). We achieved a higher prediction accuracy when PCAs were constructed with caste DEGs at the subsequent stage compared to using whole transcriptome PCAs, probably because the DEG method excluded uninformative housekeeping genes. d, Anti-body staining (VASA protein, red. RRID: AB_2893405) and in situ hybridization (vasa RNA, green) in 192-hour old embryos of M. pharaonis, showing that germline differentiation has already occurred at this stage. Among the 67 examined embryos, 18 (27%) could be documented to have no germline, indicating that it should be possible to match these presence/absence results among 192-hour embryos with BPA predictions based on 1^st instar larval transcriptomes. e, Second instar larvae of A. echinatior lack the full-body curly hairs that distinguish gynes from workers in the 3^rd larval instar, which means caste cannot be identified morphologically. We applied BPA to predict caste identities among 2^nd-instar larvae (Fig. 2a). A closer inspection showed that suspected gynes have in fact some gyne-like curly hairs, which are thicker than those in suspected workers, on their ventral thorax (arrows). These observations indicated these individuals are future gynes and were consistent with our BPA predictions. f, Among 2^nd instar larvae of M. pharaonis, PCA with whole transcriptomes showed that the overall transcriptomic difference between gynes and males was not significant (P = 0.28) while reproductive larvae of both sexes were always separated from worker larvae (P < 5e-5 and P < 5e-6 for gynes and males, respectively). This is consistent with images of 2^nd instar gyne and male larvae being indistinguishable after we used microsatellite genotyping to determine whether individuals were haploid (male) or diploid (female). Numbers of differentially expressed genes (adjusted P value < 0.05, detected with a generalized linear model that accounted for body size differences) also support this conclusion: 152 genes were differentially expressed between gyne and worker larvae, while 50 genes were differentially expressed between gyne and male larvae.

Extended Data Fig. 3. Transcriptomic canalization during caste differentiation in ants.

a, Within-stage transcriptome variation in M. pharaonis (upper panel) from 0–12 h old embryos to late pupae, plotted separately for gynes (red), workers (blue) and all individuals within each stage (black) depending on available information. The lower panel gives the same information for A. echinatior, where embryonic data were not available and 1^st instar caste phenotypes (grey) were inseparable with BPA. Transcriptome variation was quantified as 1 – r, the extent of imperfection of transcriptome-level Spearman’s correlations between a target individual and all other same-stage and same-caste individuals. Caste identities of 1^st instar individuals of M. pharaonis and 2^nd instar individuals of A. echinatior were predicted by BPA. In M. pharaonis, transcriptome variation for all individuals peaked in 36–48 h old embryos (equivalent to the gastrulation stage, 6–7, in Drosophila larvae). For both species, transcriptome variation among gynes was consistently lower than among workers. In pupal stages of M. pharaonis, transcriptome variation across all individuals exceeded transcriptome variation for the gyne and worker subsets, indicating that transcriptome differences primarily reflected realized caste differentiation, in contrast to the pattern observed across the larval stages, where the black curve was intermediate between the red and blue curves. Fourth larval instar and prepupal gyne samples of A. echinatior were excluded from this analysis, because these samples were sequenced in a different technical batch, making their transcriptome variation incomparable with the other samples. b, Developmental potential (∆) for individual gynes and workers in A. echinatior, measured as the transcriptomic distance between a focal individual and an average gyne or worker (pooling all three worker subcastes) phenotype in the next developmental stage. Developmental potential was quantified and presented as in M. pharaonis (Fig. 3a), except that all three worker subcastes were included. Caste identities for gynes and (pooled) workers in 2^nd instar larvae were predicted by BPA. As in the panel a, fourth larval instar and prepupal gyne samples were excluded to avoid a batch effect. c, PCAs of early and late pupal stage transcriptomes in M. pharaonis (early pupa gynes, n = 17; early pupa workers, n = 28; late pupa gynes, n = 30; late pupa workers, n = 22) (left) and A. echinatior (early pupa gynes, n = 18; early pupa workers, n = 47; late pupa gynes, n = 18; late pupa workers, n = 42) (right). For both species, the first PC axis (PC1) separates individual transcriptomes by developmental stage (early pupae to the left and late pupae to the right) while PC2 captures the caste-related transcriptomic variation. The overall transcriptomic difference between gynes (red) and workers (blue) increases from the early to the late pupal stage (upper panels), and the absolute values of the PC2 residuals (lower panels), representing the variation within each caste, were always lower among gynes than among workers (P < 1 × 10⁻³ for both species, two-sided t-tests). This is consistent with the mean extent of canalization being stronger in gynes than in workers. In A. echinatior, the absolute residual differences increase for the workers, consistent with A. echinatior having worker subcastes that differentiate rather late in development. Box plots show the median (centre line), 25% and 75% quartiles (boxes), outermost values (whiskers) and data points (overlapping with box and whiskers).

Extended Data Fig. 4. Early larval caste differentiation in ants.

a, Tissue-specific relative expression levels for the conserved caste-biased DEGs in early larvae, shown separately for gyne-biased (rows marked in red) and worker-biased (blue) genes. Heatmap brightness of cells reflects tissue specificity, the percentage of transcripts from targeted tissues (columns), ranging from 0% (black) to 100% (yellow). These relative abundances, based on the larval gene expression atlas of Drosophila, show that the gyne-biased DEGs in the early larval stages were mainly expressed in the midgut, fat body, and tracheal tissues, while the worker-biased DEGs were mainly expressed in the brain and central nervous system. b, Expression profiles of circadian clock-controlled protein (daywake), juvenile hormone acid O-methyltransferase-like (jhamt-like) and hexamerin among gynes and workers of the two ant species as larvae grow. All three genes are associated with the juvenile hormone signalling pathway and are significantly differentially expressed between castes in 2^nd and 3^rd instar larvae. Expression profiles are plotted against body length (log scale) to show expression dynamics as larvae grow in body length. c, DAPI staining of a representative early 3^rd instar worker larva and a representative 2^nd instar gyne larva of M. pharaonis. These animals display similar body size but wing discs (arrows) were only visible in the gyne larvae, indicating that caste determination and differentiation has already been initiated well before this early larval stage.

Extended Data Fig. 5. JH and E20 signalling pathways play key role in the regulation of canalized caste phenotypes.

a, Expression profiles of eight key regulators for insect metamorphosis that are part of the juvenile hormone and ecdysone signalling pathways (Fig. 3b), plotted against body length (log scale) of 2^nd and 3^rd instar M. pharaonis larvae. The expression levels of half of these genes (jheh2, jhamt, usp and E93) showed caste-specific body length thresholds in the 3^rd larval instar. This pattern indicates gyne and worker individuals are gated by different critical masses for entering the metamorphic molt. b, Compared to the control group (3^rd instar worker larvae fed with 10% EtOH PBS), feeding JH analogue (JHA) to 3^rd instar worker larvae delayed achieving pupation. c, JHA-fed 3^rd instar worker larvae induced inter-caste with phenotype intermediate between gyne and worker. JHA-fed workers have larger body size and developed wing buds (arrowed), however, they never developed ovaries (not show), indicating early bifurcation between colony germ–soma phenotypes. d, Compared to the control group (3^rd instar gyne larvae fed with 10% EtOH PBS), precocene I fed 3^rd instar gyne larvae were smaller and developed abnormal wings.

Extended Data Fig. 6. Canalized genes play important roles in producing adaptive caste phenotypes.

a, Tissue-specific relative expression abundances (based on Drosophila gene orthologs) among canalized genes in M. pharaonis (left panel) and A. echinatior (right panel), plotted separately for gyne-biased and worker-biased genes (M. pharaonis gyne-biased, n = 411; M. pharaonis worker-biased, n = 482; A. echinatior gyne-biased, n = 1119, and A. echinatior worker-biased, n = 899). Compared to the whole-genome background (M. pharaonis background genes, n = 8887; A. echinatior background genes, n = 7651), worker-biased canalized genes had a significantly higher relative expression in the brain, eyes, and thoracic ganglia in both ant species (one-sided t-tests; P_{M. pharaonis, brain} = 2.04 × 10⁻⁹³; P_{M. pharaonis, thoracic ganglia} = 2.29 × 10⁻⁸⁷; P_{M. pharaonis, eye} = 3.04 × 10⁻³¹; P_{A. echinatior, brain} = 9.20 × 10⁻¹²⁰; P_{A. echinatior, thoracic ganglia} = 2.73 × 10⁻¹⁸¹; P_{A. echinatior, eye} = 1.74 × 10⁻¹³⁰). In M. pharaonis, gyne-biased canalized genes had a significantly higher relative transcript abundance in the midgut, fat body and ovaries (one-sided t-tests; P_{M. pharaonis, midgut} = 7.01 × 10⁻⁵; P_{M. pharaonis, fat body} = 1.93 × 10⁻⁶; P_{M. pharaonis, ovary} = 3.04 × 10⁻¹⁰). However, in A. echinatior, gyne-biased canalized genes showed no difference with background genes for their relative transcript abundance in ovaries (one-sided t-tests; P = 0.97), consistent with the workers having retained smaller ovaries. Box plots show the median (centre line), 25% and 75% quartiles (boxes) and outermost values (whiskers); ^*P < 0.01, red for higher relative expression abundance in gyne-biased genes, blue for in worker-biased genes. b, Diagrammatic illustration of gyne-biased canalized genes being associated with traits in ovaries and wing muscles, whereas worker-biased canalized genes are associated with brain function and behaviour (see Supplementary Table 4 for full list of canalized genes). c, Canalization scores for flight related (left) and ovary specific (right) genes in the two ant species. Flight related genes were identified based on their D. melanogaster homologues associated either with flight performance itself or with striated muscle functionality (the crucial wing muscles tissue in insects). Ovary specific genes were genes having > 30% expression abundance in ovaries of D. melanogaster females compared to the sum of their expression in all tissues (see Methods). Colours of cells represents the canalization score, ranging from −3 (blue, canalized in worker-biased direction) to 3 (red, canalized in gyne-biased direction). Canalization scores in A. echinatior were calculated by comparing gyne and small worker transcriptomes. d, Developmental expression dynamics of ATP- dependent RNA helicase vasa (vas), vitellogenin receptor (yl) and serine/threonine-protein kinase Chk2 (lok) in the two ant species. All three genes are ovary specific with high expression abundance in D. melanogaster ovaries. Although these three genes showed increasing gyne-biased canalization in M. pharaonis as development proceeds, there was little expression difference between gyne and worker individuals in A. echinatior, except for yl in prepupae and to a lesser extent also in 3^rd instar larvae.

Extended Data Fig. 7. Freja’s functional role in canalizing gyne phenotypes.

a, Predicted functional domains of the protein encoded by Freja (LOC10587931), annotated with InterPro (see Methods). Freja contains a signal peptide domain at the N terminus, indicating secretion or membrane insertion. In addition, Freja contains a leucine-rich-repeat domain, suggesting its role in protein binding. b, Freja is the most strongly canalized gene in M. pharaonis, showing an increasing between-caste expression difference and a decreasing within-caste expression variance as development proceeds (n_gyne = 168 and n_worker = 161). The caste identities in 1^st instar larvae are based on BPA prediction. c, Tissue-specific RT-PCR quantification of Freja transcript abundance in adult gynes, showing Freja’s expression is restricted to the abdomen and especially highly abundant in the ovaries. d, RT-PCR quantification of the efficiency of Freja-RNAi in 3^rd instar larvae (left) and adult gynes (right). Compared to the GFP-RNAi control group, Freja-RNAi significantly reduced the expression level of Freja (p < 1e-3 in one-way ANOVAs in each age group). e, Compared with the control group, Freja-RNAi significantly reduced the number of yolky oocytes in adult gynes (p = 0.004 in one-way ANOVA). For Extended Data Figure 7b–e, box plots show the median (centre line), 25% and 75% quartiles (boxes), outermost values (whiskers) and data points (overlapping with box and whiskers).