Molecular signatures of plastic phenotypes in two eusocial insect species with simple societies

Typical Insect Genome Composition and Organization.

A single haploid male for each species was sequenced on the Illumina platform achieving 110-fold coverage. The de novo assembled P. canadensis and D. quadriceps genomes were 211 Mega-basepairs (Mbp) and 268 Mbp in size, respectively (SI Text, sections I, II.1, and II.2). These genome sequences are almost complete, with 97–99% of the conserved cluster of orthologous proteins mapped in the two genomes; 79–86% of proteins were annotated (Fig. S1 A–D and SI Text, sections II.3–II.5). The genome compositions were similar to the genome sequences of other social insects, with D. quadriceps sharing more of its predicted protein content with other ants (Formicidae), whereas P. canadensis shows more equitable levels of protein sharing with ants (Formicidae) and bees (Apidae) (Fig. 1C; Fig. S1E; SI Text, section II.6; and Dataset S1). This difference is likely to reflect the absence of any other aculeate wasp genome sequence in the public domain. Finally, the genome of P. canadensis contains more transposable elements (452,247, 12% of the genome) than D. quadriceps (217,417, 6% of the genome), most of which are simple or low-complexity repeats (Fig. 1D and SI Text, section II.7). Transposable elements were recently identified as potentially important in the evolution of social complexity in bees (6).

Low Levels of Transcriptional Differentiation Between Phenotypes.

We obtained over 100 gigabase pairs (Gbps) of brain transcriptome sequence data from 23 individual adult female brains (four to seven biological replicates each of reproductives and nonreproductives per species), generating, on average, 3.6 Mbp (20.29 ± 0.67-fold coverage) and 4.9 Mbp (17.4 ± 1.36-fold coverage) per individual for the wasp and ant, respectively (SI Text, sections III.1 and III.2 and Dataset S2). In both species, we found fewer than 1% of genes were differentially expressed (DEGs), with little evidence of functional specialization between phenotypes (5). Using the union of DEGs from EdgeR [parametric approach (16)] and NOISeq [nonparametric approach (17)] (Fig. 2; Table 1; SI Text, section III.3; and Dataset S2), we found 67 (0.4%) DEGs in P. canadensis and 147 (0.8%) DEGs in D. quadriceps. In both species, the nonparametric approach identified significantly more up-regulated genes in reproductives relative to nonreproductives (χ² = 31, P = 2.2e-08; Table 1). In P. canadensis, gene expression in nonreproductives was found to be more stochastic (noisy) than in reproductives despite similar variance of expression among the biological replicates (Fig. S2). Recent research suggests that evolution can shape noise in gene expression and that such noise can be adaptive and heritable (18–20). If noise in transcription is an indicator of phenotypic plasticity (21–23), our results would suggest that transcription in the nonreproductive phenotype is more responsive to changes in the biotic and social environment than transcription in the reproductive phenotype. Despite the small number of DEGs, significant functional enrichment of DEGs was detected in the ant reproductives, with 29 gene ontology terms significantly enriched for functions that included metabolic and ribosomal processes, regulation of expression, and an extracellular component [false discovery rate (FDR) < 0.5; SI Text, section III.4 and Dataset S2]. There was little sign of functional enrichment in the wasp (5) (although before FDR correction, oxidoreductase activity and lipid transport were overrepresented in reproductives). These data suggest there is little phenotypic specialization in the brain tissue of either species.

No Distinct Methylation Patterning Across the Genome or Between Phenotypes.

We sequenced the methylomes from three biological replicates each of individual adult brains from reproductive and nonreproductive phenotypes in P. canadensis and D. quadriceps using whole-genome bisulfite sequencing [BS-seq; 20 gigabase (GB) (>10-fold coverage) per brain] (SI Text, section IV.1 and Dataset S3).

We compared methylation patterns with the honey bee (24) to provide a reference point because the honey bee is the only close relative to our study species with comparable data on brain methylation available (SI Text, section IV.2). Global levels of methylation in the cytosine-guanine (CG) context were similar in both species, and similar to the honey bee (Table 2). P. canadensis exhibited greater methylation in the non-CG context but significantly fewer highly methylated regions than D. quadriceps (Table 2; Fig. S3 A and B; and SI Text, section IV.3). However, in comparison to the honey bee, both species showed relatively little gene body-specific methylation targeting (Fig. 3A; Table 2; Fig. S3C; and SI Text, section IV.4), together with a striking lack of consistently fully methylated cytosines (Fig. 3B). In both P. canadensis and D. quadriceps, DNA methylation is dispersed sparsely across genes (Fig. 3C), particularly in P. canadensis, whose genome lacks a DNA methyltransferase 3 (DNMT3) gene, an enzyme involved in de novo methylation (25, 26) (Fig. S4A and S5A and SI Text, section IV.5). In P. canadensis, we also found a prevalence of asymmetrical (one strand only) CG methylation, together with a variant of the DNA methylation transferase 1 (DNMT1) gene involved in the maintenance of DNA methylation, in its genome (Fig. S4 B–D and SI Text, section IV.5). As observed in A. mellifera brains (27), both study species possess and express a ten-eleven translocation (TET) hydroxylase gene and base excision repair genes involved in demethylation, and have detectable hydroxymethylation in brain tissue (Figs. S4F and S5 and SI Text, section IV.6). Together, these general features provide an epigenetic landscape that may facilitate plasticity of genome function.

