Evolution of Social Insect Polyphenism Facilitated by the Sex Differentiation Cascade

The major transition to eusociality required the evolution of a switch to canalize development into either a reproductive or a helper, the nature of which is currently unknown. Following predictions from the ‘theory of facilitated variation’, we identify sex differentiation pathways as promising candidates because of their pre-adaptation to regulating development of complex phenotypes. We show that conserved core genes, including the juvenile hormone-sensitive master sex differentiation gene doublesex (dsx) and a krüppel homolog 2 (kr-h2) with putative regulatory function, exhibit both sex and morph-specific expression across life stages in the ant Cardiocondyla obscurior. We hypothesize that genes in the sex differentiation cascade evolved perception of alternative input signals for caste differentiation (i.e. environmental or genetic cues), and that their inherent switch-like and epistatic behavior facilitated signal transfer to downstream targets, thus allowing them to control differential development into morphological castes.


Introduction
The mechanisms underlying the evolution of phenotypic novelty are hotly debated [1][2][3][4]. A fundamental question is how small genetic or epigenetic changes can produce a set of simultaneous, complementary phenotypic changes required to generate new adaptive trait combinations. A key prediction of the 'theory of facilitated variation' [5] is that regulation acts on evolutionarily conserved switch mechanisms, which then modulate expression of target loci controlling development. This process may facilitate large and complex evolutionary steps because it brings together new combinations of inputs (internal or external stimuli) and outputs (phenotypes) but does not rely on evolution of genes involved in the processes per se. Importantly, the reliance of this mode of evolution on conserved genetic and developmental processes increases the likelihood that the outputs will be functionally integrated and thus nonlethal, similar to the 'two-legged goat effect', a striking example of phenotypic accommodation in which developmental robustness allows the animal to 'adapt' to a previously unselected bipedal lifestyle [6].
Evolution by facilitated variation may be especially important to the origin of developmental polyphenisms in which organisms develop into two or more discrete forms, since polyphenisms typically result from plastic activity of regulatory genes. Additionally, it is likely that regulatory mechanisms controlling one set of polyphenism are pre-adapted to evolve control over newly evolving polyphenisms, for two reasons. Firstly, such mechanisms' pre-existing sensitivities to variable cues make it more likely that they will evolve the ability to perceive alternative gradients of novel cues, relative to constitutively expressed genes. Secondly, their downstream target genes already show inter-individual variability in expression, and the organism will thus already have evolved alternative responses to this variability. Gerhart and Kirschner [5] made predictions about the properties of the "core components" which they hypothesize to be the principal drivers of evolutionary novelty, namely that these components should display both robustness and adaptability, as well as exploratory behavior, state-dependent expression and regulatory compartmentation. The sex differentiation pathways exhibit all these properties, making them prime candidates for facilitating the evolution of new forms of polyphenism. Some components of the sex differentiation pathway (such as the doublesex-mab3 (DM) gene family; [7][8][9]) are evolutionarily ancient and conserved across diverse metazoa, and thus could potentially be involved in generating novel polyphenism in multiple distantly related taxa. In insects, the sequence of sex determination has been called hourglass-shaped [10], with highly variable input signals and downstream targets, but a small set of conserved core regulatory genes including transformer (tra) and doublesex (dsx). doublesex is alternatively spliced depending on the presence of an active TRA protein, and its sex-specific isoforms act as transcription factors causing sex-specific gene expression and development through their differential effects on multiple downstream targets [11,12].
The two social insect 'castes'-queens and workers-differ radically from one another in their developmental environment (e.g. nutritional environment) resulting in differences in size, fecundity, behavior and physiology. Ultimately, the evolution of caste polyphenism thus required concerted evolution of environmental input signals and corresponding developmental responses [13]. Eusociality has evolved at least twice within the Hymenoptera [14], but we presently lack a well-evidenced theory of the genetic mechanisms that allowed caste-specific gene expression to originate. There is increasing evidence that the evolution of polyphenism in ants, bees and wasps was achieved primarily through evolution of regulatory genes, rather than gene content or composition [15][16][17], but the core components involved are largely unknown. Here, we propose that conserved parts of the sex differentiation cascade, including the transcription factor doublesex, evolved sensitivity to new environmental input signals (e.g. nutritional signals), thereby triggering caste-specific gene expression that sends larvae on divergent developmental trajectories. To test this hypothesis, we identified dsx and its female-and malespecific isoforms, and measured their expression across life stages in the four discrete morphs (queens, workers, winged males and wingless males) of the ant Cardiocondyla obscurior. We find that dsx sex-specific isoforms are expressed both sex-specifically and morph-specifically in larvae, pupae and adults. Moreover, ninety other conserved genes with sex-biased expression showed morph-specific expression patterns during larval development, suggesting that cooption of the genes regulating sex differentiation via sex-specific alternative splicing was involved in the origin of morphologically distinct castes.

