The convergent evolution of caste in ants and honey bees is based on a shared core of ancient reproductive genes and many plastic genes

Eusociality, characterized by caste-based division of labor, has convergently evolved multiple times. However, the genomic basis of caste and degree to which independent origins of eusociality have utilized common genes is largely unknown. To elucidate these issues, we characterized caste-specific transcriptomic profiles across development and adult body segments from honey bees (Apis mellifera) and pharaoh ants (Monomorium pharaonis), representing two independent origins of eusociality. We identified a shared core of genes upregulated in the abdomens of queen honey bees and ants that is also upregulated in female flies. Outside of this shared core, few genes are differentially expressed in common. Instead, the majority of genes underlying the caste system show plastic expression, are rapidly evolving, and are relatively evolutionary young. Altogether our results show that the convergent evolution of eusociality involves the recruitment of a core reproductive groundplan along with many plastically-expressed and rapidly evolving genes.

caste-and sex-bias was not restricted to the abdomen, as there was similar highly significant 182 effect when comparing heads and thoraces, albeit with weaker effect sizes ( Figure S4). 183 Given the link between conserved caste bias and sex bias within ants and honey bees, we 184 hypothesized that caste bias is derived from ancient pathways underlying sexual dimorphism. To 185 test this hypothesis, we estimated the whole-body sex-bias of orthologs in the fruitfly Drosophila  Expression plasticity across development, caste, and tissue is correlated between species 193 While we have emphasized the conservation of abdominal differential expression between 194 queens and workers in pharaoh ants and honey bees, roughly two-thirds of abdominal caste-195 biased genes were not shared between species, and differential expression based on reproductive 196 caste or worker division of labor was largely not shared ( Figure 1A,B). Furthermore, genes were 197 often differentially expressed across many stages and tissues and sometimes in opposite 198 directions (e.g., upregulated in queen heads but downregulated in queen abdomens), while other 199 genes showed little to no expression differences ( Figure S5). To quantify the degree to which  honey bees), where we expression was averaged across the relevant comparisons, analogous to 239 our calculation of overall caste or behavior bias. However, our results are not solely driven by 240 expression levels. We performed partial correlation analysis accounting for expression to test for 241 the conditional effect of evolutionary age, connectivity, tissue specificity, and evolutionary rate 242 on caste or behavior bias. The relationship between caste/behavior bias and evolutionary age, 243 evolutionary rate, connectivity, and tissue-specificity generally all remained significant and in 244 the same direction when expression was accounted for (Table S11). The two exceptions were the 245 relationship between ant behavior bias and evolutionary age and honey bee connectivity and 246 . CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint caste bias (Table S11). However, when we removed the abdominal contribution to caste bias and 247 expression, the relationship between honey bee connectivity and caste bias remained significant 248 and negative (Table S11), revealing that the relationship between connectivity and caste bias was 249 universal except in the case of highly connected abdominally caste-biased genes, as were present 250 in our gene coexpression analysis ( Figure 2B). Together, these results indicate that each gene's 251 likelihood of recruitment to regulatory networks underlying caste or worker division of labor is In this study we present the most comprehensive transcriptomic dataset to date of caste-258 based division of labor for species representing two independent origins of eusociality. We find 259 generally low overlap between honey bees and ants in the thousands of genes that are 260 differentially expressed between queen and worker castes ( Figure 1A), and between nurses and 261 foragers within the worker caste ( Figure 1B), indicating that the independent evolution of these 262 plastic phenotypes in honey bees and ants have largely involved different genes. The one notable 263 exception is in the abdomen, where we find that about a third (~1500/4500) of genes 264 differentially expressed between queens and workers are shared between honey bees and ants. 265 We find that genes with conserved queen bias in ants and honey bees also tend to be female-266 biased in D. melanogaster, indicating that these conserved caste-biased genes are derived from 267 ancient plastically expressed genes underlying sexual dimorphism (Figure 3). Outside of this conserved core set of ancient genes associated with female reproduction, the majority of genes 269 associated with caste and worker division of labor showed plastic expression in multiple contexts 270 (e.g., between developmental stages, tissues, and castes) and were relatively young, rapidly 271 evolving, and loosely connected in regulatory networks.