Despite the general paucity of methylation patterning, we found significant conservation of methylated orthologs (Fig. 3D; SI Text, section IV.7; and Dataset S3) and a positive correlation between gene expression and CG methylated genes (Fig. S6 and SI Text, section IV.2), as seen before in other insect species (28–33). Notably, however, DEGs tended to be hypomethylated in both species (Fig. 3E), and unlike the case in brain methylomes of adult honey bees (24, 34, 35), we found no evidence that phenotypes were associated with differentially methylated genes in our two species (t test, P > 0.05; SI Text, section IV.8). Analyses of alternative splicing revealed only 28 phenotype-specific isoforms expressed in D. quadriceps and none in P. canadensis (SI Text, section IV.9 and Dataset S3). This similarity between phenotypes is likely due to the global tendency of these species to express all isoforms simultaneously (Fig. 3F). Similar to DEGs, alternatively spliced genes (ASGs) were also hypomethylated compared with non-ASGs (Fig. 3E). This result may limit the role of DNA methylation in regulating phenotype-associated gene expression or alternative splicing in our species, and contrasts with what has been described in the honey bee (26, 35–39).

MicroRNAs Are Not Preferentially Targeting Differentially Expressed Genes.

Species-specific miRNA libraries were constructed from pools of individuals to include each phenotype to determine whether large numbers of miRNAs are shared between hymenopterans to the exclusion of the other insects and to identify potential cis-regulatory elements of DEGs. From our miRNA libraries, we identified 159 miRNA families (73 in P. canadensis and 86 in D. quadriceps), including 15 previously undescribed families (Fig. S7; SI Text, section V; and Dataset S4). We identified four families that are unique to hymenopterans and an additional nine families that were shared by apocritans to the exclusion of Nasonia and other insects. We found that miRNAs (40) were not preferentially targeting phenotype-specific DEGs, because although some DEGs appeared to be highly targeted, others were not (Dataset S4). Further work is needed to investigate miRNA expression levels in large numbers of individual queens and workers to rule out a role for miRNAs in caste differentiation.

Role for Conserved Toolkit Genes and Taxon-Restricted Genes in Regulatory Networks.

Despite the low numbers of DEGs, we found evidence that DEGs were nonrandomly organized at the network level in both species. Weighted gene correlation network analyses identify groups of genes that covary significantly in expression as “modules” (41). These analyses identified 31 and 41 gene coexpression networks for the ant and wasp, respectively (SI Text, section SIII.5 and Dataset S5). DEGs were clustered nonrandomly across networks in both species (Fig. 4A). Only three (10%) and two (5%) network modules showed significant overrepresentation of DEGs in the ant and wasp, respectively, and only one network module in the ant showed evidence of functional enrichment for ribosomal terms (SI Text, section III.5). Phenotype-specific transcription in both species is therefore governed by subtle but coordinated coexpression networks.

There is a debate over the relative roles for core sets of conserved genes (42–48) and taxon-restricted genes (TRGs) (5, 44, 47, 49, 50) in the evolution of convergent phenotypes (7, 44, 46). We found evidence that both types of gene classes play peripheral roles in the molecular networks associated with phenotypic differentiation in our study species. In each species, we identified both classes of genes among DEGs, determined whether their functions were conserved, and ascertained their putative importance in the gene networks associated with phenotypic differentiation. There were significant levels of overlap in the identity of DEGs between the two species (reciprocal BLAST hits of DEGs; n = 11 genes, P < 0.003 relative to chance for both species; SI Text, section III.3 and Dataset S2), suggesting they are homologs. Some of these genes were the same as those genes that had been previously identified as conserved “toolkit” genes for alternative phenotypes in eusocial insects [e.g., cytochrome P450, vitellogenin, hexamerin-2, kruppel homolog 1 (42–48)], but others were not (e.g., fibrillin-like gene; glutaminase, esterase, and myrosinase enzymes; a gene coding for a lysozyme). Gene identity may be conserved, but not the direction of expression (5, 7): Four of 11 genes were worker-biased in the ant, whereas all 11 were queen-biased in the wasp (Dataset S2). Finally, conserved DEGs were not generally highly connected in the coexpression networks of either species (Fig. 4 C–F). This observation contrasts with eusocial insect species with phenotypes that are determined irreversibly during development, where conserved genes can play central roles in gene networks (44).

TRGs (those genes having no significant homologs in available genomic databases) were detected in DEG sets in both species (ant: 10%, n = 16; wasp: 7.5%, n = 5) (Fig. 4 C–F, Table 1, and Fig. S8) and at similar levels to TRGs across the whole genome [ant (11.6% TRGs): χ² = 0.11, P = 0.74; wasp (9.1% TRGs): χ² = 0.52, P = 0.47; Dataset S5]. Taxon-restricted DEGs are likely to be new genes (short in length relative to annotated/known genes) (49) (Fig. S8) with unknown/novel functions (“guilt-by-association” network analysis) (41), because their nearest neighbors were also taxon-restricted (unknown function) (mean = 2.3 of the 10 most connected genes had BLAST hits; Fig. S8 and Dataset S5). Finally, taxon-restricted DEGs had similar low levels of connectivity to conserved genes in the networks of both species (generalized linear model) (ant: binomial errors, P = 0.89; wasp: quasibinomial errors, P = 0.96; Fig. 4 C–F), suggesting that conserved genes and previously unidentified TRGs are similarly important in phenotypic differentiation in these two species.