Verification of haplodiploidy
Queens and workers produced from inter-population crosses were heterozygous for diagnostic microsatellite markers, whereas emerging winged and wingless males as well as one sex mosaic individual expressing both male and female characters exclusively carried the maternal alleles (S1 Table). Although single locus complementary sex determination is unlikely because the species regularly engages in inbreeding [18], C. obscurior appears to use standard haplodiploid sexual reproduction.
Identification of the functional doublesex paralog of C. obscurior The C. obscurior genome [19] has four paralogs containing the DM domain of doublesex (dsx) (pfam00751; Cobs_01393, Cobs_07724, Cobs_09254 and Cobs_18158), representing the ancestral state in holometabolous insects [20]. Sex-specific splice forms are only known from one paralog per species (e.g., in Apis [21] and Nasonia [20]), and the function of the others is unclear. In C. obscurior, only Cobs_01393 was differentially expressed in male and female larval RNAseq data (S2 Table). Moreover, Cobs_01393 had the highest sequence homology to functional dsx in other insects (S1 Fig). Finally, we found that Cobs_01393 was located within~79 kb of prospero; microsynteny of prospero and dsx is conserved across the Hymenoptera [20]. We thus conclude that Cobs_01393 is the functional paralog of dsx.

Sex-and morph-specific expression of dsx isoforms
We identified the full-length sequence and sex-specific isoforms of the functional paralog of dsx using 3' rapid amplification of cDNA ends (RACE) (S2 Fig). The first four exons are identical in both isoforms. The DM domain (pfam00751) is located in exon 2 and the dsx dimerization domain (pfam08828) in exon 4. The female-typical isoform dsx F contains one exon specific to dsx F , whereas the male-typical isoform dsx M excludes that exon but includes two others that are absent in dsx F (Fig 1A and S3 Table). This splicing pattern, with a shortened female transcript, has been inferred for the fire ant Solenopsis invicta [22], and matches dsx sex-specific isoforms in Drosophila melanogaster and Apis mellifera, but not Nasonia vitripennis [20]. While the sexsignaling function of dsx is conserved across highly divergent lineages, recent evidence shows that dsx sequence evolves rapidly [23][24][25], causing substantial inter-specific variation in dsx splicing patterns. A higher level of divergence in dsx compared to other DM domain-containing proteins in our phylogenetic analysis confirms this result (S1 Fig).
We designed primers that spanned the exon boundary of the DM domain-containing exon (to measure the overall expression of both isoforms), as well as primers specific to both isoforms for use in RT-qPCR. We found significantly higher expression of the DM domain in adult males (pooling winged males and wingless "ergatoid" males; WM and EM) compared to females (pooling queens and workers; WO and QU) (n EM = 8, n WM = 8, n WO = 10, n QU = 10; Welch two sample t-test: t 25.7 = -8.7, p<0.001, Fig 1B). Expression of the DM-domain was similar in queens and workers (t-test with Benjamini-Hochberg (BH) correction: p = 0.672), but higher in winged males compared to wingless males (p = 0.009).
We then compared expression of dsx F and dsx M across all four morphs in pupae (n EM = 10, n WM = 10, n WO = 9, n QU = 10) and adults (n EM = 7, n WM = 7, n WO = 7, n QU = 7; Fig 1C and  1D). We found morph-specific signatures of expression in both life stages for dsx F (ANOVA: pupae: F (3,35) = 42.33, p<0.001; adults: F (3,24) = 3.75, p = 0.024) as well as for dsx M (Kruskal Wallis rank sum test with df = 3: pupae: X 2 = 30.2, p<0.001; adults: X 2 = 22.6, p<0.001). Worker pupae showed significantly higher dsx F expression than queen pupae (pairwise t-test with BH correction: p = 0.013) and worker pupae and adults showed significantly lower dsx M expression than queen pupae and adults, respectively (Wilcoxon Tests with BH correction: pupae: p = 0.012; adults: p = 0.0014). Neither dsx F nor dsx M expression differed significantly between the two male morphs (dsx F : pairwise t-test with BH correction: pupae: p = 0.480, Expression of the DM domain-encoding exon in adults is higher in males than in females, and higher in winged than in wingless males. (C) dsx F is significantly higher expressed in female than in male pupae, whereas dsx M is significantly higher expressed in males than in females (D). Worker pupae express significantly more dsx F than queen pupae, and adult workers express significantly less dsx M than adult queens (C+D). Letters indicate significant differences. Sample sizes are given in parentheses, boxplots show the median, interquartile ranges (IQR) and 1.5 IQR.
doi:10.1371/journal.pgen.1005952.g001 adults: p = 0.277; dsx M : pairwise Wilcoxon tests with BH correction: pupae: p = 0.481, adults: p = 0.805). However, overall expression of both isoforms was higher in winged compared to wingless males ( Fig 1B). Our finding that dsx is differentially expressed and alternatively spliced across morphs in pupae and adults suggests that dsx might play a role in controlling polyphenic development.
Tissue specificity of dsx splicing in sex mosaics confirms function To confirm that expression of dsx isoforms corresponds with phenotypic tissue differentiation, we used qPCR to analyze dsx M and dsx F expression in male and female-typical tissues dissected from aberrant "sex mosaic" individuals that express both male and female characters. C. obscurior sex mosaics are typically laterally separated into female and male halves, indicating that intersexuality is caused by single, early developmental aberrations such as anomalous fertilization events, loss of sex locus expression or inheritance of maternal effects [26][27][28]. The expression of dsx F and dsx M was male-typical in male tissue and female-typical in female tissue for all individuals except one, which had similar levels of dsx M in both tissue types (S3 Fig). As in previous studies [29,30], we only observed individuals possessing queen and winged male traits, or worker and wingless male traits; other trait combinations were absent (S4 Table), implying that common mechanisms control morph differentiation in males and females.

