Gene expression changes during female reproductive development in a colour polymorphic insect

Pleiotropy (multiple phenotypic effects of single genes) and epistasis (gene interaction) play key roles in the development of complex phenotypes and might be especially important in polymorphic taxa. The development of discrete and heritable sympatric phenotypic polymorphisms often emerges from major-effect genes that interact with other loci and have pleiotropic effects on multiple traits and physiological functions. We quantified gene expression changes during ontogenetic colour development in an insect (damselfly: Ischnura elegans) with three heritable female colour morphs, one which is a male mimic. Using transcriptome sequencing and de novo assembly, we demonstrate that all morphs show extensive downregulation of gene expression during early colour development. Morphs then become increasingly differentiated during sexual maturation and when developing adult colouration. These different ontogenetic trajectories arise because of the heterochronic development of the male-mimicking females, compared to other two female morphs. Many loci with regulatory functions in reproductive development are uniquely regulated in male-mimicking females, including upstream and downstream regulators of ecdysone signalling and transcription factors that influence insect sexual differentiation. Our study reveals extensive epistasis in the genomic architecture of I. elegans and suggest a central role for pleiotropy in shaping the developmental trajectories of the different morphs.


Introduction
The role of epistasis (gene interaction) in shaping phenotypic variation and covariation has been discussed since the early days of the Modern Synthesis (Fisher 1930;Wright 1930;1931), and still remains a contentious subject Shao et al. 2008;Huang et al. 2012;Hemani et al. 2014;Mackay 2014;Wood et al. 2014). One way to investigate how gene interactions affect the development of phenotypic variants is to quantify the pleiotropic effects in gene expression profiles that arise as a result of allelic variance at single locus (e.g. Brem et al. 2005;Litvin et al. 2009). The number and function of differentially regulated genes between allelic variants can reveal epistatic interactions between the focal locus and other loci. Colour-polymorphic taxa are particularly suited for investigating such questions, as colour morphs typically differ in multiple physiological, developmental and life-history traits (Mckinnon and Pierotti 2010). The genetic basis of such heritable colour morphs is often quite simple and arises due to variation at one or a few loci with restricted recombination in between them (Joron et al. 2006; Thomas et al. 2008;Lamichhaney et al. 2016;Lindtke et al. 2017;Andrade et al. 2019). Pleiotropic effects of such major effect loci on the regulation of multiple other genes underlying morph differentiation and epistasis is likely to shape the development of such complex and co-adapted phenotypes.
There are two main mechanisms by which pervasive phenotypic differences can evolve in colour polymorphic systems with a relatively simple genetic architecture. Multiple loci involved in distinct pathways may be recruited into a chromosomal rearrangement, such as a chromosomal inversion with reduced recombination, effectively becoming physically linked as a 'supergene' (Kirkpatrick and Barton 2006;Schwander et al. 2014; Thompson and Jiggins 2014). Alternatively, genetic differentiation between morphs may arise from regulatory linkage (sensu West-Eberhard 2003) of a suite of co-selected traits underpinned by a polygenic architecture, whereby modularity and pleiotropy result from a few regulatory loci (Sinervo and Svensson 2002;Nijhout 2003;ten Tusscher and Hogeweg 2009). In this latter scenario, a colour-morph locus can evolve new 40 45 50 55 regulatory functions at the top of a gene-regulatory network, in addition to already present developmental effects on multiple tissues (Monteiro and Podlaha 2009). Physiological epistasisi.e. the dependency of a gene function on previously or simultaneously expressed gene products at other loci (Moore and Williams 2005) -would then underlie the pleiotropic effects of one or a few tightly linked colour-morph loci (Sinervo and Calsbeek 2003;ten Tusscher and Hogeweg 2009).
Epistatic interactions arising from such a limited number of regulatory loci will emerge during the course of development, as phenotypic differences between adult morphs increase as they approach sexual maturation.
A substantial fraction of phenotypic evolution in the history of Metazoans has resulted from regulatory changes in the time, location and volume of interactions among gene products (Gerhart and Kirschner 2007). Many evolutionary novelties have arisen through a shift in the relative timing or rate of gene expression (e.g. Gomez et al. 2008;Albertson et al. 2010;Gunter et al. 2014), rather than due to the total disruption of previous interactions or the formation of entirely novel networks (Rebeiz et al. 2011). This phenomenon, known as heterochrony was historically considered to be a source of novel phenotypic variation between independently evolving lineages (Gould 1977;Reilly et al. 1997). However, variation in the timing of developmental events can also be important at the intraspecific level (Linksvayer and Wade 2005). A small but increasing body of evidence suggests that genetic control over regulatory networks can account for substantial phenotypic variation in several colour-polymorphic taxa (Kunte et al. 2014;Timmermans et al. 2014;Yassin et al. 2016;Takahashi et al. 2018). Heterochrony can produce large phenotypic changes but with only small genetic change (West-Eberhard 2003). These findings raise questions of how complex suites of physiological traits associated with colour polymorphisms emerge, and how variation in the timing and rate of developmental processes shape such profound intraspecific differences. By investigating how and when correlated trait differences between colour morphs arise during development, we can 70 75 80 develop a better understanding of how epistasis and pleiotropy build up and shape different morphs within species.
Here, we investigated how changes in gene expression during the course of adult development differ among heritable female morphs in the trimorphic common bluetail damselfly (Ischnura elegans). This species, and several other species within the genus Ischnura have female-limited colour polymorphism that are maintained by frequency-dependent sexual conflict, in which common morphs suffer from reduced fitness due to excessive male mating harassment Gosden and Svensson 2009;Takahashi et al. 2010;Le Rouzic et al. 2015, Gering 2017.
Female morphs differ in their ontogenetic trajectories of colour change and in suites of other phenotypic traits, including fecundity Willink et al. 2019), mating rates , parasite resistance and tolerance (Willink and Svensson 2017), cold tolerance (Lancaster et al. 2017) and both larval  and adult development times (Svensson et al. submitted). Thus, the colour locus (or a set of tightly linked loci) might be involved in epistatic interactions resulting in pleiotropic effects on the expression and timing of many other developmental processes.
We asked four major questions about the developmental regulation of gene expression in colourdeveloping females of three female morphs in I. elegans. First, at which point(s) in colour development do female morphs become different in their gene expression profiles? Second, which physiological processes characterise colour development across the female morphs of I. elegans?
Third, which regulatory genes underlie these patterns? finally, are any physiological processes heterochronic between morphs with respect to colour development, and if so which ones? The answers to these questions will increase our understanding of how pleiotropy and epistasis shape the evolution of sex-limited polymorphisms.