These data support the emerging hypothesis that conserved genes, new genes, and/or new regulatory networks are important in the evolution of phenotypic diversity (5, 44, 47–51). Our analyses add to this hypothesis by identifying roles for both conserved genes and TRGs in highly plastic phenotypes.

I.1. P. canadensis.

P. canadensis wasps were collected from 10 wild colonies in Panama in June 2012 at Punta Galeta, Colon (9°21′30.877′′ N, 79°54′0.0053′′ W), under the field collection permit SE/A-20-12 from the Autoridad Nacional del Ambiente. All colonies contained a single queen (henceforth referred to as “reproductives”) that laid all of the eggs, was socially dominant, and rarely left the nest, and many workers (henceforth referred to as “nonreproductives”) that provisioned their sibling brood on their natal colony (63). Colonies were medium-sized (20–40 wasps) postemergence colonies (meaning nonreproductives are the daughters of the reproductives), and had been monitored for several weeks, during which time there had been no change of reproductives. Reproductive phenotypes were identified from observations of reproductives laying eggs, being socially dominant, and rarely leaving the colony, and nonreproductive phenotypes were identified as subordinate individuals that were observed bringing forage prey to the colony and feeding brood. Foragers tend to be the older nonreproductive individuals in their nests (1). Reproductive and nonreproductive wasps were collected directly off their nests with forceps at midday during the active foraging periods; their heads were cut off and immediately placed into RNAlater solution (Ambion), and their bodies were stored in 80% EtOH and kept at −20 °C until dissection of brains and ovaries. Ovary dissections confirmed that reproductives had mature eggs in their ovaries and had mated, whereas nonreproductives had no ovary development but some had mated, as expected in this species (63) (Dataset S2).

I.2. D. quadriceps.

Seven wild colonies of D. quadriceps were collected from the Atlantic Forest in Sergipe (11°01′23′′ S, 37°12′9′′ W) and Campo Formoso (10°2′972′′ S, 40°20′771′′ W) in Brazil from 2009 to 2011 under collection permits 10BR004553/DF and 11BR006471/DF from the Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais. Because this species nests deep in the ground, colonies had to be excavated and kept in the laboratory for several weeks to identify phenotypes. Each ant was fitted with a radio-frequency identification tag that allowed foragers (those ants that left the nest regularly to find forage) to be identified. The dominance rank of each individual was determined using multiple 30-min behavioral observation periods, during which time the type of interaction and the identities of the actor and recipient were recorded. Dominance hierarchies were then constructed from aggressive interactions for each colony (64); high-rank phenotypes were identified by frequent, high-intensity aggressive interactions, and low-rank phenotypes were identified by only rare involvement in aggressive interactions. Low-ranked individuals tend to be the older nonreproductive individuals in the nest (12). Ants were killed in liquid nitrogen; their heads were then transferred into RNAlater solution, and their bodies were transferred to 80% (vol/vol) ethanol for storage at −80 °C until dissection of brains and ovaries. Dissections recorded whether an individual had mature eggs in its ovaries and whether it had mated (Dataset S2). Gamergates (henceforth referred to as reproductives) were identified as the top-ranked dominant individual in their colony, rarely left the nest chamber (i.e., did not forage), had mature eggs in their ovaries, and had mated. Nonreproductives that foraged were identified as some of the lowest ranked (least dominant) females; they frequently left the nest (i.e., foraged), had not mated, and did not have mature eggs in their ovaries.

II.1. DNA Extraction.

DNA was extracted from a single haploid male of P. canadensis and D. quadriceps using a BioSprint 96DNA Blood Kit from QIAGEN. This kit uses the Mag Attract (QUIAGEN) magnetic particle technology for DNA purification. DNA was eluted in buffer AE (QUIAGEN). Microsatellite analysis was conducted to confirm male haploidy.

II.2. Genome Sequencing and Assembly.

DNA from whole single haploid males was sheared in a Covaris S2 instrument, and Illumina genomic libraries with a peak insert size of 513 bp (P. canadensis) and 542 bp (D. quadriceps) were sequenced on the Illumina HiSeq sequencing system using a 2 × 100-nt paired-end sequencing protocol. Mate-pair libraries with span size of 3.9 kb (P. canadensis) and 5 kb (D. quadriceps) were prepared from the same material. Mate-pair reads were obtained from 2 × 50-nt paired-end sequencing on the Illumina HiSeq platform. Data filtering and genome assembly of Illumina reads were as previously described (65). The soap de novo assembler (version 1.05) with a k-mer size of 49 was used. D. quadriceps mate-pair data contained a high percentage of read-pairs not derived from fragments encompassing the circularization junction. Such nonjunction pairs align in the same orientation and distance as regular Illumina paired-end reads. Nonjunction pairs were removed before scaffolding. To this end, 166.6 million mate-pairs were aligned against a D. quadriceps preassembly that had been calculated based on paired-end data. Read-pairs were kept that mapped uniquely with each read (81.2 million pairs). The mapping distance of read-pairs scattered around 5,000 bp, which is expected for the true mate-pairs, and around 355 bp. The mapping distance of 355 bp was the library insert size, and read-pairs with this distance represent the nonjunction pairs. Thereafter, all pairs with a mapping distance less than 1,000 bases were removed (35.0 million pairs remaining). To reduce potential errors in the scaffolding step, we further discarded duplicated read-pairs based on mapping positions. For the final D. quadriceps assembly, 3.0 million mate-pairs were used. For P. canadensis, an assembly before gap closing was deployed for downstream analysis. For D. quadriceps, we used a gap-closed assembly. Assembly metrics were calculated on all scaffolds and nonscaffold contigs sized 500 bp or larger.

