Unscrambling butterfly oogenesis

Abstract

Background

Butterflies are popular model organisms to study physiological mechanisms underlying variability in oogenesis and egg provisioning in response to environmental conditions. Nothing is known, however, about; the developmental mechanisms governing butterfly oogenesis, how polarity in the oocyte is established, or which particular maternal effect genes regulate early embryogenesis. To gain insights into these developmental mechanisms and to identify the conserved and divergent aspects of butterfly oogenesis, we analysed a de novo ovarian transcriptome of the Speckled Wood butterfly Pararge aegeria (L.), and compared the results with known model organisms such as Drosophila melanogaster and Bombyx mori.

Results

A total of 17306 contigs were annotated, with 30% possibly novel or highly divergent sequences observed. Pararge aegeria females expressed 74.5% of the genes that are known to be essential for D. melanogaster oogenesis. We discuss the genes involved in all aspects of oogenesis, including vitellogenesis and choriogenesis, plus those implicated in hormonal control of oogenesis and transgenerational hormonal effects in great detail. Compared to other insects, a number of significant differences were observed in; the genes involved in stem cell maintenance and differentiation in the germarium, establishment of oocyte polarity, and in several aspects of maternal regulation of zygotic development.

Conclusions

This study provides valuable resources to investigate a number of divergent aspects of butterfly oogenesis requiring further research. In order to fully unscramble butterfly oogenesis, we also now also have the resources to investigate expression patterns of oogenesis genes under a range of environmental conditions, and to establish their function.

Background

Successful development relies heavily on parental contribution over and above the direct effect of maternal and paternal genes. For example, maternal effect genes, which have been particularly well studied in Drosophila melanogaster, are involved in setting up; 1) the location of the germ plasm and subsequent germ cell line development in the offspring [1–3] and, 2) a basic framework of positional information, which is interpreted by the embryo’s own genetic program [4, 5]. Furthermore, insect embryos rely on nutrients for growth derived from the mother in the form of yolk deposited in the egg [6–9]. The investigation of insect egg production (i.e. oogenesis) is thus not only crucial in understanding reproductive, and consequently fitness variation [10–12], it is also a popular model system for studying epigenetic programming [13, 14], the apoptotic pathway [15, 16], stem cell behaviour [17–20], cell cycle regulation [21, 22] and developmental patterning mechanisms in general [4, 5, 23–25].

Research into the physiological mechanisms underlying insect oogenesis and egg provisioning has a rich history [26], particularly in moths and butterflies (Lepidoptera) [7, 8, 27, 28]. However, to date sufficiently detailed developmental genetic data to allow us to comprehensively understand the gene regulatory mechanisms underlying oogenesis and maternal effect gene expression controlling early embryogenesis only really exist for the model organism D. melanogaster [3–5, 15, 21]. Developmental genetic studies focussing on species other than D. melanogaster provide us with the opportunity to investigate how the Gene Regulatory Networks (GRNs) underlying insect oogenesis might have evolved [3–5, 23].

Maternal effects can have consequences that extend well beyond embryonic or juvenile development, affecting offspring fertility and longevity [28, 29]. The exact nature of the maternal effects and thus the contribution of a female to the phenotype (and fitness) of her offspring are not static, however, but to a large extent depend on her own internal state, resource availability [12, 30] and in general the environmental conditions she experienced during her life (both biotic and abiotic) [31–34]. As such maternal effects constitute a form of non-genetic transmission of environmental conditions across generations. This means that elements of the regulatory states from the oogenesis GRN of a mother can be passed on to the next generation. There is thus a developmental framework in place with mothers having the possibility to influence the fecundity and survival of their offspring in response to their own environment, thereby providing an alternative system of inheritance with profound consequences for phenotypic evolution [32, 35–38]. However, much of life history theory has been developed without regard to the actual developmental genetic basis of the variation in the traits being investigated, such as reproductive output and maternal effects [39–41]. What has been lacking is a powerful model system to study the developmental genetics of insect reproduction in an evolutionary ecological context [42]. Lepidoptera are ideal candidates to undertake such ecological evolutionary developmental (eco-evo-devo) studies given the vast amount of physiological data on oogenesis [8], as well as very detailed information, for butterflies in particular, on reproductive variability in relation to environmental variability [10, 11, 43–46].

Recently, valuable functional genomic tools have been developed for butterflies [47]; for example, for Melitaea cinxia to study life history variation [48], Bicyclus anynana to study wing colour patterning [49], the monarch butterfly Danaus plexippus to study long-distance migration [50], Heliconius species to study mimicry [51] and for both Erynnis propertius and Papilio zelicaon to study variability among populations in response to environmental heterogeneity and climate change [52]. The information that has been missing so far in butterflies is a detailed description of the ovarian transcriptome, including maternal regulation of patterning the embryo along its axes and mRNA contributed maternally to eggs. In fact, in Lepidoptera, there is a distinct lack of such developmental studies; only in the silkmoth Bombyx mori have a number of recent studies on candidate genes in maternal regulation of early embryogenesis (e.g. establishing positional information) been undertaken [53, 54].

The Speckled Wood butterfly Pararge aegeria (L.), a temperate zone species, is a popular model species for evolutionary ecology studies, for example on plasticity in female reproduction [10, 11, 55–57]. Female P. aegeria mate soon after emergence and usually mate only once [58]. At eclosion they have no or just a few [56] mature oocytes and if mated on the day of emergence, usually they start ovipositing 48 hrs later on the third day of their life [10, 11]. In female P. aegeria resources for reproduction are, to a significant degree, obtained during the larval stage and there is little opportunity to obtain more nitrogenous resources for reproduction through adult feeding [59] or nuptial gifts. Like many other butterflies [8], P. aegeria has meroistic ovaries with 8 ovarioles. Each ovariole consists of a germarium (i.e. stem cell region), previtellogenic primary oocytes, vitellogenic eggs and mature chorionated eggs [8] (Figure 1). A total of seven nurse cells transfer maternal proteins, and mRNA of maternal effect genes into developing oocytes, whilst the somatic follicle cells surrounding the oocyte are involved in choriogenesis and vitellogenesis, as well as oocyte patterning [8].

Figure 1

Overview ovarian morphology of the Speckled Wood butterfly Pararge aegeria. (A) Female P. aegeria laying an egg. (B) Complete meroistic P. aegeria ovary, consisting of a total of 8 ovarioles. Two times 4 ovarioles are attached to each other in the germarium region. Ovary in photo is still attached to the oviduct and part of the ovipositor. Only the ovaries were used for sequencing in this study. (C) Detail of previtellogenic eggs, with nurse and follicle cells visible.

In this paper, we present a comprehensive study of the genes expressed during oogenesis for the butterfly P. aegeria, using de novo transcriptome sequencing and qPCR. Given the wealth of data on reproductive physiology in Lepidoptera, the genes implicated in hormonal control of reproduction will be investigated in particular detail in this study. Furthermore, as a first step in determining the conserved and divergent elements of the butterfly oogenesis GRN (including maternal regulation of zygotic gene expression and embryonic patterning), we investigated which of the genes known to play an essential role in D. melanogaster or B. mori oogenesis were also transcribed by P. aegeria.

Although the number of ovarioles differs among D. melanogaster, P. aegeria and B. mori, these species have similar organisation of their meroistic ovaries, making for an ideal comparison. Furthermore, within Lepidoptera, the silkmoth B. mori and butterflies (including P. aegeria) belong to the more derived division Ditrysia within the infraorder Heteroneura and thus are likely to share developmental characteristics [60, 61]. Many aspects of maternal regulation of early D. melanogaster embryogenesis can be explained by the fact that it is a long germ band insect [5]. Within the order of Lepidoptera there is a transition from a short germ in the more ancestral species to something more similar to long germ in the more derived species, such as those belonging to Ditrysia [60]. This fact, again, makes for an interesting comparison between the three species.

We describe particular features of the P. aegeria ovarian transcriptome that were revealed during assembly and annotation, including orthologs of genes involved in several major conserved signaling pathways, maternal regulation of early embryogenesis, vitellogenesis and choriogenesis. We observed that P. aegeria differed most significantly from D. melanogaster (and many other insect species) in terms of stem cell maintenance in the germarium, EGF signalling in establishing oocyte polarity along anterior-posterior (AP) and dorsal-ventral (DV), and the signalling mechanisms used at the termini of the oocyte. Furthermore, we observed a high proportion of apparently unique sequences in the transcriptome, and we discuss how future exploration of the function and expression patterns of these unique sequences will undoubtedly provide valuable insights into the evolution of insect oogenesis.

Results

The main aim of this study was to identify the genes expressed in the ovaries involved in oocyte formation, establishing oocyte polarities and the RNA transcripts transferred into the eggs by the mother, which either regulate early embryogenesis or are needed during early embryogenesis. Drosophila melanogaster is arguably the best studied insect species in terms of ovarian gene expression and maternal effect gene function. Additional file 1 contains an extensively referenced list of the key essential oogenesis genes. FlyBase [62] and SilkBase [63] were used as a starting point to conduct the comprehensive literature search. The vast majority of papers thus mainly concern the model species D. melanogaster and B. mori. Furthermore, for D. melanogaster genes, a high-throughput developmental time series database was consulted for FPKM (Fragments Per Kilobase of exon per Million of fragments mapped) -based gene expression levels [64] (see also Methods), as well as an in-situ database for maternal transcript contribution to the oocyte [65]. The oogenesis genes discussed in this paper have been classified into functional groupings and were identified predominantly from D. melanogaster studies (and to a lesser extent B. mori studies). Studies on D. melanogaster oogenesis are too numerous to list exhaustively, but key relevant papers (and references therein) have been cited to enable the reader to explore the role of each particular gene during oogenesis further. It should of course be noted that quite a number of genes are expressed in different functional contexts during oogenesis, such as genes encoding the components of various signalling pathways or a gene such as cornichon, which is involved in setting up both AP and DV axis polarity as well as oocyte nucleus localisation in D. melanogaster [66]. Such genes only occur once in Additional file 1 and the tables presented in this paper, but the references to and discussion of such genes will highlight their pleiotropic functions.

Annotation and verification of expression by means of qPCR

Pararge aegeria egg and ovary RNA was sequenced using Illumina short read RNA-Seq technology. Of the 25266 contigs, 17306 contigs were of sufficient quality and length to be annotated (both automated and manually) with 30%, possibly novel or highly divergent, remaining uncharacterised (Table 1; Additional file 2; see Methods). The presence or absence of P. aegeria orthologs in the transcriptome data of 1035 essential oogenesis genes listed in Additional file 1 was verified manually; 833 were found, which is 80.5%. A total of 994 genes out of the 1035 had been identified in D. melanogaster studies. Pararge aegeria expressed 741 of these, which is 74.5%. A further 56 genes were found to be expressed for which functionality during oogenesis can be inferred, but which have not been verified experimentally. Specific genes will be discussed elsewhere in this paper. A large number of these genes are not only transcribed during oogenesis to produce an oocyte, but maternal transcripts were also found to be present in the oocyte itself (Additional file 2; Figure 2). Exceptions include genes encoding chorion proteins as well as yolk and associated proteins. Large amounts of transcripts of these genes are found in the ovaries only (Additional file 2; Table 2). A number of contigs appeared to have relatively high transcript abundance (measured by means of FPKM values; see Methods) in the oocytes compared to the ovaries, suggesting that these transcripts are important as maternal effect transcripts incorporated into the oocytes in relatively large concentrations (Table 2 and Figure 2). An example of this is the gene encoding a signal transducing adaptor molecule (STAM; Table 2 and Additional file 2), which in D. melanogaster is expressed throughout oogenesis [67], but of which transcripts are detected in very high levels in early embryogenesis [68]. On the basis of the GO terms, the 838 gene orthologs appear to be representative of the annotated genes in the transcriptome as a whole (Figures 2 and 3).