272
In the long-extinct solitary ancestors of social insects, female behaviors such as egg-  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint completely lack reproductive organs and are obligately sterile, while honey bee workers are 291 facultatively sterile and have greatly reduced ovaries relative to queens).

292
While our results point to a general genetic link between caste dimorphism and sexual 293 dimorphism, a similar link may exist between caste determination and sex determination.   The tasks performed by workers (specifically, nursing versus foraging) represent a plastic 306 phenotype that changes over the course of the worker's adult lifetime (Mikheyev & Linksvayer 307 2015). This behavioral change is known to be accompanied by a wide range of physiological 308 changes and is regulated at least in part by conserved physiological pathways, for example, those The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint proportion of shared genes was much lower in comparison to genes underlying abdominal 313 differences between queens and workers. Nonetheless, we identified a number of GO terms 314 associated with metabolism as well as development in each species (Tables S6,S7), which is 315 consistent with the notion that the transition from nurse to forager is essentially a developmental 316 process, and that common molecular pathways may provide the raw genetic material for social   Additionally, we found that orthologs with shared abdominal worker bias between honey bees 332 and pharaoh ants were on average evolutionarily younger (i.e. more taxonomically-restricted) 333 than those with shared queen bias ( Figure 1C The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint others', it is likely that layered on top of a minority of shared ancient genes associated with 336 sexual dimorphism, the bulk of the genetic architecture underlying caste-based division labor is  Similarly, tissue specificity and connectivity affect evolutionary rate because highly pleiotropic 361 (highly connected, not tissue-specific) genes experience enhanced selective constraint (Liao et al. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint to networks, highly pleiotropic, and slowly evolving. However, the emerging picture is that the 380 evolution of other traits such as social innovation may largely utilize genes on the opposite end 381 of the spectrum: genes that are peripheral elements to networks and exhibit high levels of and pharaoh ants, and these genes largely appear to be associated with queen reproduction.

388
However, it is important to note that other than the striking phenomenon of strong queen-worker 389 dimorphism that characterizes the reproductive caste system in both honey bees and ants, the The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint polygenic, such that even though only a small number of genes can be identified with large effect 423 sizes and relatively straightforward causative molecular paths, many more genes act peripherally, 424 with smaller effects and through more circuitous routes to profoundly influence trait expression   Table S1 for list of all samples). We separated adults into the three main body segments 441 (head, mesosoma, and metosoma) upon sample collection and sequenced pools of each part 442 separately. For convenience, we refer to these segments as "head", "thorax", and "abdominal" 443 tissues throughout. We sequenced whole embryos and whole bodies of larvae and pupae. Each way, we identified one-to-one, one-to-many, and many-to-many orthologous groups between A. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint one-to-one orthologs (i.e. genes for which only one gene from each species matches the given 467 OrthoDB9 ortholog group). We identified three-way 1-1-1 orthologs between A. mellifera, M. 468 pharaonis, and Drosophila melanogaster using a similar procedure based on Endopterygota 469 orthology groups from OrthoDB9.    samples and no more than three other samples total (note that honey bee queen abdomen samples 500 clustered with egg samples, while pharaoh ant queen samples did not cluster with egg samples).

