Phylogenetic analysis and trait evolution of ant cocoons

Most ant species have lost the ability to spin cocoons. To explore the evolution of cocoon loss within Formicidae, we perform an ancestral state reconstruction of cocooned pupae across a genus-level phylog-eny and use a sister clade analysis to determine the impact of cocoon evolution on ant speciation. Then, we fit models of correlated evolution between cocoon status and several other organismal traits. We find that the re-emergence of cocoons is rare and that “naked” lineages display an increased rate of speciation in 5 out of 9 sister group comparisons. Models of correlated evolution with cocoon status were favored for metapleural gland and worker polymorphism. Metapleural gland favored rates of evolution were inconclusive, while worker polymorphism displayed a higher transition rate towards polymorphism coupled with cocoon loss. These results suggest that cocoon loss may allow for other complex traits to develop and may represent a novel example of relaxed selection.


Introduction
Given the charismatic nature of ants, the scientific literature surrounding their natural history has largely focused on behaviors such as foraging and sociality (Traneillo 1989;Wheeler 2015;Reeves & Moreau 2019), and in recent years has had an increasingly molecular focus in areas such as phylogenetics (Larabee et al. 2016;Blanchard & Moreau 2017;Soloman et al. 2019), genome evolution (Williams et al. 2015;Rubin & Moreau 2016;Liu et al. 2019), and the microbiome (Lucas et al. 2017;Chua et al. 2018;Flynn & Moreau 2019;Ramalho et al. 2019). Here we address a lesser known area of research within ant life-history -cocoon evolution. Along with the rest of Endopterygota, Hymenoptera undergo holometabolous development. While the appearance of a pupal stage is maintained throughout the Formicidae, larvae pupate via two generalized approaches. In one, larvae pupate by first encasing themselves in a silken cocoon produced by larval silk glands, while in the alternate approach species develop "naked" with their cuticle fully exposed to the surrounding environment.
Here, we define a cocoon as an external structure composed of silk proteins that surrounds a pupa during development. Though silk cocoons are most well-known among Lepidoptera (most notably in the silk moth, Bombyx mori), their presence is widespread among the insects including in Hymenoptera, Coleoptera, Neuroptera, and Siphonaptera (Silverman et al. 1981;Kechington 1983;Weisman et al. 2008). While the use of cocoons is common, methods to produce the silk necessary to create them has independently evolved several times (Sutherland et. al 2010). Within the Hymenoptera alone, the ability to produce silk has independently evolved at least six times with variation in their primary use -including nest building and prey entrapment in addition to cocoon spinning (Sutherland et al. 2010(Sutherland et al. , 2012. Generally, the ant cocoon serves to protect pupae during development into the adult stage, but whether it functions to protect the pupa from predation, infection, humidity, or some combination of factors varies depending on the species (Danks 2004;Kaltenpoth et. al 2005;Danks 2006;Ellis et al. 2010;Zhang et al. 2013).
Robust conclusions for the absolute function of ant cocoons and their subsequent loss among most ant lineages has remained elusive (Wheeler 1915). It has been speculated that the presence of the insect cocoon works to protect pupae from pathogens and predation, functions as a means of climate control, and/or traps energy from light for enhanced development (Wheeler 1915;Danks 2004;Tragust et al. 2013). While it is likely that some or all of these functions have maintained cocoons in particular ant lineages, little research has focused on the ultimate causes for cocoon loss. The presence of a cocoon is the ancestral state in Formicidae, yet the trait has been lost in most extant ant species (Taylor 2007;Armitage et. al 2011). In spite of this loss, in several genera larval silk remains an integral part of their life history. For example, Oecyphylla and many species of Polyrhachis produce naked pupae even though their larvae actively produce silk for nest building (Wheeler 1915;Hölldobler & Wilson 1983;Robson et al. 2015). And interestingly in fungus growing ant species, despite not producing silk cocoons, workers actively cover their brood in a fungal sheath that appears to serve a cocoon-like role (Lopes et al. 2005;Armitage et al. 2011). These observations present a few questions: 1) When a cocoon is initially lost in a lineage, is it difficult to regain despite its ecological utility? and 2) What are the ecological factors that promote the loss or maintenance of cocoons in ant species?
Several morphological factors could reasonably correlate with and influence ant cocoon evolution. Although ant cocoons may function to guard pupae from fungal infection, metapleural glands and the sting have both been shown to produce antimicrobial compounds; the metapleural gland in particular is used for brood sanitation purposes among some species (Fernández-Marín 2006;Yek and Mueller 2011;Perito et al. 2018). Therefore, the presence of a metapleural gland and/or sting may reduce the need for a cocoon if part of its function is to protect pupae from pathogens. In addition, polymorphism in worker size may play an unexpected role. Both cocoons and monomorphic workers are the ancestral state in Formicidae (Wheeler 1915;Wilson 1953;Taylor 2007). Thus, the developmental and energetic costs required to produce polymorphic workers may impede maintenance of the cocoon phenotype over evolutionary time (Porter & Tschinkel 1985;Wheeler 1991;Hughes et al. 2003).
In addition to morphology, ecology and behavior could influence cocoon evolution. Nesting behavior and climate regime influences a colony's exposure to light, humidity, and microbes (Frouz 2000;Walker & Hughes 2011;Kadochová & Frouz 2013). As such, if a cocoon's function is to help pupae thermoregulate and guard against pathogens, these ecological factors may well correlate with cocoon presence in ant species. Diet could influence nutrient intake and thus the ability to obtain the amino acids necessary for silk production (Rudall & Kenchington 1971). Further, parasitic ant species who live among their host ants must invade nests while avoiding detection at all stages of development (egg, larvae, pupae, adult) as they would readily be removed or killed (Buschinger 2009;Chernenko 2011). Although there is no direct evidence to suggest cocoon presence/absence impacts the ability of parasitic ant species to evade host detection, we include it in this study to explore if there is a correlation between these two traits.
To address these questions in this study, we first model an ancestral state reconstruction of cocoon evolution on a recently published genus-level phylogeny (Nelson et al. 2018) of the family Formicidae. Next, we evaluate the impact of cocoon presence/ absence on diversification using sister clade analysis. Last, we measure rates of correlated evolution between cocoon presence/absence and various relevant life history traits to determine the potential influence of ecological or morphological characteristics on the lability of cocoon evolution.