II.3. Estimation of the Gene Space Completeness of the D. quadriceps and P. canadensis Genome Assemblies Using the Core Eukaryotic Genes Mapping Approach Pipeline.

We estimated the completeness of the assembly by looking at a set of “core” orthologous proteins [or “cluster of orthologous groups” (COG)] that were deemed to be highly conserved and exist in low-copy numbers in the large majority of higher eukaryotes. The latest COG set of genes includes 248 sets of “COG” proteins from six divergent eukaryotic species. Even though sequence coverage and the “N50” parameter (a measure of genome assembly quality; 536 kb and 1,359 kb in P. canadensis and D. quadriceps, respectively) are useful for characterizing a genome with regard to “base pairs,” they do not always describe the status of the gene space or whether one will be able to detect and identify genes easily and accurately in the genomic sequence. The Core Eukaryotic Genes Mapping Approach (CEGMA) pipeline has been used to determine the state of the gene space (i.e., another indicator of genome completeness) of several species in the past few years. Most (nearly) complete assemblies [using whole genome sequencing (WGS) technology] include 90% of the 248 COG proteins.

In the case of the D. quadriceps assembly, the CEGMA pipeline was able to map 99.19% (246 of 248) of the conserved COGs at least partially and 97.58% (242 of 248) over their entire length. In the P. canadensis assembly, the CEGMA pipeline mapped 98.79% (245 of 248) of the COGs at full length and 99.6% (247of 248) at least partially. The high percentage of completely detected COGs is a strong indicator that the assemblies, in terms of their “gene space,” can be considered as being of very high quality. Therefore, they were judged to be robust enough for all subsequent steps, leading to obtaining a protein-coding reference gene annotation for both species.

II.4. Obtaining a Protein-Coding Reference Gene Annotation for the P. canadensis and D. quadriceps Genomes.

Protein-coding gene annotation reference sets for the P. canadensis and D. quadriceps genomes were generated by combining the Program to Assemble Spliced Alignments (PASA r2012-06-25) and Evidence Modeler (EVM r2012-06-25) to obtain coding DNA sequences (CDS) models using three main sources of evidence: (i) aligned transcripts, (ii) aligned proteins, and (iii) gene predictions.

II.5. Functional Annotation of the P. canadensis and D. quadriceps EVM-Based Annotation.

Functional annotation was performed using InterPro, the Kyoto Encyclopedia of Genes and Genomes (KEGG), Blast2GO, and Reactome databases; information about domains, functions, and gene ontology (GO) terms was inferred from the sequences in the repositories to the proteins of interest based on sequence similarity. InterProScan (version 5) was used to scan through all available InterPro databases. BLAST alignment searches against the NCBI NR collection of protein sequences (release 2013-04) were used as input to the local Blast2GO software p2gpipe (version 2.5.0) with the databases go_201303 (March 2013). KEGG orthology (KO) groups were assigned by the KEGG Automatic Annotation Server using the bidirectional best hit method against a representative gene set from 29 different species, including D. melanogaster, A. mellifera, and Nasonia vitripennis. KO identifiers were then used for retrieval, utilizing the KEGG REST-based API service for the KEGG relevant functional annotation, such as metabolic pathways and external database references. Sets of one-two-one orthologs from PhylomeDB were used via the Reactome URL/XML query system to scan the ReactomeDB (version 46, September 2013). Additionally, SignalP was used to predict the presence and location of signal peptide cleavage sites in our proteins of interest.

A total of 86.8% of proteins predicted for D. quadriceps and 79.7% of proteins found in P. canadensis were found to have some type of annotation feature (Dataset S1). Descriptions (name) were assigned to 47.9% (D. quadriceps) and 43.4% (P. canadensis) of the proteins using Blast2GO or the KEGG. GO terms, derived from four different sources, were assigned to 62.4% (D. quadriceps) and 54.5% (P. canadensis) of the annotated proteins. GO terms were grouped by the different functional categories (Fig. S1) using CateGOrizer, utilizing GOSlim without top-level categories (molecular function, biological process, and cellular component).

Unannotated proteins (13.2% of the reference annotation for D. quadriceps and 20.3% for P. canadensis) tend to be shorter compared with the annotated ones (Fig. S8). These smaller proteins could be derived from incomplete ORFs erroneously predicted by the automatic annotation pipeline due to genome fragmentation.

