- Research article
- Open Access
Wound healing, calcium signaling, and other novel pathways are associated with the formation of butterfly eyespots
BMC Genomicsvolume 18, Article number: 788 (2017)
One hypothesis surrounding the origin of novel traits is that they originate from the co-option of pre-existing genes or larger gene regulatory networks into novel developmental contexts. Insights into a trait’s evolutionary origins can, thus, be gained via identification of the genes underlying trait development, and exploring whether those genes also function in other developmental contexts. Here we investigate the set of genes associated with the development of eyespot color patterns, a trait that originated once within the Nymphalid family of butterflies. Although several genes associated with eyespot development have been identified, the eyespot gene regulatory network remains largely unknown.
In this study, next-generation sequencing and transcriptome analyses were used to identify a large set of genes associated with eyespot development of Bicyclus anynana butterflies, at 3-6 h after pupation, prior to the differentiation of the color rings. Eyespot-associated genes were identified by comparing the transcriptomes of homologous micro-dissected wing tissues that either develop or do not develop eyespots in wild-type and a mutant line of butterflies, Spotty, with extra eyespots. Overall, 186 genes were significantly up and down-regulated in wing tissues that develop eyespots compared to wing tissues that do not. Many of the differentially expressed genes have yet to be annotated. New signaling pathways, including the Toll, Fibroblast Growth Factor (FGF), extracellular signal–regulated kinase (ERK) and/or Jun N-terminal kinase (JNK) signaling pathways are associated for the first time with eyespot development. In addition, several genes involved in wound healing and calcium signaling were also found to be associated with eyespots.
Overall, this study provides the identity of many new genes and signaling pathways associated with eyespots, and suggests that the ancient wound healing gene regulatory network may have been co-opted to cells at the center of the pattern to aid in eyespot origins. New transcription factors that may be providing different identities to distinct wing sectors, and genes with sexually dimorphic expression in the eyespots were also identified.
The evolution of novel traits and of their underlying gene regulatory networks remains a hot topic in the field in evolutionary biology. Some novel traits appear to have arisen via the qualitative modification of pre-existing traits with more primitive features [1, 2]. Examples of such novel traits include angiosperm flowers derived via the modification of the leaves into sepals and petals , the novel mammary gland derived from a skin gland that started secreting novel compounds , or of butterfly eyespots on the wing with concentric rings of color from simpler spots [5, 6]. When such evolutionary transitions take place, a fundamental question that follows pertains to the type and nature of the molecular modifications that took place to lead to the origin of the novel trait.
Proposed mechanisms for the origin of novel traits include the fusion of distinct cell types, as in the case of the insect wing [7,8,9], or the recruitment of pre-existing gene regulatory networks into novel developmental contexts [10,11,12]. In these latter cases, cells from tissues that build the novel trait generally display novel gene expression profiles due to the activation of new gene regulatory networks in those cells. Comparing the gene expression profile of these cells to those of cells serving distinct functions in the body can help identify ancestral and pre-existent gene regulatory networks that were recruited to help build the novel trait . This type of research recently pointed to a novel adaptive trait “beetle horns” having likely been co-opted from the limb gene regulatory network , and to a novel lobe in the adult genitalia of Drosophila originating from the recruitment of a gene regulatory network involved in differentiating larval breathing spiracles earlier in development . Hence, identifying the developmental origin of novel traits can begin with the identification of the underlying genes involved in building the novel trait.
In this study, we are interested in investigating the nature of the gene regulatory network underlying eyespot development. Eyespots are novel traits on the wings of butterflies that serve adaptive roles in both predator avoidance and sexual signaling [15,16,17,18,19,20,21,22]. Previous research suggested that nymphalid butterfly eyespots originated once via the co-option of pre-existing gene networks [5, 6, 23]. Network co-option was suggested because the expression of four out of five genes examined in association with eyespot development in 23 species of nymphalids and outgroups was inferred to have originated concurrently with the origin of eyespots . These data suggest that these genes may not have been recruited gradually and individually to help build the novel trait, but appeared associated with the trait in a more abrupt way via the co-option of a gene regulatory network wherein the genes were already connected to each other. Previous work, however, had already hinted at co-option events being associated with eyespot origins by the type of gene being discovered in association with eyespots. For instance, the first gene ever to be visualized in eyespots, Distal-less, pointed at a possible limb network co-option event . Later, the hedgehog regulatory circuit involved in patterning the anterior-posterior axis of insect wings was also proposed . This was followed by the wound healing gene regulatory network , and a regulatory circuit involved in wing margin development . For a clearer understanding of which, if any, of these circuits and/or networks may have been co-opted to aid in eyespot origins, as well as the discovery of new genes or networks that may have aided in eyespot origins, it is important to examine the complete set of genes involved in eyespot development.
Eyespot development in nymphalid butterflies has been studied since 1978 , but so far only 12 genes have been discovered associated with eyespots in two nymphalid model species, Junonia coenia and Bicyclus anynana (reviewed in Monteiro 2015) mostly developed against fly proteins. Thanks to innovations in gene expression profiling, such as transcriptome analysis on next-generation sequencing platforms, it is now possible to compare total gene expression across different cell types or tissues , and potentially identify the complete set of genes associated with eyespot development.
In this study, we used next-generation sequencing and transciptome analyses to identify the set of genes associated with eyespot development in the early pupal stage of B. anynana. This stage corresponds to an important signaling stage of eyespot development where cells at the center of the pattern signal to surrounding cells to differentiate the color rings [30,31,32]. In addition, this stage was chosen because of practical considerations. The fragile dorsal wing surface epidermis remains attached to the sturdier overlaying cuticle for a few hours after pupation, and this facilitates the dissection of very small pieces of epidermis containing the eyespot central cells, the pattern organizers. We used micro-dissections of wing epidermis from wild-type and from an eyespot mutant of B. anynana, Spotty, with additional eyespots on the wing, to compare the transcriptome of homologous wing sectors with and without eyespots (Fig. 1a). In addition, this approach also allowed us to identify transcription factors differentially expressed between adjacent wing sectors of the forewing, and identify differential expression of genes in male and female eyespots.
Around 400 million paired end reads (~100 bp) were generated from 15 different mRNA libraries. These libraries consisted of three biological replicates of micro-dissected tissue from two distinct wing sectors in males (Cu1 and M3) of wild-type and Spotty individuals, and three biological replicates from a single wing sector in wild-type females (Cu1; Fig. 1). A total of 135,491 contigs were assembled with a N50 value of 1061 bp. The length of the longest contig was 11,724 bp. Around 95% of reads from each library were mapped back to the transcriptome to estimate transcript abundance.
A clustering analysis that compared transcripts across all samples (15 libraries) was performed to check data quality. The analysis revealed that biological replicates and samples from the same type of tissue generally clustered together (Additional file 1: Figure S1). Differentially expressed genes were identified across the 15 libraries using a false discovery rate (FDR) p-value of <0.001 and at least a 2-fold change in gene expression. First, wild-type and Spotty mutant samples were clearly separated from each other indicating many genes were differentially expressed between these tissue samples. Second, within wild-type butterfly tissue, female samples were separated from male samples. Third, M3 and Cu1 wing sectors of wild-type male tissues were separated from each other. Lastly, all samples from Spotty wings, all containing eyespots, clustered together, but one biological replicate of the Spotty Cu1 wing sector clustered closest to tissue from the Spotty M3 wing sector. This could be due to a small number of differentially expressed genes between Cu1 and M3 wing sectors in Spotty wings (see following).
The Bicyclus anynana eyespot-associated differentially-expressed transcripts at early pupation
Genes associated with eyespot development were identified as the common set of differentially expressed genes between wing sectors that do not develop eyespots versus wing sectors that develop eyespots (Table 1). Presumably up-regulated genes represent novel gene expression patterns due to genetic co-option events associated with eyespot development, whereas down-regulated genes represent the downstream effects of these co-options on downstream targets, previously expressed at some basal level in tissues without eyespots. We calculated three sets of differentially expressed genes: between homologous M3 sectors in Wt and Spotty wings; between M3 and Cu1 sectors in Wt wings; and between M3 sectors in Wt wings and Cu1 sectors in Spotty wings (Fig. 1c). These analyses revealed the common set of transcripts that were significantly up- and down-regulated in a similar way across the three comparisons.
A total of 186 transcripts were identified as being associated with eyespot development (Fig. 1c). From these transcripts, 132 were up-regulated and 54 transcripts were down-regulated in the tissues that develop eyespots versus the tissues that do not (Table 1). The data was processed as follows: First, 4056 transcripts were found differentially expressed in M3 wing sectors between Spotty and wild-type (Wt) wings (Fig. 1c). From these transcripts, 2088 were up-regulated and 1968 were down-regulated in the M3 sector with eyespots (in Spotty). Second, 1368 transcripts were differentially expressed between Cu1 and M3 wing sectors in Wt wings (Fig. 1c). From these, 792 were up-regulated and 576 transcripts were down-regulated in the Cu1 wing sector with eyespots. Third, 2197 transcripts were differentially expressed between Cu1 wing sectors in Spotty and M3 wing sector in Wt (Fig. 1c). Of these, 1066 were up-regulated and 1131 down-regulated in the wing sector with eyespots.
To verify that this set of 186 eyespot genes was stable across the sexes, we compared the transcriptomes of Wt female tissues that develop eyespots (Cu1) with Wt male tissues that do not develop eyespots (M3). Interestingly, we found that 80% of the previously identified 186 eyespot-associated genes were differentially expressed in this comparison as well (Additional file 2: Table S1, the 20% of the genes that were not differentially expressed in this comparison are highlighted in green).
All differentially expressed genes were blasted and annotated for gene ontology (GO) terms (Additional file 3: Table S14-19). Many of the genes identified (72, 39%) were not previously annotated in the databases and this included some of the most highly up- and down-regulated genes. From the annotated set, some genes showed large fold expression differences that were highly significant (Fig. 2). Examples of such significantly (p-value <0.001) up-regulated genes in tissue containing eyespots included Alpha-aminoadipic semialdehyde mitochondrial, kal-1 protein, retrovirus-related Pol poly from transposon opus, cuticle 3-like, cysteine-rich motor neuron 1 protein, hypothetical protein KGM_03542 and chitinase-like protein. In addition, RNA-directed DNA polymerase from mobile element jockey, ubiquitin carboxyl-terminal hydrolase isozyme, kish and adenylyltransferase and sulfurtransferase MOCS3 were also four-fold up-regulated in eyespot tissues with a slightly higher p-value (p < 0.05). Other significantly up-regulated genes in eyespots, at less than two-fold expression level differences, included two Toll-like proteins. Significantly down regulated genes (four-fold) in tissues with eyespots included ryanodine receptor-like protein, larval cuticle protein lcp-22-like precursor, acid alpha-glucosidase, kinase suppressor of ras 2-like protein, hypothetical protein RR46_01763, cd63 antigen, farnesoic acid o-methyltransferase-like protein and syntaxin-1A.
In order to test whether genes associated with eyespots were enriched for any particular functional category, we performed an enrichment analysis. None of the annotated gene ontology categories were found to be enriched in eyespot tissues (Fisher’s exact test: p-value <0.05). This can be due to the fact that only a small portion of the genes had GO terms (Additional file 3: Table S14 and Table S15).
Candidate genes and candidate networks
We examined expression profiles of the genes previously associated with eyespot development in B. anynana and in other butterflies in our transcriptome data (Additional file 4: Table S3). With the exception of Ultrabithorax and wingless, the other ten genes were present in the transcriptome. However, with the exception of Distal-less and Antennapedia, most of the eyespot-associated genes were not differentially expressed between tissues with eyespots and tissues without eyespots (Fig. 3).
We further manually examined all the annotated genes (114) for a possible known function in the four gene regulatory networks previously proposed to have been co-opted to aid in eyespot origins. Some of the genes have known roles in wound healing (18 genes), limb development (3 genes), anterior-posterior axis establishment and segment polarity (4 genes), and wing margin development (4 genes), whereas most genes have no known role in these developmental contexts (Additional file 5: Table S2).
Transcription factor differences between M3 and Cu1 wing sectors
During Drosophila wing development each sector of the wing expresses a different combination of transcription factors , and veins form at boundaries of expression of some of these transcription factors [34, 35]. Butterflies may use similar mechanisms to produce their venation patterns. In addition, these transcription factors due to their presence in only specific sectors of the wing may be used as selector genes to modify the eyespot gene regulatory network in those sectors, i.e., to regulate presence/absence and/or eyespot size in a sector-specific fashion. Hence, Spotty, due to its restricted effect on only two sectors of the wing [36, 37], could be a mutation involving such a selector gene. The mutation could be either in the selector gene itself, or alternatively, in the selector’s downstream targets such as in regulatory DNA flanking eyespot network genes. This latter type of mutation is more likely, as it would have fewer pleiotropic effects such as vein disruptions. Given these assumptions we tried to identify candidate wing sector specific genes by comparing differentially expressed genes between wing sectors M3 and Cu1 in both Wt and Spotty wings, and finding the set of common differentially expressed genes, regardless of the presence and absence of eyespots. A total of 52 transcripts were differentially expressed between wing sectors M3 and Cu1 in B. anynana. Of these, 23 transcripts were up-regulated and 29 transcripts were down-regulated in wing sector M3 versus wing sector Cu1 (Table 1). The data was processed as follows: First, we found 576 transcripts up-regulated, and 792 transcripts down-regulated in wing sector M3 relative to Cu1 in Wt butterflies (Fig. 1d). Second, we found 628 transcripts up-regulated and 541 transcripts down-regulated in wing sector M3 relative to Cu1 in the Spotty mutant (Fig. 1d). We then identified the common up, and down-regulated genes between these two comparisons. From this set we identified a single annotated transcription factor, T-box 6, expressed at significantly higher levels in the M3 wing sector relative to the Cu1 sector in both Wt and Spotty wings (Wt: p-value: 1.51E-14 and logFC: 5.59; Spotty: p-value: 8.25E-011 and logFC: 6.89).
Differentially expressed transcripts between female and male Cu1 eyespots
To analyze differentially expressed transcripts between female and male eyespots, wing tissue samples were compared between Wt Cu1 male and Wt Cu1 female wing sectors (Fig. 1e). We found 2785 differentially expressed transcripts. Of these, 1714 were up-regulated (Additional file 6: Table S12) and 1071 were down-regulated (Additional file 6: Table S13) in males versus females (Table 1).
Wound healing and immune system genes are associated with eyespot development
A large fraction (16%) of the 114 annotated differentially expressed genes associated with eyespots have a role in wound healing in other systems (Additional file 5: Table S2). Previous work showed that wounding the pupal wing from 6 to 18 h after pupation induced ectopic eyespots around the area of injury [30, 38, 39]. These ectopic eyespots are similar to normal eyespots in that they have a black disc of scales surrounded by a yellow ring, but lack the central white scales . Additionally, the same transcription factors that map to the different color rings are also expressed in ectopic eyespots after wounding . Due to these similarities, Monteiro et al.  proposed that eyespot evolution involved the co-option of the wound healing gene regulatory network to aid in eyespot development. Here, the transcriptome analysis lends support to this hypothesis by identifying additional eyespot-associated genes known for their roles in wound repair.
Calcium signaling is an early trigger associated with wound repair, and this signaling pathway appears to be involved in the eyespot gene regulatory network. An increase in intracellular calcium was the proposed trigger in wound healing in worms and humans [40, 41] and wounding a different butterfly, Junonia orythia, triggers a rapid calcium wave that spreads from injured areas and leads to temporary Ca2+ increases in epidermal cells [42, 43]. One of the up-regulated genes in eyespots is endothelial-monocyte activating polypeptide ii, which functions in inducing the migration of progenitor cells into wound areas with increasing intracellular calcium mobilization in these cells . In addition, Calcium calmodulin dependent protein kinase (CAMK) is down-regulated in eyespots. CAMK is a Ca2+ dependent kinase  which functions as a negative regulator of innate immune responses to cellular damage in the worm . Its down-regulation in eyespots suggests the possible up-regulation of genes involved in wound healing. A ryanodine receptor like protein is extremely down-regulated in eyespots. Ryanodine receptor is an intracellular calcium channel  that has been implicated in the transmission of slow calcium waves in a variety of tissues. Calcium waves have been observed in Drosophila wings discs, especially associated with areas rich with Wingless signals , but, so far, this receptor is not known to be expressed in Drosophila wing imaginal discs. Wingless and phospho-MAD (a Dpp signal transducer), two other components of the wound repair pathway in Drosophila , are naturally expressed in B. anynana eyespots, shortly after the period of development examined here . These data implicate both calcium signaling and wound healing as candidate pathways involved in the natural process of eyespot development.
Other eyespot-associated genes related to wound healing are the transcription factor D-Fos/Kayak an effector of both the ERK and JNK signaling pathways , and an essential gene for initiation of epithelial repair in Drosophila ; the Ras protein, which induces wound closure in human epidermal cells ; and Chitinase-like protein which regulates immunity and wound closure in Drosophila . Another gene, however, has an opposite expression pattern to that predicted by a simple wound-healing gene regulatory network co-option scenario: Zinc finger protein is up-regulated in eyespots. This protein plays a critical role in wound healing as a negative feedback regulator of the Wnt/B-catenin pathway in mammals . Reduced expression of zinc finger protein accelerates cutaneous wound healing in mice . However, this gene is up-regulated in eyespots in our study. This could indicate that if the wound healing regulatory network was co-opted to eyespot development, it is being naturally slowed down in eyespots.
A Toll-like receptor and a Toll-like protein (also a receptor) are both up-regulated in eyespots. The Toll pathway is known for its ancestral role in activating the innate immune response , as well as a derived role in patterning the dorsal ventral axis of insect embryos . Recently, the Toll pathway has also been implicated in wound epidermal closure in Drosophila . Both in flies and in worms the immune system is activated by wounding to promote cell survival [40, 43, 55], and signals other than the presence of microbes, such as sterile wounding, can also activate the innate immune response of insects in the absence of an infection. The Toll signaling pathway is also activated in silkworms during immune response , and Toll is essential for wound closure in late Drosophila embryos . Interestingly, Toll is also expressed in the ectodermal cells that form the leading edge of dorsal closure in Drosophila embryos and at intersegmental boundaries . Most recently, the over-expression of constitutively activated Toll was sufficient to induce wound repair genes in the undamaged epidermis of Drosophila , suggesting that Toll regulates a modular network that may have been co-opted to aid in eyespot origins . Both Toll and its ligand, Spätzle, should be investigated in future regarding their possible role in eyespot development.
Other new genes and signaling pathways associated with eyespot development
A paralogue of the yellow gene, yellow b, was found to be down-regulated in eyespots. This gene is downstream of JNK signaling in Drosophila and is expressed in the leading edge of the cells involved in dorsal closure . Both JNK and Wg signaling are known to up-regulate dpp expression also in the leading edge cells  and both Wg and pMad (a Dpp effector gene) transcripts have been visualized in eyespot centers shortly after pupation . This suggests similarities between dorsal closure and eyespot gene regulatory networks in terms of gene identity, but differences at the level of gene regulation (up-regulation of yellow b in dorsal closure, down-regulation in eyespots).
This study highlighted a new signaling pathway, the fibroblast growth factor receptor (FGFR) signaling pathway, and confirmed a previously implicated pathway, the bone morphogenetic protein (BMP) pathway, associated with eyespot development at the early stages of pupation, during the focal signaling stage. Two of the differentially regulated genes in eyespots were kal-1 and Bmp-binding endothelial regulator protein. Kal-1, which was highly up-regulated in eyespots, codes for anosmin-1 protein, which is involved in the FGFR signaling pathway , a pathway with multiple roles in organogenesis, cell differentiation, and wound healing . This gene is also expressed in the embryonic cells of Drosophila that will develop into the Keilin’s organs, a larval sensory organ that develops from a sub-set of Dll expressing cells in the thorax, the remaining cells developing into the adult wings and legs [66, 67]. Kal-1 is also expressed in the ectodermal cells that form the leading edge of dorsal closure in Drosophila embryos . The second gene, Bmp-binding endothelial regulator protein, interacts and inhibits Bone morphogenetic protein (BMP) function  but it is down-regulated in eyespot centers, suggesting a role for BMP signaling in promoting eyespot development. Another BMP signaling member, the signal transducer pMAD, was previously visualized in eyespot centers in pupal wings slightly later in development .
Our data suggests that local cell division and proliferation may be happening in eyespot-containing tissues as we see the up-regulation of polymerases (RNA-directed DNA polymerase from mobile element jockey, and two other “pol-like proteins”) and up-regulation of an apoptosis inhibitor (apoptosis inhibitor 5) associated with eyespots. Scales in the eyespot centers are generally larger than in the neighboring area  and these scale cells have increased nuclear volume, suggesting DNA synthesis during the early pupal stages in butterflies . In addition, local cell division in eyespot centers was previously documented to occur during the wandering stages of development . Furthermore, tp53 regulating kinase, required to support proliferation and cell growth in Drosophila , is also up-regulated in eyespots. This suggests that cell division and/or cell growth in eyespot centers may continue from the wandering larval stage up to the early pupal stages.
New genes with a possible role in eyespot plasticity
Three genes may be playing a yet undiscovered role in eyespot plasticity of dorsal eyespots in B. anynana : alpha-aminoadipic semialdehyde, Farnesoic acid O-methyl transferase (FAMeT) and Painless. The top (annotated) up-regulated gene in eyespots is alpha-aminoadipic semialdehyde also known in Drosophila as lysine ketoglutarate reductase/saccharopine dehydrogenase (LKR/SDH), a co-factor that represses ecdysone mediated transcription of target genes . Ecdysone signaling has been implicated in the regulation of ventral and dorsal eyespot size plasticity at slightly earlier stages of eyespot development (during the wandering stage) [71, 75]. The high expression levels of this repressive co-factor in the early pupal stages may contribute to explain the relative insensitivity of dorsal eyespots to 20E injections performed at this stage of development [74, 76]. The second gene, FAMeT, is an enzyme in the biosynthetic pathway of juvenile hormone (JH) . It is down-regulated in the dorsal eyespots of the wet season form of male butterflies, used in this study. This suggests the possible involvement of JH in either eyespot development, or in the regulation of size and brightness plasticity of dorsal eyespots in males of this species across seasonal forms [73, 78]. Different expression levels of certain isoforms of this enzyme have been previously associated with caste differentiation between queens and workers of the stingless bee . The third gene, Painless, is a transient receptor potential channel that acts as a primary harmful heat detector in Drosophila . High heat (> 42.6 °C) activates Painless expression. It is possible that this gene, if sensitive to lower ambient temperatures in B. anynana, may be also involved in regulating eyespot-specific patterns of temperature-mediated plasticity in male eyespots. There is an interesting connection between this gene and calcium signaling described above. Ca2+ is known to regulate the function of the painless channel in Drosophila, and in the absence of Ca2+ painless fails to respond to heat .
Previously discovered eyespot-associated genes are poorly represented in the transcriptome data
The transcriptome data contained ten of the twelve genes previously associated with eyespot development in B. anynana, and other butterflies, at earlier and later stages of eyespot development. These genes include the hox gene Antennapedia (Antp), the earliest known up-regulated gene in larval wing discs , the additional transcription factors Distal-less (Dll) , spalt (sal) [26, 31], engrailed (en) , cubitus interruptus (ci) , the cell-surface receptors Notch (N)  and patched (ptc) , the nuclear receptor Ecdysone receptor (EcR) [82,83,84,85], the signal transducer smad , and the segment polarity gene hedgehog (hh)  (Fig. 3). Two other known eyespot-associated genes in B. anynana, Ultrabithorax (Ubx), a hox gene and wingless (wg), a ligand, were absent from our transcriptome data. The absence of Ubx is not unexpected because Ubx is expressed only in hindwings in B. anynana , and in this study tissues were dissected from forewings. The absence of wg transcripts might be related to developmental timing. The earliest wg mRNA and wg protein expression was identified in the eyespot centers at 10 h and 10.5 h after pupation, respectively [26, 87], and here our dissections were done earlier, between 3 and 6 h. Hedgehog signaling involving hedgehog (hh), and its receptor patched (ptc) were previously associated with eyespot development in Junonia coenia butterflies , but neither hh nor ptc are expressed in B. anynana eyespots , and do not seem to function in eyespot development in this species .
A surprising finding was that most of the genes previously associated with eyespot development in B. anynana were not consistently significantly differentially expressed in tissues with and without eyespots, with the exception of Antp and Dll. Several possibilities may account for this. The first is that some genes are expressed in cells other than those in the central eyespot region. This non-specificity may have lowered our power to detect significant differential expression across eyespotted and non-eyespotted samples. For instance, Antp is expressed strongly and only in eyespot center cells, unlike other genes which are expressed in other cells of the wing sector [80, 86]. For instance, en is expressed in all cells of the posterior compartment [80, 88], and N is expressed in cells along the midline of each wing sector and also in the wing margin [24, 81, 90].
A second possibility is that since most of the genes previously associated with eyespot development were identified at the protein level, with antibodies (En, Ci, N, Sal, EcR, pSmad, Wg), instead of at the mRNA level (Dll, hh, ptc) with in situ hybridizations, it is possible that the protein expression differences are the product of mechanisms of post-transcriptional regulation, rather that transcriptional regulation. Experiments in various organisms including fly and mammalian cells often show that only 40% of the variation in protein levels can be correlated with the abundance of their corresponding mRNAs, and processes involving mRNA or protein degradation can explain this lack of correlation .
A third possibility relates to gene expression patterns potentially being very dynamic over time and being absent in eyespots at the time the tissue samples were taken. Most of the genes previously associated with eyespot development were visualized either in the late larval stage and/or around 12-30 h after pupation. No gene has been directly visualized/quantified in eyespots at 3-6 h after pupation because wings are fragile and difficult to handle at this stage. It is possible, thus, that some of the genes only visualized in larval wings (ci, N) stop being differentially expressed in young pupal wings, or, alternatively, some of the genes only visualized in older pupal wings (wg, psmad) start their expression after the time interval examined here. A recent qPCR study of wg expression in whole early pupal forewings showed very low expression just after pupation, but detectable levels 6 h later .
The transcription factor Distal-less (Dll) was significantly up-regulated in the eyespotted tissues compared to non-eyespotted tissues. This is an expected result since Dll is expressed in eyespot centers of butterfly wings [24, 26, 90] and has a functional role in eyespot formation [92,93,94]. Although, two studies demonstrated that this gene is a positive regulator in eyespot development [93, 94], another study claimed a repressive regulatory function for Dll in eyespot formation . Our results support the view of a positive regulatory function of Dll in eyespot development as Dll is expressed in areas with eyespots.
Wing sector specific genes
Several genes were differentially expressed between Cu1 and M3 wing sectors, regardless of presence or absence of eyespots in these sectors. The only annotated transcription factor on our list was T-box 6. This gene was over-expressed in M3 wing sectors. A T-box 6 homolog in Drosophila is bifid, also known as optomotor-blind (omb) [95, 96]. omb is broadly expressed in the middle sectors of the fly wing imaginal disc  and is required to induce development of the L5 vein in Drosophila [35, 97]. According to the vein nomenclature of Stark et al.  the L5 vein in fly wings is homologous to the CuA1 vein (Cu1 vein here), placed between the M3 and Cu1 wing sectors. Hence, our results indicate that the B. anynana omb ortholog may have a similar expression pattern and function to omb in Drosophila. It will be interesting to determine, in future, the exact pattern of expression of this gene on a developing butterfly wing using in situ hybridizations, and the role of this gene, if any, in the control of vein patterning and eyespot development.
Genes differentially expressed in male and female eyespots
Many genes were differentially expressed between male and female Cu1 dorsal wing sectors with eyespots. Male Cu1 dorsal eyespots display different sizes and patterns of brightness plasticity relative to female eyespots [73, 78], so some of the genes below may be involved in playing a role in the control of this process. The genes that stand out are all expressed at higher levels in females and include Methoprene-tolerant, the Juvenile hormone (JH) receptor , JH-inducible protein, a gene that shows rapid induction in the presence of JH , juvenile hormone binding protein, a protein carrier of JH to the target cells [101, 102], bombyxin b-4 precursor, a mRNA precursor of bombyxin, an insulin-like peptide, which contains many family members [103,104,105,106]. The finding that a bombyxin mRNA precursor is expressed at higher levels in female eyespots is especially interesting because bombyxin family members were previously only known to be expressed in the larval brain, ovaries, fat body, and mid gut of the lepidopteran Bombyx mori . Our data suggest a novel expression domain for the B4 family members and their possible role in the regulation of eyespot sexual dimorphism in B. anynana.
The genes associated with eyespot development identified in this work may contain some of the original genes involved in eyespot origins, 90 million years ago , as well as B. anynana lineage-specific co-options that are not shared across all nymphalid species with eyespots. These lineage-specific co-options are likely to have occurred given the span of time involved since eyespot origins. Genes likely to belong to this latter category include genes regulating eyespot phenotypic plasticity, a likely derived trait within nymphalid butterflies , as well as sexual dimorphism in dorsal eyespots. Future comparative transcriptomic experiments, performed across multiple nymphalid species with eyespots, will be required to identify the core set of genes likely associated with eyespot origins and necessary for eyespot development. Currently, our data lends additional support to the hypothesis proposing that eyespots co-opted the wound healing gene regulatory network to aid in their development . The other three hypotheses were not as well supported by the current data because only fewer genes known to be associated with the proposed gene regulatory circuits and networks were also found associated with eyespots. This, however, may simply reflect our incomplete current understanding of these other developmental processes. Future studies should examine genes associated with eyespot development at earlier and later time points, although this may be experimentally challenging. Overall, our study identified many new additional genes associated with eyespot development, eyespot plasticity, and eyespot sexual dimorphism in B. anynana that can now be examined across species and also at the functional level.
Butterflies were reared in climate controlled chambers at 27 °C on a 12 L: 12D photoperiod. Humidity was kept at 80%. Larvae were fed with corn plants and adults with mashed banana.
Crossing wild-type individuals with mutant individuals
Wild-type and Spotty lines have been kept separate in the lab for many generations. As a consequence, they may have fixed distinct genetic variants that are not related to the wing pattern mutation. In order to homogenize the genetic background between these lines for the purpose of comparative transcriptomics, female Wt individuals were crossed with Spotty male individuals (Additional file 7: Figure S2). All F1 (first generation) offspring had phenotypes associated with Wt/Spotty heterozygous and displayed small intermediate eyespots in sectors M2 and M3. These F1 individuals were mated with each other to produce a F2 generation. Only wild-type and Spotty homozygous individuals were selected from the F2 generation and these were mated to individuals of the same phenotype to produce two pure-breeding F3 generation cohorts. Female and male wild-type and Spotty pupae from this generation were used for RNA extractions.
Collecting RNA from wing tissue
Samples from the eyespot centers on the dorsal pupal forewings were used to extract RNA for transcriptome analysis. Small squares of wing tissue of around 0.5 × 0.5 mm, containing the eyespot centers, were dissected between 3 and 6 h after pupation using home-made scalpels with small chips of a razor blade glued to a metal handle (see size of squares in Fig. 1a). Pupation time was estimated with time-lapse photography of pre-pupae, with 30 min intervals between photos, using a Kodak DC290 digital camera. These tissues were collected from two wing sectors, M3 and Cu1 (Fig. 1a). Dissected tissue samples were collected into an eppendorf tube containing RNAlater reagent (Qiagen) to stabilize the RNA. Around 30 pieces of tissue from the same wing sector and sex were pooled before RNA was extracted with a RNeasy Mini Kit (Qiagen), using a RNase-Free DNase (Qiagen) step. Each of these extractions had concentration between 4 ng/ul and 17 ng/ul in 70 ul of water and was used to make a separate library for RNA sequencing. Three separate libraries were made for each type of tissue to represent biological replicates. RNA concentration was measured with both a Nanodrop 1000 Spectrophotometer (NanoDrop/Thermo Scientific) and Qubit 2.0 Fluorometer (Life Technologies). RNA quality was controlled by agarose gel electrophoresis for RNA degradation and potential contamination, and Agilent 2100 bioanalyser for RNA integrity by AITbiotech. Dissections of M3 and Cu1 squares were performed in random order to avoid any gene expression biases caused by the operation.
Library preparation and RNA sequencing
RNA library construction was performed with Truseq stranded mRNA kit from Illumina by AITbiotech. mRNA was purified from equal amount of total RNA from each library using oligo(dT) beads to minimize bias for inter-library comparisons. The mRNA was fragmented randomly in fragmentation buffer, followed by cDNA synthesis using random hexamers and reverse transcriptase. After first-strand synthesis, a custom second-strand synthesis buffer (Illumina) was added, with dNTPs, RNase H and E. coli polymerase I to generate the second strand and AMPure XP beads were used to purify the cDNA. The final cDNA library was ready after a round of purification, terminal repair, A-tailing, ligation of sequencing adapters, size selection and PCR enrichment (Illumina). Library concentration was first quantified using a Qubit 2.0 fluorometer (Life Technologies), and then diluted to 1 ng/μl before checking insert size on an Agilent 2100 Bioanalyzer and quantifying to greater accuracy by quantitative PCR (Kapa Biosystems) (library activity >2 nM). Libraries were loaded into a Illumina Hiseq 2000 (next-generation sequencing machine) according to activity determined from qPCR and expected data volume of the instrument by AITbiotech. One lane was used to sequence a total 15 libraries from wild-type and Spotty lines. After sequencing, adaptors were removed from sequences and initial quality check reports for raw reads were provided by AITbiotech.
We first re-constructed transcripts from the short reads using a “de novo assembly” approach . Around 400 million paired-end reads were used for de novo assembly. Both Trinity and CLC genomics workbench were used for the transcriptome assembly using different assembly and trimming parameters. The quality of the various assembled transcriptomes was assessed according to average transcript’s length, the number of transcripts, and the coverage of aligned transcripts against the protein database of another lepidopteran, Bombyx mori. Transcriptomes were deemed to be of higher quality based on longest N50 value, fewer number of transcripts, and highest full-length coverage of all known proteins of Bombyx mori. According to these criteria, the best transcriptome was built using Trinity software using all 15 libraries and the following parameters: Quality trimming was performed with trimmomatic , which is integrated into Trinity software. Reads were trimmed at the start and the end of a read, if the threshold quality was below 5, and reads less than 30 bp were removed prior to the assembly (e.g., LEADING:5, TRAILING:5, MINLEN:30). After assembly, transcripts under 200 bp were removed from the final transcriptome.
Identifying differentially expressed genes
Differentially expressed genes were identified with the downstream applications of the Trinity software package . These methods use RSEM software for estimating transcript abundance  and EdgeR Bioconductor package for identifying differentially expressed genes across samples after sample normalization for scaling library size using a weighted trimmed mean of the log expression ratios (trimmed mean of M values (TMM)) [111, 112]. A p-value of 0.05 was chosen as the threshold to identify statistically significant up and down-regulated genes between libraries.
Gene annotations and enrichment analysis
Differentially expressed genes were blasted against a non-redundant protein database with blastx using an E-value <1e-3 (E-value (expected value) shows the number of significant alignments expected to occur by chance) and were annotated for gene functions using Blast2GO software . Afterwards, these transcripts were annotated using Blast2GO and were added into the list of previously annotated genes for further analyses. Gene ontology terms were found using non-redundant protein database in Blast2GO (E-value <1e-3). A test for functional enrichment of particular gene ontologies among the 186 eyespot associated genes was performed using Fisher’s exact test and the Drosophila database as a reference on Blast2GO PRO plugin for CLC genomic workbench.
Co-option of the gene regulatory networks
All of the annotated eyespot associated genes were searched for their functions, related to the proposed co-opted gene regulatory network using keywords; “wound healing”, “limb development”, “anterior-posterior axis”, “segment polarity” and “wing margin” in the Google Scholar search engine.
Candidate gene approach
First, the sequences of the previously identified candidate eyespot genes were downloaded from the nucleotide database of NCBI. Second, these sequences were aligned to the transcriptome with a Blastn command in BLAST+ (E-value <1e-20) and the longest transcript of the corresponding gene was selected from the Blastn results. Later, a p-value and a logFC value of the transcripts of the candidate eyespot genes were obtained from the M3 wing sectors comparison between wild-type and Spotty tissues.
Bone morphogenetic protein
Extracellular signal–regulated kinase
False discovery rate
Fibroblast growth factor
Fibroblast growth factor receptor
Jun N-terminal kinase
Logarithmic fold change
Shubin N, Tabin C, Carroll S. Deep homology and the origins of evolutionary novelty. Nature. 2009;457(7231):818–23. doi:10.1038/nature07891. PubMed PMID: 19212399
Wagner GP, Lynch VJ. Evolutionary novelties. Curr Biol. 2010;20(2):R48–52. https://doi.org/10.1016/j.cub.2009.11.010. PubMed PMID: 20129035
Strickberger MW. Evolution: Jones and Bartlett; 2000. p. 308.
Vorbach C, Capecchi MR, Penninger JM. Evolution of the mammary gland from the innate immune system? BioEssays. 2006;28(6):606–16. doi:10.1002/bies.20423. PubMed PMID: 16700061
Oliver JC, Beaulieu JM, Gall LF, Piel WH, Monteiro A. Nymphalid eyespot serial homologues originate as a few individualized modules. Proc Biol Sci. 2014;281(1787) https://doi.org/10.1098/rspb.2013.3262. PubMed PMID: 24870037; PubMed Central PMCID: PMCPMC4071533
Oliver JC, Tong XL, Gall LF, Piel WH, Monteiro AA. Single origin for nymphalid butterfly eyespots followed by widespread loss of associated gene expression. PLoS Genet. 2012;8(8):e1002893. doi:10.1371/journal.pgen.1002893. PubMed PMID: 22916033; PubMed Central PMCID: PMCPMC3420954
Clark-Hachtel CM, Linz DM, Tomoyasu Y. Insights into insect wing origin provided by functional analysis of vestigial in the red flour beetle, Tribolium Castaneum. Proc Natl Acad Sci. 2013;110(42):16951–6. doi:10.1073/pnas.1304332110.
Medved V, Marden JH, Fescemyer HW, Der JP, Liu J, Mahfooz N, et al. Origin and diversification of wings: insights from a neopteran insect. Proc Natl Acad Sci U S A. 2015;112(52):15946–51. doi:10.1073/pnas.1509517112. PubMed PMID: WOS:000367234700060
Niwa N, Akimoto-Kato A, Niimi T, Tojo K, Machida R, Hayashi S. Evolutionary origin of the insect wing via integration of two developmental modules. Evol Dev. 2010;12(2):168–76.
Gao F, Davidson EH. Transfer of a large gene regulatory apparatus to a new developmental address in echinoid evolution. Proc Natl Acad Sci U S A. 2008;105(16):6091–6. doi:10.1073/pnas.0801201105. PubMed PMID: WOS:000255356000027
Glassford WJ, Johnson WC, Dall NR, Smith SJ, Liu Y, Boll W, et al. Co-option of an ancestral Hox-regulated network underlies a recently evolved morphological novelty. Dev Cell. 2015;34(5):520–31. doi:10.1016/j.devcel.2015.08.005. PubMed PMID: 26343453; PubMed Central PMCID: PMCPMC4573913
Moczek AP, Rose DJ. Differential recruitment of limb patterning genes during development and diversification of beetle horns. Proc Natl Acad Sci U S A 2009;106(22):8992-8997. doi:10.1073/pnas.0809668106. PubMed PMID: ISI:000266580500040.
Monteiro A. Evolution and Development: Molecules. The Princeton guide to evolution. 11. Princeton, Oxford: Princeton University Press; 2014. p. 444-452.
Moczek AP, Rose DJ. Differential recruitment of limb patterning genes during development and diversification of beetle horns. Proc Natl Acad Sci U S A 2009;106(22):8992-8997. doi: 10.1073/pnas.0809668106. PubMed PMID: 19451631; PubMed Central PMCID: PMCPMC2690047.
Vallin A, Jakobsson S, Lind J, Wiklund C. Prey survival by predator intimidation: an experimental study of peacock butterfly defence against blue tits. Proc Biol Sci. 2005;272(1569):1203–7. doi:10.1098/rspb.2004.3034. PubMed PMID: 16024383; PubMed Central PMCID: PMCPMC1564111
Kodandaramaiah U. The evolutionary significance of butterfly eyespots. Behav Ecol. 2011;22(6):1264–71. https://doi.org/10.1093/beheco/arr123.
Monteiro A. Origin, development, and evolution of butterfly eyespots. Annu Rev. Entomol. 2015;60:253–71. doi:10.1146/annurev-ento-010814-020942. PubMed PMID: 25341098
Westerman EL, Hodgins-Davis A, Dinwiddie A, Monteiro A. Biased learning affects mate choice in a butterfly. Proc Natl Acad Sci U S A. 2012;109(27):10948–53. doi:10.1073/pnas.1118378109. PubMed PMID: 22689980; PubMed Central PMCID: PMCPMC3390860
Oliver JC, Robertson KA, Monteiro A. Accommodating natural and sexual selection in butterfly wing pattern evolution. Proc Biol Sci. 2009;276(1666):2369–75. doi:10.1098/rspb.2009.0182. PubMed PMID: 19364741; PubMed Central PMCID: PMCPMC2690465
Robertson KA, Monteiro A. Female Bicyclus anynana butterflies choose males on the basis of their dorsal UV-reflective eyespot pupils. Proc Biol Sci. 2005;272(1572):1541–6. https://doi.org/10.1098/rspb.2005.3142. PubMed PMID: 16048768; PubMed Central PMCID: PMCPMC1559841
Stevens M. The role of eyespots as anti-predator mechanisms, principally demonstrated in the Lepidoptera. Biol Rev. Camb Philos Soc. 2005;80(4):573–88. https://doi.org/10.1017/S1464793105006810. PubMed PMID: 16221330
Prudic KL, Jeon C, Cao H, Monteiro A. Developmental plasticity in sexual roles of butterfly species drives mutual sexual ornamentation. Science. 2011;331(6013):73–5. https://doi.org/10.1126/science.1197114. PubMed PMID: 21212355
Schachat SR, Oliver JC, Monteiro A, et al. BMC Evol Biol. 2015;15:20. doi:10.1186/s12862-015-0300-x. PubMed PMID: 25886182; PubMed Central PMCID: PMCPMC4335541
Carroll SB, Gates J, Keys DN, Paddock SW, Panganiban GEF, Selegue JE, et al. Pattern-formation and eyespot determination in butterfly wings. Science. 1994;265(5168):109–14. PubMed PMID: WOS:A1994NV30100038
Keys DN, Lewis DL, Selegue JE, Pearson BJ, Goodrich LV, Johnson RL, et al. Recruitment of a hedgehog regulatory circuit in butterfly eyespot evolution. Science. 1999;283(5401):532–4. PubMed PMID: 9915699
Monteiro A, Glaser G, Stockslager S, Glansdorp N, Ramos D. Comparative insights into questions of lepidopteran wing pattern homology. BMC Dev Biol. 2006;6:52. doi:10.1186/1471-213X-6-52. PubMed PMID: 17090321; PubMed Central PMCID: PMCPMC1654149
Held LI. Rethinking butterfly eyespots. Evol Biol. 2012;40(1):158–68. doi:10.1007/s11692-012-9198-z.
Nijhout HF. Wing pattern formation in Lepidoptera - model. J Exp Zool. 1978;206(2):119–36. PubMed PMID: WOS:A1978FT73200001
Adams J. Transcriptome: connecting the genome to gene function. Nature. Education. 2008;1(1):195.
French V, Brakefield PM. The development of eyespot patterns on butterfly wings - Morphogen sources or sinks. Development. 1992;116(1):103–9. PubMed PMID: WOS:A1992JQ86400010
Brunetti CR, Selegue JE, Monteiro A, French V, Brakefield PM, Carroll SB. The generation and diversification of butterfly eyespot color patterns. Curr Biol. 2001;11(20):1578–85. PubMed PMID: WOS:000171651700015
Nijhout HF. Pattern formation on Lepidopteran wings - determination of an eyespot. Dev Biol. 1980;80(2):267–74. PubMed PMID: WOS:A1980KS92900002
Ibrahim DM, Biehs B, Kornberg TB, Klebes A. Microarray comparison of anterior and posterior drosophila wing imaginal disc cells identifies novel wing genes. G3 (Bethesda). 2013;3(8):1353–62. doi:10.1534/g3.113.006569. PubMed PMID: 23749451; PubMed Central PMCID: PMCPMC3737175
Matsuda S, Yoshiyama N, Kunnapuu-Vulli J, Hatakeyama M, Shimmi O. Dpp/BMP transport mechanism is required for wing venation in the sawfly Athalia rosae. Insect Biochem Mol Biol. 2013;43(5):466–73. doi:10.1016/j.ibmb.2013.02.008. PubMed PMID: 23499566
Blair SS. Wing vein patterning in drosophila and the analysis of intercellular signaling. Annu Rev. Cell Dev Biol. 2007;23:293–319. doi:10.1146/annurev.cellbio.23.090506.123606. PubMed PMID: 17506700
Brakefield PM, French V, et al. Acta Biotheor. 1993;41(4):447–68. PubMed PMID: WOS:A1993NB57300012
Monteiro A, Chen B, Scott LC, Vedder L, Prijs HJ, Belicha-Villanueva A, et al. The combined effect of two mutations that alter serially homologous color pattern elements on the fore and hindwings of a butterfly. BMC Genet. 2007;8:22. doi:10.1186/1471-2156-8-22. PubMed PMID: 17498305; PubMed Central PMCID: PMCPMC1878498
French V, Brakefield PM. Eyespot development on butterfly wings: the focal signal. Dev Biol. 1995;168(1):112–23. doi:10.1006/dbio.1995.1065. PubMed PMID: 7883067
Nijhout HF. Cautery-induced colour patterns in Precis coenia (Lepidoptera: Nymphalidae). J Embryol Exp Morphol. 1985;86:191–203. PubMed PMID: 4031740
Xu S, Chisholm ADA. Galphaq-ca(2)(+) signaling pathway promotes actin-mediated epidermal wound closure in C. Elegans. Curr Biol. 2011;21(23):1960–7. doi:10.1016/j.cub.2011.10.050. PubMed PMID: 22100061; PubMed Central PMCID: PMCPMC3237753
Lansdown AB. Calcium: a potential central regulator in wound healing in the skin. Wound Repair Regen. 2002;10(5):271–85. PubMed PMID: 12406163
Ohno Y, Otaki JM. Spontaneous long-range calcium waves in developing butterfly wings. BMC Dev Biol. 2015;15:17. doi:10.1186/s12861-015-0067-8. PubMed PMID: 25888365; PubMed Central PMCID: PMCPMC4445562
Wood W. Wound healing: calcium flashes illuminate early events. Curr Biol. 2012;22(1):R14–6. doi:10.1016/j.cub.2011.11.019. PubMed PMID: 22240471
Hou Y, Plett PA, Ingram DA, Rajashekhar G, Orschell CM, Yoder MC, et al. Endothelial-monocyte-activating polypeptide II induces migration of endothelial progenitor cells via the chemokine receptor CXCR3. Exp Hematol. 2006;34(8):1125–32. doi:10.1016/j.exphem.2006.05.021. PubMed PMID: 16863920
Carafoli E, Klee C. Calcium as Cellular Regulator. New York Oxford: Oxford University Press; 1999.
Tong A, Lynn G, Ngo V, Wong D, Moseley SL, Ewbank JJ, et al. Negative regulation of Caenorhabditis Elegans epidermal damage responses by death-associated protein kinase. Proc Natl Acad Sci U S A. 2009;106(5):1457–61. doi:10.1073/pnas.0809339106. PubMed PMID: 19164535; PubMed Central PMCID: PMCPMC2629440
Lanner JT, Georgiou DK, Joshi AD, Hamilton SL. Ryanodine receptors: structure, expression, molecular details, and function in calcium release. Cold Spring Harb Perspect Biol. 2010;2(11):a003996. doi:10.1101/cshperspect.a003996. PubMed PMID: 20961976; PubMed Central PMCID: PMCPMC2964179
Balaji R, Bielmeier C, Harz H, Bates J, Stadler C, Hildebrand A, et al. Calcium spikes, waves and oscillations in a large, patterned epithelial tissue. Sci Rep 2017;7:42786. doi: 10.1038/srep42786. PubMed PMID: 28218282; PubMed Central PMCID: PMCPMC5317010.
Smith-Bolton RK, Worley MI, Kanda H, Hariharan IK. Regenerative growth in drosophila imaginal discs is regulated by wingless and Myc. Dev Cell. 2009;16(6):797–809. 10.1016/j.devcel.2009.04.015. PubMed PMID: 19531351; PubMed Central PMCID: PMCPMC2705171
Ciapponi L, Jackson DB, Mlodzik M, Bohmann D. Drosophila Fos mediates ERK and JNK signals via distinct phosphorylation sites. Genes Dev. 2001;15(12):1540–53. PubMed PMID: WOS:000169334200009
Ramet M, Lanot R, Zachary D, Manfruelli P. JNK Signaling pathway is required for efficient wound healing in drosophila. Dev Biol 2002;241(1):145-156. doi: https://doi.org/10.1006/dbio.2001.0502. PubMed PMID: 11784101.
Tscharntke M, Pofahl R, Krieg T, Haase I. Ras-induced spreading and wound closure in human epidermal keratinocytes. FASEB J 2005;19(13):1836-1838. doi: https://doi.org/10.1096/fj.04-3327fje. PubMed PMID: 16170018.
Kucerova L, Broz V, Arefin B, Maaroufi HO, Hurychova J, Strnad H, et al. The drosophila Chitinase-like protein IDGF3 is involved in protection against nematodes and in wound healing. J Innate Immun. 2015; doi:10.1159/000442351. PubMed PMID: 26694862
Lee SH, Kim MY, Kim HY, Lee YM, Kim H, Nam KA, et al. The Dishevelled-binding protein CXXC5 negatively regulates cutaneous wound healing. J Exp Med. 2015;212(7):1061–80. doi:10.1084/jem.20141601. PubMed PMID: 26056233; PubMed Central PMCID: PMCPMC4493411
De Gregorio E, Han SJ, Lee WJ, Baek MJ, Osaki T, Kawabata S, et al. An immune-responsive Serpin regulates the melanization cascade in drosophila. Dev Cell. 2002;3(4):581–92. PubMed PMID: 12408809
Lynch JA, Roth S. The evolution of dorsal-ventral patterning mechanisms in insects. Genes Dev. 2011;25(2):107–18. doi:10.1101/gad.2010711. PubMed PMID: 21245164; PubMed Central PMCID: PMCPMC3022256
Capilla A, Karachentsev D, Patterson RA, Hermann A, Juarez MT, McGinnis W. Toll pathway is required for wound-induced expression of barrier repair genes in the drosophila epidermis. Proc Natl Acad Sci U S A. 2017;114(13):E2682–E8. https://doi.org/10.1073/pnas.1613917114. PubMed PMID: 28289197; PubMed Central PMCID: PMCPMC5380074
Ma Z, Li C, Pan G, Li Z, Han B, Xu J, et al. Genome-wide transcriptional response of silkworm (Bombyx Mori) to infection by the microsporidian Nosema Bombycis. PLoS One. 2013;8(12):e84137. https://doi.org/10.1371/journal.pone.0084137. PubMed PMID: 24386341; PubMed Central PMCID: PMCPMC3875524
Carvalho L, Jacinto A, Matova N. The toll/NF-kappaB signaling pathway is required for epidermal wound repair in drosophila. Proc Natl Acad Sci U S A. 2014;111(50):E5373–82. https://doi.org/10.1073/pnas.1408224111. PubMed PMID: 25427801; PubMed Central PMCID: PMCPMC4273363
Hashimoto C, Gerttula S, Anderson KV. Plasma membrane localization of the toll protein in the syncytial drosophila embryo: importance of transmembrane signaling for dorsal-ventral pattern formation. Development. 1991;111(4):1021–8. PubMed PMID: 1879347
Monteiro A, Podlaha O. Wings, horns, and butterfly eyespots: how do complex traits evolve? PLoS Biol. 2009;7(2):e37. https://doi.org/10.1371/journal.pbio.1000037. PubMed PMID: 19243218; PubMed Central PMCID: PMCPMC2652386
Rousset R, Carballes F, Parassol N, Schaub S, Cerezo D, Noselli S. Signalling crosstalk at the leading edge controls tissue closure dynamics in the drosophila embryo. PLoS Genet. 2017;13(2):e1006640. https://doi.org/10.1371/journal.pgen.1006640. PubMed PMID: 28231245; PubMed Central PMCID: PMCPMC5344535
McEwen DG, Cox RT, Peifer M. The canonical Wg and JNK signaling cascades collaborate to promote both dorsal closure and ventral patterning. Development. 2000;127(16):3607–17. PubMed PMID: 10903184
Gonzalez-Martinez D, Kim SH, Hu Y, Guimond S, Schofield J, Winyard P, et al. Anosmin-1 modulates fibroblast growth factor receptor 1 signaling in human gonadotropin-releasing hormone olfactory neuroblasts through a heparan sulfate-dependent mechanism. J Neurosci. 2004;24(46):10384–92. https://doi.org/10.1523/JNEUROSCI.3400-04.2004. PubMed PMID: 15548653
Coumoul X, Deng CX. Roles of FGF receptors in mammalian development and congenital diseases. Birth Defects Res C Embryo Today. 2003;69(4):286–304. https://doi.org/10.1002/bdrc.10025. PubMed PMID: 14745970
Andrenacci D, Le Bras S, Rosaria Grimaldi M, Rugarli E, Graziani F. Embryonic expression pattern of the drosophila Kallmann syndrome gene kal-1. Gene Expr Patterns. 2004;5(1):67–73. https://doi.org/10.1016/j.modgep.2004.06.004. PubMed PMID: 15533820
McKay DJ, Estella C, Mann RS. The origins of the drosophila leg revealed by the cis-regulatory architecture of the Distalless gene. Development. 2009;136(1):61–71. https://doi.org/10.1242/dev.029975. PubMed PMID: 19036798; PubMed Central PMCID: PMCPMC2653810
Heinke J, Wehofsits L, Zhou Q, Zoeller C, Baar KM, Helbing T, et al. BMPER is an endothelial cell regulator and controls bone morphogenetic protein-4-dependent angiogenesis. Circ Res. 2008;103(8):804–12. https://doi.org/10.1161/CIRCRESAHA.108.178434. PubMed PMID: 18787191
Iwata M, Otaki JM. Spatial patterns of correlated scale size and scale color in relation to color pattern elements in butterfly wings. J Insect Physiol. 2016;85:32–45. https://doi.org/10.1016/j.jinsphys.2015.11.013. PubMed PMID: 26654884
Iwasaki M, Ohno Y, Otaki JM. Butterfly eyespot organiser: in vivo imaging of the prospective focal cells in pupal wing tissues. Sci Rep. 2017;7:40705. https://doi.org/10.1038/srep40705. PubMed PMID: 28094808; PubMed Central PMCID: PMCPMC5240560
Bhardwaj S, Prudic KL, Bear A, Das Gupta M, Cheong WF, Wenk MR, et al. Sex differences in 20-hydroxyecdysone hormone levels control sexual dimorphism in butterfly wing patterns. (in review).
Ibar C, Cataldo VF, Vasquez-Doorman C, Olguin P, Glavic A. Drosophila p53-related protein kinase is required for PI3K/TOR pathway-dependent growth. Development. 2013;140(6):1282–91. https://doi.org/10.1242/dev.086918. PubMed PMID: 23444356
Prudic KL, Stoehr AM, Wasik BR, Monteiro A. Eyespots deflect predator attack increasing fitness and promoting the evolution of phenotypic plasticity. Proc Biol Sci. 2015;282(1798):20141531. https://doi.org/10.1098/rspb.2014.1531. PubMed PMID: 25392465; PubMed Central PMCID: PMCPMC4262162
Cakouros D, Mills K, Denton D, Paterson A, Daish T, Kumar S. dLKR/SDH regulates hormone-mediated histone arginine methylation and transcription of cell death genes. J Cell Biol. 2008;182(3):481–95. https://doi.org/10.1083/jcb.200712169. PubMed PMID: PMC2500134
Monteiro A, Tong XL, Bear A, Liew SF, Bhardwaj S, Wasik BR, et al. Differential expression of Ecdysone receptor leads to variation in phenotypic plasticity across serial Homologs. PLoS Genet. 2015;11(9) https://doi.org/10.1371/journal.pgen.1005529. PubMed PMID: WOS:000362269000043
Monteiro A. Physiology and evolution of wing pattern plasticity in Bicyclus butterflies: a critical review of the literature. In: Sekimura T, Nijhout HF, editors. Diversity and evolution of butterfly wing patterns: an integrative approach: Springer; in press.
Vieira CU, Bonetti AM, Simoes ZL, Maranhao AQ, Costa CS, Costa MC, et al. Farnesoic acid O-methyl transferase (FAMeT) isoforms: conserved traits and gene expression patterns related to caste differentiation in the stingless bee, Melipona Scutellaris. Arch Insect Biochem Physiol. 2008;67(2):97–106. https://doi.org/10.1002/arch.20224. PubMed PMID: 18076110
Everett A, Tong X, Briscoe AD, Monteiro A. Phenotypic plasticity in opsin expression in a butterfly compound eye complements sex role reversal. BMC Evol Biol. 2012;12:232. https://doi.org/10.1186/1471-2148-12-232. PubMed PMID: 23194112; PubMed Central PMCID: PMCPMC3549281
Sokabe T, Tominaga M. A temperature-sensitive TRP ion channel, painless, functions as a noxious heat sensor in fruit flies. Commun Integr Biol. 2009;2(2):170–3. PubMed PMID: 19513273; PubMed Central PMCID: PMCPMC2686375
Saenko SV, Marialva MS, Beldade P. Involvement of the conserved Hox gene Antennapedia in the development and evolution of a novel trait. EvoDevo. 2011;2:9. https://doi.org/10.1186/2041-9139-2-9. PubMed PMID: 21504568; PubMed Central PMCID: PMCPMC3108338
Reed RD, Serfas MS. Butterfly wing pattern evolution is associated with changes in a notch/distal-less temporal pattern formation process. Curr Biol. 2004;14(13):1159–66. https://doi.org/10.1016/j.cub.2004.06.046. PubMed PMID: 15242612
Mateus AR, Marques-Pita M, Oostra V, Lafuente E, Brakefield PM, Zwaan BJ, et al. Adaptive developmental plasticity: compartmentalized responses to environmental cues and to corresponding internal signals provide phenotypic flexibility. BMC Biol. 2014;12:97. https://doi.org/10.1186/s12915-014-0097-x. PubMed PMID: 25413287; PubMed Central PMCID: PMCPMC4275937
Monteiro A, Tong X, Bear A, Liew SF, Bhardwaj S, Wasik BR, et al. Differential expression of Ecdysone receptor leads to variation in phenotypic plasticity across serial Homologs. PLoS Genet. 2015;11(9):e1005529. https://doi.org/10.1371/journal.pgen.1005529. PubMed PMID: 26405828; PubMed Central PMCID: PMCPMC4583414
Monteiro A, Prudic KM. Multiple approaches to study color pattern evolution in butterflies. Trends Evol Biol. 2010;2(1):2. https://doi.org/10.4081/eb.2010.e2.
Koch PB, Merk R, Reinhardt R, Weber P. Localization of ecdysone receptor protein during colour pattern formation in wings of the butterfly Precis coenia (Lepidoptera: Nymphalidae) and co-expression with distal-less protein. Dev Genes Evol. 2003;212(12):571–84. PubMed PMID: WOS:000181090600002
Tong X, Hrycaj S, Podlaha O, Popadic A, Monteiro A. Over-expression of Ultrabithorax alters embryonic body plan and wing patterns in the butterfly Bicyclus anynana. Dev Biol. 2014;394(2):357–66. https://doi.org/10.1016/j.ydbio.2014.08.020. PubMed PMID: 25169193
Ozsu N, Chan QY, Chen B, Gupta MD, Monteiro A. Wingless is a positive regulator of eyespot color patterns in Bicyclus anynana butterflies. Dev Biol. 2017;429(1):177–85. https://doi.org/10.1016/j.ydbio.2017.06.030. PubMed PMID: 28668322
Keys DN, Lewis DL, Goodrich LV, Selegue J, Gates J, Scott MP, et al. Butterfly eyespots and the co-option of A-P patterning genes. Dev Biol. 1997;186(2):A32–A. PubMed PMID: WOS:A1997XH77400137
Tong X, Lindemann A, Monteiro A. Differential involvement of hedgehog signaling in butterfly wing and eyespot development. PLoS One 2012;7(12):e51087. doi: https://doi.org/10.1371/journal.pone.0051087. PubMed PMID: 23227236; PubMed Central PMCID: PMCPMC3515442.
Brakefield PM, Gates J, Keys D, Kesbeke F, Wijngaarden PJ, Monteiro A, et al. Development, plasticity and evolution of butterfly eyespot patterns. Nature. 1996;384(6606):236–42. https://doi.org/10.1038/384236a0. PubMed PMID: 12809139
Vogel C, Marcotte EM. Insights into the regulation of protein abundance from proteomic and transcriptomic analyses. Nat Rev. Genet. 2012;13(4):227–32. https://doi.org/10.1038/nrg3185. PubMed PMID: 22411467; PubMed Central PMCID: PMCPMC3654667
Zhang L, Reed RD. Genome editing in butterflies reveals that spalt promotes and distal-less represses eyespot colour patterns. Nat Commun. 2016;7:11769. https://doi.org/10.1038/ncomms11769. PubMed PMID: 27302525
Dhungel B, Ohno Y, Matayoshi R, Iwasaki M, Taira W, Adhikari K, et al. Distal-less induces elemental color patterns in Junonia butterfly wings. Zoological Lett. 2016;2:4. https://doi.org/10.1186/s40851-016-0040-9. PubMed PMID: 26937287; PubMed Central PMCID: PMCPMC4774158
Monteiro A, Chen B, Ramos DM, Oliver JC, Tong X, Guo M, et al. Distal-less regulates eyespot patterns and melanization in Bicyclus butterflies. J Exp Zool B Mol Dev Evol. 2013;320(5):321–31. https://doi.org/10.1002/jez.b.22503. PubMed PMID: 23633220
Pflugfelder GO, Roth H, Poeck BA. Homology domain shared between drosophila optomotor-blind and mouse Brachyury is involved in DNA binding. Biochem Biophys Res Commun. 1992;186(2):918–25. PubMed PMID: 1497674
Smith JT. Box genes: what they do and how they do it. Trends Genet. 1999;15(4):154–8. PubMed PMID: 10203826
Cook O, Biehs B, Bier E. Brinker and optomotor-blind act coordinately to initiate development of the L5 wing vein primordium in drosophila. Development. 2004;131(9):2113–24. https://doi.org/10.1242/dev.01100. PubMed PMID: 15073155
Stark J, Bonacum J, Remsen J, DeSalle R. The evolution and development of dipteran wing veins: a systematic approach. Annu Rev. Entomol. 1999;44:97–129. https://doi.org/10.1146/annurev.ento.44.1.97. PubMed PMID: 9990717
Charles JP, Iwema T, Epa VC, Takaki K, Rynes J, Jindra M. Ligand-binding properties of a juvenile hormone receptor, Methoprene-tolerant. Proc Natl Acad Sci U S A. 2011;108(52):21128–33. https://doi.org/10.1073/pnas.1116123109. PubMed PMID: 22167806; PubMed Central PMCID: PMCPMC3248530
Dubrovsky EB, Dubrovskaya VA, Bilderback AL, Berger EM. The isolation of two juvenile hormone-inducible genes in Drosophila Melanogaster. Dev Biol. 2000;224(2):486–95. https://doi.org/10.1006/dbio.2000.9800. PubMed PMID: 10926782
Dobryszycki P, Kolodziejczyk R, Krowarsch D, Gapinski J, Ozyhar A, Kochman M. Unfolding and refolding of juvenile hormone binding protein. Biophys J. 2004;86(2):1138–48. https://doi.org/10.1016/S0006-3495(04)74188-0. PubMed PMID: 14747348; PubMed Central PMCID: PMCPMC1303906
Trowell SC. High-affinity juvenile-hormone carrier proteins in the Hemolymph of insects. Comp Biochem Phys B. 1992;103(4):795–807. PubMed PMID: WOS:A1992KC93600004
Nijhout HF, Grunert LW. Bombyxin is a growth factor for wing imaginal disks in Lepidoptera. Proc Natl Acad Sci U S A. 2002;99(24):15446–50. https://doi.org/10.1073/pnas.242548399. PubMed PMID: 12429853; PubMed Central PMCID: PMCPMC137736
Aslam AF, Kiya T, Mita K, Iwami M. Identification of novel bombyxin genes from the genome of the silkmoth Bombyx Mori and analysis of their expression. Zool Sci. 2011;28(8):609–16. https://doi.org/10.2108/zsj.28.609. PubMed PMID: 21801003
Nagasawa H, Kataoka H, Isogai A, Tamura S, Suzuki A, Ishizaki H, et al. Amino-terminal amino acid sequence of the silkworm prothoracicotropic hormone: homology with insulin. Science. 1984;226(4680):1344–5. https://doi.org/10.1126/science.226.4680.1344. PubMed PMID: 17832633
Satake S, Masumura M, Ishizaki H, Nagata K, Kataoka H, Suzuki A, et al. Bombyxin, an insulin-related peptide of insects, reduces the major storage carbohydrates in the silkworm Bombyx Mori. Comp Biochem Physiol B Biochem Mol Biol. 1997;118(2):349–57. PubMed PMID: 9440228
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. https://doi.org/10.1038/nbt.1883. PubMed PMID: 21572440; PubMed Central PMCID: PMCPMC3571712
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170. PubMed PMID: 24695404; PubMed Central PMCID: PMCPMC4103590
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512. https://doi.org/10.1038/nprot.2013.084. PubMed PMID: 23845962; PubMed Central PMCID: PMCPMC3875132
Li B, Dewey CNRSEM. Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. https://doi.org/10.1186/1471-2105-12-323. PubMed PMID: 21816040; PubMed Central PMCID: PMCPMC3163565
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. https://doi.org/10.1093/bioinformatics/btp616. PubMed PMID: 19910308; PubMed Central PMCID: PMCPMC2796818
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3) PubMed PMID: WOS:000277309100013
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6. https://doi.org/10.1093/bioinformatics/bti610. PubMed PMID: 16081474
We thank Laurent Del Basso for training on the command line in an UNIX-like environment and scripting, Heidi Connahs for discussion, and Firefly farms for providing corn for larvae.
This work was supported by the Singapore International Graduate Award (SINGA) fellowship from A*STAR to NÖ and the Singapore Ministry of Education Awards MOE2014-T2-1-146, and MOE R-154-000-602-112 to AM and NUS award R154-000-587-133 to AM.
Availability of data and materials
Raw reads and the transcriptome were submitted to the National Center for Biotechnology Information (NCBI) database under BioProject PRJNA398486 (http://www.ncbi.nlm.nih.gov/bioproject/398486). Raw reads were deposited in the Sequence Read Archive under the accession number SRP116259. The transcriptome was deposited in the Transcriptome Shotgun Assembly (TSA) database under the accession number GFVN00000000 and the version described in this paper is the first version, GFVN01000000.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Hierarchical clustering of similarities in gene expression across the fifteen libraries. Color bars at the top of the figure represent libraries from the following tissue samples: Pink (female wild-type Cu1), blue (wild-type Cu1), green (wild-type M3), yellow (spotty Cu1), red (Spotty M3). Rows represent significantly differentially expressed genes across libraries clustered by their similar expression patterns. Color key indicates up (red) and down (green) regulated genes according to their log fold changes. (TIF 10.2 mb)
P-values and logFC (log2-fold change) of differentially expressed genes associated with eyespot development across all three comparisons. Negative logFC values represents up-regulated genes in eyespots, because the values from tissues without eyespots were divided by those tissues with eyespots. (XLSX 152 kb)
Blast2GO annotations for up-regulated genes associated with eyespots. Table S15. Blast2GO annotations for down-regulated genes associated with eyespots. Table S16. Blast2GO annotations for up-regulated sector-specific transcripts in M3 versus Cu1 wing sector. Table S17. Blast2GO annotations for down-regulated sector-specific transcripts in M3 versus Cu1 wing sector. Table S18. Blast2GO annotations for up-regulated genes in female versus male eyespots. Table S19. Blast2GO annotations for down-regulated genes in female versus male eyespots. Each table is in a different tab. (XLSX 82 kb)
Total transcript raw read counts and TMM values of genes previously associated with butterfly eyespot development across the four types of tissue sample (each with three replicates) used to identify the genes associated with eyespot development (see Fig. 1c). Multiple rows for each gene represent different non-overlapping transcripts obtained for the same gene. (XLSX 45 kb)
Eyespot associated genes in the previously proposed co-opted gene regulatory networks. (DOCX 122 kb)
Up-regulated genes in M3 wing sector in Spotty versus M3 wing sector in wild-type. Table S5. Down-regulated genes in M3 wing sector in Spotty versus M3 wing sector in wild-type. Table S6. Up-regulated genes in Cu1 wing sector in wild-type versus M3 wing sector in wild-type. Table S7. Down-regulated genes in Cu1 wing sector in wild-type versus M3 wing sector in wild-type. Table S8. Up-regulated genes in Cu1 wing sector in Spotty versus M3 wing sector in wild-type. Table S9. Down-regulated genes in Cu1 wing sector in Spotty versus M3 wing sector in wild-type. Table S10. Up-regulated genes in M3 wing sector in Spotty versus Cu1 wing sector in Spotty. Table S11. Down-regulated genes in M3 wing sector in Spotty versus Cu1 wing sector in Spotty. Table S12. Up-regulated genes in Cu1 wing sector in male (wild-type) versus Cu1 wing sector in female (wild-type). Table S13. Down-regulated genes in Cu1 wing sector in male (wild-type) versus Cu1 wing sector in female (wild-type). Each table is in a different tab. (XLSX 790 kb)
Crosses between wild-type and Spotty butterflies. Female Wt individuals were crossed with Spotty male individuals. All F1 (first generation) offspring had phenotypes associated with wt/Spotty heterozygous and displayed small intermediate eyespots in sectors M2 and M3. These F1 individuals were mated with each other to produce a F2 generation. Only wild-type and Spotty homozygous individuals were selected from the F2 generation and these were mated to individuals of the same phenotype to produce two pure-breeding F3 generation cohorts. Female and male wild-type and Spotty pupae from this generation were used for RNA extractions. (TIFF 860 kb)