Ecological and morphological data collection
We collected ecological and morphological data on 268 ant genera from the primary literature, online databases (AntWiki.org), previously compiled data from Blanchard & Moreau (2017), and personal observations by authors (Supplementary Table S1). For each genus, data on cocoon presence and/or absence were recorded. In genera where a subset of species had unknown cocoon status and all known species within the genus displayed a singular character state, we assumed the genus was monomorphic for that state. For those genera that are polymorphic for the cocooned state (both absent and present) they were recorded as having both states (0&1). Genera with completely unknown cocoon status were recorded as missing information and removed from the analysis. For each genus, morphological data were collected on the presence of a metapleural gland, sting, and size polymorphism among workers. Ecological data were collected on nesting preference, diet, parasitism behavior, and geographic distribution (tropical/subtropical/temperate).
All geographic information was gathered from Antmaps.org (Janicki et al. 2016;Guénard et al. 2017). Upper and lower limits for a geographic range were based on the borders of the most northward and southward countries/regions in which they are found. For a genus to be included in a geographic zone, it must have at least one country/region that entirely exists within that zone. We used the Tropic of Cancer and the Tropic of Capricorn to designate the upper and lower bounds of the tropical zone (23.5˚ N/S). The area between 23.5˚ N/S and 35˚ N/S we designated as the subtropical zone and all areas beyond 35˚ N/S as the temperate zone.
For all morphological and ecological data collected, characters were initially recorded as binary, multistate, and/or polymorphic depending on the trait in question. However, for ecological correlation analysis all states were converted into binary characters (see Ecological Correlations, below).