II.6. Comparison of Genomes.

The P. canadensis and D. quadriceps genomes differed in the composition of genes and proteins, as expected from phylogenetically distinct species. We constructed a phylome resource for each species. From this phylome resource, we identified orthologous genes and proteins with genomic databases (P. canadensis: www.phylomedb.org/phylome_134, D. quadriceps: www.phylomedb.org/phylome_135). D. quadriceps showed 70% similarity with proteins in the other ants (Formicidae) and 46% of all proteins shared with the only other sequenced ponerine ant (H. saltator) (Fig. 1C). Conversely, P. canadensis shares similar levels of protein similarity with the bees (Apidae, sharing 34%) and ants (Formicidae, sharing 22%). Less than 25% of proteins were the same in the two species. However, the top 10 GO terms and the GOSlim terms were similar between the two species (Fig. S1).

II.7. Phylogenomic Relationships in the Aculeate Hymenoptera.

We reconstructed the phylogeny of the aculeate Hymenoptera using only the available data on published genome sequences for social insects, including P. canadensis and D. quadriceps. The discussion on whether more taxa or more genes are better is not settled, and there are arguments for both. Because recent multigene analyses of Hymenoptera have challenged the phylogenetic relationships based on a few genes (67), we chose to maximize the opportunity of analyzing phylogeny across hundreds of genes. A maximum likelihood tree was produced based on the combined analysis of a concatenated alignment of 897 single-copy genes present in all species. This tree placed ants as the basal group to wasps and bees (Fig. S1E). This result suggests that the wasp and ant do not belong to the same historically named family Vespoidea (5, 67), and that, instead, wasps share more ancestry with bees than with ants. In support of this conclusion, P. canadensis shared more similarity with Apoidea (bee) proteins than Formicidae (ant) proteins. The tree was reconstructed using randomized accelerated maximum likelihood utilizing the PROTein, GAMMA distributed rates, LG (Lee and Gascuel) model (PROTGAMMALG). Bootstraps lower than 100 are indicated in the nodes of the tree. We reran the analysis with a smaller dataset of 93 genes previously obtained for P. canadensis (5) and found the wasp to be basal to the ants and bees, as also suggested by a recent transcriptome-based analysis of the hymenopteran phylogeny (67). This analysis indicates that topology changes depending on the number of genes used, suggesting that the choice of genes strongly influences resulting topology.

II.8. Transposon Annotation.

Repeats were detected using RepeatMasker (version 4.0.1) with repeat library 20120418. Searches were limited to repeats deriving from Apocrita species. Hits were grouped according to the annotated repeat class, and the grouped classes were analyzed separately (Fig. 1D).

III.1. Transcriptome Sequencing.

Total RNA was extracted from single brains using the QIAGEN All Prep DNA/RNA Mini Kit according to the manufacturer’s instructions. Where possible, samples of reproductives and nonreproductives were matched within nests. For RNA-seq libraries, 50 to 200 ng of total RNA was enriched for mRNA using Dynabeads Oligo(dT)₂₅ from Invitrogen in two subsequent steps of purification with fresh beads. Fragmentation was done by incubation of mRNAs for 5 min at 94 °C in the First-Strand Buffer (Invitrogen), and directly followed by cDNA synthesis using a SuperScript III Kit (Invitrogen) according to the manufacturer’s instructions. dUTPs were incorporated for second-strand synthesis for library orientation. cDNA was end-repaired, A-tailed, and ligated with a methylated Adaptor Oligo Kit (Illumina) using the NEB Next Kit (New England Biolabs) according to the manufacturer’s instructions. dUTP excision was done before amplification using USER mix (New England Biolabs). Libraries were amplified with 16 cycles using 2× Phusion HF buffer (New England Biolabs). Size selection and cleaning between steps were performed with the AMPure XP system (Agencourt) to select DNA fragments between 250 bp and 500 bp. Paired-end libraries were sequenced on either an Illumina GAIIx or Illumina HiSeq system.

III.2. EVM and de Novo Assembly.

The EVM/PASA approach (SI Text, section II.5) provides a comprehensive protein-based annotation from multiple sources of data but can miss transcripts due to tissue-specific transcription of some genes, low-coverage genes, and noncoding elements. To obtain the most comprehensive transcriptome annotation, we performed additional de novo assemblies with two different and widely used methods (Cufflinks and Trinity) and merged these assemblies with the EVM/PASA annotation. Using Cufflinks, we started with the independent assembly for each sample, without providing any annotation and with the default parameters. In the next step, we merged all of the individual annotations plus the EVM annotation (described in SI Text, section II.5) by using the Cuffmerge program. We then proceeded by removing duplicate transcripts from this merged annotation. We used the program Gffread from the Cufflinks package (−M and −d) to collapse matching transcripts and obtain a final set of transcripts without duplicates.