Co-option of sex-specific isoforms in larval morph differentiation
We analyzed published RNAseq data [31] from individual early 3 rd instar larvae (QU, EM, WM, WO; n = 7 each) on an exon-level with DEXSeq [32]. We found morph-biased expression in each of the seven dsx exons, and confirmed sex-specific expression of the DM domain, dsx F , and dsx M in the early 3 rd larval stage (S4 Fig and S5 Table). Overall, dsx expression was higher in males than in females, and higher in wingless morphs compared to winged morphs (EM > WM, WO > QU).
We hypothesized that other genes with sex-specific alternative splicing have been similarly co-opted for morph differentiation. Using a conservative false discovery rate of 0.005, DEX-Seq analysis identified 179 exons of 91 genes with sex-biased expression (S6 Table). Dsx exon 5 (= dsx F ) is ranked 5 th among the top 10 differentially expressed exons and exons 6 and 7 (= dsx M ) are the two most significant differentially expressed exons across all samples. To test for co-option of this set of exons into morph differentiation, we performed a hierarchical clustering analysis based on log-transformed exon counts. Queens and workers, as well as winged and wingless males, were clearly separated by the set of sex-biased exons, with the exception of two male samples that clustered with the wrong male morph (bootstrap node support: QU/WO = 75, WM/EM = 68) (Figs 2, S5 and S6 for bootstrap support for all nodes). Because terminal switch points for morph differentiation in male and female larvae may differ [31], misclassification of two male samples (WM34 & EM29) in hierarchical clustering may reflect higher plasticity in males compared to females at this particular developmental stage. Accordingly, in C. obscurior 3 rd instar larvae, more genes are differentially expressed between queens and workers than between winged and wingless males [31].
To identify the sex-biased exons that most strongly affect separation between sexes and morphs, we performed a principal component analysis (PCA) of the 179 normalized exon counts. PC 1 separated sexes (29.9% explained variation), PC 2 (15.3%) and PC 4 (6.8%) separated female and male morphs, respectively (Fig 3; linear discriminant analysis using Wilk's test on PCs 1, 2 and 4; factor sex: F (1,28) = 95.81, p < 0.001; factor morph: F (3,28) = 27.70, p < 0.001), while PC 3 (7.7%) did not separate between sexes or morphs (linear discriminant analysis using Wilk's test on PC 3; factor sex: F (1,28) = 0.06, p = 0.80; factor morph: F (3,28) = 1.81, p = 0.17). From the 179 exons, we identified those with the strongest influence on sex (PC 1), female morph (PC 2), and male morph (PC 4) by extracting the exon loadings that fell in either the 10% or 90% quantiles for each PC (S6 Table). Using these lists, we identified dsx (replicating the RT-qPCR results) and seven other genes that showed both sex-specific and morph-specific alternative splicing, of which kr-h2 has a putative transcription factor function  (Table 1). All eight genes are conserved across the Insecta, and a Gene Ontology (GO) term enrichment analysis with topGO [33] suggests that they serve basic metabolic and other core functions (S7 Table).