Ancestral state reconstruction
We used the ant phylogeny reconstructed by Nelsen et al. (2018) for all subsequent analyses. We modified this tree to the genus level, which contained 317 of the 334 currently known extant ant genera and was based on 1,730 ant species. We selected species to represent genera based on those for which relevant ecological data was available. Within the R analysis platform (R Core Team, 2018), we used the drop. tip function in package ape v5.1 (Paradis 2004) and pruned fourteen outgroup taxa from the tree as well as an additional eighteen genus-level tips for which either cocoon status was unknown or there was a severe lack of morphological/ecological data available. Further, we dropped additional tips from the thirteen polyphyletic genera from the Nelsen et al. (2018) tree, leaving one tip to represent each genus. We retained the tip that most accurately mirrored the phylogenetic relationships found in Blanchard & Moreau (2017), allowing for a tree represented by solely monophyletic groups.
Ancestral state reconstruction (ASR) of cocoon evolution was estimated and visualized at each internal node of our phylogeny using the function rayDisc in the R package corHMMv1.22 (Beaulieu et al. 2012(Beaulieu et al. , 2017. This package is useful as it generates a maximum likelihood (ML) model of evolution that allows us to include polymorphic character states (both cocoon absent and present species within a genus). Ancestral states at each internal node were calculated using marginal reconstruction. To determine the best model for cocoon evolution, we calculated and compared the Akaike information criterion values corrected for small sample size (AICc) across the equal rates (ER), symmetrical rates (SR), and all rates different (ARD) models. A lower AICc indicates a more supported model; when comparing models, we chose a difference of four as a sufficient cut-off to determine significant increase in model fit (Guthery 2003;Alfaro et. al 2009).

Sister clade analysis
To assess the impact of cocoon presence/absence on speciation, we implemented a sister clade analysis on our ancestral state reconstruction exemplar tree using the richness. yule.test function in the R package ape (Paradis et al. 2004). Although state-dependent speciation and extinction (SSE) models such as BISSE have shown promise in addressing many of the issues of earlier models, criticism has fallen upon these methods of diversification analysis to accurately detect relationships between character state evolution and phylogenetic clades (Ng & Smith 2014;Maddison & Fitzjohn 2015;Rabosky & Goldberg 2015). The hidden state speciation and extinction model (HISSE) has sought to address many of the concerns brought towards the BISSE model, however the most recent version of HISSE (1.9.4) does not support genus-level phylogenies as it does not allow for the inclusion of clade specific species numbers or fractions (Beaulieu & O'Meara 2016). Therefore, we determined that a model-based sister clade analysis is the most robust analysis for the effects of character states on our exemplar genus level phylogeny as it incorporates the number of species in each clade as well as the date since divergence between sisters as noted in our ASR tree (Paradis 2011).
A sister clade was included if one clade of exclusively cocooned species shared a common ancestor with a second clade of exclusively naked species. This was the case for nine sister clades. The clade including the polymorphic genus Polyrachis could also form a 10th sister clade pair with the genus Colobopsis (Fig. 1). As such, we ran two analyses; one with Polyrhachis clade included (ten sister clade pairs) and one without (nine sister clade pairs). For the analysis with Polyrhachis included, the genus was coded as cocooned as this is the state for the majority of known Polyrhachis species.

Ecological correlations
First, we visualized potential patterns of paired evolution with cocoon presence and a given morphological or ecological trait of interest (e.g. Fig. 2). Then, we used the corDISC function in the R package corHMMv1.22 (Beaulieu et al. 2012(Beaulieu et al. , 2017 to fit models of correlated evolution between cocoon status and the remaining morphological and ecological traits collected. corDISC is useful as it allows for missing data in the traits being compared. For this test, we compared a four-parameter independent model to an eight-parameter dependent model of correlated evolution. For traits in which correlated evolution with cocoons was favored, we then performed a parameter restriction test, comparing a seven-parameter model to the eightparameter model. This helps determine which rate correlation is favored when the direction of evolution of one paired trait state is compared to the direction of evolution for a separate paired trait state (e.g. whether correlated rate of evolution favors cocooned genera then evolving size polymorphism as opposed to genera with naked pupae evolving size polymorphism.) As certain rate correlations may not have biological relevance, we deemed a rate correlation as having potential biological significance if it met one of three criteria: 1) From a given trait pair, the evolution of trait A is favored over the evolution of trait B or vice versa (Fig. 3A). 2) the transition of trait A1 to A2 or A2 to A1 is favored, where B remains constant (Fig. 3B). 3) A significant difference in the rate of evolution of trait B dependent on the state of trait A (Fig. 3C).
corDISC requires binary states, all multistate and polymorphic traits had to be to converted to fit the required format. For cocoon status to be placed into only binary states, all polymorphic genera (0&1) had to be reclassified as present or absent. This was the case for seven out of the 262 genera included in the analysis. We chose cocoon presence as state "0" to remain consistent with the ant taxonomic literature (Urbani et al. 1992). We implemented a "majority rule" such that if over 50% of the known species within a polymorphic genus are of a particular state, the entire genus was reclassified into the majority state for the analysis. Six of the seven polymorphic genera were thus classified as cocoon present (state 0): Ambylopone, Formica, Hypoponera, Lasius, Polyrhachis, and Ponera. The remaining genus, Proceratium, was classified as absent (state 1). Cocoon state was then correlated with the presence of a metaplural gland, presence of a sting, diet, nesting habitat, social parasitism, worker polymorphism, and geographic zone. The morphological characters metapleural gland and sting were both coded as 0 = absent, 1 = present. For the metapleural gland, polymorphic genera for this organ were classified by a majority rule (i.e. if over 50% of species within the genus possess a metapleural gland, the genus was labeled as 1). As the vast majority of ants possess a metapleural gland, if the existence of the gland was unknown for a genus it was assumed present. Sting characters were already a binary state and therefore did not require a transformation. In addition, social parasitism was recorded as a binary character as  In each case red arrows depict one of the two possible paths evolution could take. Red letters illustrate where the evolutionary change takes place. A) From a given trait pair, the evolution of trait A is favored over the evolution of trait B or vice versa. b) the transition of trait A1 to A2 or A2 to A1 is favored, where B remains constant. c) A significant difference in the rate of evolution of trait B dependent on the state of trait A. either (0) no species in the genus parasitize other ants or (1) some to all species in a genus are social parasites.
In contrast, size polymorphism data were originally collected as a multistate character as either monomorphic (0), some size variation (1) or true polymorphism (2). To correlate size polymorphism against cocoon presence, it was binned and measured in two separate schemes: "strict polymorphism" meaning monomorphic/some variation (0) vs. polymorphic (1) and "variable polymorphism" meaning monomorphic (0) vs. any size variation (1).