Table 1 Transcript abundance

Figure 2

Gene Ontology manually annotated genes. The presence or absence of orthologs of essential oogenesis genes listed in Additional file 1 has been manually verified. The Gene Ontologies (GO) of genes that were present were determined by BLAST2GO and GO terms were subsequently condensed using the generic GO Slim subset. The histogram details the number of Pararge aegeria manually verified contigs (note, as has been observed for many de novo assemblies, for some genes multiple contigs were present in the transcriptome) for each GO term. FPKM estimates were used to compare the levels of transcripts found in the ovaries and as maternal transcripts in the egg. Using a Log2 fold change threshold of 1, genes were classified in the histogram as present in similar amounts in the egg and ovarian transcriptome (labelled Ubiquitous), used predominantly during oogenesis to make an egg, but not or hardly used as a maternal transcript (labelled Ovary), or highly concentrated in the egg as maternal transcripts (labelled Egg).

Table 2 Sequencing and annotation summary

Figure 3

Gene Ontology total transcriptome. The Gene Ontologies (GO) of succesfully annotated genes in the total transcriptome were determined by BLAST2GO and GO terms were subsequently condensed using the generic GO Slim subset. The histogram details the number of Pararge aegeria contigs (note, for some genes multiple contigs were present in the transcriptome) for each GO term. FPKM estimates were used to compare the levels of transcripts found in the ovaries and as maternal transcripts in the egg. Using a Log2 fold change threshold of 1, genes were classified in the histogram as present in similar amounts in the egg and ovarian transcriptome (labelled Ubiquitous), used predominantly during oogenesis to make an egg, but not or hardly used as a maternal transcript (labelled Ovary), or highly concentrated in the egg as maternal transcripts (labelled Egg).

For of a subset of 17 genes, sampled across the functional groups identified in Additional file 1, the expression in the ovarioles and the presence of transcripts in the oocyte were confirmed further by means of RT-qPCR. These genes were: argonaute 2 (AGO2), caudal (cad), decapentaplegic (dpp), egalitarian (egl), exuperantia (exu), Fragile X mental retardation 1 (Fmr1), nanos-like (nos-like), nanos-M (nos-M), nanos-O (nos-O), ornithine decarboxylase antizyme (Oda), anterior open (aop), par-1, piwi, chorion b-ZIP transcription factor (CbZ), staufen (stau), vitellogenin receptor yolkless (yl; VgR) and vitellogenin (Vtg/Vg). Two further genes, which have not been explicitly studied in the context of oogenesis (references in Additional file 1), were investigated: embryonic lethal abnormal vision (elav) and minibrain (mnb). Furthermore, 3 housekeeping genes were selected to be used as reference genes: RNA polymerase II 215 KD subunit (RPII215), TATA binding protein (Tbp) and zwischenferment (zw, G6PDH) (Additional file 3).

The qPCR results were used to confirm the presence of expression as well as the levels of expression (as indicated by means of FPKM values) in the transcriptome dataset (Figure 4; Additional files 4, 5, and 6). Transcripts of vitellogenin were not transferred into the oocytes and very few dpp transcripts were transferred into the egg (Figure 4). All of the other oogenesis genes investigated by means of qPCR were included as maternal effect gene transcripts in the oocytes (see also Additional file 2). Specific qPCR results will be discussed in the remainder of the paper.

Discussion

Germ-line and ovarian stem cells

In D. melanogaster three major signalling pathways play a significant role in cystoblast differentiation, and the maintenance and division of germ-line and ovarian stem cells; TGF-beta, Wnt and hedgehog signalling [69–71]. Components of all three signalling pathways have been identified for P. aegeria (Table 3 and Additional file 1). However, it is not clear, to what extent these signalling pathways are essential in the Lepidopteran germarium, as they were not identified as such in B. mori using SAGE analyses [72]. Rather than signalling, for example, a previously unidentified non-coding RNA appears to regulate cystoblast differentiation in B. mori [72].

Figure 4

qPCR results. Normalised relative abundance of transcripts for 19 genes of interest. Data above the midline (the median gene expression level set at 1) indicate a relatively high number of transcripts in the oocyte compared with the ovary. Boxes represent the interquartile range. Whiskers represent the minimum and maximum observations. Note Vtg/Vg transcripts were not found in the oocyte.

Table 3 Maintenance and division of germ-line and ovarian somatic stem cells

The TGF-beta ligands glass bottom boat (gbb) and dpp were expressed in P. aegeria ovarioles (qPCR results; Table 3). The type I TGF-beta receptors used were thickveins (tkv) and an activin type 1 receptor similar to baboon (ATR1) (Additional files 1 and 2), the latter of which is present in the D. melanogaster oocyte as a maternal transcript necessary for early embryogenesis [73]. No evidence, however, could be found for an ortholog of activin type I receptor saxophone (sax) (Table 3). No ortholog of the activin type II receptor punt (pnt) was found, although PACG16964 was found to be a type II BMP receptor (Additional file 2). The P. aegeria transcriptome contained orthologs of two SMAD family genes; Mothers against dpp (Mad) and Smad on X (Smox), but not of medea nor of the anti-SMAD Daughters against decapentaplegic (Dad), which have been shown to be of importance in D. melanogaster germline stemcell maintenance [71]. Furthermore, the negative regulator of Dpp signalling dullard (dd) was found to be expressed in P. aegeria ovaries. In D. melanogaster this gene plays a role in wing vein formation [74], and although it has been found to be maternally deposited [65], its role in oogenesis has not been verified. Another negative regulator of Dpp signalling, brinker (brk), which plays a role in eggshell patterning in D. melanogaster [75, 76], was also expressed by P. aegeria. In D. melanogaster, bag of marbles (bam) interacts with Dpp signalling to regulate stem cell maintenance and differentiation in the germarium [77]. However, bam is a Drosophila unique gene and is not found in P. aegeria.

During oogenesis P. aegeria females express two Wnt receptors, which show orthology to frizzled-2 and frizzled-7 (Table 4 and Additional file 1). Furthermore, they express the Wnt receptor l(2)43Ea (boca), which plays a role in D. melanogaster vitellogenesis [78], as well as dishevelled (dsh), which is part of the Wnt receptor complex (Table 3 and Additional files 1 and 2). Other components of the Wnt pathway expressed include armadillo (arm), pangolin (Tcf/LEF), groucho (gro), axin (axn), sugarless (sgl), legless (lgs), pygopus (pygo) and shaggy (sgg; Zw3), as well as wntless (wls)(Table 3 and Additional file 2; references in Additional file 1). Maternal transcripts of each of these genes were found in the oocyte (Table 3; Additional files 1 and 2), with the exception of sgl. Asymmetric localisation of maternal axn RNA has been shown to be involved in AP formation in Tribolium castaneum [79]. Rather interestingly, the ligand wingless (wg) was not found in the assembled transcriptome (Table 3 and Additional file 2). However, 201 ovary and 100 oocyte raw RNA-seq reads mapped against the complete wg CDS from our unpublished P. aegeria genome (approximately between 3.2× and 6.5× coverage, displaying a discontinuous transcript with a number of gaps not covered by reads; Additional file 7). In D. melanogaster, transcripts of wg are not found in the oocyte [65] and although Wnt signaling has been established as present during oogenesis [69], expression levels of wg are extremely low [64], making it hard to detect the transcripts. It is clear that in P. aegeria there is strong maternal contribution to zygotic Wnt signaling (Additional file 2), but whether Wnt signaling plays a role during oogenesis needs to be further investigated.

Table 4 Cytoskeleton and actomyosin contractile ring assembly

No ortholog of Drosophila wnt4 (a vertebrate wnt9 ortholog) was found (Table 3), which in D. melanogaster is involved in regulating cell movement during ovarian morphogenesis [80]. Finally, transcripts of an ortholog of shifted (shf) were present both in the ovary and oocyte in P. aegeria (Table 3 and Additional file 2). This gene encodes an EGF-like protein acting as a Wnt inhibitory factor 1, which in D. melanogaster stabilises hedgehog signalling and transcripts of which are deposited in the oocyte [81]. Hedgehog (hh) itself, as well as components of the pathway including smoothened (smo), fused (fu), Suppressor of fused (Sufu), and cubitus interruptus (ci) were all found to be expressed and maternal transcripts of all were present in the oocyte (Table 3; Additional files 1 and 2). Both costa (cos2) and the receptor patched (ptc) were not expressed during oogenesis by P. aegeria (Table 3; Additional file 1). Although Ptc protein has been detected in the D. melanogaster germarium [70], detecting ptc transcripts may prove more difficult because ptc appears to be transcribed in very low amounts [64], and it is possible that this is why ptc transcripts were also not found in P. aegeria. As has been observed for Wnt signalling, there is a maternal contribution to zygotic Hh signalling, but presently it is not clear whether this signalling pathway plays a significant role during P. aegeria oogenesis.

Cytoskeleton and actomyosin contractile ring assembly

Orthologs of the vast majority of genes that have been described as affecting the cytoskeleton and actomyosin contractile ring during D. melanogaster oogenesis were expressed in P. aegeria (Table 4). One of the genes not found is ovarian tumor (otu), which plays a crucial role during D. melanogaster oogenesis. Otu is involved in cytoskeletal formation, cyst formation in germ-line cells, nurse cell chromosome dispersion and gurken (grk) mRNA localisation [82]. For 14 genes no P. aegeria orthologs could be found in the dataset (Table 4). For a number of these, this is not surprising, as in general it has proven to be difficult to find orthologs outside the genus Drosophila; for example dicephalic (dic), mushroom body defect (mud), hold up (hup) and stand still (still)(references in Additional file 1).

Pararge aegeria females were found to express E-Cadherin (Table 4). E-Cadherin-dependent adhesion underlies the positioning of the oocyte at the posterior of the cyst, which in turn plays a role in establishing the AP polarity in D. melanogaster during very early oogenesis [83].

Oocyte determination (including fusome formation) and formation of the anterior-posterior polarity during the early stages of oogenesis

Three genes have been described in the literature as important in D. melanogaster follicle ring canal formation; visgun (vsg), nasrat (fs(1)N) and scraps (scra)[84, 85]. Only fs(1)N was not transcribed by P. aegeria females (Additional file 1). Fusomes, regions of spectrin-rich cytoplasm, are essential in D. melanogaster to establish a system of directional transport between cystocytes underpinning oocyte determination and subsequent oocyte polarity [86]. The majority of genes that are expressed early in D. melanogaster oogenesis regulating the formation of the fusome (e.g. alpha and beta spectrin and hu-li tai shao) were also transcribed by P. aegeria, as well as the genes involved in establishing initial AP polarity, including par-1 and egalitarian (egl) (Figure 4 qPCR results and Table 5; references in Additional file 1). Par-1 in particular is essential in D. melanogaster for both oocyte determination and for establishing AP polarity through its effects on the organisation of the microtubule cytoskeleton in conjunction with a number of other proteins [87]. Among the proteins with which Par-1 interacts in establishing AP polarity are Bazooka (Baz/Par3), Bicaudal D (BicD), Lkb1/Par4, Egl, 14-3-3epsilon, and Dynein proteins (references in Additional file 1). The genes encoding these proteins were all expressed by P. aegeria (Table 5). Transcripts of both par-1 and egl were also present in the oocyte (Figure 4 qPCR results and Additional file 2).

Table 5 Oocyte determination, fusome and AP polarity

