The dynamic landscape of gene regulation during Bombyx mori oogenesis
BMC Genomics volume 18, Article number: 714 (2017)
Oogenesis in the domestic silkworm (Bombyx mori) is a complex process involving previtellogenesis, vitellogenesis and choriogenesis. During this process, follicles show drastic morphological and physiological changes. However, the genome-wide regulatory profiles of gene expression during oogenesis remain to be determined.
In this study, we obtained time-series transcriptome data and used these data to reveal the dynamic landscape of gene regulation during oogenesis. A total of 1932 genes were identified to be differentially expressed among different stages, most of which occurred during the transition from late vitellogenesis to early choriogenesis. Using weighted gene co-expression network analysis, we identified six stage-specific gene modules that correspond to multiple regulatory pathways. Strikingly, the biosynthesis pathway of the molting hormone 20-hydroxyecdysone (20E) was enriched in one of the modules. Further analysis showed that the ecdysteroid 20-hydroxylase gene (CYP314A1) of steroidgenesis genes was mainly expressed in previtellogenesis and early vitellogenesis. However, the 20E–inactivated genes, particularly the ecdysteroid 26-hydroxylase encoding gene (Cyp18a1), were highly expressed in late vitellogenesis. These distinct expression patterns between 20E synthesis and catabolism-related genes might ensure the rapid decline of the hormone titer at the transition point from vitellogenesis to choriogenesis. In addition, we compared landscapes of gene regulation between silkworm (Lepidoptera) and fruit fly (Diptera) oogeneses. Our results show that there is some consensus in the modules of gene co-expression during oogenesis in these insects.
The data presented in this study provide new insights into the regulatory mechanisms underlying oogenesis in insects with polytrophic meroistic ovaries. The results also provide clues for further investigating the roles of epigenetic reconfiguration and circadian rhythm in insect oogenesis.
In eukaryotes, differentially fine regulation of various pathways at the level of transcription determines correct ontogenesis. Recently, relatively simple processes such as the oogenesis of the domestic silkworm Bombyx mori (Lepidoptera) and that of the fruit fly Drosophila melanogaster (Diptera) have been widely used to reveal the molecular mechanisms of developmental regulation . Consequently, much research progress has been made toward understanding the changes in gene expression and the role of hormones in insect oogenesis [1,2,3].
The domestic silkworm B. mori and the fruit fly D. melanogaster have similar ovary structures, that is, polytrophic meroistic ovaries, which represent one type of insect ovaries. Oogenesis in polytrophic meroistic ovaries can sequentially be subdivided into three developmental periods: previtellogenesis, vitellogenesis and choriogenesis. Previous studies have suggested that the rise in the molting hormone 20-hydroxyecdysone (20E) titer in early pupae could trigger previtellogenic development and vitellogenesis in the domestic silkworm [3, 4]. The cause is the induction, by 20E, of the expression of a conserved cascade of regulatory genes such as nuclear receptors E75, HR3, and Br-C . However, the transition from vitellogenesis to choriogenesis requires another cascade of gene expression that is controlled by a decline in the 20E titer. While Drosophila oogenesis is regulated by both juvenile hormone (JH) and 20E , it is found that the oogenesis of the domestic silkworm depends only on 20E signaling . This finding suggests a conserved role for 20E signaling in the oogenesis of both the domestic silkworm and fruit fly. Thus, the role of 20E in insect oogenesis has become a focus of research. A recent study used genetic mosaic analysis to screen putative 20E–responsive genes for novel roles in the control of the earliest stages of Drosophila oogenesis .
Although the roles of several transcription factors in the oogeneses of the domestic silkworm and fruit fly have been well known, the dynamic landscape of gene regulation at the genome scale during the oogenesis of polytrophic meroistic ovaries remains to be determined. In the domestic silkworm, each ovary consists of 4 ovarioles, and each ovariole can contain up to 75–80 follicles or eggs, which are organized into a single array (Fig. 1a). Each developing follicle is separated by 2–2.5 h of developmental time from its immediate neighbor . This unique feature of the system makes it possible to simultaneously isolate all stages of follicle development from a single animal for physiological, biochemical and gene expression studies. Thus, the developing ovariole of the domestic silkworm provides an excellent model for studies on the changes in gene expression during the execution of long-term developmental programs [3, 8]. To exploit this advantage, in this study we obtained multiple transcriptomes at different time points of ovariole development by next-generation RNA sequencing. The goal of the study was to determine the dynamic landscape of gene regulation during domestic silkworm oogenesis by analyzing time-series transcriptome data.
Drastic morphological changes during oogenesis
In polytrophic meroistic ovarioles, the germarium is the area at the tip of the ovarioles where egg formation is initiated, containing germ line stem cells, somatic stem cells, and their niches . The microscopic structure of the B. mori germarium showed that only one oocyte at the center was enveloped by follicle cells in disordered fashion during previtellogenesis (Time 1 in Fig. 1a). Proceeding to the vitellogenesis stage (Time 3 in Fig. 1a), the follicles were entirely covered by the vitelline membrane. Meanwhile, the interconnected nurse cells began to shrink and gradually degenerated. In contrast, follicle cells close to the oocyte-nurse cell interface commenced to migrate centripetally between the oocyte and the nurse cells, and the oocyte grew and then completely occupied the follicular volume during late vitellogenesis (Time 8 in Fig. 1a). Furthermore, a loose chorion layer formed at the onset of choriogenesis (Time 9 in Fig. 1a). The translucent chorion layer then became dense, and cavities formed between the vitelline membrane and the follicle cells layer during middle choriogenesis (Time 13 in Fig. 1a). During late choriogenesis, the nuclei of follicle cells became sparse and decreased in number during late choriogenesis (Time 15 in Fig. 1a).
Transcriptome sequencing and confirmation of developmental transition
We investigated the profiles of gene expression during follicle maturation in B. mori by analyzing RNA-seq data of transcriptome at various development stages. The overall mapping rates for every sample ranged from approximately 60% to 80% (Additional file 1), and the expression profiles of typical genes validated by qRT-PCR were consistent with the transcriptomic changes of RNA-seq data (Additional file 2). Principal component analysis (PCA) results showed that all samples were divided into two distinct groups (Fig. 1b), separated at the transition from vitellogenesis to choriogenesis. Hierarchical clustering also divided the samples into two clusters (Fig. 1c). In addition, the expression patterns of marker genes showed distinct, stage-specific features. Genes encoding Foxo were exclusively expressed during vitellogenesis, and the members of the chorion family were divided into early, middle and late subgroups during choriogenesis according to the temporal expression patterns of typical chorion genes, like the early chorion genes (ErB.2 and ErB.4), the middle chorion genes (BmCho-1) and the late chorion genes (HcB 12) (Additional file 3). Taken together, five groups could be referred to different oogenesis stages. The first (Time 1 to Time 6) and second groups (Time 7 and Time 8) included the stages from previtellogenesis to middle vitellogenesis and late vitellogenesis. The third (Time 9 and Time 10), fourth (Time 11 and Time 12) and fifth groups (Time 13 to Time 15) represented the early, middle and late choriogenesis stages, respectively. These results are consistent with the morphological and physiological changes that occur during oogenesis.
Drastic changes in gene expression at developmental transition points
To identify the differentially expressed genes during follicle maturation, we compared the gene expression levels between adjacent stages. A total of 1932 unique genes were identified to be differentially expressed between all consecutive stages. Similarly to the observed morphological changes, differential expressions of most genes mainly occurred at the two transition points from previtellogenesis to vitellogenesis and from vitellogenesis to choriogenesis (Fig. 2a). The functional GO terms of those genes were mainly enriched in reproduction, cell differentiation, morphogenesis and transmembrane transport (Fig. 2b), highlighting the typical biological progress of germ cell differentiation. Using KEGG enrichment analysis, we found that differently expressed genes were enriched in the classifications of metabolism of energy and cofactors, and environmental information processing, including ABC transporters and ECM-receptor interaction pathways (Fig. 2c). Additionally, 686 and 542 genes were up-regulated and down-regulated at the transition between late vitellogenesis (Time 8 in Fig. 1a) and early choriogenesis (Time 9 in Fig. 1a), respectively. The down-regulated genes are specifically involved in enzyme regulator activity, isomerase activity, proteinaceous extracellular matrix, and the lipid metabolic process. Conversely, the genes involved in cellular localization, death, cell localization, macromolecule localization and regulation of biological quality were specifically up-regulated (Additional file 4), showing drastic changes in biological progress around the transition from vitellogenesis to choriogenesis.
We also identified six modules corresponding to distinct development stages during oogenesis (Fig. 3). Each module was enriched in specific functional classification and pathway. For example, the genes in Module 1 were exclusively and highly expressed at the onset of oogenesis (Time 1 and Time 2 in Fig. 1a) and significantly enriched in immune system process (Additional file 5). Furthermore, the dorso-ventral axis formation and insect hormone biosynthesis pathways were highlighted as key regulatory pathways (p-value <0.05, FDR < 0.15) (Fig. 4a). In Module 2, the genes were significantly enriched in methyltransferase activity, ncRNA metabolic process and spliceosome pathway (Fig. 4b), suggesting specific regulatory processes. In Module 3, the genes were significantly enriched in biosynthetic process, transport, localization, modification process and oxidative phosphorylation pathway, indicating the main preparation and protein storage phases (Fig. 4c). Module 4 mainly included genes involved in cell adhesion, transmembrane transporter activity and oxidoreductase activity and pathways of fatty acid metabolism, fatty acid biosynthesis and biosynthesis of unsaturated fatty acids (Fig. 4d), suggesting that late vitellogenesis was an important period for energy storage. Moreover, the genes in Module 5 were specifically expressed in early choriogenesis and were enriched in the GO categories of reproduction and cytoskeleton and the fatty acid degradation pathway (Fig. 4e). Module 6 nearly covered the entire choriogenesis period, mainly including the genes involved in the Wnt, Hippo, Hedgehog, TGF-beta signaling pathways as well as autophagy regulation (Fig. 4f); biological processes such as cell communication, signal transduction, stimulus response, autophagy, cell proliferation, cell division, mitotic cell cycling and cell death were also enriched in this module.
To identify the key regulatory factors in each module, we first identified 761 transcription factors (TFs) and 90 cofactor genes by genome-wide screening. Total TFs were classified into 7 superclasses and 60 TF families based on the characteristics of their DNA-binding domains (Additional file 6), including the Forkhead and HMG_box families [10, 11], which are known to play a regulatory role in oogenesis. In Module 1, many key factors’ responses to 20E signaling were observed, for example, Br-c, Eip78C and Notch [12, 13]. For Module 2, SIX4, HMGA and eggless in particular were also highly expressed throughout nearly the entire vitellogenesis period, suggesting that TFs and epigenetic factors might play a joint role in a tight regulatory system. Module 3 included MafK, E74A and c-myc, which have been reported to play a role in D. melanogaster oogenesis . CbZ and FTZ-F1 in Modules 4 and 5, respectively, were previously identified as the main regulatory genes for choriogenesis [15, 16]. The diapause-associated genes timeless, Ets and Tret1 were mainly expressed in Module 6, suggesting that choriogenesis may be the critical period for egg diapause .
Hormone regulation of B. mori oogenesis
Pathway analysis showed that the insect hormone biosynthesis pathways were enriched during the initial stage of oogenesis. Furthermore, nearly every module harbored many TF responses to 20E signaling, suggesting that 20E signaling plays an important regulatory function in oogenesis. To confirm the regulatory function of 20E in oogenesis, an experiment in which exogenous 20E and 20E agonist tebufenozide (RH-5992) were injected into ligated pupa abdomens was performed. The ovaries of individuals treated by being injected with exogenous 20E developed normally, whereas the development of the ovaries in control individuals was arrested (Additional file 7A, B), indicating that 20E is essential for the first transition from previtellogenesis to vitellogenesis. However, although the ovaries in individuals injected with RH-5992 could develop, the development of follicles was soon blocked . Thus, the second transition from vitellogenesis to choriogenesis could not be observed (Fig. 5c). These observations suggested that 20E was the key regulator for the formation of the two transitions during oogenesis. Furthermore, we examined the changes in the 20E titer in female ovarioles during pupal-adult development. As shown in Fig. 3, the 20E titer reached its maximum level during the early stage of vitellogenesis, then decreased drastically before the transition from vitellogenesis to choriogenesis. The changes in the 20E titer were consistent with the expression profiles of the genes involved in 20E synthesis and metabolism (Fig. 5a). The expression levels of ecdysteroidogenic enzyme genes showed distinct scheduling (Fig. 5b). The gene Shd (CYP314A1) encoding ecdysteroid 20-hydroxylase, the key enzyme that transforms ecdysone into 20E, was highly expressed at the initial stage of oogenesis. This finding corresponded with the high 20E titer (Fig. 3). In addition, the expression level of the gene CYP18A1 regulating steroid hormone inactivation specifically increased before the transition point from vitellogenesis to choriogenesis (Fig. 5b, Additional file 2), suggesting that this metabolic enzyme may result in a decrease in 20E during this transition. Furthermore, several intermediate receptors (EcR, DHR and InR) for hormone signaling exhibited stage-specific expression patterns during the early stage of oogenesis (Fig. 3), suggesting that this stage is highly important for the functions of hormone regulatory signaling. Thus, we performed injection experiments involving several inhibitors of the insulin signaling pathway. It was observed that the development of follicles was clearly delayed and the follicles showed abnormalities (Additional file 7), and the expressions of marker genes, like Cyp18a1 and early chorion gene, were up-regulated at choriogenic stages. Probably, insulin signaling may be related to the normal progression of oogenesis.
Consensus oogenesis regulatory network between B. mori and D. melanogaster
Because B. mori and D. melanogaster have similar ovary structures and the data on D. melanogaster oogenesis have been made available, we performed a comparative analysis of the oogenesis regulatory network. Similarly, the genes encoding the york proteins and vitelline membrane proteins were specifically expressed during the early period of oogenesis from stage 1 to stage 10 of D. melanogaster oogenesis, while genes encoding chorion proteins were mainly expressed from stage 11 to stage 14 (Additional file 8A). Thus, D. melanogaster oogenesis was also divided into two distinct phases, previtellogenesis with vitellogenesis and choriogenesis. The shd gene was expressed during previtellogenesis and vitellogenesis before the transition, while the expression level of Cyp18a1 increased during the late stages (from stages 11 to 13, Additional file 8B), suggesting that a low 20E titer was required for the transition from vitellogenesis to choriogenesis. This result is similar to that obtained for B. mori. In addition, we identified four consensus modules that harbor 176 TFs and chromatin remodeling factors. In these modules, 233 genes were known to play an important role in oogenesis based on the mutation phenotypes in D. melanogaster (Additional file 9). Strikingly, the ecdysteroid synthesis-related genes, such as nvd, dib, sro, phm and sad, and the known factors of oogenesis responding to 20E signaling, such as Pxt and Iswi, are enriched in the consensus modules [19, 20]. This indicated that the ecdysteroid synthesis and some hormone signaling pathways in ovaries were conserved in both B. mori and D. melanogaster [19, 21].
Each follicle is a functional unit containing a single oocyte in insects, which is essential for the generation of healthy oocytes. However, the mechanisms that regulate the oogenesis have not been well known. Two typical morphological changes occur during the oogenesis of insects, the accumulation of the yolk protein and the formation of chorion, resulting in two main transitions from previtellogenesis to vitellogenesis and from vitellogenesis to choriogenesis. Our time-series transcriptome data obtained during silkworm oogenesis do reflect that the molecular process of follicle development is much different between vitellogenesis and choriogenesis. Furthermore, it can be inferred from our results that the synthesis of amino acids and ribosomes in vitellogenesis occurs in preparation for macromolecular modification in choriogenesis (Additional file 5). It should be realized that the developing follicles consist of two types of cells: germline cells (oocyte/nurse cells) and somatic cells (follicular epithelium). Due to the small size of follicles, cells of different types were difficult to dissect, respectively. This limited the investigation of separate developmental dynamics of different cell types. However, due to the molecular substance and signal exchange between different cell types, the homogenized developing follicles could be considered as a functional unit. In addition, direct identification of maternal genes was difficult using this data, because of the lack of the cell type information and the change of protein levels. To a certain extent, morphologic observation could supplement the missing information for dynamics of different cell types. For example, the genes expressed continually during choriogenesis were likely to be maternal genes, because at that time the follicle cells gradually experienced apoptosis and oocyte was uniquely reserved in the mature egg . Previous studies have also suggested that maternal mRNAs are mainly polyadenylated during late oogenesis in D. melanogaster and that regulation of poly (A) tails in oocytes shapes the translatomic landscape of embryos . Indeed, we found that Wispy, a typical marker gene for studying maternal factors, was highly expressed during choriogenesis. This result is consistent with the accumulation of maternal mRNAs during B. mori choriogenesis. In addition, we found that several signaling pathways important for embryogenesis were also enriched in choriogenesis (Fig. 4f). This result implies that these molecular processes may occur in preparation for late embryo development. Energy metabolism is generally believed to be important for molecular processes. Our data showed that fatty acid biosynthesis and unsaturated fatty acid biosynthesis mainly occurred during the late stage of vitellogenesis, whereas their metabolism occurred during early stage of choriogenesis (Fig. 4d, e). This finding may reflect that fatty acid metabolism provides the energy needed for the transition from vitellogenesis to choriogenesis.
Recent studies have shown that various hormones may be responsible for the expression regulation of the genes involved in the two transitions during oogenesis in insects. In Tribolium castaneum, Blattella germanica and Locusta migratoria, the initiation of oogenesis depends on JH [24,25,26]. In D. melanogaster, the appropriate balance between JH and 20E is crucial for the normal progress of oogenesis . Our ligation experiments on prepupae suggested that 20E played a central role in B. mori oogenesis. Thus, a high 20E titer is essential for the transition from previtellogenesis to vitellogenesis, whereas the decline in the 20E titer results in the transition from vitellogenesis to choriogenesis. This result is in accord with the 20E titer changes we observed during various stages of oogenesis. Furthermore, we investigated the expression patterns of 20E synthesis- and metabolism-related genes. Our results showed that the stage-specific expressions of these related genes reflected the 20E titer changes at the transition point during oogenesis. In addition, we compared the regulatory network of genes between B. mori oogenesis and D. melanogaster oogenesis because these two insects have similar ovary structures. The source of ecdysone controlling oogenesis in B. mori ovaries is different from that in D. melanogaster ovaries. The ecdysteroid 20-hydroxylaseare-encoding gene (Shd) was expressed at the onset of oogenesis in B. mori. At this stage, ecdysteroid is not completely synthesized in the ovaries, which is induced by 20E [28, 29]. Thus, the ecdysone controlling oogenesis is mainly released from the prothoracic glands, indicating that the ecdysteroids produced by the ovaries have no autocrine/paracrine or endocrine function to regulate the progression of oogenesis in the silkworm. However, in D. melanogaster, the Shd gene was expressed until late vitellogenesis, in accord with sequential expressions of other ecdysteroid synthesis-related genes in follicles. This result implies that the ecdysone produced in ovarian follicular cells, which is stimulated by JH, may play a paracrine or autocrine regulatory role in oogenesis [30,31,32,33]. Although B. mori oogenesis is only related to 20E and D. melanogaster oogenesis is controlled by both 20E and JH, the mechanisms of ecdysteroid synthesis and metabolism in the two insects are similar to each other. For the two insects, increase and decrease in the 20E titer are required for the transitions from previtellogenesis to vitellogenesis and from vitellogenesis to choriogenesis, respectively. This finding suggests that the gene regulatory subnetwork stimulated by 20E in ovaries is evolutionarily conserved. In addition to finding consistency in 20E regulation, the different subnetworks may shed some light on the unique characteristics of hormonal regulation of these two insects oogenesis.
Previous studies put forward a model to address the question of how promoter architecture facilitates developmentally accurate chorion gene regulation, which described both cis-elements and their corresponding transcription factors (TFs) were important for the developmental regulation of chorion genes . The BmC/EBP transcription factor was considered both as an activator and a repressor during choriogenesis in B. mori. BmGATAβ-binding was possibly related to early repression of early and middle chorion gene expression and late activation of late chorion gene pairs. BmHMGA-factor-binding to middle chorion gene promoters recruits other TFs and causes bending of the promoter [1, 35]. The expressions of chorion genes were correlated with the expression of these factors of oogenesis. Furthermore, other factors, like CbZ and FTZ-F1, were expressed around the transition point from vitellogenesis to choriogenesis. Thus, we identified the TFs and cofactors co-expressed with chorion genes and the known regulators in each module. Many zinc finger proteins, like ush and fork head domain transcription factors were in the same module with the known factors (Additional file 10A). The subnetwork of chorion and their co-expressed TFs also gave a hint to find the putative regulators, like Sox15, FOXG1 and HLF (Fig. 6b). In this study, a co-expression analysis of the genes involved in 20E and insulin signaling regulation was also performed. Similar to the previous findings that USP and FoxO have physical interactions , we found that InR, FoxO, USP and HR96 are co-expressed during oogenesis (Fig. 6c). An experiment involving the injection of insulin signaling inhibitors also confirmed that insulin is an important factor for maintaining the normal development of follicles. Although the morphological changes are similar by RH-5992 and inhibitors of insulin signaling pathway, follicles developed slowly and the yolk accumulation decreased, the inhibitory effects were different. RH-5992 blocked the 20E pathways, and the developmental arrest by RH-5992 cannot be rescued by culturing the follicles in vitro. However, the developmental arrest under small doses of inhibitors of insulin signaling pathway can be rescued. Indeed, we observed that only the duration of pupa kept longer while the ovary developed normally (data not shown). Whether insulin signaling pathway is related to oogenesis or whether there are dosage effects of insulin signaling inhibitors on oogenesis remains to be further investigated in future. In addition, recent studies have shown that the central circadian clock genes period and timeless, which do not encode actual TFs but encode regulatory proteins that modulate transcription, are rhythmically expressed in the PG and play a role in controlling development rhythms. Physical interactions among CDK8, CyclinC, EcR, and USP have also been detected . Similar interactions between JH signaling and circadian clock genes that act at the organ-autonomous level in Pyrrhocoris apterus and Aedes aegypti have been well documented [38, 39]. Thus, the gene expression levels of a certain set of genes in the follicles might be transcriptionally regulated under the control of the central clock network, which may elucidate the autonomous regulatory mechanism of oogenesis in B. mori. Our results showed that ecdysteroid biosynthesis enzymes such as Sro, Dib and Sad and the circadian clock components Cyca and Cycb were co-expressed during oogenesis (Fig. 6a). This finding suggests that the autonomous temporal expression of biosynthesis enzymes we observed may be the result of strictly circadian regulation. Certain circadian rhythm genes, such as CRY1 and CycE, are also co-expressed with USP and InR, implying that circadian genes may function downstream of multiple regulatory hormone signaling.
Additionally, histone methylation is likely the most extensively studied epigenetic modification during oogenesis. In D. melanogaster, Piwi negatively regulates Polycomb group proteins mediated by histone 3 lysine 27 trimethylation (H3K27m3) in maintaining germline stem cells and oogenesis. Indeed, our data showed that the histone methylation-related gene eggless was highly expressed during vitellogenesis, whereas histone demethylase (lysine-specific demethylase lid) was highly expressed during choriogenesis (Fig. 3, Additional file 10B), suggesting a molecular switch in the dynamics of histone lysine methylation at the transition from vitellogenesis to choriogenesis. Further dedicated experiments are required to determine the contribution of epigenetic reconfiguration to the regulatory reorganization of the transcriptional program with 20E [40, 41].
In summary, our findings reveal the important role of 20E in inducing the related cascades of molecular processes during B. mori oogenesis and accordingly elucidate the dynamic landscape of gene regulation. These results will expand our understanding of the details of insect oogenesis. Strikingly, the data presented in this study also provide clues for further investigating the roles of epigenetic reconfiguration and circadian rhythm in insect oogenesis.
Experimental sample collection and next-generation RNA sequencing
The polyvoltine strain 305 of the domestic silkworm was reared on fresh mulberry leaves at 25 °C under a 12 h-light:12 h-dark photoperiod. The ovarioles in pupae on day 7 (P7D) were dissected in physiological saline, and the dissected ovarioles were then infiltrated with 90% glycerol for 10 min at room temperature to distinguish the −1/+1 transition point from vitellogenesis to choriogenesis. Every three sequential follicles from 3 individuals were homogenized as one time sample from the transition point to both sides, the last part of ovarioles was discarded and the retained follicles at the onset were pooled as the first sample. In total, 15 samples labeled Time 1 to Time 15, from previtellogenesis to choriogenesis, were obtained (Fig. 1a). Total RNA was extracted from two biological repetitions using TRIzol (Invitrogen, USA) according to the manufacturer’s instructions and treated with DNase I. The mRNA was purified, and used to library construction. The library was sequenced using an Illumina HiSeq™ 2000 platform with a 125 paired-end module at Novogene (Beijing, China).
Assembly of transcriptome and annotation
The raw reads were trimmed and filtered using Trimmomatic software (v0.33) to remove adapters and low-quality bases. The clean reads were then realigned and assembled into transcripts using TopHat (v2.0.13) and Cufflinks (v2.2.1) with the B. mori genome from SilkDBv2.0 as a reference [42, 43]. Cuffmerge and Cuffnorm were used to calculate expression at the gene and transcript levels in fragments per kilobase of exon model per million mapped fragments (FPKM) and reads per kilobase of exon model per million mapped fragments (RPKM). All assembled transcript sequences were annotated by the UniProt databases, NCBI NR databases and Kyoto Encyclopedia of Genes and Genomes (KEGG) database using BLAST with an E-value threshold of 10−5; the Gene Ontology (GO) terms were inferred from the best hit using Blast2GO . The replicates of all samples were analyzed by the Pearson correlation method. Principal component analysis (PCA) and hierarchical clustering were performed using the sample profiles to obtain an overview of the sample relationships in R (v3.3.1).
Identification of transcription factors, cofactors and typical gene families
All protein sequences predicted by TransDecoder were searched against the Pfam database by hmmscan in the HMMER package (v3.1b1), with 10−4 as the E-value . Then, sequences with hidden Markov model (HMM) profiles from the Animal Transcription Factor DataBase (AnimalTFDB) and REGULATOR database were identified as TFs, except for the known receptors, enzymes and proteins associated with histone and transposon [46, 47]. The genes with GO terms related to chromatin remodeling were categorized as cofactors. The typical gene families involved in follicle maturation, including nutrition proteins (30 K and ESP) and structure proteins (chorion and vitelline membrane protein), were also identified.
Differential expression analysis and construction of gene co-expression network
The gene expression profiles of adjacent stages were compared by DEseq2 . The genes that had a two-fold difference in expression level between two stages and a false discovery rate (FDR) less than 0.05 were identified as differentially expressed genes. A co-expression network was constructed using the weighted gene co-expression network analysis (WGCNA) package . Total gene profiles were imported to generate the modules with a soft threshold parameter β = 10, and modules were merged when their eigengene correlation was higher than 0.7. The 20E titers of different follicles were added as phenotypic features to evaluate the modules’ relationships with regulatory hormone. After subjecting the data to strict quality control assessments, the top six modules with good consistency were selected for downstream analysis. The co-expressed interactions among regulators and genes involved in signaling pathways were visualized using Cytoscape (3.2) . The GO function and pathway enrichment analyses were tested using the FDR (Benjamini and Hochberg correction method, FDR < 0.05). Only genes having gene symbols in top terms are shown in Additional file 5.
Consensus network construction of cross-species analysis
The oogenesis transcriptomes from D. melanogaster at different development stages were collected from the NCBI Sequence Read Archive (SRA) and Gene Expression Omnibus (GEO) (Accession: PRJNA257462, PRJNA326885, GSE83616) [23, 51]. We obtain the RPKM expressions of total coding genes in each samples using the same procedure in B. mori. We then detected the consensus modules between D. melanogaster and B. mori oogenesis using a procedure similar to that reported in previous studies . The powers for transforming the co-expression matrices were scaled by topology model fit value. Module preservation statistical tests were used to analyze module preservation across different datasets. We calculated 200 permutations of the preservation statistics and generated a Z-summary value by averaging them. The Z-summary indicates whether a module is strongly preserved (Z-summary score > 10), moderately preserved (5 < Z-summary score < 10), or not preserved (Z-summary score < 5). To visualize module structures, we extracted the genes with the highest module membership and the strongest gene-gene connections among these from the signed topology overlay matrix to visualize the networks.
Real-time quantitative PCR analysis
Ovary samples were dissected from B. mori as described above and stored at −80 °C. A TransZol Up Plus RNA Kit (Transgen, China) was used to isolate total RNA, and 1.5 μg of total RNA free of genomic DNA contamination was used to synthesize cDNA with oligo d (T) primers. The relative expression levels of 12 genes that serve as biomarkers for oogenesis and typical pathways, such as the ecdysteroid biosynthesis pathway and the insulin signaling pathway, were measured using a Bio-Rad CFX System (Bio-Rad, USA) with SYBR Premix Ex Taq II Supermix (Takara, Japan) according to the manufacturer’s instructions. PCR amplification was conducted under the following conditions: 95 °C for 30 s, followed by 40 cycles at 95 °C for 10 s and at 60 °C for 30 s. The gene expression levels were calculated by the 2-△Ct method and normalized against an internal reference gene, BmRpl3 (GeneBank accession: HQ615689). For each sample, three biological replicates were taken based on the average value, and the standard error was calculated. The primers used in this study were designed by Primer Premier 5 and are shown in Additional file 11
Exogenous hormone and inhibitors in vivo injection
We made the ligation in the chest in the prepupa stage and then injected the exogenous 20E and 20E agonist tebufenozide (RH-5992) into the pupa abdomen on day 1 at a concentration of 2 μg/μl. The control groups were injected with 30% alcohol. The ovary was dissected from pupa on day 6. To examine the role of the insulin signaling pathway in follicle development, we used several inhibitors (LY294002, U0126 and Rapamycin, AbMole BioScience, USA) to block the pathway of insulin signaling. LY294002 is a highly selective inhibitor of phosphatidylinositol 3 (PI3) kinase. U0126 is a highly selective inhibitor of both MEK1 and MEK2, a type of MAPK/ERK kinase. Rapamycin is an inhibitor of mTOR. All inhibitors were used in previous researches of the ecdysone biosynthesis in the prothoracic gland [53,54,55]. The inhibitors were injected into the hemolymph of pupae on day 2 and day 4, respectively. The concentration of the injection was twice the dosage used in the prothoracic gland [53,54,55]. The control groups were injected with DMSO. The ovary was dissected to be observed 48 h after injection. Ovarioles of pupa on day 6 were cultured in Grace’s medium (Gibco) at 25 °C for 48 h in the absence or in continuous presence of inhibitors (LY294002, U0126 and Rapamycin), then dissected follicles at choriogenic stages (+1 − +6). The gene expression of typical genes (Cyp18a1 and Era) by qRT-PCR. The method of qRT-PCR was described above.
Paraffin sectioning and 20E titer detection
Follicle samples of B. mori were collected as described above. The 20E titer at every developmental timepoint was detected by B. mori 20E ELISA Kit (Huijia BioScience, China) according to the accompanying manual. The same samples were fixed with 4% paraformaldehyde at 4 °C and then cut using the normal paraffin sectioning method and stained by the typical hematoxylin-eosin method.
false discovery rate
fragments per kilobase of exon model per million mapped fragments
hidden Markov model
Kyoto Encyclopedia of Genes and Genomes
pupae on day 7
principal component analysis
20E agonist tebufenozide
reads per kilobase of exon model per million mapped fragments
weighted gene co-expression network analysis
Papantonis A, Swevers L, Iatrou K. Chorion genes: a landscape of their evolution, structure, and regulation. Annu Rev Entomol. 2015;60:177–94.
Belles X, Piulachs M-D. Ecdysone signalling and ovarian development in insects: from stem cells to ovarian follicle formation. Biochim Biophys Acta. 1849;2015:181–6.
Swevers L, Iatrou K. The ecdysone regulatory cascade and ovarian development in lepidopteran insects: insights from the silkmoth paradigm. Insect Biochem Mol Biol. 2003;33:1285–97.
Swevers L, Eystathioy T, Iatrou K. The orphan nuclear receptors BmE75A and BmE75C of the silkmoth Bombyx Mori: hornmonal control and ovarian expression. Insect Biochem Mol Biol. 2002;32:1643–52.
Schonbaum CP, Perrino JJ, Mahowald AP. Regulation of the vitellogenin receptor during Drosophila Melanogaster oogenesis. Mol Biol Cell. 2000;11:511–21.
Tsuchida K, Nagata M, Suzuki A. Hormonal control of ovarian development in the silkworm,Bombyx mori. Arch Insect Biochem Physiol. 1987;5:167–77.
Ables ET, Hwang GH, Finger DS, Hinnant TD, Drummond-Barbosa D. A genetic mosaic screen reveals Ecdysone-responsive genes regulating Drosophila Oogenesis. G3amp58 GenesGenomesGenetics. 2016;6:2629–42.
Niwa YS, Niwa R. Transcriptional regulation of insect steroid hormone biosynthesis and its role in controlling timing of molting and metamorphosis. Dev Growth Differ. 2016;58:94–105.
Yamauchi H, Yoshitake N. Developmental stages of ovarian follicles of the silkworm,Bombyx Mori L. J Morphol. 1984;179:21–31.
Song J, Li Z, Tong X, Chen C, Chen M, Meng G, et al. Genome-wide identification and characterization of fox genes in the silkworm, Bombyx mori. Funct Integr Genomics. 2015;15:511–22.
Wei L, Cheng D, Li D, Meng M, Peng L, Tang L, et al. Identification and characterization of sox genes in the silkworm, Bombyx mori. Mol Biol Rep. 2011;38:3573–84.
Xu J, Gridley T. Notch signaling during oogenesis in Drosophila Melanogaster. Genet Res Int. 2012;2012:1–10.
Yang C, Lin Y, Liu H, Shen G, Luo J, Zhang H, et al. The broad complex isoform 2 (BrC-Z2) transcriptional factor plays a critical role in vitellogenin transcription in the silkworm Bombyx Mori. Biochim Biophys Acta. 1840;2014:2674–84.
Sun G, Zhu J, Li C, Tu Z, Raikhel AS. Two isoforms of the early E74 gene, an Ets transcription factor homologue, are implicated in the ecdysteroid hierarchy governing vitellogenesis of the mosquito, Aedes aegypti. Mol Cell Endocrinol. 2002;190:147–57.
Sourmeli S, Papantonis A, Lecanidou R. BmCbZ, an insect-specific factor featuring a composite DNA-binding domain, interacts with BmC/EBPγ. Biochem Biophys Res Commun. 2005;338:1957–65.
Lecanidou R, Papantonis A. Silkmoth chorion gene regulation revisited: promoter architecture as a key player. Insect Mol Biol. 2010;19:141–51.
Ikeno T, Tanaka SI, Numata H, Goto SG. Photoperiodic diapause under the control of circadian clock genes in an insect. BMC Biol. 2010;8:116.
Swevers L. The ecdysone agonist tebufenozide (RH-5992) blocks the progression into the ecdysteroid-induced regulatory cascade and arrests silkmoth oogenesis at mid-vitellogenesis. Insect Biochem Mol Biol. 1999;29:955–63.
Tootle TL, Spradling AC. Drosophila Pxt: a cyclooxygenase-like facilitator of follicle maturation. Development. 2008;135:839–47.
Ables ET, Drummond-Barbosa D. The steroid hormone Ecdysone functions with intrinsic chromatin remodeling factors to control female Germline stem cells in drosophila. Cell Stem Cell. 2010;7:581–92.
Beckstead RB, Lam G, Thummel CS. Specific transcriptional responses to juvenile hormone and ecdysone in drosophila. Insect Biochem Mol Biol. 2007;37:570–8.
Carter J-M, Baker SC, Pink R, Carter DR, Collins A, Tomlin J, et al. Unscrambling butterfly oogenesis. BMC Genomics. 2013;14:283.
Lim J, Lee M, Son A, Chang H, Kim VN. mTAIL-seq reveals dynamic poly(a) tail regulation in oocyte-to-embryo development. Genes Dev. 2016;30:1671–82.
Sheng Z, Xu J, Bai H, Zhu F, Palli SR. Juvenile hormone regulates vitellogenin gene expression through insulin-like peptide signaling pathway in the red flour beetle, Tribolium castaneum. J Biol Chem. 2011;286:41924–36.
Comas D, Piulachs MD, Bellés X. Induction of vitellogenin gene transcription in vitro by juvenile hormone in Blattella Germanica. Mol Cell Endocrinol. 2001;183:93–100.
Minakuchi C, Zhou X, Riddiford LM. Krüppel Homolog 1 (Kr-h1) mediates juvenile hormone action during metamorphosis of Drosophila Melanogaster. Mech Dev. 2008;125:91–105.
Dubrovsky EB, Dubrovskaya VA, Berger EM. Juvenile hormone signaling during oogenesis in Drosophila Melanogaster. Insect Biochem Mol Biol. 2002;32:1555–65.
Ohnishi E, Chatani F. Biosynthesis of ecdysone in the isolated abdomen of the silkworm, Bombyx mori. Dev Growth Differ. 1977;19:67–70.
Hanaoka K, Onishi E. Changes in ecdysone titre during pupal-adult development in the silkworm, Bombyx mori. J Insect Physiol. 1974;20:2375–84.
Petryk A, Warren JT, Marqués G, Jarcho MP, Gilbert LI, Kahler J, et al. Shade is the drosophila P450 enzyme that mediates the hydroxylation of ecdysone to the steroid insect molting hormone 20-hydroxyecdysone. Proc Natl Acad Sci U S A. 2003;100:13773–8.
Sonobe H, Yamada R. Ecdysteroids during early embryonic development in silkworm Bombyx Mori: metabolism and functions. Zoolog Sci. 2004;21:503–16.
Ramaswamy SB, Shu S, Park YI, Zeng F. Dynamics of juvenile hormone-mediated gonadotropism in the lepidoptera. Arch Insect Biochem Physiol. 1997;35:539–58.
Swevers L, Iatrou K. Ecdysteroids and ecdysteroid signaling pathways during insect oogenesis. In: Ecdysone, structures and fucntions. Part II in the postgenomic era, Ecdysteroid genetic hierarchies in insect growth and reproduction. Dordrecht: SpringerVerlag; 2009. p. 127–64.
Papantonis A, Lecanidou R. A modified chromatin-immunoprecipitation protocol for silkmoth ovarian follicular cells reveals C/EBP and GATA binding modes on an early chorion gene promoter. Mol Biol Rep. 2009;36:733–6.
Tsatsarounos SP, Rodakis GC, Lecanidou R. Analysis of developmentally regulated chorion gene promoter architecture via electroporation of silk moth follicles: promoter analysis via electroporation. Insect Mol Biol. 2015;24:71–81.
Koyama T, Rodrigues MA, Athanasiadis A, Shingleton AW, Mirth CK. Nutritional control of body size through FoxO-Ultraspiracle mediated ecdysone biosynthesis. Elife. 2014;3:e03091.
Xie X-J, Hsu F-N, Gao X, Xu W, Ni J-Q, Xing Y, et al. CDK8-Cyclin C mediates nutritional regulation of developmental transitions through the ecdysone receptor in drosophila. PLoS Biol. 2015;13:e1002207. Schneider DS, editor
Bajgar A, Jindra M, Dolezel D. Autonomous regulation of the insect gut by circadian genes acting downstream of juvenile hormone signaling. Proc Natl Acad Sci U S A. 2013;110:4416–21.
Shin SW, Zou Z, Saha TT, Raikhel AS. bHLH-PAS heterodimer of methoprene-tolerant and cycle mediates circadian expression of juvenile hormone-induced mosquito genes. Proc Natl Acad Sci. 2012;109:16576–81.
Czech B, Preall JB, McGinn J, Hannon GJ. A transcriptome-wide RNAi screen in the drosophila ovary reveals factors of the germline piRNA pathway. Mol Cell. 2013;50:749–61.
Kronja I, Yuan B, Eichhorn SW, Dzeyk K, Krijgsveld J, Bartel DP, et al. Widespread changes in the posttranscriptional landscape at the drosophila oocyte-to-embryo transition. Cell Rep. 2014;7:1495–508.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78.
Duan J, Li R, Cheng D, Fan W, Zha X, Cheng T, et al. SilkDB v2.0: A platform for silkworm (Bombyx Mori ) genome biology. Nucleic Acids Res. 2010;38:D453–6.
Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinforma Oxf Engl. 2005;21:3674–6.
Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, et al. Pfam: the protein families database. Nucleic Acids Res. 2014;42:D222–30.
Zhang H-M, Liu T, Liu C-J, Song S, Zhang X, Liu W, et al. AnimalTFDB 2.0: A resource for expression, prediction and functional study of animal transcription factors. Nucleic Acids Res. 2015;43:D76–81.
Wang K, Nishida H. REGULATOR: a database of metazoan transcription factors and maternal factors for developmental studies. BMC Bioinformatics. 2015;16:114.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Eichhorn SW, Subtelny AO, Kronja I, Kwasnieski JC, Orr-Weaver TL, Bartel DP. mRNA poly(a)-tail changes specified by deadenylation broadly reshape translation in drosophila oocytes and early embryos. Elife. 2016;5
Talukdar HA, Foroughi Asl H, Jain RK, Ermel R, Ruusalepp A, Franzén O, et al. Cross-tissue regulatory gene networks in coronary artery disease. Cell Syst. 2016;2:196–208.
Gu S-H, Lin J-L, Lin P-L, Chen C-H. Insulin stimulates ecdysteroideogenesis by prothoracic glands in the silkworm, Bombyx mori. Insect Biochem Mol Biol. 2009;39:171–9.
Gu S-H, Young S-C, Lin J-L, Lin P-L. Involvement of PI3K/Akt signaling in PTTH-stimulated ecdysteroidogenesis by prothoracic glands of the silkworm, Bombyx mori. Insect Biochem Mol Biol. 2011;41:197–202.
Gu S-H, Yeh W-L, Young S-C, Lin P-L, Li S. TOR signaling is involved in PTTH-stimulated ecdysteroidogenesis by prothoracic glands in the silkworm, Bombyx mori. Insect Biochem Mol Biol. 2012;42:296–303.
We thank Dr. Yi-Hong Shen at Southwest University, China, for her help in collecting domesticated silkworm samples and Dr. Hirohisa Kishino at the Graduate School of Agricultural and Life Sciences, The University of Tokyo, for discussion on the network construction. We would like to thank all other members of Prof. Zhang’s lab for their laboratory assistance.
This work was supported by the National Natural Science Foundation of China (No. 31272363, 31471197) and the Hi-Tech Research and Development (863) Program of China (2013AA102507–2).
Availability of data and materials
The RAW data we used can be found in the FTP (ftp://126.96.36.199/) (Google browser is a better choice). The statistics of sequencing samples were provided in Additional file 1. We also provided the detailed datasets of the differentially expressed genes (DEGs), co-expressed Modules genes and conModules genes in the website (http://188.8.131.52:3838/silkoodb).
Ethics approval and consent to participate
The sample strain 305 is owned by Sericulture & Agri-food Research Institute, Guangdong Academy of Agriculture Science, and provided by YX. We have got the permissions of the research institution.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The statistics of sequencing samples. (XLS 29 kb)
The quantitative PCR expression results of typical marker genes. The explanation of gene identities is in the Additional file 12. (PDF 444 kb)
The expression profiles for typical protein families during oogenesis in B. mori. (PDF 1351 kb)
The GO classifications of differentially expressed genes between Time 8 and Time 9. The up-regulated and down-regulated genes are marked in red and blue, respectively. The classifications of molecular function, cell component and biological process are represented in different colors. (PDF 600 kb)
The key factors in co-expressed modules during oogenesis in B. mori. (XLS 27 kb)
The statistics of transcription factors. (A) The classification of TFs. (B) TF families with the top 10 numbers. (PDF 301 kb)
The morphological changes of ovaries following hormone and inhibitor treatments. (A-C) The morphological changes observed following the 20E, 30% alcohol and RH-5992 treatments, respectively. (D-G) The morphological changes in the DMSO control and various inhibitors of insulin pathways, LY294002, U0126 and Rapamycin, respectively. (H-I) The quantitative PCR expression results of typical marker genes. The explanation of gene identities is in the Additional file 12. (PDF 14262 kb)
The expression profiles of typical genes during oogenesis in D. melanogaster. (A) The expression profiles for typical protein families during oogenesis in D. melanogaster. (B) The expression profiles for ecdysteroid synthesis and metabolism pathways during oogenesis in D. melanogaster. The explanation of gene identities is in the Additional file 12. (PDF 655 kb)
The key factors in consensus modules between B. mori and D. melanogaster. (XLS 29 kb)
The quantitative PCR primers of typical marker genes. (XLS 28 kb)
The expression profiles for typical TFs and cofactors during oogenesis in B. mori. (A) The expression profiles for typical TFs during oogenesis in B. mori. (B) The expression profiles for typical cofactors during oogenesis in B. mori. (PDF 603 kb)
The explanation of gene identities. (XLS 284 kb)
About this article
Cite this article
Zhang, Q., Sun, W., Sun, BY. et al. The dynamic landscape of gene regulation during Bombyx mori oogenesis. BMC Genomics 18, 714 (2017). https://doi.org/10.1186/s12864-017-4123-6