Ancestral state reconstruction
After pruning, our phylogeny consisted of 262 genera representing 15 extant ant subfamilies-Martialinae and Agroecomyrmicinae were excluded due to lack of relevant data (Fig. 1). Genera with a cocooned state comprised 33.5% of the tree (88 Genera; representing 3,422 species), while naked pupa consisted of 63.7% (167 Genera; representing 8418 species), and the remaining 2.6% being polymorphic for the trait (7 Genera; representing 1,300 species). When modeling cocoon evolution, all rates different (ARD) had the lowest AICc value and was thus supported as the most appropriate rate model for the ancestral state reconstruction (Supplementary Table S2). The ancestral state reconstruction (ASR) analysis inferred the presence of a cocoon as the ancestral state in the ants in concordance with previous work (Wheeler 1915;Taylor 2007;Armitage et al. 2011). Analysis of the internal nodes suggests that when the cocoon state is lost, the trait rarely evolves again; across all 262 extant tips, there are only two possible instances of the re-emergence of a cocoon after it has been lost. Further, genera that are polymorphic for the cocoon trait (n = 7 independent lineages) exclusively emerge from clades with a cocooned ancestral state.

Sister clade analysis
Five out of the nine sister clade pairs had greater speciation in the naked clade (Table 2). All sister clades in which the most recent common ancestor was over 100 Mya (n = 3) had greater speciation in the naked character. The naked state was associated with a significant increase in diversification, but the rate differed depending on the number of sister pairs analyzed. In our nine-sister pair analyses, there was a strong signal of Table 1. Binning schema for all traits measured for ecological correlation analysis. All traits were converted into binary characters for use in the ecological correlations analyses.