Soon after the posterior localisation of the oocyte in the D. melanogaster cyst, EGF signalling takes place in the posterior between the oocyte (Grk ligand) and the overlying follicle cells (Torpedo receptor) [88, 89], further consolidating AP polarity. Orthologs of the fast-evolving grk are difficult to find outside the genus Drosophila [24]. Two genes encoding EGF ligands and likely to be paralogs of grk, spitz (spi) and keren (krn), are involved in the regulation of border cell migration in D. melanogaster [90]. A single spi/krn-like EGF ligand has been found in the genomes of N. vitripennis and T. castaneum, and has been argued to be functionally similar to grk in DV patterning in these species [24]. Pararge aegeria females expressed an ortholog of this single spi/krn-like EGF ligand, with the sequence displaying significant similarity to Harpegnathos saltator spi (Additional file 2; Table 6). Large amounts of these transcripts were detected in the P. aegeria oocyte (Additional file 2), suggesting a significant role for its use during early embryogenesis as observed in D. melanogaster [65]. Given the expression of a spi/krn in P. aegeria and the significance of EGF signalling in insect oogenesis in general, and establishing oocyte polarity in particular [24], it is very surprising that only weak evidence was found for expression of egfr, the gene encoding the EGF receptor, in P. aegeria ovaries (Table 6). None of the contigs in our de novo assembly could be clearly identified as an egfr transcript. However, 780 raw RNA-seq reads did map against the complete efgr CDS from our unpublished P. aegeria genome (approximately 7.1× coverage, displaying a discontinuous transcript with a number of gaps not covered by reads; Additional file 7). Intriguingly, all of the raw reads that mapped successfully came from the ovariole transcriptome, not the oocyte transcriptome, consistent with the importance of EGF signalling during oogenesis itself. Transcript levels of egfr are low to moderate in D. melanogaster ovaries [64], and thus there is always the possibility, as was suggested for the absence of ptc transcripts in our study, that P. aegeria egfr transcript levels were not high enough to be accurately detected. However, it is intriguing that as for a number of other components of the EGF pathway involved in DV patterning in D. melanogaster, P. aegeria also did not transcribe, for example, rho during oogenesis (Table 6). Spatial restriction dorsally of rhomboid (rho), encoding a ligand-processing protease in the EGFR pathway, is necessary in D. melanogaster both for DV axis formation as well as for correct patterning of the eggshell [89] (further references in Additional file 1). Although further study is required, at present it thus seems that EGF signalling either does not play a significant role in P. aegeria during oogenesis or a highly divergent one. This will be discussed further in the next section.

Table 6 Follicle cell gene expression and border cell migration

Genes acting early in the ovariole to establish dorsal-ventral polarity and genes promoting follicle cell motility such as border cell migration

Quite a number of genes involved in establishing DV polarity in the oocyte are also important for choriogenesis and dorsal appendage formation in D. melanogaster (references in Additional file 1). Apart from aforementioned grk, pipe was also not expressed by P. aegeria. Pipe plays an essential role in establishing DV polarity in D. melanogaster oocytes, with its expression being confined to ventral follicle cells as a result of localised EGF signalling [91]. Recently, however, it has been proposed that pipe is not necessary in a number of insect species studied [4] and even in D. melanogaster there appears to be a second mechanism in establishing DV [92] that may involve delayed induction by graded maternal Dpp signalling in the perivitelline space [93]. Whatever the mechanism employed by Lepidoptera, it is clear from B. mori research that the factors determining DV polarity are associated with the egg cortex [94].

Despite significant differences found in expression patterns of genes involved in EGF signalling in a number of insects, this pathway has been argued to be the ancient mechanism for establishing DV polarity in insect eggs [4]. Transcription factors that have been discussed as mediators of EGF signalling include pointed (pnt), aop and capicua (cic) [91]. Only the latter two were expressed by P. aegeria and present as maternal transcripts, but whether they play a role in establishing DV polarity remains to be investigated (Tables 6 and 7, and Additional file 2; qPCR results). The ETS transcription factor Aop also plays a role in border cell migration and does not receive input exclusively from EGF, but from a number of signalling pathways including Notch [95]. All components of the Notch signalling pathway were expressed in the ovarioles, with only Notch (N) itself not being present as maternal transcripts in the oocyte (Table 6 and Additional file 2). Maternal N transcripts are also not found in D. melanogaster.

Table 7 Dorsal ventral polarity

The Notch pathway interacts with the EGF pathway in establishing oocyte polarity in D. melanogaster, in particular through its effects on follicle cell differentiation at both termini of the oocyte [96]. As has been established in this study, there is only weak evidence at present for the use of the EGF pathway during P. aegeria oogenesis, and it is striking that the iroquois-class homeodomain protein Mirror is not expressed by P. aegeria (Table 7). This protein appears essential in D. melanogaster in integrating EGF and Notch signalling in follicle differentiation and thus establishing AP and DV polarity [97]. Apart from the EGF pathway, Notch interacts with a number of other proteins in patterning the follicle cells surrounding the oocyte, including Toucan and Daughterless (references in Additional file 1). These were expressed by P. aegeria (Table 6), suggesting that the Notch pathway is essential for correct patterning of the follicle cells and possibly oocyte polarity, but in P. aegeria it may not require an interaction with the EGF pathway. Further studies are required to establish whether butterflies have dispensed with EGF signalling and localised pipe expression in establishing oocyte polarity and instead rely on, for example, the Notch and Dpp pathway.

Anterior and posterior system genes

The Lepidopteran Bombyx mori displays features of both short and long germ band type insects, in which orthodenticle (otd) and cad maternal mRNA are localised to establish the embryonic AP-axis [53]. Both were expressed during P. aegeria oogenesis (Table 8) and indeed were present as mRNA in the oocytes (Additional file 2; Figure 4 qPCR results for cad). Bicoid (bcd) is Drosophila-specific and although no ortholog was found to be expressed, the genes that are involved in bcd localisation were, including exu and stau, but not swallow (swa) (Table 8; Figure 4 qPCR results). As observed in D. melanogaster, transcripts for both exu and stau were also present in significant amounts in P. aegeria oocytes (Figure 4 qPCR results; Additional file 2) [65]. The use of bcd in translational repression of cad is unique to Drosophila. It is very likely that the ancestral mechanism for translational repression of cad is by means of the KH-domain containing protein encoded for by mex-3 [98]. Pararge aegeria females expressed an ortholog of mex-3 (Table 8). Furthermore, in D. melanogaster, bcd interacts with genes such as bicoid interacting protein 3 (bin3), eIF4E, larp1, polyA binding protein (pAbp) and AGO2 in order to repress cad translation [99]. All of these were found to be expressed in P. aegeria, and similarly to D. melanogaster [64, 65], present as maternal transcripts in the oocytes (Tables 8 and 9, and Additional file 2; Figure 4 qPCR results for AGO2).

Table 8 Maternal specification of embryonic anterior-posterior axis

Table 9 Maternal specification of embryonic posterior

Drosophila melanogaster includes maternal hunchback (hb) transcripts into the egg, the protein of which will form an AP gradient during early embryogenesis and cooperate with Bcd to specify the anterior of the embryo, whilst being repressed at the posterior by Nos [100]. Although there is variation between insect species as to whether maternal hb RNA or protein is transferred to the egg, as well as in the significance of the maternal contribution to the Hb gradient for AP patterning, the transcription of hb during oogenesis appears conserved [5, 101]. For example, although only zygotic Hb is necessary for AP patterning in the grasshopper Schistocerca americana embryo, maternal hb transcripts appear to be involved in distinguishing embryonic from extra-embryonic cells along the AP axis, whilst in D. melanogaster maternal and zygotic Hb are redundant for AP patterning of the embryo [101]. In B. mori, the hb transcripts detected appear to be transcribed by the zygote, not the mother [53, 101]. Pararge aegeria also did not express hb during oogenesis (Table 8), suggesting that Lepidoptera, or at least Ditrysia, may have dispensed with a maternal contribution to the Hb gradient in the embryo.

Nanos is involved in both the differentiation of the germ plasm and posterior patterning in D. melanogaster [102], although these two functions can be mechanistically uncoupled [103]. Lepidopteran primordial germ cells (PGCs) develop in a midventral position and in the germ disk after blastoderm formation, not posteriorly before the blastoderm is formed as in D. melanogaster [54]. It is therefore unlikely in Lepidoptera that the genes involved in setting up the embryonic posterior will interact with and be dependent on the genes involved in the localisation of germline determinants, as shown to occur in D. melanogaster [54, 60]. Bombyx mori contains a number of nos paralogs (nos-M, -O, -P and –like (also called –N)), which indeed appear to have divided up these functions [54]. Although it has been argued that B. mori does not have a germ plasm, the location of maternal B. mori nos-O transcripts in the embryo seems to correspond with where the PGCs will form [54]. These nos paralogs, with the exception of nos-P, are expressed during oogenesis in both B. mori and P. aegeria, with maternal transcripts detectable in P. aegeria eggs (Figure 4 qPCR results; Additional file 2 and Table 9) [53]. Nanos-P is primarily zygotically expressed during embryogenesis in B. mori and may be implicated in stabilising the embryonic AP-axis [53]. The nos paralogs have also been found in the monarch butterfly (D. plexippus) genome [50] and phylogenetic analysis of nos sequences shows nos-P to be quite different from the other paralogs (Additional file 8), suggesting it may have a different functional role.

Translational repression of D. melanogaster nos RNA is accomplished during oogenesis by proteins encoded by glorund (glo) and in the early embryo by smaug (smg) [104]. Transcripts of both are found in D. melanogaster oocytes [65]. A P. aegeria ortholog of smg was found, which was present as RNA in the oocyte, but not of glo (Table 9 and Additional file 2). Furthermore, Smg protein bound to the nos 3’ UTR recruits the deadenylation complex CCR4-NOT in D. melanogaster [105]. Rapid deadenylation leads to decay of nos RNA, which is essential in establishing the AP gradient of nos RNA [105]. Although it has been argued above that Lepidoptera in all likelihood do not use nos paralogs during oogenesis in establishing the posterior, P. aegeria does express all the genes that encode proteins that form this complex, despite the absence of an obvious ortholog for twin/CCR4 (Table 9). In D. melanogaster it is the germ plasm protein Oskar (Osk) that prevents rapid deadenylation at the posterior pole, establishing nos as a posterior defining gene [105]. Ditrysia appear not to possess an osk ortholog [3], which could be another reason why the identified nos paralogs may not being involved in AP axis formation during oogenesis. Indeed, P. aegeria also does not possess an ortholog of osk (Table 9; unpublished P. aegeria genome).

Germ plasm, polar granules, nuage and p-bodies

Although a germ plasm type structure has been identified cytologically in the moth Pectinophora gossypiella [2], it is not clear whether Lepidoptera possess a proper germ plasm as they lack osk, which has been argued to have been co-opted as the essential gene in germ plasm formation in holometabolous insects [1, 3]. Pararge aegeria may not possess an osk ortholog, but it does express two genes, which in D. melanogaster silence osk translationally during oogenesis; bruno [106] and cup [107] (Table 9 and Additional file 1). It should be noted, however, that these genes are expressed in a number of functional contexts during oogenesis in D. melanogaster (e.g. cell cycle regulation; references in Additional file 1). As part of the germ plasm, Oskar induces polar (or germ) granule formation and in doing so interacts with a number of genes that characterise these polar granules, in particular tudor (tud), vasa (vas) and valois (vls) [3, 103]. Only valois (vls) could not be found in the P. aegeria transcriptome (Tables 9 and 10).

Table 10 Ovarian nuage and piRNA pathway