501
Because the same genes were not always present in such a bicluster, we tabulated the 502 number of queen abdomen biclusters each gene was found in and retained genes present in a 503 higher proportion of biclusters than a given cut-off, determined by inspection of frequency 504 distributions of bicluster presence. In pharaoh ants, we found a large set of genes present in 505 greater than 90% of queen abdomen biclusters, and we retained these genes for further analysis 506 (N = 1039 genes; Figure S10A, i.e. the same set of genes was repeatedly found). In contrast, 507 honey bee queen abdomen biclusters tended to contain one of two groups of genes, as the 508 frequency of presence in the bicluster peaks at 60% and 30% ( Figure S10B). Out of 1245 genes 509 present in greater than 60% of the identified biclusters, 877 were differentially expressed and The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint upregulated in queen abdomens relative to worker (also note that this set of genes exhibited 511 much higher expression in eggs than the latter set). In contrast, out of 1057 genes present in 25-512 35% of biclusters, 611 out of were differentially expressed and upregulated in worker abdomens, 513 compared to 47 upregulated in queen abdomens. Therefore, it is clear that the more common 514 bicluster represents genes associated with queen abdomens, so we retained this set of genes for 515 further analysis (N = 1245 genes). 516 We proceeded with our analysis using these identified sets of genes, which we term  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint decided to focus on the difference between "ancient" genes, which we defined as displaying 534 orthology with non-insect animals, and a number of hierarchical younger categories: "insect", 535 "Hymenoptera", "Aculeata", "ant", "bee", and "novel" (where "ant" refers to genes found in M. 536 pharaonis and other ants but not in any other species, "bee" refers to genes found in A. mellifera 537 and other bees but not in other species, and "novel" refers to a gene found only in A. mellifera or annotated hymenopteran genomes (48 total). We added to this group ten non-hymenopteran 543 insect genomes and ten non-arthropod genomes (see Table S12 for a full list of included 544 genomes). Therefore, a gene labeled as "ancient" displayed a significant BLASTp hit to one of 545 the ten non-arthropod genomes. While phylostratigraphy typically employs an extremely large 546 database containing all available representative taxa, we reasoned that for our study resolution 547 between categories such as "Bilateria" and "Eukaryota" was unnecessary. Furthermore, adding 548 extraneous genomes effectively dilutes the database, such that more similarity is needed to pass 549 an E-value threshold. Because we included only a sample of non-hymenopteran genomes, we 550 were therefore able to stringently identify orthologs (E-value 10 -10 in comparison to a typical 551 value of 10 -5 (Drost et al. 2015)) and accurately place them along the hymenoptera phylogeny. for each sex). While much recent work has detailed tissue-specific (and even cell-specific) 559 expression in D. melanogaster, we reasoned that whole bodies (lacking the same body parts as 560 we used in this study) would give us reasonable overall levels of sex bias to compare our data to. 561 We aligned reads to the current release of the D. melanogaster genome available on NCBI 562 (assembly release 6 plus ISO) using the same pipeline as described for A. mellifera and M. 563 pharaonis, and we performed differential expression analysis between males and females as 564 described above.  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Estimation of tissue specificity 578 While comparing expression between castes across the entirety of development in many tissues 579 was unfeasible for this study, a previous study performed RNA-sequencing on twelve tissues in 580 A. mellifera in worker nurses and foragers (Jasper et al. 2015). We downloaded the available data 581 and aligned reads to the current version of the A. mellifera genome, using the pipeline described Estimation of dN/dS 587 We estimated evolutionary rate using dN/dS, the ratio of non-synonymous to synonymous 588 nucleotide changes. We estimated pairwise dN/dS between each focal species and a second 589 closely related species. For A. mellifera, we compared protein-coding sequences to A. cerana, 590 and for M. pharaonis, we compared protein-coding sequences to Solenopsis invicta. We 591 identified 1-1 orthologs using OrthoDB9, as described above. For each 1-1 ortholog pair, we  orthologs. We identified GO terms associated caste-biased or behavior-biased genes using the P-604 value of differential expression between queens and workers or nurses and foragers. We 605 identified GO terms associated with overall caste or behavior bias using the value of overall 606 caste bias or behavior bias as input. We identified enriched terms with P-value < 0.05.