Discussion
The role of dsx in mediating development Our study suggests provides evidence that the sex differentiation pathway has been co-opted to control morph-specific development, as we predicted from the theory of facilitated variation. The major candidate gene dsx was alternatively spliced in males and females, and differentially expressed between queens and workers and between winged and wingless males. We independently replicated these results using qRT-PCR and RNAseq data from different individuals and life stages. Strikingly, we found that exons showing sex-biased expression were also differentially expressed between morphs, suggesting that dsx and other sex-biased genes mediate polyphenism within each of the sexes. The RNAseq analysis conservatively identified eight genes that have sex-specific and morph-specific alternative splicing; all of these genes were evolutionarily conserved and had GO terms associated with basic cellular functions. While dsx encodes sex-specific transcription factors and co-ordinates expression of a large number of downstream genes [34], except for a putative role of kr-h2 (see below) the other genes exhibit no transcription factor function. We confirmed that the sex-specific isoforms of dsx correlated with tissue type by analyzing male and female-typical tissue dissected from aberrant sex mosaic individuals. Finally, we reaffirmed that sex mosaics are always either hybrids of a queen and a winged male, or a worker and a wingless male, implying common morph differentiation control mechanisms in both sexes, especially regarding winglessness.
Interestingly, dsx has been shown to be a central hub gene involved in generating evolutionary novelty and polyphenism in other taxa. In a butterfly, genetic variation in dsx is associated with a heritable female-limited wing color/shape polymorphism, suggesting that dsx has been co-opted to control a novel, female-limited trait as well as maintaining its function in sex differentiation [24]. In the genus Drosophila, new localizations of dsx are thought to have facilitated the evolution of a novel male-limited trait (the sex combs), highlighting how the preexisting sex determination system was co-opted to produce a new polyphenism [35]. In the dung beetle Onthophagus taurus, RNAi experiments suggested that variation in dsx splicing mediates the difference in the presence of horns between males and females, and also controls a nutritionally dependent, male-limited polyphenism between large-horned and small-horned males [36]. A subsequent study of another horned beetle showed that different dsx isoforms control the sensitivity of the mandibles to juvenile hormone (JH), such that male mandibles are stimulated to grow by JH while those of females are not [37]. Thus it appears that dsx first evolved to mediate male-limited expression of horns by elevating the sensitivity of male horn tissue to JH [37] and perhaps also the IGF signaling pathway [38], and was then secondarily co-opted to control a nutrition-sensitive, male-limited polyphenism. The beetle dsx data are thus highly congruent with the theory of facilitated variation: the male polyphenism evolved using pre-existing genetic switches and developmental mechanisms to link a novel combination of stimuli and outputs (here, larval nutrition and horn phenotype).
Pre-and posttranscriptional genetic tools are not yet well established in ants but there is circumstantial evidence for similar links between dsx and JH in C. obscurior. A previous experiment showed that JH is involved in the development of larvae of both sexes into winged morphs [39], and the present study found differences in dsx splicing and expression between winged and wingless morphs. Thus, we speculate that the isoforms of dsx may mediate the responsiveness of developing tissues to JH, as hypothesized for beetles [37]. Significant differences in feminizer expression between queens and workers in the stingless bee Melipona [40], an upstream signal of dsx in bees [41], likewise suggests co-option of sex differentiation genes into caste differentiation in bees. There are no homologues of csd and feminizer in ants, because csd evolved in the Apis lineage by duplication of feminizer [42]. In ants, the closest homologue to feminizer is transformer. In C. obscurior, we could not detect morph-specific expression of the two transformer paralogues (tra1: Cobs_03145 and tra2: Cobs_18309), although they were expressed in a sex-specific manner.
In addition to dsx, we found a second sex-biased transcript with putative regulatory function. This ortholog to kr-h2 was alternatively spliced in queen and worker larvae, rendering the Kruppel homolog family a promising candidate for modulating plastic responses to the environment. kr-h2 has structural similarity to the JH-inducible transcription factor kr-h1 [43], which is involved in the initiation of metamorphosis in other insects [44,45]. kr-h2-induced differences in developmental timing may explain why metamorphosis is delayed in queens compared to workers [39], and further points to a link between sex-specific transcription, function in transcriptional regulation, sensitivity to JH, and evolutionary co-option into within-sex polyphenism.