Study system
Females of the Common Bluetail damselfly (Ischnura elegans) occur in three heritable colour morphs whereas males are monomorphic (Fig. 1a). Previous studies using breeding experiments in controlled laboratory environments across multiple generations have revealed that colour morph development is governed by a single autosomal locus (or a set of tightly linked loci) with sexlimited expression to females (Cordero 1990;Sánchez-Guillén et al. 2005). These female morphs also undergo unique ontogenetic colouration changes, which are correlated with the onset of reproductive capacity and are not apparent in males (Cordero et al. 1998;Svensson et al. 2009;Willink et al. 2019). Female morphs differ in both their sexually mature (final) colour patterns as well as in the qualitative colour changes that are expressed during sexual development (Fig. 1). In I. elegans, these female morphs are typically referred to as 'Androchrome', 'Infuscans' and 'Infuscans-obsoleta'. To be consistent with previous publications (Willink and Svensson 2017), we hereafter refer to these three heritable female morphs, including both their immature and mature developmental phases, as A-, I-and O-females.
All three female morphs in I. elegans undergo pronounced ontogenetic changes in their colour pattern (Fig. 1b). Shortly after emergence from their last nymphal moult, A-and I-females exhibit a thoracic violet background colour with black antehumeral stripes. In A-females, the violet colour becomes turquoise-blue over post-emergence development, and therefore the final thoracic colour pattern of sexually mature A-females strikingly resembles the colour of males (Fig. 1). This adaptive colour change in A-females is consistent with male mimicry, for which there is experimental and observational evidence in species of Ischnura (Robertson 1985;Gosden and Svensson 2009;Gering 2017). In contrast, the violet thorax colour of I-females becomes greenishbrown in association with the onset of reproductive potential (Willink et al. 2019;Fig. 1b). Finally, O-females are also sexually dimorphic when sexually mature, but they do not exhibit distinct 120 125 130 135 antehumeral stripes at any point in development and their thoracic colour changes from salmonpink to brown (Cordero et al. 1998;Willink et al. 2019;Fig. 1b).
The two sexually dimorphic I-and O-females exhibit dramatic ontogenetic changes in the phenotypic expression of a blue patch on the eighth abdominal segment (Fig. 1b). This blue abdominal patch is initially present in sexually immature individuals of all three female morphs ( Fig. 1b). However, in I-and O-females the blue patch subsequently becomes concealed by dark pigment over sexual development, whereas in males and the male-mimicking A-females, the blue colouration is instead retained during the entire course of sexual development (Henze et al. 2019;Fig. 1). All three female morphs therefore undergo ontogenetic colour changes over sexual development, but there are pronounced differences between A-females and the two heterochrome female morphs (I-and O-females). These colour changes occur continuously over the span of a few days (Svensson et al. submitted,Fig. 1b and S1). During this transitional period, individual females with intermediate phenotypes between the immature and mature developmental phases can be visually distinguished (Fig. 1b).

Fieldwork and sampling design
We sampled three females of each heritable colour morph at three stages of colour development, using the detailed description and classification of these colour stages by Cordero et al. (1998) and Svensson et al. (2009). The immature developmental phase was defined in all morphs by the presence of a contrasting blue colour patch in the distal portion of the abdomen, and for A-females by the simultaneous expression of the violet thoracic colouration. Females were classified to the mature developmental phase if they expressed the final colour phenotype of sexually mature individuals of their respective morph. Finally, transitional females expressed an intermediate phenotype between these two extremes. For I-and O-females these transitional individuals were identified by partial concealment of the blue patch with pigment, which results in a brownish-blue 145 150 155 160 abdominal patch (Fig. 1b). Transitional A-females were characterised by a violet-turquoise thoracic colouration, which precedes the development of the final turquoise-blue background colour (Fig.   1b). All individuals in this study were collected in the field in southern Sweden during the summer of 2015 and from populations that form part of a long-term longitudinal study, in an area of approximately 40 x 40 Km 2 (TableS1; Le Rouzic et al. 2015;Willink and Svensson 2017).