Both the ovarian nuage, an electron-dense perinuclear structure found predominantly in nurse cells [108], and polar granules are characterised by a number of the same genes, including tud, vas and vls (references in Additional file 1). The nuage appears not only to play a role in protecting germline cells against the expression of selfish genetic elements in the majority of animals, but also in establishing the polar granules in D. melanogaster [108, 109]. It is therefore not surprising that PIWI proteins and their bound PIWI-interacting RNAs (piRNAs) have been identified as important for both nuage and polar granule formation [109, 110]. Many of these genes encode TUDOR-domain containing proteins and seem to evolve rapidly making it difficult to find orthologs outside Drosophila; e.g. vreteno (vret), Brother of Yb (BoYb) and Sister of Yb (SoYb) [110]. Indeed, no orthologs of these genes could be found in the P. aegeria transcriptome (Table 10). Other genes encoding TUDOR-domain containing proteins seem more conserved, such as TDRD1, tejas (TDRD5), TDRD7 and spindle E/homeless (TDRD9) [3, 110] and these were expressed by P. aegeria (Table 10). What is interesting about TDRD7 is that it shares the OST-HTH/LOTUS functional domain with osk [1, 3]. It is likely that this domain is involved in RNA binding and thus for regulating mRNA translation and/or localisation in germ cell development [111].

There are three genes that encode PIWI proteins; piwi, aubergine (aub) and argonaute 3 (AGO3) [112]. All three were expressed during oogenesis by P. aegeria (Figure 4 qPCR results; Tables 1 and 10). Piwi also plays an essential role in the D. melanogaster germarium and is thus involved in the establishment, maintainance and renewal of germline stem cells [113]. Furthermore, mutations in D. melanogaster piRNA (Piwi-interacting RNA) pathway genes often disrupt the axes of the developing oocyte, through their effects on the microtubule cytoskeleton; for example maelstrom (mael), zucchini (zuc) and squash (squ) affect DV polarity [114, 115]. The latter two also interact with aub in D. melanogaster in silencing osk translation during oogenesis [115]. Similarly, the RNAi pathway gene armitage (armi) affects axis formation and is involved in osk translational silencing in D. melanogaster [107]. Neither zuc nor squ was found in the P. aegeria transcriptome, but mael and armi were (Tables 7 and 10).

Ovarian processing bodies (i.e. P-bodies) are aggregates of translationally inactive ribonucleoproteins (RNPs). In D. melanogaster these can be found in nurse cells, but also appear to be involved in compartmentalisation of mRNA decay and translation repression, for example of osk [116, 117]. With the exception of EDC4/Ge-1 and pacman (pcm), genes that encode the essential components of P-bodies were expressed in P. aegeria (described in the context of oogenesis or otherwise, Table 11 and references in Additional file 1). RNA of P-body components, for example Dcp1, are also transferred to oocytes during D. melanogaster oogenesis and are necessary for early embryogenesis [116]. This was also observed in P. aegeria (Additional file 2).

Table 11 Ovarian processing bodies

Once the germ plasm has been established at the posterior in D. melanogaster, a number of (late-acting) maternal-effect genes are essential in germline formation during early embryogenesis ([118]; further references in Additional file 1). Pararge aegeria females do express similar genes to the fruit fly, including genes associated traditionally with D. melanogaster pole plasm, such as arrest/bruno (aret) and imp [119]. However, there are some notable exceptions, the most significant of which are germ cell-less (gcl) and polar granule component (pgc) (Tables 12, and 13, and Additional file 1). These genes are essential in D. melanogaster, but there are no known pgc orthologs outside the genus Drosophila. Although orthologs can be found for gcl even in vertebrates, none can be found in genomic databases for the Lepidoptera, including the new data presented here. The gene wunen (wun) is involved in germ cell migration in D. melanogaster embryos (references in Additional file 1). Pararge aegeria females also include wun transcripts in the oocyte (Table 13 and Additional file 1).

Table 12 Germ plasm formation and germline viability

Table 13 Maternal effect genes

Maternal transcripts involved in regulating early embryogenesis – dorsal-ventral patterning of the embryo and early neurogenesis

Drosophila melanogaster uses an elaborate network of genes to pattern the DV axis during embryogenesis on the basis of the oocyte polarity established during oogenesis (discussed in [89, 120]; further references in Additional file 1). As discussed elsewhere in this paper, the two genes essential for establishing DV polarity in D. melanogaster oocytes, grk and pipe (the latter of which is repressed dorsally [120]), were absent from the P. aegeria transcriptome. The genes that are subsequently involved in establishing the ventral side of the D. melanogaster embryo are co-opted from the Toll innate immune defense pathway (including a serine protease cascade [121]). A similar cascade has been described in T. castaneum, but at present it is not known whether it is restricted to the ventral perivitelline space [4]. This protease cascade and associated (ventral) genes were also expressed in P. aegeria, but at present it is unclear in which functional context they are used. These genes include; windbeutel (wind), nudel (ndl), gastrulation defective (gd), snake (snk), easter (ea), spn27A, spz, tube (tub) and pelle (pll) (Tables 7 and 13; Additional files 1 and 2). No orthologs for the zinc-finger gene weckle (wek) have yet been found outside Drosophila, and wek was also not found in P. aegeria (Table 13). In D. melanogaster, Toll receptor protein accumulates during the embryonic syncytial stage prior to nuclear migration, and is activated ventrally as the result of a serine/protease cascade (references in Additional file 1). The Toll-like receptor expressed by P. aegeria during oogenesis was found to be an ortholog of 18 wheeler (18w), rather than toll (tl) (Tables 6 and 13). In D. melanogaster 18w is involved in dorsal appendage formation and follicle cell migration [122], and DV patterning [89]. While P. aegeria eggs do not have dorsal appendages, 18w may be involved in DV patterning. In D. melanogaster 18w expression in relation to eggshell patterning, and thus DV polarity, is dependent on input from Dpp and EGF signalling pathways [89]. As discussed elsewhere in the paper, there is not much evidence for EGF signalling in P. aegeria oogenesis, but there is for Dpp signalling (e.g. Figure 4 qPCR results). Furthermore, analyses of Toll receptors have shown that B. mori tl and 18w sequences were more similar to each other, than to D. melanogaster toll [123]. It thus remains to be investigated exactly which functional role 18w fulfils during oogenesis in Lepidoptera.

Pararge aegeria did express cactus (cact) and dorsal (dl) (Table 13). Dorsal protein is distributed evenly in a D. melanogaster embryo, but a gradient in the uptake of Dorsal protein into the nucleus (high on the ventral side) is essential for subsequent DV patterning in the D. melanogaster embryo. Dorsal protein activates some genes, whilst repressing others along the DV axis [120, 124]. While there are some differences in detail, the gene regulatory network underlying embryonic DV patterning is largely conserved in all insects [4]. The Dorsal protein represses dpp ventrally and the protein encoded by grainyhead (NTF-1/grh) acts as co-repressor [124]. RNA of grh is deposited maternally into the oocyte to be translated and used ventrally during embryogenesis [124]. Repression of dpp by a Dorsal gradient does not, however, occur in T. casteneum [4]. A high concentration of Dpp will eventually be restricted to the dorsal side of the D. melanogaster embryo and its concentration is further restricted ventro-laterally by Short gastrulation (Sog), which in D. melanogaster may also be maternally provided [120]. Rather interestingly, this antagonistic interaction between Dpp and Sog may already be employed during oogenesis for the establishment of DV polarity in the oocyte [125]. The vrille (vri) gene encodes a Bzip transcription factor that interacts in D. melanogaster with Dpp signalling, acting as dominant maternal enhancers of embryonic DV patterning defects caused by ea and dpp mutations [126]. Two P24 proteins encoded by eclair (eca) and baiser (bai) are essential for the activity of maternal Tkv, a type I Dpp receptor [127]. Pararge aegeria females did transfer maternal transcripts of grh, dpp, tkv, eca, bai and vri into the oocyte, but did not express sog maternally (Figure 4 qPCR results; Tables 3 and 13; Additional files 1 and 2).

Drosophila melanogaster females express a group of genes called the yema genes (yema 2.8, 3.4, 3a, 3b, 3c, 4 and 9.5) during oogenesis, with most of them displaying strict maternal expression. This may be of importance in the development of the central nervous system of the embryo [128]. However, the exact functional roles of the yema genes are not known and there are no orthologs outside Drosophila [128]. No orthologs were found for these genes in the P. aegeria transcriptome (Table 13 and Additional file 1). Pararge aegeria females did, however, express a number of other genes that are implicated in embryonic brain development or in general in the nervous system; e.g. neuralized (neu), elav, brainiac (brn), Fmr1, brain tumor (brat), mnb, and terribly reduced optic lobes (trol) (Tables 3, 6 and 13; Additional file 1). Of these, mnb and elav have not been explicitly studied in the context of oogenesis (references in Additional file 1). Although maternal transcripts of these genes may play a role in embryonic neural development in D. melanogaster, these genes appear to be important in establishing polarity of the oocyte and its differentiation during oogenesis (references in Additional file 1). The expressions of three of these were further investigated by means of qPCR: elav, Fmr1 and the serine/protease encoding mnb (Figure 4 qPCR results). To date, of these three, only Fmr1 has been described as present in D. melanogaster oocytes, but elav, Fmr1 and mnb were all found in P. aegeria oocytes (Figure 4 qPCR results) [129]. Compared to the ovaries, the amount of elav and Fmr1 transcripts in the oocytes was quite low (Figure 4 qPCR results; Additional file 2), suggesting they are important during oogenesis. Whether these genes play a role of significance in establishing oocyte polarity in P. aegeria needs to be investigated.

Terminal genes

The Torso receptor tyrosine kinase (RTK) pathway has been implicated in a number of different processes during D. melanogaster oogenesis, including vitelline membrane (or envelope) biogenesis [130] and in particular terminal region specification [131]. The maternal-effect gene torso (tor) encodes a receptor whose ligand is most probably encoded for by trunk (trk). Furthermore, the protein encoded by torsolike (tsl) plays a role upstream of trk in activating the Tor receptor in a localised manner, and is thought to be essential for terminal specification [132]. Although both tor and tsl are involved in terminal specification in T. castaneum, different tissues are patterned and Torso signalling plays a role in defining the posterior growth zone during embryogenesis in this short germband insect [133]. Torso signalling is by no means the default mechanism for terminal specification, as the honey bee (Apis mellifera) has the gene tsl, but not tor and trk in its genome [134]. The honey bee seems to rely on other mechanisms for terminal specification [135]. Pararge aegeria does not express clear orthologs of either tor or trk during oogenesis, but does express tsl (Table 14). Bombyx mori does have a RTK in its genome (BGIBMGA003976), which shows similarity to torso, as well as to tie-like and Cad96Ca. Pararge aegeria did not express tie-like (Table 6), but did express Cad96Ca (PACG18092; Additional file 2). This transcript was not present in oocytes and was found only in the ovarioles (Additional file 2). Furthermore, a TBLASTN of the putative B. mori tor against the P. aegeria transcriptome showed that transcript PACG7078 (complete CDS; Additional file 2) was similar (E-value= 5.0 E^-50), although it had greater similarity to the receptor tyrosine kinase Fps85D than to tor. This transcript is present in both P. aegeria oocytes and ovarioles, but its role in oogenesis has not been described in the literature. It is clear that P. aegeria uses RTK signalling during oogenesis and that the sequences of its ligands and receptors have diverged from those of other insects. However, at present it is unclear in which functional context RTK signalling takes place.

Table 14 Terminal specification

Chromatin regulation during oogenesis, DNA replication, general transcription and maternal regulation of zygotic transcription in general