An extended evo-devo framework for social insect polyphenism
We believe that the hypothesis advanced here, i.e. co-option of sex differentiation pathways into social insect caste polyphenism, is complementary to a previous theory regarding the proximate mechanisms underlying the origin of eusociality, termed the reproductive groundplan hypothesis (RGPH). Based on the ovarian ground plan hypothesis [46], the RGPH posits that eusociality arose via changes in the regulation of pre-existing gene sets relating to reproductive physiology and behavior, for example when genes involved in nest provisioning and brood care began to be expressed in unmated, non-reproductive individuals [47]. Research on the RGPH has stressed the importance of genes with nutrition-sensitive expression in delimiting the queen and worker "genetic toolkits", in light of evidence that caste fate is nutrition-sensitive [48], that diet preference, reproduction and behavior are pleiotropically linked [49], and that some nutrition-related genes such as IRS and TOR influence caste fate [50]. Juvenile hormone, which is involved in regulatory feedback loops with some nutrition-related gene networks, has also been linked to caste differences [48,50]-including in our study species C. obscurior [31,39]-as well as to within-caste polymorphisms (e.g. [51]). Our hypothesis and the RGPH both argue that regulatory evolution caused conserved genes to acquire caste-specific expression. Our hypothesis is distinct in that it explicitly proposes that this regulatory evolution takes place in sex differentiation genes, but leaves the targets of these genes unspecified. By contrast, the RPGH makes predictions about which gene networks produce caste-biased phenotypes (e.g. ovary development, [52]), but makes no prediction regarding the identity of the regulatory sequences controlling these networks. Thus, the hypotheses do not overlap, and both may be correct. Analyses of potential regulatory links between the pathways presented here and those implicated with the RGPH will reveal to what extent they are connected.

