Conserved roles of Osiris genes in insect development, polymorphism and protection

Much of the variation among insects is derived from the different ways that chitin has been moulded to form rigid structures, both internal and external. In this study, we identify a highly conserved expression pattern in an insect‐only gene family, the Osiris genes, that is essential for development, but also plays a significant role in phenotypic plasticity and in immunity/toxicity responses. The majority of Osiris genes exist in a highly syntenic cluster, and the cluster itself appears to have arisen very early in the evolution of insects. We used developmental gene expression in the fruit fly, Drosophila melanogaster, the bumble bee, Bombus terrestris, the harvester ant, Pogonomyrmex barbatus, and the wood ant, Formica exsecta, to compare patterns of Osiris gene expression both during development and between alternate caste phenotypes in the polymorphic social insects. Developmental gene expression of Osiris genes is highly conserved across species and correlated with gene location and evolutionary history. The social insect castes are highly divergent in pupal Osiris gene expression. Sets of co‐expressed genes that include Osiris genes are enriched in gene ontology terms related to chitin/cuticle and peptidase activity. Osiris genes are essential for cuticle formation in both embryos and pupae, and genes co‐expressed with Osiris genes affect wing development. Additionally, Osiris genes and those co‐expressed seem to play a conserved role in insect toxicology defences and digestion. Given their role in development, plasticity, and protection, we propose that the Osiris genes play a central role in insect adaptive evolution.