In general, the genes that encode proteins involved in chromatin remodelling, DNA replication and transcription are highly conserved across insects and often across the Metazoa in general (references in Additional file 1). A large number of these genes have been studied specifically in the context of oogenesis in D. melanogaster (Table 15; references in Additional 1). Pararge aegeria was found to express orthologs of a number of these genes (Table 15 and Additional file 1). The genes not expressed by P. aegeria seem to either have no clear insect orthologs outside Drosophila, or no such orthologs have been reported in Lepidoptera, such as B. mori. Genes not expressed by P. aegeria, but for which Lepidopteran orthologs exist include TATA box binding protein-related factor 2 (Trf2), sex combs on midleg (scm), and Arginine methyltransferase 1 and 8 (DART1 and DART8, Table 15 and Additional file 1). The gene scm is a member of the polycomb group (PcG) and similar to D. melanogaster polyhomeotic (ph-p) gene. Both play versatile and important roles in D. melanogaster oogenesis, particularly in ovarian follicle formation [136, 137]. Pararge aegeria females did express and transfer orthologs of other PcG genes into the oocyte. These include the polycomb repressive complex 1 (PRC1) genes sex combs extra (sce), polycomb (ph), posterior sex combs (psc), the PRC2 genes extra sex combs (esc), Enhancer of zeste (E(z)) and the polycomb related genes Enhancer of polycomb (E(ph)) and additional sex combs (asx) (Table 15, Additional files 1 and 2; references therein). Recently these genes have also been identified in B. mori embryogenesis [138]. These genes encode proteins that regulate DNA and histone methylation patterns and general chromatin remodelling. However, they also appear to be important specifically during oogenesis and embryogenesis and may be implicated in transferring gene regulatory states from one generation to the next, being regarded as candidate genes in epigenetic processes [139], with possible involvement in transgenerational effects in relation to environmental heterogeneity.

Table 15 Regulation of transcription and chromatin structure

Genes influencing the cell cycle regulators of mitosis and meiosis

A large number of genes that regulate mitosis have been studied in a reproductive context in D. melanogaster. These genes are not only involved in stem cell maintenance and differentiation in the germarium, but also in relation to endocycling in nurse cells and selective amplication of genes (such as chorion genes) important in oocyte production (further references in Additional file 1). As before, the genes that were not expressed by P. aegeria in a mitotic context seemed either to have no clear insect orthologs outside Drosophila, or no such orthologs have been reported in Lepidoptera such as B. mori (Table 16). Among these are dacapo (dap), matrimony (mtrm), microcephalin (MCPH1) and chiffon (chif) (Additional file 1). The full list of genes in Table 16 contains a large number of cyclins, which regulate cyclin dependent kinases (CDKs). Orthologs of two common cyclins could not be found in the P. aegeria transcriptome: cyclin E and J (see the discussion on choriogenesis elsewhere in this paper).

Table 16 Cell cycle tregulation during mitosis and meiosis

The cell cycle becomes arrested in meiotic prophase I in the majority of Metazoans oocytes. This is initiated during the first stages of oogenesis in region 2 of the D. melanogaster germarium [140]. The intriguing fact is that the gene bruno is not only essential in regulating the translation of a number of genes during oocyte differentiation, but it also appears to be involved in regulating the silencing of Cdk1 activity in order to achieve primary meiotic arrest [140]. It should be noted that oocyte AP and DV polarity is established during primary meiotic arrest and only once the oocyte is properly patterned by stage 14 is this arrest broken [140]. As indicated before, bruno was expressed by P. aegeria females (Table 9).

Meiosis during butterfly and moth oogenesis is characterised by the absence of crossing over and the formation of chiasmata [141, 142]. Cytological studies have established that female Lepidoptera may form synaptonemal complexes (SC) in early meiotic prophase I, but no recombination nodules (RN) are formed subsequently. Instead, a structure called elimination chromatin is formed [143]. Usually chiasmata are formed from retained pieces of the SC in which a RN is, or has been, present [144]. The formation of the chiasmata takes place in the cell destined to become the oocyte in the D. melanogaster germarium [140]. Four genes appear essential in D. melanogaster for SC formation and thus possibly chiasmata formation: crossover suppressor on 2 of Manheim (c(2)M); crossover suppressor on 3 of Gowen (c(3)G); corona (cona) and nipped-B (references in Additional file 1). No genes specific for RN alone could be identified on FlyBase [62]. Pararge aegeria females only express nipped-B (Table 16 and Additional file 1), which is involved in a number of cellular processes in D. melanogaster including mitosis [145]. It is also the only one of the four SC genes for which orthologs outside Drosophila can be identified. Rather interestingly, a large proportion of the genes involved in D. melanogaster meiotic chromosome cohesion and segregation also appeared to be Drosophila or Diptera specific and were not identified in the P. aegeria transcriptome. These include grauzone (grau), corona (cona), orientation disrupter (ord) and mei-S332 (Table 16; references in Additional file 1). A number of genes are, however, highly conserved and orthologs have been found in Lepidoptera as males do display crossing-over [141, 142]. These include both mei-W68 and mei-218 but in particular includes the essential meiotic checkpoint gene pch2 (references in Additional file 1). Female P. aegeria did not express any of these genes (Table 16 and Additional file 1). The P. aegeria oogenesis transcriptome described here is thus in accordance with the previous observations made during cytological studies on female Lepidoptera [141–143].

Vitellogenesis and lipid storage

Not only is cell cycle regulation coordinated with oocyte differentiation in D. melanogaster [140], but also with resource provisioning of the oocyte [22]. The gene greatwall (gwl), for example, is both essential in D. melanogaster for maternal provisioning of the egg during vitellogenesis and to ensure secondary meiotic arrest by stage 14 of oogenesis in metaphase I [22]. It is a highly conserved gene in Metazoa and P. aegeria females did express this gene during oogenesis (Table 16 and Additional file 1). Furthermore, gwl (antagonistically) interacts with polo kinase (polo) in mitotic regulation particularly during early embryogenesis, and is maternally provided (references in Additional file 1). Transcripts of both were detected in P. aegeria oocytes (Table 16; Additional files 1 and 2).

Vitellogenesis during insect oogenesis is characterised by the accumulation in the developing oocytes of large lipid transfer proteins (LLTPs; i.e. yolk protein precursors), such as Vitellogenin (Vtg/Vg) and Apolipophorins (ApoLPs) [8, 9]. Predominantly, LLTPs are produced in the fat bodies and secreted into the hemolymph [8, 9], but not all yolk proteins are extraovarian [146]. Follicle cells not only allow extraovarian yolk protein to reach the oocytes, they also produce significant amounts of LLTPs themselves in a number of insect species, including D. melanogaster [146]. Vitellogenic behaviour of follicle cells is under hormonal control [146]. LLTPs are transported into the oocytes via clathrin-dependent endocytosis mediated by the receptors VgR (in D. melanogaster Yolkless, Yl) and LpR [9, 147]. Nurse cells transport yl/VgR RNA into previtellogenic oocytes, thus preparing the oocyte for Vtg uptake [148]. Pararge aegeria females expressed not only Vtg/Vg, apoLp-III, apoLp, their receptors yl/VgR and LpR, but also the genes described in D. melanogaster vitellogenic endocytosis (references in Additional file 1). These genes include clathrin heavy and light chain (chc and clc), sec5, sec6, garnet (G) and jagunal (jagn) (Figure 4 qPCR results; Tables 2 and 17; further references in Additional file 1).

Table 17 Reproductive physiology and vitellogenesis

The major yolk proteins, such as vitellogenins, share sequence similarities with lipases. Although not catalytically active, the vitellogenin region with sequence similarity to lipases is argued to be involved in steroid hormone binding, thus providing a possibility for a direct interaction with the hormones that regulate their production [149]. For example, maternal ecdysteroids are bound as ecdysteroid-phosphates to the Vtg cleaved product Vitellin (Vn) in yolk granules in B. mori and released as ecdysteroids during yolk uptake in the embryo as a result of dephosphorylation by ecdysteroid-phosphate phosphatase (EPPase)[150]. Pararge aegeria did express EPPase (Table 18). Furthermore, a significant component of yolk in a B. mori egg is the ovarian egg-specific protein ESP, a minor yolk protein [151]. The gene encoding ESP is intriguing, as convincing orthologs for minor yolk proteins outside the moths Galleria mellonella (yolk protein/yolk polypeptide 2) and Samia cynthia (ESP) had not been found [149]. More recently, however, a further two sequences with strong sequence similarity to G. mellonella yolk protein 2 have been discovered in D. plexippus and Plodia interpunctella, whilst ESP does show significant sequence similarity with genes encoding the KK-42 binding proteins in Antheraea moth species [152] (Additional file 9). Sharing the same ABhydrolase lipase region, The KK-42 binding proteins and the minor yolk proteins also show strong sequence similarity to lipases identified in species such as D. melanogaster, in particular lipase-1 and 3 (lip-1 and 3) [149]. Lepidoptera may have evolved to use paralogs of these genes in yolk formation. Rather interestingly, although not functioning as a yolk protein, lip-1, but not lip-3, is expressed in vitellogenic follicles in D. melanogaster [149]. An orthologs of lip-1, and possibly lip-3 (very short partial contig), was expressed by P. aegeria, whilst no clear ortholog of a minor yolk protein was found (Table 17; Additional files 2 and 9).

Table 18 Yolk consumption

Among the most highly transcribed genes in P. aegeria ovarioles is an ortholog of the slime mold Physarum polycephalum gene spherulin-2A. No transcripts were found for this gene in eggs (Table 2 and Additional file 2). Lepidopteran orthologs of the protein encoded by this gene have been shown to function as a subunit Yp4 of follicular epithelium yolk protein produced by follicle cells [153].

Yolk is a food source for the developing embryo and a number of genes encoding Cathepsins and Vacuolar Proton ATP-ases are maternally expressed during oogenesis to facilitate yolk uptake in the embryos (references in Additional file 1). Pararge aegeria females were found to express all described yolk uptake genes, with the exception of the acid phosphatase 1 gene (acph-1) (Table 18 and Additional file 1).

Physiology of oogenesis

Reproductive output depends on female nutritional status which not only affects the rate and duration of oogenesis significantly, but also whether previtellogenic egg chambers will enter the vitellogenic stage or apoptose [154]. Two signalling systems are involved; insulin and hormone signalling [155]. In D. melanogaster, for example, absence of the insulin receptor substrate (IRS) Chico precludes vitellogenesis, whilst a sharp increase in 20-hydroxy-ecdysone (20E) relative to juvenile hormone (JH) results in apoptosis of the egg chamber before vitellogenesis is initiated or completed [16, 155]. Although the two signalling systems operate simultaneously and interact, both have been shown to be able to independently terminate egg chamber progression before vitellogenesis takes place in D. melanogaster [155]. Furthermore, the Lepidoptera express a set of unique genes encoding insulin-like peptides, the Bombyxins (Bbx) [156]. The bbx genes are expressed predominantly in the brain, but some may also be expressed in ovaries [156]. Moths, in particular B. mori, possess a large number of bbx-like genes in their genome [156], but the genome of the butterfly D. plexippus appears to have only three such genes [50]. Orthologs of 2 of these 3 (bbxA1-like and bbxA3-like) were transcribed in P. aegeria ovarioles, whilst a third partial IRS transcript showed more sequence similarity to chico than to any bbx-like gene (Table 17 and Additional file 1). The insulin-like receptor (InR) was also expressed by P. aegeria during oogenesis (Table 17 and Additional file 1). Furthermore, P. aegeria expressed a large number of downstream target genes of insulin signalling including genes encoding the serine/threonine protein kinase Akt, the various protein phosphatase 2A subunits (PP2A, e.g. Widerborst) and the lipid storage droplet proteins 1 and 2 (Lsd1 and Lsd2). Please refer to Table 17 and references in Additional file 1 for additional details.