Ecological correlations
Of the 17 ecological traits we measured against cocoon status, the dependent model of correlated evolution was favored for three conditions: both bins of size polymorphism (strict and variable) and metapleural gland presence (Table 3; data on correlation models for all 17 traits located in Supplementary Table S2). Within strict size polymorphism, the parameter restriction test found thirteen favored rate correlations with cocoon state; four of which may be biologically relevant given our criteria (Fig. 4). We found that the transition rate for the loss of a cocoon is favored in the presence of polymorphism compared to in the presence of monomorphism (Fig. 4A). However, once in the "naked" state polymorphism loss is favored relative to the gain of polymorphism in the "naked" state ( Fig. 4C). This is further supported as the naked-polymorphic state has a significantly higher rate of change to monomorphism as opposed to the evolution of a cocoon (Fig. 4D).
Within variable size polymorphism our parameter restriction test found fifteen favored rate correlations with cocoon status of which six could have biological relevance for our analysis. Four of these rates displayed the same directionality as those found in the strict size polymorphism bin (Fig. 4). Of the additional rates found, one correlation demonstrates that when cocooned, there is a higher rate toward polymorphism loss rather than its gain, supporting the central results from the strict polymorphism bin (Fig. 4E). The second additional rate shows a somewhat contrary result, with a higher rate of polymorphism loss when pupae are naked (Fig. 4F).
Our parameter restriction test for metapleural gland presence found nine favored rate correlations with cocoon status, of which two could have biological relevance given our parameters (Supplementary Fig. S1). Our analysis displays an increased rate Table 2. Sister clade analysis results, showing the clades being compared (Figure 1), the number of known species in each clade, and the estimated time since divergence in millions of years (mya). Sister pair 2 included Polyrachis, and was included in a seperate sister clade anaylsis.   towards the gain of a metapleural gland as opposed to its loss both when pupae are cocooned and naked. This result suggests that the presence of a metapleural gland is not associated with cocoon status, and is likely due to ant genera without a metapleural gland being underrepresented in our phylogeny. Details on the parameter restriction tests for biologically relevant correlations of cocoons with worker polymorphism and metapleural gland presence can be found in Supplementary Table S3.

Discussion
Cocoon formation is an important feature in the life history of many insect orders, however its evolution has been understudied in several taxa, including the ants. Within Formicidae, hypothesized roles for cocoons include, among others, shielding pupae from predators and pathogens and functioning as a means of thermoregulation (Wheeler 1915;Frouz 2000;Walker & Hughes 2011). However, despite the supposed utility of cocoons, most ants have lost this feature from their life history. This is even more surprising given that ancestral lineages are hypothesized to have been cocooned (Taylor 2007;Armitage et al. 2011) and few investigations to date have addressed possible ecological factors contributing to cocoon loss. In this study, we conduct a phylogenetic analysis of cocoon evolution within the Formicidae and model the potential for correlations between cocoon presence/loss and several morphological and ecological traits.

Ancestral state reconstruction
Research examining the emergence of novel character states is widespread, while less attention is given to trait loss (Ellers et al. 2012). Here, we find that while cocoon presence is the ancestral state within Formicidae, the majority of extant ant genera included in the analysis have naked pupa (roughly 64% of extant ant species), representing at least 12 independent losses of cocoons. Once the ancestral cocoon state is lost, possible reemergence only occurs twice -in the genus Lioponera and the clade containing Eciton, Nomamyrmex, and Labidus. These losses may be the result of relaxed selection by which the appearance of some other ecological/morphological factor through evolutionary time impedes cocoon maintenance or compensates for its disappearance (Coss 1999;Lahti et al. 2009). Loss of fundamental traits through relaxed selection has been identified in a diverse array of taxa, including eyesight in P. andruzzii cave fish, auditory sensitivity in noctuid moths, and photosynthesis in early eukaryotes (Fullard et al. 2007;De Castro et al. 2009;Calderoni et al. 2016). In each case, a change in the environment allowed for trait loss and the reallocation of energy towards other phenotypes. Spinning cocoons within ant species may follow similar trends. For example, the possible redundancy of function with fungal sheath use in leaf-cutter ants or the intricate use of silk for nest weaving in Oecophylla species provide insights into ecological factors that may result in relaxed selection for cocoon presence (Lopes et al. 2005;Robson et al. 2015).
Alternatively, there may be an unknown utility in pupae remaining naked. Prior research has shown that naked pupae develop faster, likely due to a smaller energetic burden (Tragust et al. 2013). In species which produce naked and cocooned pupae within a single nest, this feature could be a strategy to alternate between less investment in vulnerable but fast-developing naked pupae growing alongside protected, cocooned pupae that develop more slowly (Purcell & Chapuisat 2012). Further research into specific environmental variables that either reduce cocoon function or lead to the utility of naked pupae should be explored.