608
General statistical analyses 609 We performed all statistical analyses and constructed plots in R, version 3.4.0 (R Core Team        A. mellifera (right). "Head", "thorax", and "abdomen" refer to body segments of adults, while pupa and larva refers to whole bodies. Genes are divided by color into categories based on gene-gene orthology between the two species. Insets show the proportion of each category of gene out of all differentially expressed genes at that stage. C) Genes with conserved queen bias (upregulated in queen abdomens of both species) are comprised of proportionally ancient genes than those with conserved worker bias or genes with non-conserved bias *: the category "larva" represents differential expression across larvae of all stages for which caste can be identified (second to fifth larval stage).
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure 2. Abdominal caste bias (log2 fold change queen versus worker) is correlated with connectivity within the queen-abdomen specific bicluster (i.e. module) in A) ants (rho = 0.536, P < 0.001) and B) honey bees (rho = 0.617, P < 0.001), indicating that centrally located genes in the bicluster are highly caste-biased. Genes upregulated in queens are in red, while genes upregulated in workers are in blue. Connectivity is proportional to the most highly connected gene in the module. Furthermore, connectivity within the queen abdominal module is higher in genes found in the module for both species (shared) versus genes found in the module for only one species (not shared) in C) ants and D) honey bees. ***P < 0.001 (Wilcoxon test) . CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure 3. Caste bias is derived from sex bias. Abdominal caste bias (queen vs worker log2 fold-change) is correlated to abdominal sex bias (queen vs male log2 fold-change) in A) M. pharaonis (rho = 0.715, P < 0.001) and B) A. mellifera (rho = 0.774, P < 0.001) and abdominal sex bias is correlated between the two species (rho = 0.280, P < 0.001) (C). Genes are colored by conservation of differential expression in abdomens, with red indicating genes upregulated in queens in both species, blue indicating genes upregulated in workers of both species, and grey indicating genes that exhibited non-conserved expression patterns. Finally, genes with conserved queen bias tend to be femalebiased in D. melanogaster based on whole-body adult samples while genes with conserved worker bias tend to be male-biased in D. melanogaster (likely reflecting downregulation in females)  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S2. Pearson correlation of log-fold change between nurses and foragers as measured in each tissue in M. pharaonis and A. mellifera for each 1-1 ortholog (N = 7640). Error bars indicate pearson correlation 95% confidence intervals.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S3. Log2 fold-change at each stage/tissue in each phylostrata category. Positive values indicating higher expression in queens compared to workers. Log2 fold-change has been adjusted relative to the median value at that stage/tissue, in order to compare across tests. "Ancient" genes indicate any genes shared beyond insects (i.e. with vertebrates).
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S4. Pearson correlation of caste (queen/worker) and sex (queen/male) expression bias in ants and honey bees. Error bars represent pearson correlation 95% confidence intervals. Correlations are significant in all cases (P < 0.001), but abdominal correlations are strongest.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S5. Number of times each gene is upregulated in queen and workers across all comparisons (larva, pupa, and adult head, thorax, and abdomen). Color brightness is logarithmically proportional to the number of genes in each cell. N = 10804 genes (M. pharaonis); 11775 genes (A. mellifera).
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S6. Overall caste bias (A) and overall behavior bias (B) is correlated between ants and honey bees. "Overall" bias refers to the euclidean distance of all log2 fold-change values (queens/worker for caste, nurses/foragers for behavior). Red line is trendline of linear model; Spearman correlation P < 0.001 in all cases.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S7. Overall caste bias and overall behavior bias were correlated within A) ants and B) honey bees. "Overall" bias refers to the euclidean distance of all log2 fold-change values (queens/worker for caste, nurses/foragers for behavior). Red line is trendline of linear model; Spearman correlation P < 0.001 in all cases.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S8. Genes that exhibit more behavior bias across body segments have younger estimated evolutionary ages (A,B) and tend to be loosely connected (C,D; ant: rho = -0.099, P < 0.001; honey bee: rho = -0.157, P < 0.001) and rapidly evolving (E,F; ant: rho = 0.079, P < 0.001; honey bee: rho = 0.226, P < 0.001). "Overall behavior bias" combines nurse forager log2 foldchange values across all adult body segments. Connectivity is calculated using all samples and genes and scaled proportionally to the highest value.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S9. Genes exhibiting more behavior bias tend to be tissue-specific. There was a positive correlation (Spearman correlation, P < 0.001 in each case) between caste/behavior bias and tissue specificity, where tissue specificity ( ) is estimated using data from 12 honey bee tissues. = 1 indicates a genes is expressed in only one tissue, while lower values indicate genes are more ubiquitously (i.e. evenly) expressed across tissues.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Figure S10. Histogram of the frequency with which genes were placed in the queen abdomen bicluster (out of 1000 runs). Plaid biclustering is a non-deterministic process, so different sets of genes can be present in each run. In ants (A), 1039 genes were present in >90% of queen abdomen biclusters and retained for further analysis. There are two peaks in the frequency distribution for honey bees (B). The lower frequency peak is made up of worker-associated genes (downregulated in queen abdomens) while the higher frequency peak (~60%) is made up of queen-associated genes. We retained genes with >60% frequency for further analysis.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S1. Full listing of sample types and number of each sample collected. "L1" and "L2" refers to larvae of the first and second stage, etc. We began caste-specific sampling at stage two because caste is determined and regulated in M. pharaonis by the end of the first larval instar (Warner et al. 2016). After the first larval instar in M. pharaonis, worker-destined larvae can be distinguished from reproductive-destined larvae, which include male-destined and queendestined larvae (Peacock & Baxter 1950). As such, our "queendestined" ant larvae samples likely contain some male-destined larvae, but the proportion is expected to be low, as the sex ratio is known to be heavily queen-biased (Schmidt et al. 2011 Table S2. Number of differentially expressed genes (DEGs) between queens and workers for each comparison (FDR < 0.1). "L2", "L3", etc refer to the 2nd and 3rd larval stage, respectively, while "larva_overall" is the result of differential expression with caste as main effect across all larval samples. Differentially expressed genes are divided into "queen associated", which exhibited higher expression in queens, and "worker associated", which exhibited higher expression in workers. Differential expression analysis performed with N = 10804 (ant) and 11775 (honey bee) genes.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S3. Enriched gene ontology terms based on Gene Set Enrichment Analysis (GSEA) of differential expression between queens and workers in ants. P-value derived from Kolmogorov-Smirnov tests.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S4. Enriched gene ontology terms based on Gene Set Enrichment Analysis (GSEA) of differential expression between queens and workers in honey bees. P-value derived from Kolmogorov-Smirnov tests.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S5. Number of differentially expressed genes (DEGs) between nurses and foragers for each comparison (FDR < 0.1). Differentially expressed genes are divided into "nurse associated", which exhibited higher expression in nurses, and "forager associated", which exhibited higher expression in foragers. Differential expression analysis performed with N = 10804 (ant) and 11775 (honey bee) genes.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S6. Enriched gene ontology terms based on Gene Set Enrichment Analysis (GSEA) of differential expression between nurses and foragers in ants. P-value derived from Kolmogorov-Smirnov tests.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S7. Enriched gene ontology terms based on Gene Set Enrichment Analysis (GSEA) of differential expression between nurses and foragers in honey bees. P-value derived from Kolmogorov-Smirnov tests.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S8. Hub genes of the queen abdominal module in ants. Hub genes were defined as genes with intra-modular connectivity in at least the 90th percentile, and log2 fold-change (queen/worker) greater than 2.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S9. Hub genes of the queen abdominal module in ants. Hub genes were defined as genes with intra-modular connectivity in at least the 90th percentile, and log2 fold-change (queen/worker) greater than 2.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S10. Enriched gene ontology terms based on overall caste or behavior bias in ants and honey bees. GO terms are derived from D. melanogaster orthologs. P-value is from gene set enrichment analysis (Kolmogorov-Smirnov test) . CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S11. Partial correlation between connectivity, evolutionary rate (dN/dS), evolutionary age (phylostrata), and tissue-specificity (tau) and caste or behavior bias while accounting for expression. Analysis was performed separately for each species and comparison (i.e. separately for caste bias and behavior bias), as well as while including or excluding abdomen in calculations of caste bias and expression. Connectivity is total connectivity measured across all samples and genes. Phylostrata is a measure of estimated evolutionary age, with higher values indicating younger genes. Tau is the degree to which genes exhibit tissue-specific expression across 12 honey bee tissues (results presented only for honey bees). N = 10520 genes (ants), 10011 genes (honey bees).
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint Table S12. List of species used for phylostratigraphy analysis, with the NCBI Taxonomy ID.
. CC-BY 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/454645 doi: bioRxiv preprint