Next, we conducted a de novo transcript assembly using Trinity. This procedure recovered additional transcripts that were not found in the previous analyses. Additional transcripts were added to the final assembly, where no existing sequence was previously described. The process involved the use of “in-house” scripts to automate raw read filtering with a “fastq_quality_filter” (Hannon laboratory FASTX-Toolkit, hannonlab.cshl.edu/fastx_toolkit/) to remove 5′ and/or 3′ regions with Phred scores less than 30 and ribosomal RNA removal using the SortMeRna (version 1.99 beta) pipeline, using Insecta 28S and 18S databases downloaded from SILVA (www.arb-silva.de/). Quality checking was performed by FastQC. Next, Trinity was used with default settings, except for −SS_lib_type FR [forward/reverse (mate-pair data)]. This step was performed on a high-memory node of the Advanced Computing Research Centre (University of Bristol, www.bris.ac.uk/acrc/).

For P. canadensis, this procedure resulted in 15,595 genes with 53,181 transcripts. For D. quadriceps, 12,765 genes and 57,730 transcripts were found.

RSEM (RNA-Seq by Expectation-Maximization) was used to calculate transcript abundance quantifications for each individual by mapping individual RNA-seq reads (the same reads that were filtered for Trinity) against the combined EVM, Cufflinks, and Trinity merged assembly. This process created gene and isoform raw counts and fragments per kilobase per million (FPKM) reads.

III.3 Identifying Differentially Expressed Genes (DEGs).

III.4. GO Enrichment Analyses of DEGs.

For the GO analysis, we used the Blast2GO Fisher exact test (FDR < 0.5) and selected enrichments with at least five genes in the reference set. We used the assembled list of GO terms associated with each gene and tested for enrichment in all DEGs (SI Text, section III.3.1.1) and in specific coexpression modules (SI Text, section III.5) against a background of all GO terms for each annotated gene.

III.5. Network Analyses of Coexpressed Modules of Genes.

Genes alone give a very limited understanding of transcriptome dynamics. Grouping genes that covary in expression into modules can be more informative than analyzing individual genes. To obtain systems-view of gene expression, we performed coexpression network analysis by using Weighted Gene Coexpression Network Analysis (WGCNA) (41), which allowed us to group genes together that have similar expression patterns (i.e., covary) into “modules.” Genes with correlated expression (in a module) are assumed to have some common role in a particular pathway/functionality. To identify modules of importance in caste differentiation, we then correlated the connectivity of the genes within the module with a gene significance measure (log₁₀ q-value) for differential expression between castes. Genes that are differentially expressed between castes and are highly connected within a module (i.e., positive correlation) are predicted to be drivers of expression in that module, and may reflect an important molecular pathway in caste differentiation. We used this prediction to identify modules of putative importance in caste differentiation in each species and then compared them between species to look for evidence of conserved networks underlying caste differentiation. We tested for nonrandom distribution of DEGs across modules by determining whether the proportion of DEGs is drawn from the same binomial distribution or whether DEGs are clustered within modules. To do this test, we fitted to general linear models, one with a single parameter (the global proportion) and the other with a separate parameter for each of the modules (the saturated model). We compared the best fit of the two models for each species separately using ANOVA, and then also compared the two species.