Apart from nutritional status, environmental factors such as temperature can affect hormone concentrations, providing a possibility for environmental control of reproductive output [7, 26]. The interplay between 20E and JH is dynamic and complex, as both 20E and JH also play a role in regulating choriogenesis [157]. Both hormones have a range of pleiotropic effects during oogenesis and their exact developmental role is not only titre related, but also dependent on the dynamic spatio-temporal expression patterns of the receptors and modulators of hormone signalling [157].

There has been extensive investigation of JH signalling [7, 26], but the signal transduction pathway, including the JH receptor, remains poorly understood [158–160]. The most likely candidate gene for the JH receptor proposed to date is the basic helix–loop–helix (bHLH)/Per-Arnt-Sim (PAS) domain gene methoprene-tolerant (met) [158–160]. It may form a homodimer, or possibly may form a JH-dependent transcriptionally active complex with another member of the bHLH-PAS family. The most likely candidate for the complex is the steroid co-activator NCoA-1/p160 FISC, encoded by the gene taiman (tai) in D. melanogaster [158, 160]. The tai gene was originally discovered as a gene that was expressed in follicle cells in the functional context of border cell migration and was described as an ecdysone co-receptor (Table 6; references in Additional file 1). Pararge aegeria females expressed both met and tai (Tables 6 and 17 and S2; contigs for tai PACG7006 and PACG13674 in Additional file 2). An ortholog for tai (UNIPROT: G6DPV9) can also been found in the genome of D. plexippus [50].

Not much is known about which genes are transcriptionally regulated by the JH activated receptor complex [161]. The gene kruppel-homolog 1 (krh1) has been described as a JH response gene, inhibiting 20E induced broad (br) expression in D. melanogaster, but not in the specific context of oogenesis [159]. Both khr1 and br were expressed by P. aegeria females (Additional file 1). Furthermore, JH may either directly or indirectly upregulate ornithine decarboxylase (odc), which regulates polyamine biosynthesis and appears to be essential for vitellogenesis [162]. Both odc and its antagonist gutfeeling (oda), also a mitotic cell-cycle regulator, were expressed in P. aegeria. Maternal transcripts of odc and oda were found in eggs (Figure 4 qPCR results; Table 17, Additional files 1 and 2).

In order to regulate the precise amount of JH in both hemolymph and organs, two sets of enzymes are involved in JH degradation; the JH epoxide hydrolases (JHEHs) and the JH esterases (JHEs) [163]. JHEs function predominantly in the hemolymph and degradation is reversible, whilst JHEHs regulate the amount of JH in organs and degradation is irreversible [163]. Apart from JHEH, five recently discovered JHEH-like protein genes have been characterised in B. mori [163] and in addition to JHEH, P. aegeria expressed orthologs of three of these; jheh-lp1, jheh-lp3 and jheh-lp5 (Table 17 and Additional file 1). With the exception of jheh-lp5, moderate amounts of transcripts of JHEHs were found in the eggs (Additional file 2). The females did not express a clear ortholog of jhe, but did express an ortholog of a gene encoding an intracellular binding protein of JHE presumed to be involved in its transport (JHEbp or DmP29, Drosophila mitochondrial protein 29, Table 17). Significant amounts of maternal JHEbp transcripts were found in P. aegeria eggs (Additional file 2).

Juvenile hormone itself may be bound by JH binding proteins (JHbp) to enable immobilisation, regulate degradation or enable transport [28]. Four complete JHbp CDSs were identified in P. aegeria ovaries; JHbp, cytosolic JHbp (cJHbp), hemolymph JHbp (hJHbp) and a sequence showing strong orthology to takeout (to) identified in D. melanogaster as involved in JH binding (Table 17). Transcripts of both cJHbp and to were transferred to the eggs by P. aegeria (Additional file 2). Given that JH itself can be transferred maternally into eggs in Lepidoptera, it has been argued that JH binding proteins such as cJHbp will protect the developing embryo against the teratogenic effects of any excess JH transferred from the mother [28].

There is a significant amount of life-history variation among insects and consequently in the relative importance of 20E and JH on oogenesis [26], even within Lepidoptera [8]. Lepidoptera have been categorised into four (physiological) groups based on the hormones used to initiate vitellogenesis, choriogenesis and thus the timing of mature egg production [7]. Nymphalids, like P. aegeria, have been argued to best match the criteria for group 4 [7] where JH is the essential gonadotropic hormone. Juvenile hormone in this group is necessary for: a) synthesis of Vtg in the fat body and possibly the ovary (results supporting the latter in this study); b) inducing patency of ovarioles; c) uptake of Vtg by the oocyte (follicle cells deform to facilitate this uptake and this deformation is under JH control) and d) choriogenesis by the follicle cells. Whilst 20E modulates JH signalling in Nymphalids, it plays a more significant role in vitellogenesis and choriogenesis regulation in B. mori and D. melanogaster [7, 146].

Ecdysone signalling, including its target genes, is in general better understood than JH signalling [164]. Bombyx mori appears to be capable of producing ecdysteroids in the ovaries [8], as does D. melanogaster [165]. Drosophila melanogaster expresses start1 during oogenesis in significant amounts in nurse cells, most likely in response to ecdysone signalling. The cholesterol transporter Start1 may in turn facilitate ecdysteroid production from cholesterol-based precursors [165]. Another gene expressed in the nurse cells essential during D. melanogaster cholesterol conversion in the ovaries is defective in the avoidance of repellents (dare), which encodes an Adrenodoxin reductase [166]. Furthermore, in D. melanogaster the SGT1 protein homolog ecdysoneless (ecd) and disembodied (dib) have been described as essential for ecdysone, both for functionality and its production in the ovaries [165, 167]. Maternal transcripts of D. melanogaster start1 are hypothesised to be deposited into the egg to facilitate ecdysteroid signalling in the developing embryo [165]. Rather intriguingly P. aegeria females did not express dib, but did express ecd, start1, and dare. We observed the transfer of transcripts of all three genes into the oocytes (Table 17 and Additional file 2). Start1 has been implicated in ecdysteroid synthesis in the prothoracic gland in B. mori [168]. Further investigation is needed to determine whether ecdysteroids can be produced in P. aegeria ovaries and if the transfer of maternal start1 and dare transcripts is involved in ecdysteroid signalling in early embryos. In common with the majority of insects [8, 157], P. aegeria females did express ecdysone receptor (EcR) and its partner ultraspiracle (usp; labelled chorion factor 1 (cf1) in B. mori) in the ovaries (Table 17). Although JH may be the gonadotropic hormone in P. aegeria, it is clear from the expression results presented here that 20E signalling does play a significant role in vitellogenesis and that there may be maternal regulation of ecdysteroid signalling in early embryos.

Among the so-called early genes in the hierarchy of genes up-regulated in response to activation of EcR in B. mori ovaries are the orphan nuclear receptor genes hr3 and E75(a,b, c and d), the transcription factor gene E74 and the Broad-Complex gene Br-C [151]. The genes encoding the two receptors Hepatocyte nuclear factor 4a and 4B (HNF4A and HNF4B) are up-regulated with a delay in B. mori and their expression increases during vitellogenesis [169]. With the exception of E74, all of these genes were expressed in P. aegeria (Tables 6, 17 and Additional file 1). In B. mori Hr3 regulates the expression of ESP during vitellogenesis, and it regulates the expression of GATAbeta (i.e. transcription factor BCFI) during choriogenesis [151]. As discussed before, P. aegeria females did not express ESP, but did express the related gene lip-3 (Table 17). Furthermore, they also expressed GATAbeta (Table 19 and Additional file 1).

Table 19 Eggshell formation

Vitelline membrane formation and choriogenesis

Vitellogenesis and choriogenesis are carefully coordinated, primarily by hormone signalling. The vitelline membrane (i.e. the inner eggshell layer) is formed halfway through vitellogenesis [170], for which RTK signalling is necessary as discussed elsewhere in this paper. The formation of the vitelline membrane is of significance in maternal regulation of embryonic AP and DV patterning, as some maternal factors become localised in the perivitelline space in D. melanogaster and interact with localised factors inside the oocyte [170]. This also appears to be the case in B. mori [94], although the genes involved remain uncharacterised. As discussed before, Ndl protein (also tellingly called ovarian serine protease in B. mori) is expressed in all follicle cells and is essential for DV patterning of the embryo in D. melanogaster [171]. Ndl is an unusual protein in that not only is its structure reminiscent of an extracellular matrix protein, but that it also has a catalytically active serine/protease domain [171]. As such, it is involved in both vitelline membrane formation as well as acting as the basis of the serine/protease cascade ventrally, essential for the maternally regulated DV patterning of the D. melanogaster embryo [170]. Pararge aegeria females expressed ndl and as in D. melanogaster, no transcripts were found in the oocyte (Table 6 and Additional file 2). It remains to be seen whether Ndl plays a similar dual role in P. aegeria.

Insect vitelline membrane protein (VMP) genes show tremendous sequence diversity. For example, no clear orthologs can be found for D. melanogaster VMP genes outside the genus Drosophila. The best-characterised VMP gene in Lepidoptera is VMP30 [172], for which orthologs can be found in both moths and butterflies and which was also expressed in P. aegeria ovarioles. Once again, no transcripts were found in the oocyte (Table 19 and Additional file 2).

After the follicle cells have secreted proteins to form the vitelline membrane, endocycling takes place in D. melanogaster and clusters of chorion genes are selectively amplified or expressed at very high levels [170, 173]. Perhaps rather surprisingly, P. aegeria did not express an ortholog of G1/S specific cycE, which in D. melanogaster is essential for chorion gene amplification and endocycling in general ([173]; Table 16; further references in Additional file 1). There is a possibility that Lepidoptera do not selectively amplify the chorion genes prior to the onset of choriogenesis, as no evidence was found for this in B. mori [174]. However, nurse cells do become polyploid during B. mori oogenesis [8]. Pararge aegeria females did express the G1/S specific genes cycC and cycD, as well as the S-phase regulators E2f1 and dp (Table 16; further references in Additional file 1).

Choriogenesis as a whole is coordinated by genes such as chorion peroxidase (pxt) in D. melanogaster [170], which was also expressed by P. aegeria (Table 19). Furthermore, apart from aforementioned GATAbeta, a number of specific transcription factors are involved in the critical regulation of the spatio-temporal expression patterns of the various chorion genes in the later stages of oogenesis in Lepidoptera. All chorion genes in B. mori have multiple cis-regulatory binding sites for CCAAT/enhancer binding protein (C/EBP) transcription factors and their expression levels are C/EBP concentration dependent [175]. The D. melanogaster ortholog of C/EBP is slbo, which is also expressed in follicle cells though predominantly involved in border cell migration (references in Additional file 1). High mobility group protein A (HMGA) is essential for B. mori choriogenesis as it induces chorion gene promoter bending and recruits C/EBP and GATAbeta [176]. Pararge aegeria expressed C/EBP (i.e. slbo), its negative regulator tribbles (trbl) and HMGa (Tables 6, 16 and 19), but it is not known in which functional context slbo is used. Another transcription factor for which cis-regulatory binding sites have been identified for chorion genes, in both D. melanogaster and B. mori, is the C₂H₂ zinc finger protein Chorion factor 2 (Cf2) [177]. Furthermore, a chorion-specific b-ZIP transcription factor (CbZ) has been described in B. mori [175] and orthologs can be found in butterfly genomes, such as that of D. plexippus [50]. However, the exact function of CbZ during choriogenesis has not been characterised. Both cf1 and CbZ were transcribed by P. aegeria, with transcripts of the latter rather intriguingly found to be present in the oocyte (Figure 4 qPCR results; Table 19).