Conclusion
Co-option of conserved genes involved primarily with sex differentiation in novel contexts allows functionally integrated gene networks to produce discrete phenotypes. Together with the horned beetle data reviewed above, our study suggests that core components of the sex differentiation pathway such as dsx can produce evolutionary novelty by acting as a switch for nutrition and JH-sensitive growth and development. Although many mechanisms of gene regulation have been implicated in controlling caste-specific development in social insects (e.g. methylation [53], transcription factors [31], small RNAs acting post-transcription [17], RNA editing [54] or structural chromatin modification [55]), all of these depend on some higherlevel genetic switch to trigger differential activity. We propose that highly conserved hub genes such as dsx, which can translate variable input signals into large transcription differences using intermediate-level regulators, were the most basic mechanism responsible for the repeated evolutionary transition to eusociality and caste polyphenism.

Verification of haplodiploid genetic sex determination
Crosses between five queens of a C. obscurior population from Japan (JP) and five wingless males of a C. obscurior population from Brazil (BR) were set up by placing sexual pupae together with some brood and~20 workers in plaster-filled Petri dishes. Nests were checked twice a week, provided with water, honey and pieces of dead insects and kept at constant conditions (12h 28°C light, 12h 23°C dark). We sampled emerging F1 hybrid QU, WO, EM and WM pupae and extracted DNA from the 10 parental and 71 F1 individuals (23 EM, 3 WM, 22 QU, 22 WO, 1 GY = gynandromorph, for sample sizes per family see S1 Table). Each individual was analyzed at three variable microsatellite loci (Cobs_1.1, Cobs_8.3, Cobs_8.4; for primer sequences see S8 Table). PCRs were performed using the BIO-X-ACT Short Mix (Bioline) and microsatellite analyses were carried out on an ABI PRISM (Applied Biosystems).

Verification of functional dsx and its sex-specific isoforms in C. obscurior
To find the functional dsx ortholog of C. obscurior, we identified DM domain-containing proteins of Drosophila melanogaster, Nasonia vitripennis, Apis mellifera, Pogonomyrmex barbatus, Acromyrmex echinatior and C. obscurior by BLASTp and tBLASTn analyses (S9 Table) and aligned them with MUSCLE [56]. We extracted the DM domain region from the manually corrected alignment (S7 Fig) and built a phylogenetic tree in MEGA [57], applying a WAG+G+I phylogenetic model and bootstrap resampling with 1,000 replicates (S1 Fig).
We reanalyzed previously published RNAseq data of larvae [31]. After removing adapter sequences with cutadapt and performing quality filtration with Trimmomatic, the reads were mapped against the reference genome with tophat2 (v2.0.8) and bowtie2 (v2.1.0) in sensitive mode. We generated count tables with HTseq based on the Cobs1.4 official gene set and used DESeq2 [58] to assess sex-specific expression of the four dsx paralogs following size factor normalization.
We applied RACE (Rapid Amplification of cDNA Ends) for identification of dsx isoforms. Total RNA was extracted from three females (QU adult, QU pupa, WO pupa) and three wingless males (one pupa, two adults) using the peqGOLD MicroSpin Total RNA Kit (peqlab). Transcription to cDNA was performed with the AffinityScript Multiple Temperature cDNA Synthesis Kit (Agilent Technologies), using the 3' RACE Adapter GCGAGCACAGAATTAA TACGACTCACTATAGGTTTTTTTTTTTTVN. 3' RACE was performed in a nested PCR using two gene-specific 3' primers (dsx4_for4, Co_dsx_p3_for, for primer sequences see S8 Table) and the 5' primer provided in the First Choice RLM-RACE Kit (Ambion). PCRs were performed using the BIO-X-ACT Short Mix (Bioline) with the following protocol: 94°C (3 min), followed by 35 cycles 94°C (30 sec), 60°C (30 sec), 72°C (2 min) and a final elongation of 72°C (7 min). The products were purified with the NucleoSpin Gel and PCR Clean-up (Macherey-Nagel) and Sanger sequenced at LGC Berlin.