RNA sample preparation and sequencing
In addition to the 27 females sampled for the differential gene expression (DGE) analysis, an additional 58 mature-coloured females were sampled and used for the de novo transcriptome assembly below (Table S1). Field-caught individuals were bisected at the second abdominal segment and the two sections were preserved in RNAlater immediately following bisection and stored at -80 C within 1-4 h of capture. Samples were then shipped to the Beijing Genomics Institute (1/F, 16th Dai Fu Street, Tai Po Industrial Estate, Tai Po, Hong Kong) for extraction and sequencing. Total RNA was extracted using the Trizol LS Reagent (Life Technologies) following the manufacturer's instructions. RNAseq libraries were prepared for individual samples using the Illumina TruSeq kit with mRNA enrichment. Each sample was sequenced using 100 bp paired-end reads at a depth of 2 GB per sample on an Illumina HiSeq 2000. Raw data will be deposited at the National Center for Biotechnology Information (NBCI). Sample information and reads will be available through Biosample links and sequence read archives.

Transcriptome assembly and annotation
Raw reads were trimmed using Trimmomatic v.036 (Bolger et al. 2014) to remove adapter sequences and ambiguous nucleotides. We cropped the first 10 nucleotides of each read and trimmed regions of four nucleotides that had an average quality score below 20. Reads shorter than 25 nucleotides after trimming were discarded. All trimmed reads were pooled to assemble a 170 175 180 185 reference transcriptome using Trinity v. 2.2.0 (Grabherr et al. 2011). Assembly quality was assessed by generating Ex50 estimates using scripts from the Trinity package (Haas et al. 2013), and by quantifying the ortholog hit ratio (OHR) (Vera et al. 2008;O'Neil et al. 2010;O'Neil and Emrich 2013). Ex50 estimates are generated by plotting the N50 of contigs for different expression quantiles, searching for the maximal N50 for regions capturing 70 -90% of the expression data.
OHR estimates were generated using the predicted protein set from the genome of the banded demoiselle damselfly, Calopteryx splendens (Csple_OGS_v1.0.faa; Ioannidis et al. 2017), then using blast (tblastn; Camacho et al. 2009) to search for the best hit of each protein in the I. elegans transcriptome assembly. OHR was then estimated by dividing the aligned length of an I. elegans transcript contig by the length of its orthologous C. splendens protein, wherein ratios of 1 indicate the full coverage of the contig for a given protein. The result was a single OHR estimate for each non-redundant protein from a related damselfly, the ortholog of which was found in the I. elegans transcriptome. The OHR distribution provides a detailed assessment of a transcriptome, indicating how well 1000's of genes have been assembled. Prior to analysis, redundancy due to isoforms and duplications in the C. splendens protein set was controlled by collapsing it by the longest protein isoform of each protein cluster group. This was based upon >90% amino acid identity using CD-HIT v4.5.4 (Li and Godzik 2006) with parameters used in making the UniRef90 database (Suzek et al. 2014). Annotation of the reference transcriptome was conducted using blastx v. 2.6.0+ with an evalue of 10e-5 against the National Center for Biotechnology Information (NCBI) non-redundant database (nr). Blast hits were then mapped to Gene Ontology (GO) terms using Blast2GO v. 4.1.9 (Conesa et al. 2005). If the blasted sequences returned more than a single hit, only the highest scoring hit was used in subsequent gene ontology (GO) enrichment analysis.

Data Pre-processing
Prior to the DGE analysis we used the CORSET algorithm (Davidson and Oshlack 2014) to reduce the number of spurious isoforms, which often arise during de novo transcriptome assemblies of non-195 200 205 210 model organisms due to sequencing artefacts. CORSET hierarchically clusters contigs based on shared reads and expression levels. In this case, clusters were based on the 27 samples to be used in our DGE analysis. The read counts are then summarised by these clusters, providing gene-level transcript abundances. Lowly expressed genes (i.e. with fewer than one read per million reads in at least three samples) were discarded at this stage and gene expression distributions were normalised via the weighted trimmed mean of M-values (TMM) method (Robinson and Oshlack 2010).