Introduction
The chitin-based insect cuticle is a key ingredient in the spectacular ecological and evolutionary success of the insects. Chitin forms parts of both the external and internal anatomy of an insect, forming a hard skeleton that has been modified in nearly endless ways and facilitated anatomical, physiological and functional divergence in insect form (Emlen & Nijhout, 2000). On the outside, chitin is a barrier to the external environment, but also light and flexible enough to form the skeletal basis of locomotary appendages (legs and wings). On the inside, rigid invaginations of the exoskeleton (apodemes) are muscle attachment points that have allowed insects to evolve amazing strength and speed (Larabee & Suarez, 2014), but also form the respiratory tracheal system and the fore and hindgut (Terra, 1990). Chitin is also secreted in the midgut of insects as part of digestion and immunitythe peritrophic matrix/membrane is a chitinaceous structure that surrounds ingested material and is a barrier between the digestive tract and surrounding tissues (Lehane, 1997). The peritrophic matrix thus enables the encapsulation and isolation of toxins and infectious agents.
Where as a great deal is known about the structure and development of insect cuticle (Wigglesworth, 1957;Andersen, 1979;Hopkins & Kramer, 1992;Merzendorfer, 2011;Tomoyasu & Fujiwara, 2017), some studies using next-generation transcriptomic techniques, such as RNA sequencing (RNAseq), find genes of unknown function coincident with the timing of cuticle deposition (Ren et al., 2005;Sobala & Adler, 2016). Some of these genes are the most highly expressed genes during particular developmental stages, such as several Osiris family genes during stages of pupal development (Graveley et al., 2011;Gelbart & Emmert, 2013). Intriguingly, very little is known about many of these highly expressed genes, including those in the Osiris family.
The Osiris genes are a family of genes that have a high degree of synteny across all studied insects and appear to have radiated via gene duplications near the origin of the Insecta, approximately 400MYA (Shah et al., 2012;Misof et al., 2014). The majority of Osiris genes,~20, lie in an 160-kbp region (on chromosome 3R in D. melanogaster). The locus containing these genes has been known for several decades in D. melanogaster due to its embryonic lethality when present in haploid or triploid copy (the triplo-lethal locus, Tpl; Lindsley et al., 1972;Dorer et al., 1995Dorer et al., , 2003. Although this result clearly suggests dosage-dependent lethality, it is unclear how or why the lethality occurs. Previous work established that embryos with an abnormal copy number of tpl fail to form midgut and die as late embryos or upon hatching from the egg (Smoyer et al., 2003). The function of only one Osiris gene is known, Osi21, and this particular gene is not in the main cluster. Osi21 in D. melanogaster, also known as diehard4, inhibits the recycling of rhodopsin in the retina as part of the endocytic pathway (Lee et al., 2013).
Despite a longstanding knowledge of their existence, their precise role remains unclear. These genes are primarily characterized by a protein domain of unknown function (DUF1676) and have no gene ontology (GO) terms associated with their molecular function. In addition to DUF1676, Osiris genes contain a putative signal peptide, a hydrophobic region that is likely a transmembrane region, and an AQXLAY motif at their 3 0 end (Shah et al., 2012). This structure is compatible with the known role of Osi21 in endosomal trafficking (Lee et al., 2013).
We used comparative transcriptomics to shed light on this enigmatic family of genes. We predicted that the conservation of synteny in the Osiris cluster has a functional basis through gene regulation and that genes with physical proximity would share similar expression patterns across species. Due to their high degree of sequence conservation, we also hypothesized that the timing of expression of these genes would be consistent even across clades that diverged hundreds of millions of years ago and that their expression would coincide with suites of other genes, forming coregulated expression modules. By leveraging species producing multiple divergent phenotypes, the social insect castes (queen, worker, and male), we were also able to test for differences in Osiris gene expression responsiveness to the developmental environment (Anderson et al., 2008). In the context of studies on insect chemical protection, our data suggest that Osiris genes may be involved in the regulation of insect cuticle, both in the generation of external morphology as well as in forming a defensive barrier.

Ortholog detection
The D. melanogaster set of Osiris genes is the most curated set and was used as a basis for detection of orthologs and paralogs. The methods used in this research was similar to that of Shah et al. (2012), but for consistency was repeated for all genomes included in this paper. A combination of approaches was used to detect Osiris genes. First, reciprocal best hit BLAST (rbhBLAST), with an e-value threshold of e < 0.0001, was performed between each genome and that of D. melanogaster (using the 24 D. melanogaster protein models as queries). To widen this search to paralogs, all identified Osiris genes from rbhBLAST were used as BLAST queries into their own genome, and hits with e < 1E-5 were then used as reciprocal queries into D. melanogaster (kept if their best hit was an Osiris gene). As an independent, and complementary, approach, hidden Markov models were used to search proteins. A model was built, function hmmbuild, using all 24 D. melanogaster Osiris proteins and used as a template to search all other genomes using the function hmmsearch in HMMER3.1b2 (Wheeler & Eddy, 2013). Evidence from all of these approaches was combined to determine likely Osiris genes.
For a simplified tree, we used a subset of the above taxa, but still representative of major insect clades, only 1:1:1 orthologs, and the same parameters as above. To help resolve internal nodes, we forced monophyly of Osiris ortholog groups based on the results of the larger tree. The resulting tree is Fig. 1.

DNA and RNA source and processing
RNA-sequencing data from D. melanogaster were from Graveley et al. (2011, Table S9) for 12 embryonic stages, six larval stages, six pupal stages and six adult stages. Data from P. barbatus were from Smith et al. (2015) (NCBI project PRJDB3493), and consisted of 16 libraries spread across four developmental stages. These were as follows: early last (4th-instar) instar larvae (prior to differences in mass), late last instar larvae (queens and workers are distinguishable by mass), pupae and adults. In each developmental stage, there were two worker and two queen libraries. The raw RNA-sequencing reads were mapped to the newer NCBI annotation (Thibaud-Nissen et al., 2013) of the P. barbatus genome, Pbar_-UMD_V03, at NCBI (annotation release 100), using tophat2 (Kim et al., 2013), and bowtie2 for genome indexing (Langmead & Salzberg, 2012). Feature counting was performed using HTseq (Anders et al., 2015). B. terrestris RNA-sequencing data were from Harrison et al. (2015;NCBI project PRJEB9366), and remapped as above to the assembly Bter_1.0 and annotation release 102 performed at NCBI (Thibaud-Nissen et al., 2013). The B. terrestris libraries consisted of larval, pupal and adult stages spread over queen, male and worker castes. We only used nonreproductive adult workers when comparing adults. Each stage by caste library set had three biological replicates, except for queens, where there was only a single library per stage. See the supplement for sample code used in RNA mapping.
RNA-sequencing data from Formica exsecta were from Morandin et al.
(GenBank Biosample SAMN02046301-SAMN02046306) and consisted of 12 libraries spread across three developmental stages (pupae, emerging adult and old adult). In each developmental stage, there were two biological replicates of queen and worker. The raw RNA-sequencing reads were mapped to the Formica exsecta genome (NCBI project) using tophat2 (Kim et al., 2013) and Cufflinks (Trapnell et al., 2010); RSEM was used for transcript quantification (Li & Dewey, 2011). Gene expression data, as RPKM, are included in Data S1.

Gene expression correlation over space and time
Osiris gene expression data in each species, including all stages and castes, was used to calculate expression similarity. A correlation matrix of Osiris gene expression was performed in R using the cor function in the base package (R Core Team 2017). Visualization of correlation was performed with the corrplot package (Wei et al., 2016), using AOE clustering. Gene expression over space and evolutionary history was examined only in D. melanogaster. To examine correlation over space, the distance between pairwise gene starting positions was calculated and was correlated with correlation coefficients of Osiris developmental gene expression. Similarly, pairwise branch length distances between all genes in D. melanogaster were calculated using the cophenetic.phylo function from the ape package (Paradis et al., 2004) using a tree with all other species omitted. These pairwise distances were then correlated with Osiris developmental gene expression correlation coefficients. Similarity in Osiris gene expression among species was calculated for each species using all librariesas above, (dis)similarity distance matrices were calculated using Bray-Curtis distance. Pairwise Mantel tests compared the similarity of Osiris gene expression between species pairs.

Analysis of expression and co-expression
Only features with at least one count per million (CPM) in at least four or three libraries (i.e. the number of developmental stages in each data set), for P. barbatus and B. terrestris, respectively, were kept for downstream analysis. Downstream analysis and graphing were performed using the R statistical environment (R Core Team 2017), with graphics produced in ggplot2 (Wickham, 2016). Differential expression analysis between castes, within developmental stages, was performed using edgeR , using SVA (Leek & Storey, 2007;Leek et al., 2017) first to remove possible unwanted variation; caste was used as a factor in the model using the function svseq. Differences among libraries were adjusted using TMM normalization (Robinson & Oshlack, 2010), and comparisons between groups were performed using a linear model using the sveq result as a covariate and caste 9 developmental stage as a factor.
Co-expression modules were calculated using WGCNA on CPM expression count (Langfelder & Horvath, 2008) with one-step network construction and module detection. All modules were correlated against developmental stages (in the case of D. melanogaster) or developmental stage 9 caste (for the social insects). The input data set was filtered using WGCNA cleaning step that removes genes with too many missing values or zero variance. Due to the interest of this study, we mapped Osiris genes to modules and looked most closely at modules loaded with Osiris genes, which were also positively correlated with pupae of the worker caste. Gene overlap lists and Venn diagrams were calculated with gplots (Warnes et al., 2009) and VennDiagram (Chen & Boutros, 2011).
WGCNA recommend the use of filtered normalized expression count or FPKM, RPKM, without it altering the final results. To verify this assumption, the above WGCNA analysis was repeated with RPKM values as input data set, and similar results were obtained. Genemodule associations are available in Data S2, and for Osiris genes (with their Osiris ortholog names) in Data S3.

Gene ontology (GO) analysis
InterProScan5 (Jones et al., 2014) was used to search all proteins in each annotation for matching protein domains Pfam database, (Finn et al., 2016) and GO terms. The package topGO (Alexa et al., 2006;Alexa & Rahnenfuhrer, 2016) was used to detect enrichment of terms in co-expression modules using the weight01 algorithm to take into account GO term relatedness.

Quantitative real-time PCR
Colonies of German cockroaches, Blatella germanica, were housed in plastic containers and fed ad-lib commercial cat food and water. Individuals were anaesthetized with wet ice, and heads were measured using a calibrated stereo-microscope camera to assign instars (Tanaka & Hasegawa, 1979). After head measurements, the samples were immediately deep-frozen. RNA was extracted by homogenizing whole bodies in TRIzol (Life Technologies) using a pestel, and then, we followed the manufacturer's protocol for the Direct-zol TM RNA Miniprep kit (Zymo Research). Purified RNA (0.25-1 lg) was reverse-transcribed using the QuantiTect RT kit (Qiagen). Quantitative real-time PCR (qrtPCR) was performed using the QuantiTect SYBR Green PCR Kit (Qiagen) with a total volume of 10 lL with 10 lM each forward and reverse primer and 1 ng cDNA. Primers were designed using NCBI primer designing tool (Ye et al., 2012) for four Osiris genes (Osi7, two paralogs of Osi9a and b, and Osi20) as well as two control genes, GAPDH and Actin. Protein sequences were obtained from Baylor College of Medicine Human Genome Sequencing Center (NCBI project PRJNA203136). The primers used were ( curves were run to confirm the appropriate product. All primer pairs had efficiency~0.95, and it was verified that GAPDH and Actin transcript abundance did not vary predictably with development. All reactions were performed in triplicate. The difference in cycle threshold (DC t ) between each gene and each control gene was calculated and averaged. Results were graphed in ggplot2 (Wickham, 2016).

Results and discussion
Ortholog identification and phylogeny We used methods similar to Shah et al. (2012) to find putative orthologs and establish orthology among Osiris genes across divergent insect species. We used a combination of protein BLAST (Altschul et al., 1990) and hidden Markov model searches with HMMER 3.1.2b (Eddy & Wheeler, 2015) using all 24 annotated D. melanogaster Osiris genes. Once putative Osiris genes were found, we performed BLAST searches for additional paralogs within each genome. To establish orthology relationships, we used reciprocal best hit BLAST (rbhBLAST) between putative Osiris genes and D. melanogaster. Finally, once lists were compiled for each species, we constructed phylogenies of both all potential orthologs and 1:1:1 orthologs across insect taxa using BEAST2 (Drummond et al., 2012;Bouckaert et al., 2014); to help resolve internal nodes, we forced monophyly of groups with posterior probability support > 0.9. The phylogeny of representative orthologs across insects (Fig. S1) supports monophyly of nearly all Osiris ortholog groups found in D. melanogaster, in agreement with the results of Shah et al. (2012). Figure 1 shows the relationships of a subset of 1:1:1 orthologs across the major insect clades. The major radiation of Osiris genes dates at least to the origin of insects,~400 MYA (Misof et al., 2014). Many internal nodes in the phylogeny lack robust support. There appear to be two major ancient lineages of Osiris genes, the red cluster and the blue-green cluster of Fig. 1. As the noninsect Osiris orthologs group with Osi 17/24, this group is the putative ancestral Osiris gene. The green cluster (Osi 18-20) resulted from an internal duplication within the blue cluster. The duplications in the blue-green group likely gave rise to the genes at the beginning and end of the current cluster (Figs 1 and S2). The red group is duplications that filled in the middle of the cluster (Fig. 1).

Spatial/phylogenetic regulation of osiris genes
To explore the dynamics of Osiris gene expression, we used three primary gene expression data sets, the D. melanogaster modENCODE developmental transcriptome that spans 12 embryonic, six larval, six pupal and 3 1 ( 2 0 1 8 ) 5 1 6 -5 2 9 six adult developmental stages (Graveley et al., 2011), a data set from the bumble bee, Bombus terrestris, with larvae, pupae and adults of each males, workers and queens (Harrison et al., 2015) and a data set from the red harvester ant, Pogonomyrmex barbatus, with two larval, pupal and adult stages of both workers and queens (Smith et al., 2015). A fourth data set on the wood ant, Formica exsecta (Morandin et al., 2015), was included for comparative purposes, but only includes pupae and adult stages of queens and workers. Downstream analysis was conducted on normalized gene expression after filtering out genes with low transcript abundance; all analyses were performed using the R statistical framework version 3.3.3 (R Core Team 2017).
Transcript abundances of Osiris genes tend to be positively correlated within each species, though there are clear 'blocks' (Figs 2-3, Figs S3-S6, Data S1-S3) where expression among genes of a block is very highly correlated (r = 0.8-1). Among pairwise correlations of Osiris genes in D. melanogaster, 82% were positive and 51% were statistically significant (P < 0.05, Fig. S6). To assess similarity of expression among species, a similarity matrix of expression was made for each species using the vegan package in R (Oksanen et al., 2013). Genes had similar patterns of developmental expression across species. The gene expression matrices were all highly correlated (Mantel test, P < 0.05 in all comparisons, Fig. S7). Despite these species diverging > 300 MYA (Misof et al., 2014), the Osiris genes have maintained similar within-cluster patterns of gene expression, suggesting that they are coregulated. Precise temporal regulation of these genes is likely necessary given their putative function in the regulation of chitin because chitin deposition is periodic during the development of insects.
We calculated both pairwise linear chromosomal distances and pairwise phylogenetic distances for D. melanogaster to examine potential mechanisms of regulation. Osiris gene expression similarity (pairwise gene correlations) was correlated with physical gene distance (P < 0.0001, Fig. 3), and phylogenetic distance (P = 0.004, Figs S8-S11), suggesting either co-expression due to locality or related ancestral regulatory mechanisms. Co-regulation of neighbouring genes is a pervasive pattern in genomes (Boutanaev et al., 2002;Spellman & Rubin, 2002;Michalak, 2008). Intriguingly, there are Hymenoptera-specific ultraconserved elements (UCEs) in the 5 0 UTR of six Osiris genes (Osi7,8,9,16,20;Davies et al., 2015). These elements do not appear to have a conserved RNA secondary structure, but may play a role in translational regulation. Because we detect both a pattern of phylogenetic conservation of expression, as well as spatial co-expression, it is unclear whether the observed co-regulation of genes in this cluster is due to neighbourhood effects, inherited sequence patterns or a combination of the two. A previous analysis (Morandin et al., 2016) examining signatures of selection in ants found that all of Osi15 the ten Osiris genes examined have a ratio of nonsynonymous to synonymous nucleotide changes less than one, suggesting that they have likely evolved (at least in the ants) under purifying selection. In summary, Osiris genes are highly conserved at the coding and noncoding sequence level, and conservation extends to both synteny and developmental gene expressionover hundreds of millions of years.

Temporal expression patterns
There are three clearly distinct times of Osiris upregulation in D. melanogaster development, 14-to 18-h embryos, second-instar larvae and 48-h-old pupae (Fig. 4). During these times, some Osiris genes are among the most expressed genes in the transcriptome.
Osi6 is the sixth highest in 14-to 16-h embryos and Osi3, 7, 9 and 12 are first, third, ninth and second during 48-h-old pupae, at well over 1009 mean gene expression in those developmental stages. Similarly, Osiris gene expression peaks in pupal development in both B. terrestris and P. barbatus. There is a global pattern across data sets, ubiquitously higher expression in workers (and males) compared to queens -Osiris gene expression in worker pupae was > 1000-fold higher than in queens for some genes (e.g. Osi3, 5, 6). Due to low power to detect significant differences in expression between castes, only a few of the Osiris genes (Osi 3, 6, and 9) were statistically different after correction for false discovery rate (FDR < 0.05, Table S1). The pattern of greater expression in workers relative to queens suggests a possible role of these genes in generating, or regulating, caste differences. Differences in dispersal and reproduction are among the most defining differences among social insect castes. In B. terrestris and both ant species, there are differences in the reproductive capacities of queens and workers, as well as body size differences. There are also wing polymorphisms in all three species. In B. terrestris, in which all castes have functional wings, worker and male wings are more similar to each other than either are to queen wings in both shape and size (Medler, 1962;Gerard et al., 2015). Worker castes in both ant species lack wings. Interestingly, Osiris genes are among the most highly expressed genes in the early stages of D. melanogaster wing development (Sobala & Adler, 2016), when their expression coincides with the timing of cuticle envelope formation. Whereas temporal peaks of most Osiris genes are relatively consistent, Osi18-20 (the green group, Fig. 1) are temporally different relative to the others, peaking later in both embryonic and pupal development in D. melanogaster (Fig. 4, Figs S8-S11), and peaking in larvae of male and worker B. terrestris (Fig. 4). This difference in developmental timing, coupled with their physically tight linkage, suggests differential regulation. Interestingly, despite both minor and major inversions in the Osiris cluster of some insects, Osi18-20 remained physically proximate. The largest rearrangement of the Osiris cluster can be seen in the flour beetle, Tribolium castaneum, where the entire middle of the cluster was inverted and translocated to the 5 0 of Osi1the tail end of this inversion, Osi14-16, was further translocated over 2Mbp away (Fig. S2).

Co-expression modules
To gain insight into the function of the Osiris genes, we examined which other genes have correlated expression using weighted gene co-expression analysis using the WGCNA package in R (Langfelder & Horvath, 2008). This analysis was performed for all four species, and in each, we correlated module eigengenes with pooled developmental stages, or caste 9 development for the social insects. There were five co-expression modules that were positively correlated with the pupal stage of D. melanogaster, two of these had Osiris genes and one had the majority of the cluster (11 Osiris genes of 220 total genes) and included those with the highest levels of expression (Fig. S12). There was only a single coexpression module significantly correlated with worker pupae in both B. terrestris (Fig. S13) and P. barbatus (Fig. S14). Osiris genes were abundant in each of these worker-pupa-correlated modules; there were 11 Osiris genes of 88 total genes in the B. terrestris worker-pupa module, and seven Osiris genes of 182 total, with another seven Osiris genes (of 243 total) in a second module weakly correlated with worker pupae in P. barbatus. No modules were significantly correlated with worker pupae in F. exsecta (Fig. S15), likely due to limited power with few distinct developmental stages, and fewer samples. Nonetheless, F. exsecta Osiris genes were divided into two modules (seven and 11 genes in each) with both modules having a weak association with worker pupae.
We annotated gene ontology (GO) terms for all genes in the above modules, concentrating on those modules enriched with Osiris genes. The 'structural constituent of cuticle' gene ontology term was enriched in both ants and D. melanogaster, but was not statistically significant in B. terrestris (P = 0.17; Fig. 5). Gene ontology terms that were enriched in at least two species were 'chitin binding' and 'serine endopeptidase activity.' The Osiris genes do not have annotated molecular function GO terms, and thus did not bias this analysis. Overall, the enriched GO terms strongly suggest that these modules play a role in chitin and cuticle formation. The term for serine endopeptidase activity is also consistent with a role in chitin processes as these enzymes are known to cleave chitin synthases, which are secreted as zymogens, and increase the rate of chitin synthesis (Merzendorfer & Zimoch, 2003;Broehan et al., 2007Broehan et al., , 2010. Among the primary co-expression modules for all species, there was a four-species overlap of three genes, Osi8-9 and CG6055 (Fig. 5). Three-species overlaps included Osi3, 6, 7, as well as five genes of unknown function (CG7031, CG12964, CG14237, CG14636 and CG32816), and multiple genes implicated in cuticle/ chitin formation and assembly (Fig. 4). The function of dyl has been established as the deposition of chitin, both around bristles (Nagaraj & Adler, 2012) and in the expansion of wings (Ren et al., 2005). The relevance of St2 sulfotransferase is not clear, although sulfotransferase activity is required for proper epidermis formation in nematodes (Kim et al., 2005). Genes CG6055 and CG14866 encode C-type lectin-like proteins, which in other species are known to bind chitin as part of antifungal immune responses (Bueter et al., 2013). The genes dusky (dy), Mind-the-Gap (mtg), Gasp and obstructor A (obst-A) all bind chitin, with the latter two known to effect epidermal cuticle integrity (Tiklov a et al., 2013).
Interestingly, a butterfly with a seasonal wing colour polymorphism upregulates several Osiris genes under long-day conditions (Vilcinskas & Vogel, 2016). Upregulated along with the Osiris genes is dyl, but upregulated under the alternative short-day condition are multiple cuticle proteins, C-type lectins, a carboxypeptidase and Gasp. This result (and see Osiris gene expression in toxicology and immunity, below) suggests that genes within these modules can be decoupled for different purposes, though acting on the same phenotypes (in this case, wings).

Osiris gene expression in the hemimetabola
The only currently available nonholometabolous insect gene expression developmental series is in the termite, Zootermopsis nevadensis (Terrapon et al., 2014). Osiris gene expression is high in the egg sample, but is otherwise very low. To gain a better understanding of the evolution of gene expression in Osiris genes, we conducted quantitative real-time PCR (qrtPCR) on the hemimetabolous cockroach, Blatella germanica. We designed primers to amplify Osi7, Osi20 and two copies of Osi9 (named a and b), along with control genes Actin and GAPDH. We sampled 56 individuals representing all nymphal instars from a laboratory colony. Osiris gene expression was undetectable, or nearly so, in all but the last larval instar, along with some Osi9a expression in one adult sample (Fig. 6). When Osiris gene expression peaked, it had very high abundance relative to GAPDH and Actin, ubiquitously expressed genes. These data suggest a conserved timing of expression for Osiris genes in the penultimate instarthe last nymphal instar for the hemimetabola, and the pupa for the holometabola. In both insect groups, the penultimate stage coincides with  Fig. 1. The x-axis scale starts at the location of Osi1, and scales are different for each species; the identity and total length of the genome region pictured is on the left.
the final transition to adulthoodthe maturation/acquisition of traits associated with dispersal and reproduction.

Osiris gene knock-down
Despite their unknown function, Osiris genes are known to be essential for proper development. Abnormal copy number of the tpl locus (encompassing the Osiris cluster) in D. melanogaster results in embryonic or post-hatching lethality. The iBeetle knock-down screen (Schmitt-Engel et al., 2015) has embryonic and pupal phenotypes for ten Osiris genes. Most of these phenotypes are lethal in late-larval injections, with lethality up to 100% for Osi7, 17 and 24. The larval knock-down of Osi7 fails to metamorphose and arrests as a prepupa, while others have increased lethality in pupae or early adults. Embryonic knock-down phenotypes indicate major developmental errors, from segment loss to appendage malformation, and of course death (Schmitt-Engel et al., 2015).  Osiris gene expression in toxicology and immunity Given the chitinous nature of the peritrophic matrix and its role as an infection barrier by isolating food and toxins, and the possible role of Osiris genes in chitin regulation, it is not surprising that the Osiris genes are upregulated in response to immune challenge. When fed supernatant from cultured Bacillus thuringiensis (Bt), the flour beetle, T. castaneum, upregulated seven Osiris genes (Greenwood et al., 2017). Interestingly, these same genes were down-regulated when challenged with Bt (regardless of whether they were immune primed). Also of interest from this study is a decoupling of molecular function, as noted from GO term enrichment, in Bt primed and Bt challenged individuals. The structural constituent of cuticle GO term is enriched with Bt priming, and the serine endopeptidase activity term is under-represented with priming, although these terms reverse upon Bt challenge. This is in contrast to our study where these terms were enriched in pupae/ worker pupae, again suggesting a possible decoupling of developmental gene co-expression modules in immunity. Several studies in honeybees, Apis mellifera, have made associations with Osiris genes and immunity.
When challenged with chalkbrood fungus, Ascosphaera apis, Osi6 was down-regulated (Aronstein et al., 2010). On the other hand, a study infecting A. mellifera larvae with foulbrood bacteria, Paenibacillus, led to upregulation of 17 Osiris genes (Cornman et al., 2013). The authors of this last study suggested that this effect could be confounded by development, as other genes like dusky (dy) and dusky-like (dyl) were also upregulated. However, in the light of the frequency with which Osiris genes are associated with infection protection, at different life stages, it is likely that the observed pattern is biologically relevant. Interestingly, this effect of Osiris genes may not be localized to the midgut, as much of their expression is in other tissues, especially in the pupal fat body (Gelbart & Emmert, 2013). A proteomic study of the sexually transmitted honeybee fungal pathogen, Nosema apis, found Osi7 protein upregulated in drone sperm infected with the pathogen (Grassl et al., 2017). Also among the upregulated proteins were chitinases, which act to break down chitin in fungi, suggesting a logical connection to the upregulation of Osi7 in this study.
The Osiris genes are also upregulated upon challenge with plant secondary chemical defences. In response to feeding on plants producing glucosinolates, the leaf-rolling fly, Scaptomyza flava, upregulated eight Osiris genes along with various genes involved in chitin binding and that are structural constituents of cuticle (Whiteman et al., 2012). Both D. sechellia and some populations of D. yakuba have independently evolved host specificity for the noni fruit, Morinda citrifolia, which produces a potent insecticide, octanoic acid. Both species show evidence of strong selection on the Osiris cluster (Yassin et al., 2016). Moreover, several Osiris knock-downs (Osi6-8) lose resistance to octanoic acid (Andrade L opez et al., 2017). Interestingly, the opposite effect was found in adult tissue-specific knock-downs of these same genes, particularly in fat body and salivary gland. Additionally, D. melanogaster upregulated several Osiris genes upon chemical challenges, including at several life stages. For example, Osi6 has high to very high expression when exposed to heavy metals, caffeine and the Sindbis virus (Gelbart & Emmert, 2013).

Conclusions
With this study, we aim to address a major knowledge gap regarding the importance of Osiris genes for insect development, ecology and evolution. Osiris genes have seemingly contrary function. They are upregulated in social insect castes that have no wings (or smaller wings), yet are also upregulated in wing development of nonsocial species. They are upregulated and downregulated in immune responses, and similarly in the presence of a toxin. However, the essential nature of these genes in insect development is apparent, with lethality at multiple developmental time-points in their absence/knock-down. Furthermore, their role in toxin resistance is also becoming clear, as is their role in the generation of alternative phenotypes. Given the breadth of their action, these genes could have very important roles in insect adaptation and success, from host-parasite and herbivore-plant coevolution to morphological novelty and phenotypic plasticity. The expansion of these genes at the root of the insect tree is further suggestive that they may have played an important role in wing evolution and/or the radiation of early insects.

Acknowledgments
The authors are very grateful to J. Gadau, the Bornberg-Bauer and Kurtz laboratory groups, and the entire Institute of Evolutionary Biology at U. M€ unster, for being very generous with their ideas and time. H. Lerner advised on phylogenetic methods, thank you. Joe Coolon provided helpful discussions and input on the manuscript. Funding for this project was provided as a fellowship by the Evolution Think Tank (U. M€ unster) to CRS and grant support from the Indiana Academy of

Supporting information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1 Phylogeny of insect Osiris orthologs. Figure S2 Synteny of Osiris genes within the Holometabola. Figure S3 Correlation matrix of developmental Osiris expression, B. terrestris. Figure S4 Correlation matrix of developmental Osiris expression, P. barbatus. Figure S5 Correlation matrix of developmental Osiris expression, F. exsecta. Figure S6 Plot of correlation statistical significance by strength for all four species. Figure S7 Comparative NMDS plots of Osiris developmental expression. Figure S8 Guide phylogeny of D. melanogaster Osiris genes. Figure S9 Phylo-heat-map of expression: 16 h embryo vs. 48 h pupa. Figure S10 Phylo-heat-map of expression: 18 h embryo vs. 2nd instar larva. Figure S11 Phylo-heat-map of expression: 48 h vs. 72 h pupa. Figure S12 Co-expression module development correlation: D. melanogaster. Figure S13 Co-expression module development correlation: B. terrestris. Figure S14 Co-expression module development correlation: P. barbatus. Figure S15 Co-expression module development correlation: F. exsecta. Table S1 Between caste differences in Osiris gene expression. Data S1 Relative gene expression tables (RPKM, reads per kilobase of transcript per million mapped reads) for species (except fruit fly) used in gene expression analyses. Data S2 Mapping of genes to co-expression modules for all species used in co-expression analyses. Data S3 Osiris gene name and co-expression module grouping for all species used in co-expresion analyses.