Sister clade
The results of our sister clade analysis found that naked clades overall experienced a higher rate of speciation (Table 2). Although there was a decrease in significance once we included Polyrhachis into the analysis, we still detected a signal towards increase in species diversity when pupae develop naked for five out of nine clades with naked pupae development. It should be noted that our sister-clade analysis did not account for extinction events, so it is feasible that the observed higher speciation in naked lineages may be in part driven by lower extinction in these same lineages, though this was not explicitly tested.
This measured increase in diversification among naked taxa could be explained by relaxed selection in particular ecological contexts. Ecological factors such as brood sanitation through increased brood care or in some cases the antimicrobial activity of metapleural glands may remove some of the primary functions of cocoons, leading to their redundancy. Redundancy in cocoon function can lead to loss, which in turn allows for energy that was previously allocated to cocoon formation to be used for other means such as faster development (Purcell & Chapuisat 2012;Tragust et al. 2013). An increased developmental rate may afford a competitive advantage when competing against similar-sized cocooned ants or other insects in similar niches (Sebens 1982;Hill et al. 1993;Reitz & Trumble 2002). In addition, this conserved energy could allow for development of novel traits such as extreme polymorphism or the use of silk in other capacities -for example, Ooechypylla with their heavy use of silk in nest construction and Melissotarsus species which are the only ant species known to produce and use silk as adults (Fisher & Robertson 1999;Crozier et a1. 2010).

Ecological correlations
Based on our analyses, we observed an overall trend toward cocoon loss as worker polymorphism evolves (Fig. 4). In addition, when ants are in a naked/monomorphic state the evolution of polymorphism is favored over the evolution of cocoons (Fig.  4B). Interestingly, once in a naked/polymorphic state, there is an evolutionary trajectory toward monomorphism, albeit with a continued resistance toward the gain of a cocoon (Fig. 4C, D). In no rate change scenario is the gain of a cocoon favored. This finding is reflected in our phylogeny, as the re-emergence of cocoons after its initial loss is rare. The loss of the cocoons is likely far easier than the evolution of cocoon spinning as it represents a complex trait with both physiological and behavioral components. In many cocoon spinning species, partial burial of late instar larvae by nurses is required for silk spinning to commence (Gobin et al. 1998).
Correlations with metapleural gland presence gave mixed results ( Supplementary  Fig. S1). This inconsistency is likely caused by the lack of several independent evolutionary events of metapleural gland loss among ant genera (Bolton 2003;Yek & Mueller 2011). Robust analysis of ecological correlation decreases the fewer independent events which have occurred, and based on our binning procedure only 2-3 independent losses of the metapleural gland occurred ( Supplementary  Fig. S2).
While we find support for an ecological link between cocoon absence and worker polymorphism, it must be noted that the emergence of worker polymorphism does not appear within all genera with naked pupae. Other ecological components not investigated in this study likely also contribute. Since we ran our analyses on a genus level scale across all of Formicidae, we may lose important ecological signals within particular genera. Certain lineages may be pushed toward cocoon loss due to diversification of silk use, while others forwent cocoons for faster development. Analysis of ecological correlations at these finer scales can elucidate these additional patterns. However, it is notable that out of the ecological factors included in this study, worker polymorphism stands out despite the great variation in genus-level ecology and life history and thus likely exerts a strong generalized influence on ant cocoon evolution.
We also note that that we did not find support for ecological correlations of cocoon evolution with the evolution our established parameters for a stinger, social parasitism, diet, nesting habitat, or climate regime. This could be a result of under sampling, and as more ecological data is collected from various ant genera we encourage re-analysis of these traits to determine if these results are maintained. In addition, other traits not tested in this analysis, such as colony size or foraging behavior, may also contribute to cocoon evolution but should be investigated in further studies.

Conclusions
Taken together, our study sheds new light on our understanding of the evolution of cocoons among ants and provides evidence for ecological mechanisms that may that have led to their disappearance among various genera. The correlation between worker polymorphism and naked pupae was a surprising but interesting discovery. The underlying driving factors for this correlation are not addressed in this study, but developmental or energetic constraints could be contributing factors. Our results suggest that cocoon loss may be a result of relaxed selection, however further inquiry is necessary to confirm this hypothesis. We do find that ant genera which lack cocoons display a higher speciation rate, and relaxed selection may be in part responsible. This preliminary study acts as a starting point for continued research into how the evolution of differing metamorphosis strategies among ants is shaped by their ecology and influences their evolutionary trajectories.