From the combined assembly, only genes with expression above 1 FPKM across all samples were selected for network analyses (11,427 in P. canadensis and 11,758 in D. quadriceps). As a gene significance measure, we used the −log₁₀ q-value obtained from the q-value of the NOISeq differential expression test between castes for each species. Default settings were used, except “BETA = 8” in P. canadensis and “BETA = 12” in D. quadriceps. Thirty-one modules were recovered in D. quadriceps [mean number of genes per module is 12 (range: 12–3,715)], and 41 modules were recovered in P. canadensis [mean number of genes is 68 (range: 16–2,268) (Dataset S5). By using the log₁₀ P value (for connectivity measurement in the module, Dataset S5) as a significance score and the log₂ fold change to indicate caste-biased differential expression, we identified three (light yellow, magenta, and dark turquoise) and two (yellow, red) modules positively correlated in the ant and wasp, respectively. We then looked into the functionality of the modules, by checking for GO enrichment in the gene terms (as described SI Text, section III.4). The magenta module in the ant was the only one to be significantly functionally enriched, with multiple ribosomal terms. We were also able to use this analysis to pick out key hub (highly connected) genes for phenotypic differentiation.

III.6. Variance Analysis.

FPKM values (SI Text, section III.3) were log-transformed, and statistics were performed using the MatrixStat Package in R. The coefficient of variation (CV), variation, and SD values were calculated for each transcript (Fig. S2 A, C, and D). For Fig. S2D, nonlinear regression lines (RLs) for reproductive and nonreproductive states were calculated based on an exponential decay model curve. CV distances from RLs were plotted, and values below 0.05 of √ (reproductive^2 + nonreproductive^2) were excluded. Transcripts with a distance greater than 0.05 for a specific reproductive state and less than −0.05 for the other reproductive state were marked as specific and highly variable transcripts. Statistical tests (χ²) between variable and nonvariable groups were performed between each reproductive state.

IV.1. DNA Extraction and BS-Seq Library Preparation for Methylomes.

Genomic DNA was extracted from single brains using the QIAGEN All Prep DNA/RNA Mini Kit according to the manufacturer’s instructions. For BS-seq libraries, synthetic genomes of P. canadensis and D. quadriceps were prepared to serve as an unmethylated control genome in the analysis of global 5-methyl-cytosine (5mC). Synthetic genomes for each species were generated from an individual brain extraction using the REPLI-g Ultra-fast Mini Kit from QIAGEN. One microliter of these synthetic genomes was used to generate BS-seq libraries as outlined below. Between 200 ng and 500 ng of input genomic DNA was used per library. One library per species was additionally spiked with unmethylated lambda DNA at 1:1,000 (wt/wt) to provide an estimation of BS conversion efficiency (Bismark reports in Dataset S4). After sonication (Covaris), DNA was end-repaired, A-tailed, and ligated with a methylation Adaptor Oligo Kit (Illumina) using the NEB Next kit according to the manufacturer’s instructions. The adaptor-ligated DNA was treated with sodium-BS using an Imprint DNA Modification Kit from Sigma–Aldrich according to the manufacturer’s instructions for the two-step protocol. BS-treated DNA was amplified using KAPA HiFi Uracil + DNA Polymerase (KAPA Biosystems) with 15 cycles. Size selection and cleaning between steps were performed with an AMPure XP system (Agencourt) to select DNA fragments between 250 bp and 500 bp.

IV.2. BS-Seq Analysis and Quantification.

IV.3. Non-CG Methylation Validation.

To eliminate the possibility that non-CG methylation was false positively called due to SNPs in the genome assembly, we analyzed the downstream base of the reads from which the non-CG calls originated (Fig. S3A). The first 1 million reads of each library were used to count the methylated and nonmethylated cytosines that align to WT or variable sequences to estimate the levels of the non-CG methylation wrongly assigned due to genomic CpG SNPs. Both contexts of non-CGs (CHG and CHH where H signifies adenine, thymine, or cytosine) were analyzed.

IV.4. CpG Island Analysis.

CG methylation drops around the transcription start site (TSS) more drastically in A. mellifera compared with D. quadriceps and P. canadensis; both of our study species seem to maintain a permissive transcriptional environment at the TSS (i.e., with no methylation; Fig. 3A), a feature that often indicates genomic CpG islands (CGIs) in vertebrates (70). To detect CGIs in our species, we first applied a widely used algorithm developed for the human genome (71). We found only 256 CGIs in P. canadensis, whereas more than 26,000 were detected in D. quadriceps, comparable to the number of CGIs found in other ant species (6). Conventional algorithms for the detection of CGIs have been reported to underpredict CGIs in nonhuman species (72). Indeed, by plotting the CG density around the TSS in P. canadensis, we found significant enrichment of CGs around the TSS compared with the basal intergenic level, although P. canadensis has shorter stretches of CG enrichment compared with D. quadriceps, as expected from the CGI prediction numbers (Fig. S6A). Our findings show that such genomic features are conserved well beyond vertebrate genomes (72) and indicate a lack of maintenance of CGI boundaries and density specifically in insects.

IV.5. Phylogenetic Analyses of Methylation Machinery.

The NCBI NR database was searched for homologs of Dnmt1, Dnmt2, Dnmt3a/3b, Uhrf1, Tdg, Smug, Ung, and Tet enzymes using the BLAST algorithm. The collection was extended by acquiring the sequences of additional proteins identified in the Transcriptome Shotgun Assembly Database (tsa_nr) and high-throughput genomic sequences databases. Next, the sequence collection was sorted and verified for the completeness of protein sequences, proper assignment to the protein groups, and presence of characteristic protein domains and catalytic motifs. When possible, the protein sequences were also retrieved via BLAST searches from annotated protein databases [official gene sets (OGSs) are listed in Dataset S3] available for hymenopteran genomes. The phylogenetic tree was based on the NCBI taxonomy and generated with the interactive tree of life (Fig. S5A).

IV.6. Analysis of Global 5mC and 5-Hydroxymethyl-Cytosine Levels by Liquid Chromatography/MS.

Five hundred nanograms of genomic DNA was incubated with 5 units of DNA Degradase Plus (Zymo Research) at 37 °C for 3 h. The resulting mixture of 2′-deoxynucleosides was analyzed on a Triple Quad 6500 mass spectrometer (AB Sciex) fitted with an Infinity 1290 LC system (Agilent) and an Acquity UPLC HSS T3 column (Waters), using a gradient of water and acetonitrile with 0.1% formic acid. External calibration was performed using synthetic standards, and for accurate quantification, all samples and standards were spiked with isotopically labeled nucleosides as described by Bachman et al. (74) (Fig. S5C).

IV.7. Functional Enrichment Analyses of Methylated Genes (GO) and Ortholog Identification.

Functional enrichment of methylated genes (Dataset S4) was assessed by generating a NR list of genes from an initial transcript list and analyzing this list with Blast2GO, using a background of all genes. Groups with more than 10 genes with a significance level of P < 0.05 (Fisher test) were taken to be significant and used to generate functional enrichment tables. To find methylated orthologs between species, two BLAST databases were created using protein sequences. Each protein then underwent sequence searching by BLAST against the database from the other species. When a pair of proteins both matched the other’s database with ≥40% identity and ≥30 amino acids, they were classed as orthologs.

IV.8. Identification of Differentially Methylated Genes (DMGs).

To identify DMGs associated with phenotypes in P. canadensis and D. quadriceps, we used the 1-kb window probe quantification method (SI Text, section IV.2). For each probe, the methylation levels of the three reproductive replicates were compared with the three nonreproductive replicates using a Student’s t test (P < 0.05) and applying Benjamini and Hochberg multiple testing correction. For both species, the outcome of this test failed to detect any candidate genes (overlapping with at least one significant probe) having a phenotype-associated methylation signature.

IV.9. Alternative Splicing (AS) Analysis.

For the alternative splicing analysis, we used only protein-coding genes from the EVM-based annotation.

Reads from all RNA-seq samples were aligned to the reference assembly and EVM-based annotation using gemtools RNAseq pipeline (version 1.6.2), which uses Graphical Environment Manager (GEM; algorithms.cnag.cat/wiki/The_GEM_library) as a mapper. GEM is “split-aware” and finds gapped matches; that is, it can correctly align reads originating from exon junctions into the genome. Additional filtering steps were performed to allow fewer than two mismatches and five multimappings. On average, 74 ± 4% (D. quadriceps) and 82 ± 2% (P. canadensis) of the reads were mapped across samples, with 72 ± 4% (D. quadriceps) and 81 ± 2% (P. canadensis) of the reads mapping uniquely. Flux Capacitor (version 1.2.4, sammeth.net/confluence/display/FLUX/Home) was used to quantify genes, transcripts, exons, and splice junctions in each sample separately. Expression levels were given in pure read counts and in reads per kilobase per million (RPKM).

V.1. Library Construction and Analysis.

miRNAs are important regulators of phenotype differentiation in social insects (40, 76). miRNA libraries were constructed for pools of individuals of each phenotype for both species to determine whether large numbers of miRNAs are shared between hymenopterans to the exclusion of the other insects. Total RNA was extracted from multiple individuals, developmental stages, and castes for both taxa using a TRIzol (THERMO) extraction. The Illumina TruSeq Small RNA Sample Preparation Kit was used for the library construction, with 27.3 million (D. quadriceps) and 31.5 million (P. canadensis) reads sequenced. These reads were processed using Galaxy and analyzed using miRDeep2 in conjunction with a variety of customs scripts to reduce the number of misannotated miRNAs. Previously unidentified miRNA families were identified following established criteria, where miRNA loci were considered distinct families if there were mismatches in the seed region (nucleotides 2–8) (77, 78) so that only high-confidence miRNA families were annotated (i.e., the presence of both 5′ and 3′ products and clear homogeneity in the processing). In addition, BLAST searches of the genomes were conducted to identify any miRNAs that were not present in the small RNA libraries but have been annotated in other species.

This combined approach allowed the identification of 73 miRNA gene families, including four previously unidentified families in P. canadensis and 86 miRNA gene families, including 11 previously unidentified families in D. quadriceps (submitted to miRBase). The predictive nature of the miRNA complement (79) allows an assessment of the completeness of sequencing, which is high, with only 13 miRNA families missing (four in D. quadriceps and nine in P. canadensis), and equates to 8% missing data (genome screens for Acyrthosiphon indicate ∼24% missing data), although some of those data missing will be genuine secondary losses rather than incomplete sampling. These data allow us to reconstruct the evolutionary history of many eusocial insect-specific miRNAs; previous experimental data have been confined solely to Apis. Thus, we identify five miRNA families shared by all apocritans, and nine shared by aculeatans (Fig. S7).

V.2. miRNA Identification and Targeting.

miRNA can play a role in post-transcriptional gene regulation either by precisely binding to target sites in the 3′ UTR of the candidate mRNA, which is then cleaved, leading to degradation of the mRNA, which is common in plants, or by imprecise binding of the 5′ end of the miRNA onto the target site on the 3′ UTR of the mRNA, leading to down-regulation and repression of the mRNA, which is more common in animals. Thus, we investigated whether miRNAs are preferentially targeting DEGs between phenotypes in comparison to those genes that are nondifferentially expressed (Dataset S4).

TargetSpy was used to identify miRNA/target interactions, and it uses two input files: a list of miRNAs and a list of candidate 3′ UTRs. Custom scripts were used to look for the last CDS of a gene, and this process was then extended on the 3′ side until it reached 500 bp, until it reached the end of the genome sequence, or until another gene was hit. This 3′ UTR was then extracted for all 35 DEGs and 3,121 non-DEGs in D. quadriceps and for all 33 DEGs and 2,223 non-DEGs in P. canadensis. Against these 3′ UTR sequences, 195 miRNA products (known 5′ and 3′ sequences) and 141 miRNA products were investigated for D. quadriceps and P. canadensis, respectively. Because 3′ UTR length was variable (D. quadriceps: non-DEGs = 319 nt, DEGs = 372 nt; P. canadensis: non-DEGs = 313 nt, DEGs = 429 nt), and because the longer the 3′ UTR, the greater was the likelihood of target sites being identified, the data were normalized per 100 nucleotides.

Our results (Fig. S7) showed that there was no difference in the number of miRNA target sites between the non-DEGs (2.8 hits per 100 nucleotides, SD = 1.5) or the DEGs (3.0 hits per 100 nucleotides, SD = 1.4) in P. canadensis or the non-DEGs (3.4 hits per 100 nucleotides, SD = 2.1) and the DEGs (3.8 hits per 100 nucleotides, SD = 2.1) in D. quadriceps.