Differential gene expression analysis
DGE analyses were conducted using the packages edgeR McCarthy et al. 2012) and limma (Ritchie et al. 2015) in R v. 3.4.4 (R Core Team 2018, and were based on the workflow describe by (Law et al. 2016). In order to apply linear-based statistical modelling to the expression data, we estimated the mean-variance relationship in the log transcript counts per million (log-cpm), and incorporated this mean-variance trend into a precision weight for each normalised observation, using the voom method (Law et al. 2014;Liu et al. 2015). The advantage of this statistical approach, which estimates the mean-variance relationship instead of specifying a generating probabilistic distribution, is that it more appropriately controls for type I error rates when sample sizes are small and when sequencing depth varies between samples. In the limma pipeline these precision weights are then taken into account when fitting the linear models to the log-cpm gene expression data (Law et al. 2016). Empirical Bayes moderation was subsequently applied to these fitted models to more accurately estimate expression variability among genes and the 'robust' procedure was employed to allow variance outliers (Smyth 2004;Law et al. 2014;Phipson et al. 2016). Multiple testing across statistical contrasts was conducted using the "global" method, which is recommended when the number of differentially expressed (DE) genes is interpreted as a measure of the strength of a contrast (Smyth et al. 2018). The global method performs multiple testing adjustment across the entire set of contrasts and genes being tested. The P-value cut-off was 220 225 230 235 adjusted according to Benjamini and Hochberg's false discovery rate (Benjamini and Hochberg 1995).
Our goal was to compare DGE across developmental stages and between the different colour morphs. Thus, we fitted a fully-factorial linear model including the fixed effects of female colour morphs, developmental phases and their interactions (~0 + Morph + Phase + Morph:Phase). With this model structure, the intercept is removed from the first factor (Morph), but kept in the second factor (Phase). Therefore, within each morph the immature phase is the intercept, and expression differences between groups can be compared by specifying a contrast matrix. Three types of comparisons are possible with the interaction term: (1) contrasts among morphs within each phase, (2) contrasts between subsequent phases within morphs and (3) and contrasts between developmental phases among colour morphs (full contrast matrix in Table S2). Because we sampled three females of each of three morphs and at each of three phases in colour development (27 females in total), developmental comparisons in (2) and (3) span one of two windows of colour development: early colour development, between the immature and transitional phases, and late colour development, between the transitional and mature phases. Our main focus in this study was on (3), because these contrasts represent the aspects of development that are regulated in opposite directions or with different magnitude among morphs. We used (2) to address the question of which genes are differentially expressed during the course of colour development and which are shared, either among morphs or between different developmental phases. We also used contrasts in (2)  A-females ( Fig. 1b and S1). Given our classification of developmental colour phases, this differential rate of colour change occurs only during early colour development, whereas late colour development is of a similar duration for all female morphs ( Fig. 1b and S1). In spite of these differences in the time to transition among developmental phases, the early and late developmental windows represent a similar fraction of the full extent of colour changes that occur over development in each morph.

Correlations among subsets of differentially expressed genes
We used vector correlations to determine if developmental expression changes shared similar patterns across morphs for each colour transition (early and late colour development). Vector correlations have been more commonly used in studies of multivariate phenotypic evolution (Zelditch et al. 2012), but they are also useful to assess the similarity between estimated contrasts of gene expression across a set of multiple genes (Zinna et al. 2018). The correlation between two mean-centred vectors, each containing the estimated magnitude of expression differences between two experimental categories, is: (1) where a⋅b is the dot product between the two vectors of statistical contrasts and ‖a ‖ and ‖b ‖ represent the magnitudes of vectors a and b respectively, and are given by the square root of the sum of squares of all contrast estimates. The gene vectors a and b in this study correspond to logtransformed coefficients estimated from the linear model and contrast matrix described above (Table S2)  We followed the approach of Zinna et al. (2018) to determine if vector correlations between contrasts that have been considered significantly different from zero were particularly strong, compared to correlations between vectors of a random sample of genes, irrespective of the degree of statistical significance. This would indicate that those genes identified as significantly regulated across a period of colour development are also more similar (or dissimilar) in the direction of their regulation than a random gene vector of equal length. We defined in each case an experimental vector, for example, a vector containing all differentially expressed genes during early colour development that are shared between two female morphs. We then obtained an empirical distribution of vector correlations, drawn from 10 000 randomly sampled vectors of the same length as the experimental vector. We considered a correlation for the experimental vector as extreme when it fell outside the 95 percentile of correlations in equal-length random vectors.

Gene ontology enrichment analysis
To uncover the physiological changes that follow colour transitions in the different female morphs, we tested for significantly enriched GO terms in groups of developmentally regulated genes in comparison to the reference transcriptome, using the bioconductor package TopGO (Alexa and Rahnenfuhrer 2016) in R. For these analyses, we focused on contrasts between developmental colour phases within morphs (contrast group (2); Table S2). In the cases where multiple transcripts associated with a single gene (i.e. a Corset cluster) differed by at least one GO term (2 231 out of 30 134 mapped transcripts), all GO terms mapped to the different transcripts were pooled, and every unique term was assigned to the gene. Given that our interest was to explore how physiological processes change over development in the three female morphs, here we focused on biological process ontologies. The significance of enriched terms was tested using Fisher's exact test and the 'elim' algorithm, in order to account for the hierarchical structure of gene ontologies (Alexa et al. 2006). With this algorithm, more general annotations are excluded for genes already counted as 300 305 310 315 significantly enriched for a more specific GO term. GO terms with P-values < 0.01 were considered significantly enriched in the gene set of interest, compared to the reference transcriptome.

Colour development, heterochrony and female fecundity
Our recent findings of faster colour development in A-females (Svensson et al. submitted,Fig. 1b) and a higher number of differentially expressed genes during late colour development in this morph only (see Results) indicates that A-females reach their final developmental colour phase earlier than I-and O-females. This also suggests that this colour change might be decoupled from other developmental processes, including reproductive maturation. To test the idea that colour development could be faster (i.e. heterochronic) than most other developmental changes in Afemales, we first compared early and late expression changes within each female morph, using the vector correlation approach described above. A-females reach their mature colour at about the same absolute age as I-and O-females develop their transitional phase (Svensson et al. submitted,Fig. 1b). If most other developmental changes proceed at similar rates in all female morphs, we would expect expression contrasts in early and late colour development to be more correlated in A-females than in the other two morphs, as in A-females these two developmental periods span a shorter absolute time. In contrast, if the fast colour changes of A-females indicates overall faster development, we would expect early and late expression contrasts to be as different in A-females as in I-and O-females, despite the shorter duration of early colour development in A-females.
Furthermore, if A-females accelerate their colour maturation in relationship to reproductive development, we expect that field-caught A-females in their final developmental colour phase would more often exhibit reproductive failure than I-and O-females, because at any given time a fraction of the mature-coloured A-females would still be sexually immature. Our fecundity data to test this prediction come from our surveys of 18 natural populations in Southern Sweden, which were visited frequently during the mating season (June to August) for 2-17 years, between 2001 and 2017 Le Rouzic et al. 2015;Willink and Svensson 2017). During these visits, mating couples were collected in the field. Mated females were placed in individual plastic containers with moist filter paper for oviposition. After 72 h eggs were scanned at 1 200 dpi and subsequently counted using ImageJ (Schindelin et al. 2012).
In these surveys, mated females were classified by morph and developmental phase. However, transitional females were not classified as such, and before 2015 we did not distinguish between Aand I-females in the immature developmental phase (Fig. 1b). Therefore, we modelled the probability of reproductive failure only for the immature and mature developmental phases, and all data on immature A-and I-females come from the last three years of our population survey (2015)(2016)(2017). We used a mixed-effect model fitted by MCMC in the R package MCMCglmm (Hadfield 2010) to test for a difference between morphs in the probability of reproductive failure, which we defined here as mated females that produced no viable eggs. We specified a binomial model with morph, developmental phase and their interaction as fixed effects. We also allowed for a random interaction between the population and season terms to affect the variance around the fixed effect estimates. All code to reproduce these analyses as well as the reproductive failure data will be uploaded to Dryad.