Chorion protein (cp) genes evolve possibly even faster than vitelline membrane protein genes [178] and sequence similarity between D. melanogaster cp genes with those identified in Lepidoptera, including P. aegeria, is very low indeed (Table 19; further references in Additional file 1). The infraorder Heteroneura, to which B. mori and butterflies belong, possess unique helicoidal lamellar chorions, which may provide additional strength [61]. Furthermore, the two species for which chorion genes have been characterised and studied in some detail, Lymantria dispar and B. mori, have an extensively derived chorion in which the helicoidal lamellar framework is modified by expansion and densification [61]. Expression patterns of these chorion genes are also dynamically very complex. Gene families in Lepidoptera encoding the structural chorion proteins are characterised by numerous gene duplications, occasional subsequent gene loss, gene conversion, and in general rapid sequence divergence [61, 179]. As a result, determining orthology between individual chorion genes of different species is very difficult and chorion protein phylogenetic trees are characterised by species-specific clusters (i.e. families) of genes [179]. Automatic annotation of butterfly chorion genes in the D. plexippus genome and from our P. aegeria ovarian transcriptome was performed on the basis of the most significant BLAST hit to available moth chorion gene sequences (Additional file 2 and Table 19). It is very doubtful, however, that true orthology has been uncovered in this way, as chorion genes within a species tend to be more similar to each other than to those found in other species. The phylogenetic tree of Lepidopteran chorion genes in Additional file 9 shows distinct clustering between moths and butterflies for each of the chorion gene families. Pararge aegeria chorion genes were highly transcribed during oogenesis (Table 2 and Additional file 1). As well as expressing these chorion gene families, Bombyx mori expresses a gene encoding protein 80 (BmEP80), which forms part of the eggshell and is produced by the follicle cells [180]. BmEP80 is also highly transcribed during P. aegeria oogenesis (Tables 2 and 19; Additional data file 1).

Apoptosis and autophagy

Programmed cell death is an essential process during oogenesis in D. melanogaster and B. mori, with nurse and follicle cells undergoing apoptosis as oogenesis progresses, while complete egg chambers may apoptose in response to environmentally induced hormonal signals such as starvation [15, 16, 154, 181]. Often, apoptosis and autophagy operate synergistally [181] and are to some extent integrated in D. melanogaster ovaries, where the effector caspase Dcp-1 and the inhibitor of apoptosis protein BIR-superfamily domain protein Bruce (also called survivin in B. mori) regulate both autophagy and starvation-induced cell death [182]. Recently, all apoptosis-related genes have been characterised in B. mori, and the results of the study by Zhang and co-workers showed that most of these genes are highly conserved [183]. Furthermore they demonstrated that a number of gene duplications have occurred in the Lepidoptera (e.g. genes ecoding BIR-superfamily domain proteins)[183]. Many of the known genes involved in autophagy and apoptosis have been studied in a reproductive context in D. melanogaster (references in Additional file 1) and the majority of these were expressed during oogenesis by P. aegeria (Table 20). In particular, P. aegeria expressed buffy, three orthologs of bruce (Additional file 2) and the Lepidopteran ortholog of D. melanogaster dcp1, caspase-1 (Table 20).

Table 20 Growth regulation, apoptosis and autophagy

General growth regulators (including the Hippo Pathway)

Hippo is a highly conserved serine-threonine kinase 3-like signalling protein (also called STE20). It is essential for regulating tissue size and growth [184]. Hippo signalling interacts with various other cellular processes in this functional context, including programmed cell death and cell cycling [184]. Hippo signalling is, however, required in a wide variety of developmental contexts, not just tissue growth [184]. In D. melanogaster oogenesis, for example, it is essential for establishing AP polarity in the oocyte as it regulates the expression of the downstream effector of Notch signalling, the gene hindsight/pebbled (hnt), which is required for posterior follicle cell maturation [184]. Orthologs of all the Hippo signalling related genes (i.e. Hippo signalling components, as well as up- and downstream factors) have been identified as being essential in D. melanogaster oogenesis (references in Additional file 1) and were transcribed by P. aegeria, with possibly two exceptions: merlin (mer; ERM2) and mob as tumor suppressor (mats, mob1) (Table 21). Merlin/ERM2 is a member of the band 4.1 protein superfamily and is characterised by a highly conserved FERM (Four.1 protein, Ezrin, Radixin, Moesin) domain involved in crosslinking the cell membrane and the actin cytoskeleton and so is thus important in localising proteins [184]. Pararge aegeria expressed a highly similar gene, ERM1 (Table 9), which in P. aegeria shows a highly significant sequence similarity to ERM2 (Table 9). In D. melanogaster ERM1 is important for Osk localisation [185], but clearly it cannot function in this way in P. aegeria, which lacks Osk. Likewise, P. aegeria appeared to express paralogs that are significantly similar to mob1; mob2 and mob4-like (i.e. preimplantation protein in B. mori) (Table 21). The latter is most likely the Lepidopteran ortholog of D. melanogaster mob1.

Table 21 Growth regulation and Hippo pathway

Heat shock proteins and their control of protein abundance during oogenesis

Heat shock proteins (Hsps) provide a possible mechanism for environmental control of development in ovaries and as maternal effects. The transcription of genes encoding Hsps, or molecular chaperones in general, is not only regulated in response to various environmental factors (e.g. temperature), but is also essential during many developmental processes, including oogenesis. It is thought that Hsps are important for both developmental buffering and differentiation [72, 186](further references in Additional file 1). The functional contexts in which Hsps operate are incredibly varied [186]. In D. melanogaster, for example, Hsp60C is essential in organising and maintaining cytoskeletal and cell adhesion components and thus for establishing AP and DV oocyte polarity [186], whilst Hsp70 affects border cell migration through its effects on the actin cytoskeleton [187]. A large number of genes encoding Hsps and related proteins have been described in a functional context during D. melanogaster oogenesis (references in Additional file 1) and orthologs of all of these were transcribed during P. aegeria ovarioles, often very abundantly (e.g. heat shock protein cognate 3, hsc3) (Tables 2 and 22; Additional file 2).

Table 22 Heat shock proteins

Ribosomal machinery needed for increased ovarian protein synthesis and early embryogenesis

Genes encoding ribosomal proteins, rRNA and other proteins involved in translation (e.g. RpA1) are among the most highly transcribed genes during Metazoan oogenesis, as large amounts of the translation machinery are needed both during oogenesis and by the developing embryo [188]. Just like Hsps, specific ribosomal proteins have been studied in a wide variety of functional contexts during D. melanogaster oogenesis and early embryogenesis (Tables 12 and 18; further references in Additional file 1). Ribosomal genes were also among the most highly transcribed in P. aegeria oogenesis (Table 2; Additional file 2).

Immune defense and Wolbachiainfection

Orthologs of the majority of the genes identified from the literature as being involved in immune response during oogenesis were also found to be expressed by P. aegeria and present as maternal transcripts in the oocytes (Table 23; Additional files 1 and 2). Apart from the aforementioned Toll innate immune defense pathway, which may have been co-opted for DV patterning of the embryo (Table 13), these include a large number of genes encoding Serpins (Table 23). Drosophila melanogaster spn27A (the ortholog of which is called serpin-3 in B. mori), has been implicated in DV axis formation [120].

Table 23 Immune defense

The facultative reproductive parasite Wolbachia sp. is an endocytosymbiont in many arthropod species affecting oogenesis in a multitude of ways and the Bacterium is maternally transmitted [189–191]. In D. mauritiana, Wolbachia increases egg production by affecting the maintenance and division of germ-line stem cells [20], while in the wasp Asobara tabida, Wolbachia confers a reproductive advantage to the females by properly regulating apoptosis during oogenesis via its regulation of iron metabolism and ferritin expression [190, 192]. However, in D. melanogaster highly infected females suffer from a range of oogenesis defects mediated via grk signalling [193]. Pararge aegeria females were also found to be infected with Wolbachia, but how this affects oogenesis in this species is at present not known. However, we did observe that the gene encoding an ortholog of the Ferritin 2 light chain protein (FER2-LCH) was amongst the most highly transcribed genes during P. aegeria oogenesis (Tables 2 and 23), but at present it is unknown whether this effect is due to Wolbachia or whether elevated expression levels are a normal part of female P. aegeria reproduction.

Egg activation, ovulation, gene regulation in oviduct upon mating and maternal effect genes involved in fertilisation

As discussed elsewhere in this paper, after vitellogenesis both the D. melanogaster and the Lepidopteran oocyte are in a secondary meiotic arrest in metaphase I [60, 194]. Unlike in Lepidoptera [60], egg activation in D. melanogaster is not triggered by the act of fertilisation, but due to the mechanical pressure experienced by the oocyte when moving from the ovary into the small and tight oviducts [194]. Egg activation involves eggshell modifications, resumption of meiosis, translation and subsequent degradation of maternal mRNAs, and cytoskeletal changes [194]. A small number of genes have been described as important in D. melanogaster in the latter stages of oogenesis in the general functional context of egg activation (references in Additional file 1). Orthologs for only around half of these were found in the P. aegeria transcriptome (Table 24), which may indicate observed differences in the mechanism of egg activation between the Lepidoptera and D. melanogaster. Among the genes found in the P. aegeria transcriptome is wispy (fs(1)M19/wisp) (Table 24). In D. melanogaster it is a maternal effect gene, encoding a GLD-2 family protein with polynucleotide adenylyltransferase activity and is essential for the oocyte-to-embryo transition [195]. The D. melanogaster Wisp protein is required for poly(A) tail elongation of bcd, toll, and tor transcripts upon egg activation. It is thus important for proper patterning of the embryo [195], but is also required to maintain a high level of active (phospho-) mitogen-activated protein kinases (MAPKs)[195]. Given that P. aegeria females did not express bcd and tor, it remains to be investigated whether wisp is of any importance in patterning of the embryo.

Table 24 Egg activation

Conclusions

A large proportion of the genes currently described in the literature as being essential during insect oogenesis (in particular D. melanogaster oogenesis) were transcribed by P. aegeria and transcripts were transferred to the oocytes. As this was an ovarian transcriptome study, the precise functional context in which these genes were transcribed has not been identified. Differences in the functional context in which particular genes are expressed are to be expected compared to model organisms such as D. melanogaster and even B. mori. What is perhaps more revealing, however, is the absence of certain transcripts in the database, in particular where these transcripts concern paradigms of maternal regulation for various aspects of early insect embryogenesis [3–5, 24]. Pararge aegeria differed most significantly from D. melanogaster (and quite a number of other insect species), both in terms of stem cell maintenance or differentiation in the germarium and in establishing (and maintaining) polarity along AP, DV and at the termini of the oocyte. In particular, although Pararge aegeria females expressed an ortholog of a spi/krn-like EGF ligand and possibly its receptor, many components of the EGF pathway involved in patterning of the axes in D. melanogaster embryos, as well as pipe and mirror, were not expressed. This may either suggest that there is not much evidence for a significant role of EGF signalling in establishing P. aegeria oocyte polarity, or that its functional role and genes involved is divergent from other insects. This requires further study, as well as the functional role and significance of Dpp and Notch signalling in this context.

Although the more derived species such as B. mori within the Ditrysia are argued to be long germ band-like [94], it is more appropriate to describe them as intermediate germ band [53, 54], as they have a very unusual preblastoderm stage. Like D. melanogaster, cleavage in B. mori and the butterfly Pieris rapae is superficial but nuclear migration to the periphery of the oocyte and subsequent cellularisation occurs in an anterior to posterior gradient, after which they display long germ band characteristics [60]. It is very likely that this has a bearing on maternal effect gene expression regulating axes patterning after oocyte polarity has been established during the pre-vitellogenic stages in Ditrysia compared to D. melanogaster, and this could be reflected in the gene expression data presented in this study (e.g. the absence of maternal expression of hb). Although progress has been made in investigating B. mori embryonic patterning [53, 54], how polarity is established during oogenesis in Ditrysia and in the Lepidoptera as a whole is not known. This needs further investigation, and P. aegeria may prove an ideal model these future studies.