Expression of dsx and its isoforms in morphs and sex mosaics with realtime quantitative PCR (qPCR)
Total RNA was extracted from adults (8 EM, 8 WM, 10 QU, 10 WO) using the RNeasy Plus Mini Kit (Qiagen) and transcribed to cDNA using the AffinityScript Multiple Temperature cDNA Synthesis Kit (Agilent Technologies). Expression of the DM domain was quantified by qPCR using the primer pair dsx4_for4/dsx4_rev1 and normalized with two housekeeping genes (RPS2_new, RPL32; see S8 Table for primer sequences).
We further used qPCR to measure isoform-specific dsx expression (dsx M and dsx F ) in pupae and adults of all four morphs, and in tissue from four sex mosaic pupae. We dissected the head and thoraces of the sex mosaics (for morphological descriptions see S4 Table) laterally into male and female halves and stored male and female tissue parts separately in RNAlater-ICE (Ambion), resulting in one female and one male sample per individual. We extracted total RNA from 9-10 pupae and seven adults of each of the four morphs, and from the sex mosaic tissue using the peqGOLD MicroSpin Total RNA Kit (peqlab) including a DNA digestion step with the peqGOLD DNase I Digest Kit (peqlab). After cDNA synthesis with iScript cDNA Synthesis Kit (Bio-Rad) we quantified gene expression of dsx F and dsx M using isoform specific, intron-spanning primers (dsx F : 4for/F5rev, dsx M : 4for/M5rev; see Fig 1A for position of primers) and two housekeeping genes (RPS2_new, Y45F10D_JO1). All qPCR reactions were performed in triplicates (repeatability was uniformly high, so we took the mean of the three replicates prior to analysis). Data analysis was carried out according to [59], using the geometric mean of the two housekeeping genes for normalization.

Testing whether sex-biased exons differ between same-sex morphs
We analyzed published RNAseq data [31] from 3 rd star instar larval QU, WO, WM and EM (n = 7 each) and assessed differential exon-specific expression with DEXSeq [32]. Raw reads were trimmed and passed through quality filtration as described in [31] and mapped to the reference genome Cobs1.4 [19] using STAR [60]. We corrected the dsx and tra gene model using the RACE results for dsx, and split the tra gene model into two paralogs (tra1 and tra2), as observed in other ants [42,61,62]. For all other genes we used gene models of the Cobs1.4 official gene set. We followed the default workflow of DEXSeq and tested for differential exon usage between males and females based on a false discovery rate of 0.005. In the resulting 179 sex-specific exons, we tested for morph-specific exon profiles using hierarchical clustering (implemented by the R function hclust using the ward.D2 method [63]) of pairwise Manhattan distances between log-transformed normalized exon counts. We assessed the support for each node in the cluster analysis using bootstrap resampling with 10,000 replicates using the pvclust package in R 3.1.2.

Identifying genes with sex-and morph-specific exons
We conducted a PCA with normalized exons counts. We visually identified principal components that best separated between sexes (PC 1), female morphs (PC 2) and male morphs (PC 4) and confirmed that these components suffice to separate among sexes and morphs with linear discriminant analysis and subsequent Wilk's tests in R 3.1.2. Based on loadings of exons on each component, we identified exons that fell in either 10% or 90% quantiles (S6 Table) as those with the strongest influence on PC 1, PC 2 and PC 4. From this list, we extracted only those genes that contained multiple exons with strong influence on both sex (PC 1) and morph (PC 2 and/or PC 4). This yielded a list of eight candidate genes showing alternative splicing between sexes as well as morphs (Table 1).
Supporting Information S1 Table. Results of microsatellite analyses of F1 individuals from interpopulational crosses. For each family, a Japanese (JP) Cardiocondyla obscurior queen was mated with a Brazilian (BR) male. Emerging F1 individuals were genotyped using three population-specific microsatellite markers. This showed that all F1 males (EM = ergatoid, wingless males, WM = winged males) and one gynandromorph (GY) exclusively carried the maternal (JP) allele, whereas emerging females (QU = queens, WO = workers) carried both parental alleles (JP+BR). Sample sizes are given in parenthesis.