Reference transcriptome
Using an extensive collection of RNA-Seq data, we assembled a transcriptome of 889,001 transcripts. Grouping transcripts by their expression levels and looking at those with the highest levels of expression in relation to their assembly lengths, a peak of transcript assembly length (N50) of nearly 3kb was found at expression levels containing 70-80% of transcripts (i.e. ExN50 ~ 3kb), suggesting sufficient data for a high quality assembly. We also estimated the OHR of our assembly by querying our assembly with the full protein set of the only published odonate genome, the 350 355 360 365 banded demoiselle damselfly (Calopteryx splendens), which shared a common ancestor with I. elegans approximately 85 million years ago (the average from two dated phylogenies, 92 Ma in Thomas et al. 2013 and 79 million years in Waller and Svensson 2017). Of the 22,507 C. splendens proteins that were non-redundant, we found good hits (e-value less than 10e-5) for roughly 19,962, and of these, we found OHR of 0.7 and higher for nearly 70% (13,497; see Fig. S2 and Table S3).
Given the extensive divergence between the two taxa, the number of putative orthologs between species recovered is high and the strong skew towards high OHR values is indicative of a high quality transcriptome assembly.

General trends in early colour development
In all female morphs, early colour development, between the immature and transitional phases was characterised by extensive downregulation of gene expression. This was evidenced by comparisons between consecutive developmental phases and within morphs (contrast group (2); Table 1; Table   S2). Across-the-board downregulation exceeded upregulation by a factor of more than 4, in Ofemales, to more than 25, in A-females (Table 1). Of all downregulated genes during early colour development, 3 158 were shared by the three female morphs, whereas there were no upregulated genes shared by all morphs. The downregulated genes shared by all female morphs corresponded to ~ 62%, 75% and 67% of all downregulated genes in A-, I-and O-females, respectively. These genes were significantly enriched for biological processes involved in protein synthesis, transport and catabolism, energy metabolism, signal transduction, and cathecolamine biosynthesis ( Table 2, S4). Broadly shared downregulation of gene expression also resulted in highly correlated vectors of significantly regulated genes between all morph pairs (Fig. 2a).