Unfortunately, maternal effect gene expression and regulation have received significantly less research attention in Lepidoptera compared to vitellogenesis, choriogenesis and reproductive physiology [8]. This is reflected in the discussion of the results in this paper. Although the latter aspects of oogenesis are well suited to studies of reproductive output under a variety of environmental conditions, many of the genes discussed in this study highlight the interconnectedness of all stages during oogenesis, for example eggshell production and oocyte polarity. Furthermore, key candidate genes that have the potential to play an important role in transgenerational maternal effects have been identified. Among these are genes encoding heat shock proteins and proteins involved in chromatin remodelling.

This study has taken a much-needed first step in determining the conserved and divergent elements of the butterfly oogenesis GRN (including maternal regulation of embryonic patterning) and establishes P. aegeria as an eco-evo-devo model system for the study of butterfly oogenesis. In order to fully unscramble butterfly oogenesis, an investigation of the spatio-temporal expression patterns of the genes discussed in this study, as well as establishment of their function, is required. Further studies are also required to establish the function and expression patterns of the uncharacterised contigs identified in this study, which make up 30% of the total contigs found, and are undoubtedly composed of genes that are of high importance in butterfly oogenesis.

Methods

Butterfly rearing and sample collection

As butterflies were used in this study, no ethical approval was required. Eggs were collected from a large outbred laboratory population of P. aegeria (kept at 300–400 individuals per generation). This population originated from a woodland population from the south of Belgium (St. Hubert; established from 50 eggs) and by the time of the experiment, the butterflies had been reared in the laboratory for 10 generations. Newly hatched larvae were placed on potted host plants (4 larvae per plant) of Poa trivialis L. with access to ad libitum food and were reared until eclosion in a climate room under a regime (24±0.3°C, LD 16:8) that promotes direct development (i.e. no diapause). On the day of eclosion (i.e. day −1, between 9 and 12 h) females from this laboratory stock placed individually in netted cages (0.5 m³) along with a potted P. trivialis plant for oviposition and an artificial flower containing a 10% honey solution [55]. Later the same day (between 13.00 and 16.00 h) a virgin male was introduced to the cage and the mating pair was left undisturbed for 24 h.

Eggs from 50 mated 4-day old females were collected within 20 minutes of being laid, which is well before the onset of cleavage and thus early embryogenesis in butterflies [60]. The eggs were placed immediately in 1ml TRI-Reagent (Sigma-Aldrich, Dorset, UK) and homogenised thoroughly. Furthermore, 2 mated females aged 4 days were sacrificed by severing the nerve cord, after which the abdomen was removed and the ovaries dissected out in ice-cold PBS (1×), with dissection taking no longer than 15 minutes to avoid RNA degradation. The ovaries were pooled and likewise homogenised immediately in 1ml TRI-Reagent.

RNA extraction and quality control

The homogenate (both of eggs and ovarioles/ovary) was first centrifuged at 13000g for 10min primarily to remove the yolk, after which the supernatant was vortexed with 200μl of chloroform. Phases were separated at 13000g for 15min at room temperature. The aqueous phase was removed and precipitated in 0.5ml isopropanol [196]. The RNA samples were further purified using the RNeasy Mini Kit and re-eluted in 30μl nuclease-free water, following the manufacturer’s instructions (Qiagen, Hilden, Germany). Preliminary yield and quality for each RNA extraction were assayed using a Nanodrop, while RNA integrity was verified using the Agilent BioAnalyzer 2100 PicoRNA Chip (Agilent Technologies, Winnersh, UK) (Additional file 10).

De novotranscriptome assembly

Pararge aegeria egg and ovary RNA was sequenced by Source BioScience (Nottingham, UK) using Illumina short read RNA-Seq technology. Both total RNA samples went through polyA selection, fragmentation and double stranded cDNA conversion to produce two separate libraries (300bp insert size) in accordance with the Illumina mRNA-seq library preparation protocol (Illumina, San Diego, USA). Sequencing was performed on the Illumina Genome Analyzer IIx platform with one flowcell lane allocated to each library. A total of 61,400,070 single-reads of 38 base pairs (bp) in length were obtained from the ovary and egg flowcell lanes (31,836,256 and 29,563,814 reads for ovary and egg samples respectively) which were pooled to produce a de novo assembly in CLC Genomics Workbench v4.0 (CLC bio, Aarhus, Denmark) using the default settings for short read data (automatic word and bubble size) [197]. The assembly generated 25266 contigs (Additional file 2) of an average length of 535bp (N50=671bp), 41.06% GC content and an estimated average coverage of 124× per nucleotide.

The RNA-seq data was analysed by FASTQC on the Galaxy platform [198, 199]. Adaptor dimer or overruns in the reads (stretches of sequence matching the library preparation primers/adaptors) were trimmed from both egg and ovary data sets using CLC Genomics Workbench. Furthermore, the sequences were trimmed down to 25 bp from the 5’ end and sequencing artefacts discarded using the FASTX-Toolkit on Galaxy. Subsequently, the trimmed reads were mapped using default parameters against the de novo assembly using TopHat on the Galaxy server [200]. FPKM values were estimated from the TopHat output using Cufflinks [201] with quartile normalisation and multi read correct enabled. The estimates were limited to a reference general feature format file containing locations of the predicted coding regions from the automated annotation if available.

Annotation

The 25,266 contigs generated by the de novo assembly (Additional file 2) were processed through a similarity-based annotation workflow. Open reading frames (ORF) over 200 bp were identified and extracted with the EMBOSS tool “getorf” in Galaxy. The GC content increased to 42.23% when limited to possible coding regions. The predicted ORF and contig sequences were then processed through different BLAST strategies to provide the most suitable annotation possible (Additional files 11 and 12). The alpha group compared the predicted ORF sequences against protein databases to identify complete or highly conserved transcripts. The beta group compared the full contigs against protein databases to identify incomplete or out of frame transcripts. Sequences not identified in the alpha and beta group were compared further against nucleic acid coding sequences (delta) and finally the whole nucleotide database (zeta). Each search strategy was attributed a different rank, ranging from A to I. Identity was inferred based on similarity to the top ranking hit. Similarity scores (SS) were assigned to each hit based on the bitscore (S’), number of positives in each alignment (P) and original contig length (L). Similarity score was calculated using the formula:

$$\mathit{SS}=S^{\prime }\frac{P}{L}$$

Effectively this required hits with higher bitscores to also have good query coverage and positive matches. Any hit attaining an SS below 18 (lower SS threshold) was discarded from each rank, using the next best hit (which may be in a lower rank or group) (Additional file 11). Hits were sorted based on group, positives, rank and SS to determine the top hit that would be used to infer the nature of each sequence. Similarity scores also allowed an initial indication of possible homology; SS above the upper threshold (>/=40) were considered High, those above the lower SS threshold (>/=18) were considered Mild and any others were considered Low. Any hit with a bitscore below 40 was excluded from inferring any possible identity or homology (Additional files 12 and 13).

The output from the automated annotation was checked manually for any errors (Additional file 2). Furthermore, using FlyBase [62] and SilkBase [63] as a starting point, a comprehensive literature search was conducted to identify those genes that have been studied in the context of insect oogenesis and maternal regulation of early embryogenesis (1035 genes, of which 994 have been studied in D. melanogaster; fully referenced in Additional file 1). For a further 56 genes functionality during oogenesis can be inferred, but their expression during oogenesis has not always been verified experimentally. The presence or absence of orthologous P. aegeria transcripts in both the oocyte and the ovarioles was verified for each of the 1091 genes and these transcripts were further annotated manually (indicated as such in Additional file 2).

The final BLAST results (1 top hit per sequence) used for annotation, including those genes annotated manually, were used as input in the BLAST2GO software [202] and assigned with Gene Ontology (GO) terms where possible. To help provide an overview of the GO based on the BLAST results, the GO terms were condensed using the generic GO Slim subset.

Transcript abundance and qPCR of genes involved in oogenesis and maternal regulation of early embryogenesis

For of a subset of 19 genes the expression in the ovarioles and the presence of transcripts in the oocyte were confirmed further by means of RT-qPCR (Additional file 3). For both ovary and oocyte, cDNA was generated from 500 – 1000 ng of RNA using the Verso RT Kit (Thermo Fisher, Surrey, UK). The reverse transcriptions were primed by a 3:1 mix of random hexamers:oligo-dT taking place in 20μl total volume reactions at 42°C for 30 min after an initial 5 min denaturation step at 70°C. Negative reverse transcription (NRT) controls were run in parallel without both Verso RT enzyme mix and primers. A final heat deactivation at 95°C for 2 min was also implemented to deactivate the RT enhancer. The resulting cDNA was stored at −20°C.

For the qPCR stage, suitable primer pairs were selected automatically using the online Primer3+ primer design service and tested in-silico via the Integrated DNA Technologies online structure prediction package (Oligo Analyzer). Only those primers exhibiting the best stability were selected. Each primer pair was tested on a 3-step 5-fold dilution series of the ovary cDNA in triplicate, which enabled the primer pair efficiencies to be determined using the CFX Manager software (Bio-Rad Laboratories, California, USA). Primers with adequate efficiency (>65%) were then used for investigating the transcript abundance in the egg and ovary cDNA (Additional file 3).

All qPCR runs were performed on the CFX96 Real-Time PCR Detection System (Bio-Rad) on white 96-well plates in ABsolute Blue qPCR SYBR Green Mastermix (Thermo Fisher, Surrey, UK) with the recommended amount of ROX reference dye (Additional file 14). Test samples were measured in triplicate, while no template controls (NTC) and NRTs were present in duplicate on each plate. The CFX96 data generated was recorded by the CFX manager program using automatic threshold determination. The quantification cycle (Cq) values are listed in Additional file 4.

Relative transcript abundance (i.e. ovary versus egg) was used to reveal whether any individual transcript was used as a maternal effect gene transcript or was merely necessary for oocyte production. Relative transcript abundance in the ovaries and eggs were obtained using the relative expression software tool REST v2.0.13.0 software package [203], which used the 3 available reference genes to normalise the measurements obtained from the egg and ovary derived cDNA (Additional file 5).

The number of reads mapping to a transcript of a particular gene in RNA-seq data was argued to be correlated linearly with the number of transcripts of that gene [204]. Rather than using read counts, it is considered to be more appropriate to use a corrected relative value, taking transcript length and total number of mapped reads into account [204]. Cufflinks generated such corrected values, the FPKM values, which can be used for the reliable determination of transcript abundance for each of the genes discussed in this study (Additional file 2). In fact, for the 22 genes in the P. aegeria transcriptome investigated by means of qPCR, transcript abundance calculated on the basis of Cq values by means of the methods described in [205] showed significant positive correlation with FPKM values in the combined oocyte and ovary transcriptome (Pearson regression, with null hypothesis that correlation is >0: t₄₁ = 2.37, P = 0.011; Additional file 6).

Annotated contigs and accession numbers of raw data

The sequence read data reported in this manuscript have been deposited in the NCBI Sequence Read Archive and are available under the accession numbers SRR771147 (ovarian reads) and SRR772253 (oocyte reads). Additional file 15 provides the fasta format sequences of the assembled contigs, including the suggested annotated names (top BLAST results as well as information on the manual annotation listed in Additional file 2). Additional file 2 provides information on the start and end of the coding regions in the contigs.