General trends in late colour development
No single gene was either downregulated or upregulated in all female morphs during late colour development. This absence of shared regulatory patterns was caused by the distinctiveness of A-375 380 385 390 395 females. There was a perfect negative correlation between A-and I-females in the regulation of genes that were differentially expressed in both morphs (Fig. 2a). Of 269 such genes, 164 were upregulated in I-females and downregulated in A-females and 105 exhibited the opposite regulatory pattern. There was also a negative correlation in the direction of regulation for the genes differentially expressed in A-and O-females, but this correlation was not particularly extreme (P = 0.14), as most genes, whether significant or not, exhibited opposite developmental contrasts during late colour development. In contrast, there were 66 genes with significant developmental regulation in both I and O-females and all of these genes were regulated in the same direction in both morphs (Fig. 2a). The biological processes enriched in the DE genes during late colour development were also markedly different between A-females and both I-and O-females. A-females upregulated genes enriched for functions in cell cycle progression, meiosis and the regulation of transcription, while downregulating multiple processes associated with energy metabolism (Table S6, (Table S6, S7).

Colour development heterochrony and female fecundity
Unlike I-and O-females, which underwent a qualitative change in the direction of gene regulation at the onset of late colour development, regulation of gene expression was positively, although not extremely, correlated between early and late colour development in A-females (Fig. 2b). Such correlation between developmental periods in A-females was not extreme when compared to a random sample of expression contrasts drawn from equal-length gene vectors regardless of their significance level (P = 0.163; Fig. 2b). Nonetheless, downregulated genes during both early and late colour transitions in A-females were enriched for two shared biological processes, both associated with ATP biosynthesis (Table S4,  Across these three female morphs, immature developmental phases exhibited a nearly four-fold probability of reproductive failure compared to their mature-coloured counterparts (PMCMC < 0.001; Fig. 3). There were no differences in the rate of reproductive failure among the three female morphs in the immature phase (all PMCMC > 0.05; Fig. 3a). Among mature-colour individuals, reproductive failure was somewhat higher in A-females compared to I-females (PMCMC = 0.019; Fig. 3b), yet this difference was small relative to differences between mature and immaturecoloured females (Fig 3).

Regulatory differentiation during colour development
There were more substantial morph differences in developmental regulation of gene expression during the second phase of colour development (contrast group (3) Table S2; Fig. 4). During this period, A-females became increasingly differentiated from both I-and O-females as revealed by the regulatory pattern of 856 genes (Fig. 4). Statistical contrasts for all these genes were of the same direction between A-and I-females and between A-and O-females.
To examine how epistatic effects of the colour-morph locus could mediate the differentiation of Afemales, we queried the two subsets of genes which were differentially regulated between A-and both I-and O-females (contrast group (3); Table S2) for regulatory functions. First, we used the search term 'juvenile hormone' (JH), and 'ecdysone' (E), two major hormones involved in regulating insect reproduction and development (Simonet et al. 2004;Flatt et al. 2005). While we found no differentially regulated genes directly involved in JH metabolism, four genes related to ecdysone metabolism and signalling were upregulated in A-females compared I-and O-females during late colour development. To discern the potential functions of these genes, we used blastx 430 435 440 445 against annotated proteins in the Drosophila melanogaster genome in FlyBase (Thurmond et al. 201). Three of these ecdysone-related genes were annotated to D. melanogaster genes in the Halloween group, which code for Cytochrome P450 enzymes that participate in ecdysone biosynthesis (Table 3, Fig. 5). The fourth gene was annotated as coding for the ecdysone-induced protein E74, a transcription factor that responds to concentrations of 20E in the ovary of D. melanogaster (Table 4; Fig. 6; Ables and Drummond-Barbosa 2010).
We further queried the gene sets differentiating A-females from both I-and O-females for regulatory functions, by identifying all genes in such subsets annotated with the GO term "regulation of transcription". This query revealed 18 genes, all differentially regulated during late colour development (Table 4). Most of these genes (~89%) increased in expression in A-females, while remaining relatively constant or decreasing in I-and O-females. We again used blastx against the D. melanogaster genome to examine the potential functions of these regulatory genes. We focus on the 13 genes that had good hits (e-value < 10e-5) against annotated proteins in SwissProt/TrEMBL (UniProt Consortium 2019). These genes comprised several transcription factors, including E74, a ribonucleoprotein and two proteins involved in histone acetylation (Table   4, Fig. 6).

Discussion
In this study, we investigated the regulatory changes that follow colour development in adult females of the trimorphic damselfly I. elegans. During the early phase of colour development, thousands of genes, underlying multiple metabolic processes, were significantly downregulated in all three female morphs (Table 1-2, S4). In striking contrast, upregulated genes were fewer and more idiosyncratic in function between morphs (Table 1, S5). This widespread downregulation of developmental gene expression resulted in strong correlations in the expression contrasts between the three different morphs (Fig. 2a). Soon after emergence from the final nymphal moult, gene expression is largely downregulated, potentially causing a slowdown in basic metabolism. The same pattern is present in holometabolous insects where more genes are downregulated than upregulated after pupation (e.g. Wang et al. 2010;Zheng et al. 2012), particularly in females (Graveley et al. 2011). Once energetically expensive tissues such as flight muscle have developed, resting metabolic rates also decline with age in adult insects (Hack 1997;Piiroinen et al. 2010;Niitepõld and Hanski 2013). Such metabolic cutback may contribute to delaying senescence (Sohal and Weindruch 1996).
Regulatory patterns during early colour development were also consistent with heterochronic colour development in A-females compared to the other two female morphs. Our recent work has shown that colour development is accelerated in the male-mimicking A-females compared to the two other female morphs (Svensson et al. submitted , Fig. S1). However, it it is still unclear whether the accelerated colour changes indicate that adult A-females develop overall faster or if they instead decouple their colour changes from other developmental processes. Disentangling between these two different alternatives is important to understand if ontogenetic colour changes in females act as signals of reproductive suitability, and to assess if the reliability of such signals varies among the different female morphs, which have alternative reproductive strategies (Willink et al. 2019).
We found that the expression contrasts during early and late colour development were negatively correlated in only I-and O-females (Fig. 2b). Some of the biological processes downregulated in all female morphs during early colour development, were also downregulated during late colour development in A-females (Table S4, S6). Moreover, mature A-females were more likely to exhibit reproductive failure than mature I-females, albeit this difference is relatively small (Fig. 3). A possible explanation for this morph difference is that A-females accelerate their reproductive development relative to other developmental changes, but such heterochrony comes at a fecundity cost, which sometimes results in reproductive failure. Across Ischnura species, developmental colour changes in heterochrome females (such as the I-and O-morphs) correlate with the onset of 480 485 490 495 reproductive capacity and are used by males to identify suitable partners (Fincke 1987;Langenbach 1993;Takahashi and Watanabe 2011;Willink et al. 2019). Our data show that the accelerated ontogenetic colour changes in A-females of I.elegans might also affect their reproductive potential.
In fact, A-females do also have lower average fecundity than both I-and O-females in the field Willink and Svensson 2017).
In adult females of I. elegans, most regulatory differences between morphs arise late in colour development (Fig. 4), when the colour differences between them are most apparent (Fig. 1b). The colour-morph locus in I. elegans likely interacts with numerous other genes via hormonal regulation. We found that three I. elegans genes with inferred roles in the ecdysone biosynthesis pathway were among the uniquely upregulated genes of A-females during late colour development.
One of these genes was inferred to code for 20EMO, the enzyme that catalyses the conversion of E to 20E in the ovarian follicle and nurse cells of D. melanogaster (Table 3 The other two genes code for enzymes associated with upstream points in the ecdysone biosynthesis pathway (Table 3; Fig. 5b,c). Finally, an ecdysone responsive gene essential for oogenesis also showed upregulation during late colour development in only A-females (Table 4; Fig. 6a). In D.
melanogaster, E74 responds to ecdysone titres to promote germline proliferation during early oogenesis, and is key for preventing germline degeneration during later stages of oogenesis (Ables and Drummond-Barbosa 2010).
The fact that both ecdysone biosynthesis and ecdysone-responsive genes are uniquely upregulated in A-females (Fig. 5, 6a) suggests that the colour morph locus might have pleiotropic effects via ecdysone signalling, regulating the timing of female reproductive development. Accelerated reproductive development in A-females during late colour development might partially compensate for their shorter period of colour transition, explaining why mature-coloured A-females are reproductively mature, just like I-and O-females, but have, on average, lower fecundity than these two other morphs Abbott 2005, Willink andWillink et al. 2019).
This interpretation is also supported by the expression pattern of several other regulatory genes presumably involved in reproductive development (Table 4) (Fig. 6b-f). For instance in D. melanogaster, a heterogeneous ribonucleoprotein (hnRNP) binds to poly(ADP-ribose) and this complex induces the expression of cadherin, which is in turn essential to maintain cell-cell adhesion during germline stem cell proliferation (Ji and Tulin 2012). hnRNPs can also be responsive to ecdysone levels through association with the protein on ecdysone puffs (PEP) (Amero et al. 1993). In I. elegans, both hnRNP and cadherin genes increase in expression during the late colour development of adult Afemales, compared to both I-and O-females (Fig. 6b, S4a). Consistently, a gene coding for poly(ADP-ribose) glycohydrolase, the enzyme that degrades poly(ADP-ribose), was instead downregulated during late colour development in A-females compared to O-females (Fig. S4b), whereas PEP was upregulated (Fig. S4c).
In insects, changes in hormone titres or their intracellular response machinery can occur at critical periods causing qualitative developmental responses, such as those underlying the development of environmentally induced polyphenisms (Nijhout and Wheeler 1982;Nijhout 1999). Our results suggest a similar mechanism in the development in some of the pervasive phenotypic differences that distinguish A-females from the two other female morphs (heterochrome females). Instead of an environmental cue, the colour morph locus in I. elegans, might act at the top of a ecdysonemediated hierarchy and gene cascade that governs the differential development of the three discrete adult female morphs.
The inferred functions of the regulatory genes uniquely expressed in A-females compared to I-and O-females are also consistent with the biological processes that characterised late colour development of A-females (Table 4, (Table S7). These results suggest that during adult development in I. elegans the colour-morph locus controls a regulatory cascade with pervasive effects on cell proliferation. These morph differences in the regulation of cell proliferation might occur primarily in the ovaries, as suggested by two of the biological processes enriched in the upregulated genes of A-females ("meoitic cell cycle" and "piRNA metabolic process"). While meiosis is a necessary step to produce gametes, piRNA protects the fidelity of genome transmission, by interacting with PIWI proteins to silence transpons in the germline (Malone et al. 2009). Accordingly, three out of four putative PIWI proteins in the I. elegans transcriptome were also upregulated in A-females during the late colour transition (Fig. S5).
The regulatory genes with unique expression patterns in A-females during late colour development also revealed a potential role of the colour-morph locus in inducing morph differentiation in somatic tissues of I. elegans. One of these genes mapped to Dmrt (Table 4) which is a gene family involved in sex determination and differentiation across animal taxa (e.g. Shen and Hodgkin 1988;Kato et al. 2011;Matson et al. 2011). This gene was upregulated in A-females compared to both I-and O-females. In somatic tissues of D. melanogaster the Dmrt gene, doublesex (dsx) integrates sexual identity, patterning and hormonal signals to induce the development of sexually dimorphic structures (Kopp 2012). Alternative splicing of dsx also underpins the development of male-mimicry in the female polymorphic butterfly Papilio polytes (Kunte et al. 2014). In the I. elegans congener I. senegalensis, heterochrome females are distinguished from both A-females and males in the expression pattern of a dsx isoform (Takahashi et al. 2018). dsx is therefore also a plausible regulator of the development of male mimicry in I. elegans, as it is differentially expressed between adult males and females (Chauhan et al. unpublished manuscript). However, we found no evidence of differential expression of dsx in I. elegans females, during the period of developmental colour changes (Fig. S6a) et al. 2005), but as all other Dmrt genes it contains a DM DNA binding motif, which confers regulatory capacity over pseudopalindromic DNA targets (Zhu et al. 2000).
Holometabolous insects have multiple Dmrt genes (Wexler et al. 2014), and I. elegans has at least two (Fig. 6h, S6a). However, the functions of the Dmrt genes other than dsx are poorly known in insects, although these genes might also integrate sexual identity, spatial and temporal signals during sexual differentiation (Picard et al. 2015). While sex determination mechanisms are diverse across animal taxa, their downstream actions which result in sexual differentiation are typically more conserved and often rely on transcriptional regulation by the DM domain DNA binding motif (Zhu et al. 2000). In line with the possibility that this Dmrt gene may drive morph-specific development in I. elegans, we found that both a homeotic gene (Table 4, Fig. 6g) and two known targets of dsx (Fig. S6b-e) in D. melanogaster were also differentially expressed between malemimicking A-females and both I-and O-females during late colour development.
In A-females, the expression of a homeotic gene annotated to Sex combs reduced (Scr) decreased as the expression of the Dmrt gene increased (Fig. 6g-h). Such epistatic interactions between Dmrt and homeotic genes have been described in D. melanogaster where both dsx and Scr are required for the proper localization of sexually dimorphic sex combs in male legs (Kopp 2011). Interestingly, we also found that three genes with similarity to bric-á-brac (bab1/2), a direct target of dsx in D.
melanogaster (Christiansen et al. 2002) were also upregulated during late colour development in Afemales ( Fig. S6b-d). In D. melanogaster, bab1/2 integrates dsx and homeotic inputs during the development of sexually dimorphic pigmentation (Kopp et al. 2000). Finally, the signalling peptide FMRFamide, another direct target of dsx in D. melanogaster (Luo et al. 2011), was also differentially regulated between A-and heterochrome females during late colour development. In D. melanogaster, FMRFamides modulate locomotory escape behaviour (Klose et al. 2010). In I. elegans, expression of the FMRFamide declines as the expression of the Dmrt gene increases (Fig.  S6e). Together, these results suggests that conserved regulatory hierarchies downstream of sex determination might also partly regulate the developmental differentiation between the malemimicking A-females and the two heterochrome I-and O-female morphs.

Conclusion
Here, we have investigated developmental regulation of gene expression during ontogenetic colour changes among heritable female colour morphs in I. elegans. Many regulatory changes occur shortly after the last larval moult (early colour development) and involve extensive downregulation of core metabolic processes in all female morphs. However, key regulatory differences between morphs emerge at a later ontogenetic stage (late colour development) when females develop their final colour patterns. These morph differences reflect heterochrony in colour development as malemimicking females develop their final colouration faster, albeit at a fecundity cost. The regulatory genes that differentiate male mimics during late colour development suggested a unique temporal pattern of germ cell proliferation, presumably driven by ecdysone signalling that influence morphspecific differences in reproduction. Male-mimicking A-females also differ from the two other morphs in their expression of a Dmrt gene and its potential downstream targets. This suggests that the mechanisms behind somatic sexual differentiation and male mimicry might have a common underlying genomic and developmental basis. The colour-morph locus in I. elegans is therefore involved in many epistatic relationships with numerous other genes, presumably via hormonal pleiotropy. Our results provide insights on how pronounced phenotypic differences between these female morphs emerge during the course of development and reveal how epistasis and pleiotropy play key roles in the generation of complex and co-adapted phenotypes within this colour polymorphic species.      their immature phase until sexual maturity (see Fig. S1). In A-females, the thorax colour changes during sexual development from violet to dull violet-turquoise, and finally to bright turquoise-blue.
In I-and O-females, both thorax and distal abdomen undergo noticeable colour changes presumably involving extensive pigmentation. While the thoracic colour patterns of these two morphs are different, they share developmental darkening as the blue patch in the distal abdomen portion becomes completely concealed at sexual maturity. A-and I-females express very similar immature phases. However, these two morphs can be distinguished by the black "arrow-like" pattern in the eighth abdominal segment, which is thicker and marked by a central notch in I-females compared to A-females.    Table   S2). (a) early colour development, (b) late colour development. The number of genes with no evidence of morph differences in differential expression between developmental colour phases is shown in the bottom right corner of each plot.  melanogaster. Expression data (log counts per million) have been weighed by within-sample variability and by the mean-variance relationship among genes, using 'voom' (see Methods).
Corresponding gene products and inferred modes of regulation and functions are shown in Table 4.   Table S2. Contrast matrices for the three types of comparisons of differential gene expression in females of Ischnura elegans, given the interaction term in the model formula: ~0 + Morph + Colour Phase + Morph: Colour Phase. There are three morphs (Androchrome, Infuscans, and Obsoleta) and three colour development phases: immature (Imm), transitional (Tra) and mature (Mat). The immature colour phase is treated as the intercept within each morph (i.e. 'Androchrome' refers to immature Androchrome females). Three types of contrasts were estimated: (1) contrasts among morphs within each colour phase, (2) contrasts between subsequent colour phases within morphs and (3) and contrasts between developmental changes among morphs. For (3), developmental changes corresponded to one of two developmental windows: early colour development (E), between the immature and transitional phase and late colour development (L), between the transitional and mature colour phase.  splendess genome that is covered by a single (longest) I. elegans transcript. Around 70% of the sequences are covered to 70% or more of their full expected length, indicating the transcriptome quality is acceptable.    Figure S4. Expression pattern of genes co-regulated and potentially in epistasis (see Discussion) with a heterogeneous ribonucleoprotein during adult development of I. elegans females. Transcripts were annotated against the non-redundant database of the National Center for Biotechnology Information (NCBI).
Expression data (log counts per million) have been weighed by within-sample variability and by the meanvariance relationship among genes, using 'voom' (see Methods). Expression data (log counts per million) have been weighed by within-sample variability and by the meanvariance relationship among genes, using 'voom' (see Methods). Abbreviations: PARG = poly(ADP-ribose) glycohydrolase, PEP = protein on ecdysone puffs. Figure S6. Expression patterns of genes in the I. elegans transcriptome with significant blast hits against doublesex (dsx) and two of its direct targets, bric-á-brac 1 and 2 (bab1/2) and FMRFamide, in D.
melanogaster, in the non-redundant database of the National Center for Biotechnology Information (NCBI).
Expression data (log counts per million) have been weighed by within-sample variability and by the meanvariance relationship among genes, using 'voom' (see Methods).