20-Hydroxyecdysone-responsive microRNAs of insects

ABSTRACT

20-Hydroxyecdysone (20-HE) plays essential roles in coordinating developmental transitions of insects through responsive protein-coding genes and microRNAs (miRNAs). The involvement of single miRNAs in the ecdysone-signalling pathways has been extensively explored, but the interplay between ecdysone and the majority of miRNAs still remains largely unknown. Here, by small RNA sequencing, we systematically investigated the genome-wide responses of miRNAs to 20-HE in the embryogenic cell lines of Bombyx mori and Drosophila melanogaster. Over 60 and 70 20-HE-responsive miRNAs were identified in the BmE cell line and S2 cell line, respectively. The response of miRNAs to ecdysone exhibited a time-dependent pattern, and the response intensity increased with extending exposure to 20-HE. The relationship between ecdysone and the miRNAs was further explored through knockdown of ecdysone-signalling pathway genes. Specifically, ecdysone regulated the cluster miR-275 and miR-305 through the coordination of BmEcR-B and downstream BmE75B, and the interaction between BmEcR and miR-275 cluster was strengthened by the feedback regulation of BmE75B. Ecdysone induced miR-275-3p and miR-305-5p through the ecdysone response effectors (EcREs) at the upstream of the pre-miR-275 cluster. Overall, the results might help us further understand the relationship between ecdysone signalling pathways and small RNAs in the development and metamorphosis of insects.

KEYWORDS:

• 20-hydroxyecdysone
• ecdysone-signalling pathway
• microRNA
• insect
• Bombyx mori
• Drosophila melanogaster

Introduction

The steroid hormone 20-hydroxyecdysone (20-HE), commonly known as ecdysone, functions as a crucial coordinator of cell proliferation, differentiation and apoptosis throughout the life cycle of insects and other arthropods [1–3]. Ecdysone pulses guarantee normal metamorphosis at larval-pupal transformation and almost influence all aspects of insects’ living activities [4–9]. Ecdysone triggers the cascade of gene transcription by binding to a heterodimeric nuclear hormone receptor complex formed by the nuclear ecdysone receptor (EcR) and the heterodimer partner ultraspiracle protein (USP) [3,10–13]. EcR and USP transcripts are expressed in numerous tissues of Drosophila adults, but the activity of receptor varies depending on tissue type and adult age [8,14,15]. The 20-HE-EcR-USP dynamic heterodimer acts as a ligand-dependent transcription factor and activates a small set of early responsible genes including E74, E75, Broad complex (BR-C) and E93 by binding to ecdysone response elements (EcREs) which are mostly identified distantly from promoters of these early genes, typically within 1–3 kb of their transcription start site [1,4,11,16,17]. In the presence of ecdysone, the early responsive genes function at the top of the signalling hierarchy and coordinate the activation of a large group of downstream effector genes via sequence-specific DNA binding [17,18]. The EcR gene is directly and immediately induced by ecdysone in a tissue- and stage-specific manner and plays a dual role as a repressor or activator of ecdysone-responsive genes depending on ecdysone binding [3,13,19,20]. Through the use of two promoters and as a result of alternate splicing, the EcR gene encodes three major protein isoforms (EcR-A, EcR-B1, EcR-B2) which have a conserved carboxy-terminal region and unique amino-terminal regions responsible for different functions in insects [8,13,19,21–25]. In Bombyx mori, EcR-A is highly expressed in the silk gland but EcR-B is mainly expressed in the midgut, epidermis and fat body during larval-pupal transformation [26]. 20-HE-EcR-USP induces E75 but the three members of nuclear receptor superfamily (designated E75A, E75B and E75C) encoded by E75 exhibit different responses to 20-HE [27]. Overexpression of E75B represses ecdysone biosynthesis and slows the metamorphosis of silkworm [28]. E74A and E74B show different sensitivities to 20-HE and have different effects on the regulation of programmed cell death in the control of 20-HE in larval-specific tissues of silkworm [29,30]. In summary, the steroid hormone ecdysone functions in the developmental transitions of insects through a series of signalling hierarchy pathways which involve a group of early target genes functioning at the top of the signalling hierarchy and coordinating the activation of a larger group of late genes.

In addition to these conventional pathway genes, emerging studies show that ecdysone pathway also involves some noncoding RNAs. Especially, microRNAs (miRNAs) are involved in the context of ecdysone and juvenile hormone regulation and even act as the nodes in the hormone-signalling pathway. These endogenous non-coding RNAs play a wide range of roles in processes such as metabolism, life span, cell proliferation and cell death, stress resistance and adaptive immunity [31–37]. Some miRNAs are components of additional systemic signalling pathways of insects, including ecdysone, insulin, stress, immune and circadian pathways [38]. The induction of let-7-Complex (let-7, miR-100 and miR-125) was clearly observed in cultured Drosophila S2 cells experimentally exposed to 20-HE for over 24 hours (hrs) [39–41]. Both EcR and BR-C are required in S2 cells for 20-HE to induce let-7 and miR-125 [42]. As a classic antagonist of 20-HE, juvenile hormone (JH) represses the induction of let-7 in cells treated with 20-HE [40]. 20-HE induces the upregulation of let-7 in Bactrocera dorsalis larvae while let-7 regulates the ecdysone-signalling pathway through the exact dose of E75 gene [43]. Ecdysone sensitizes its own signalling pathway and acts in a positive autoregulatory loop involved in EcR-dependent induction of EcR and miR-14-mediated repression of EcR [44–46]. Ecdysone signalling downregulates miR-965 at the onset of Drosophila pupariation to increase histoblast proliferation and migration for normal tissue morphogenesis [47]. Mosquito-specific miR-1890 controls JHA15 mRNA stability in a stage- and tissue-specific manner and is activated by the ecdysone pathway likely through an interaction with the ecdysone receptor and the early genes, E74B and Broad Z2 and even the miR-1890 promoter [48]. Ecdysone-signalling pathway activates expression of ovary-enriched miR-318 in epithelial development during Drosophila oogenesis [49]. The miRNA bantam promotes systemic growth largely by repressing the production of ecdysone [50,51]. Taken together, an increasing number of miRNAs have been confirmed to act in a variety of ways in the ecdysone-signalling pathway to control the post-transcriptional regulation of insect metamorphosis development.

Despite those advances, however, the molecular cues linking ecdysone-signalling and most miRNAs in many insects still remain largely unknown. To date, over 500 miRNAs have been identified in the silkworm, and most of them show clear temporal and spatial expression patterns [36,52–57], but only a few of them have been reported to interact with genes in the ecdysone-signalling pathway in this species. For example, bmo-let-7 exhibits a clear association with ecdysone pulse in its temporal expression [52], and targets genes, FTZ-F1 and Eip74EF (E74), key regulatory factors in the ecdysone pathway [58]; bmo-miR-281 regulates the development of Malpighian tubules through suppression of EcR-B expression [59]; two mature miRNAs of the precursor bmo-mir-14 jointly function in the ecdysone-signalling pathway: bmo-miR-14-3p targets E75 and EcR-B to maintain ecdysone homoeostasis for normal development and metamorphosis [60], and the non-dominant bmo-miR-14-5p suppresses at least nine genes to control ecdysone biosynthesis [61]. In conclusion, although increasing evidence shows a direct or indirect relationship between miRNA regulation and hormone pathway in silkworm and other insects, the complex relationship between miRNA and ecdysone is far from clear. Here we addressed this issue through systematic analysis of 20-HE upon miRNAs in different cell lines, tissues and individuals of silkworm, and compared the responses of miRNAs at the whole-genome level between the embryonic cell lines of B. mori and D. melanogaster. We further showed that miR-275 and miR-305 were regulated by 20-HE through the binding site of EcR at the upstream of pre-miR-275 cluster locus. This study might lay a basis for further investigation of the molecular mechanisms involved in the biogenesis and function of ecdysone-mediated miRNAs.

Results

Screening for 20-hydroxyecdysone-responsive miRNAs in the embryonic cells

By using Hiseq 2500 SE50, we conducted small RNA sequencing of silkworm embryonic cell line (BmE), which was incubated with 20-HE for 6 hours (6 h PIE) and 6 days (6 d PIE) (see Methods; Additional file 1). The raw data of sequencing have been deposited in the National Center for Biotechnology Information Sequence Read Achieve database (Accession NO. PRJNA632996). The average number of total reads from the 12 libraries was 16,951,010 and the average total clean reads with 18 to 30 nt were 14,315,141, accounting for 84.05% of the total reads (Additional file 2). The genome-matched reads were screened out from the total reads and the average genome-matched reads of Ecd_6 h and Ecd_6d were 7,783,007 (49.81%) and 6,039,351 (46.48%), respectively (Additional file 3). Two peaks were shown at 22 nt and 28 nt for all genome-matched reads, as are normally considered to be the classical length of miRNAs and piRNA, respectively, averagely 69.54% of the clean reads for the 22 nt were predicted to be miRNAs and 0.3% of the clean reads for the 28 nt were piRNAs (Additional file 4). The total reads of miRNAs at 6 h PIE and 6d PIE were 12,100,670 and 9,050,041, only accounting for 12.76% and 11.76% of the total clean reads, respectively (Additional file 4). For the 12 libraries, the first nucleotide of miRNA sequences with 18 nt displayed clear C bias, but those with longer sequences exhibited clear U bias at the first nucleotide; the nucleotides at 2, 4, 11, 13, 17 showed G bias, the nucleotides at 3, 5, 8, 14,15, 16, 20, 21 and 27 presented clear A bias, the nucleotides at 6, 9, 10, 12, 19, 22 through 26 displayed clear U bias (Additional file 5). In all, 416 known miRNAs and 88 novel miRNAs were identified in this work, as summarized in Additional file 6 and Additional file 7, respectively. Two hundred one out of the 416 known miRNAs were derived from two arms of one precursor (named 5p/3p miRNA), and most of these 5p/3p miRNAs showed an asymmetric expression pattern in the non-treated cells (Additional file 8). For example, in the control cells of 6 h PIE, bmo-bantam-3p showed at least 2,900 times higher in the expression than bantam-5p. However, there was another case where nine pairs of 5p/3p miRNAs were almost equally expressed in these samples, including miR-190-5p/3p, miR-2757-5p/3p, miR-2778b-5p/3p, miR-281-5p/3p, miR-33-5p/3p, miR-3389-5p/3p, miR-6497-5p/3p, miR-965-5p/3p and miR-9b-5p/3p (Additional file 8). Totally, 46 clusters of miRNAs were revealed based on the genome locations of miRNA precursors (<10 kb) (Additional file 9) and 30 groups of paralogs or families were found according to the seed regions of mature sequences of miRNAs (Additional file 10).

All miRNAs identified in each sample were then submitted to overlap analysis in the three biological replicates of each group and the overlapped members in each group accounted for the majority of known and novel miRNAs (Fig. 1A). Hierarchical cluster analysis showed that samples treated with 20-HE were clustered in the same clade and were separated from respective controls (Fig. 1B, C). In all, 22 miRNAs showed significant responses (P < 0.05) at 6 h PIE, including 17 up-regulated ones (e.g. let-7-3p, miR-13b-5p, miR-13b-3p, miR-2738, miR-275-5p, miR-275-3p, miR-305-5p, miR-305-3p, miR-8) and 5 down-regulated ones (e.g. miR-14-5p, miR-34-3p, miR-317-5p, miR-3350, miR-3405) (Fig. 1B, Additional file 11). Six pairs of 5p/3p miRNAs showed significant response to 20-HE at 6 h PIE, including miR-13b-5p/miR-13b-3p, miR-275-5p/miR-275-3p, miR-2a-5p/miR-2a-3p, miR-2b-5p/miR-2b-3p, miR-305-5p/miR-305-3p and miR-71-5p/miR-71-3p (Additional file 8). Most 5p/3p miRNAs exhibited obvious asymmetry in their response to 20-HE (Additional file 8). For example, in the DMSO-treated BmE, miR-275-3p was about 5.09 times as high as miR-275-5p, but miR-275-5p was 2.46 times as high as miR-275-3p in the 20-HE-treated BmE. Furthermore, the one with lower expression level in each 5p/3p pair showed a stronger response to 20-HE (Additional file 8). For example, let-7-5p was 11.8 times as high as let-7-3p in the control cells but only let-7-3p was significantly up-regulated by 20-HE; miR-14-3p is 71.31 times as high as miR-14-5p in the control cell but only miR-14-5p was significantly down-regulated by 20-HE.

Figure 1. 20-HE-responsive miRNAs revealed by small RNA sequencing in the BmE cell line. (A) Overlap analysis of known miRNAs confirmed in all examined samples. (B) Hierarchical cluster analysis of miRNAs responding to 20-HE at 6 h PIE. (C) Hierarchical cluster analysis of miRNAs responding to 20-HE at 6 d PIE. S01-S03 and S07-S09 are Ecd-treated samples, S04-S06 and S10-S12 are controls incubated with DMSO.

At 6 d PIE, 66 miRNAs showed a significant response to 20-HE (p < 0.05), 39 of which were significantly up-regulated and 27 showed significant down-regulation (Fig. 1C, Additional file 12). There were more miRNAs responding to 20-HE at 6 d PIE than at 6 h PIE likely due to the accumulating effect of 20-HE. Top three of the most strongly up-regulated miRNAs at 6 d PIE were miR-13-b-3p, miR-2a-2-5p and miR-2798, which were increased to 15.08 folds, 12.72 folds and 12.14 folds, respectively. Top three of the most strongly down-regulated miRNAs were miR-34-5p, miR-317-5p and miR-317-3p, which were decreased to 104.46 folds, 65.08 folds and 25.92 folds, respectively. However, the accumulating effect of 20-HE was not observed for all responding miRNAs. Especially, miR-3350 and miR-3405 were significantly down-regulated by 20-HE at 6 h PIE but showed no obvious response at 6 d PIE. The expression asymmetry of 5p/3p miRNAs also existed at 6 d PIE (Additional file 8). For example, the expression level of bantam-3p was over 2,000 times higher than that of bantam-5p. Totally, 12 pairs of 5p/3p miRNAs displayed upregulation response and two pairs (miR-317-5p/3p and miR-34-5p/3p) showed downregulation, whereas either expression level or response degree of 5p miRNA was different from that of 3p miRNA.

miRNAs respond to ecdysone in a time-dependent pattern

To further confirm the response of miRNAs to 20-HE revealed by sequencing, we incubated BmE cells with 20-HE and collected cells at different time points for qRT-PCR assay of some up-regulated, down-regulated and non-responding miRNAs (Additional Fig. S1). As the up-regulated miRNAs revealed by sequencing, let-7-3p and miR-8-3p remained up-regulated from 3 h PIE to 6 d PIE, miR-2795 was up-regulated from 12 h PIE to 48 h PIE, miR-100 and miR-87 began to be up-regulated at 24 h PIE and maintained this response to 20-HE until 6d PIE, let-7-5p did not show up-regulation until 48 h PIE. On the contrary, miR-34-3p and miR-14-5p showed significant down-regulation at 6 h PIE, which was well consistent with the response revealed by sequencing. The response of miR-34-5p began to be observed 3 hrs later than those of let-7-3p and miR-8-3p and remained down-regulated until 6 d PIE. Although coexisting in a cluster on the genome, bmo-let-7, miR-100 and miR-2795 did not respond to 20-HE synchronously, which was consistent with the responses revealed by sequencing. Another cluster of miRNAs, bmo-miR-275 and miR-305, was significantly up-regulated at 6 h PIE and 6 d PIE (Fig. 2A), showing the same responses to 20-HE as revealed by sequencing (Fig. 1B, C). Although the time-dependent response was clearly observed for these miRNAs, bantam and miR-281-5p did not respond to 20-HE through 3 h PIE to 6 d PIE, showing that their expressions are unrelated to 20-HE in the embryonic cell line of silkworm.

Figure 2. Time-course responses of miRNAs to 20-HE and JHA in the BmE cell line. (A) Response of miR-275 and miR-305 to 20-HE at 6 h PIE and 6 d PIE. (B) Time-course upregulation of miRNAs by 20-HE. (C) Time-course downregulation of miRNAs by 20-HE. (D) Responses of miRNAs to JHA. (E) Buffering effect of 20-HE and JHA on the responses of miRNAs and EcR/USP to 20-HE and JHA. All cells were collected at 24 h PIE for (E). The same letters in (E) represent no significant differences, while different letters indicate significant differences. The data are shown as mean ± SEM. *P <0.05; **P <0.01.

A larger number of responding miRNAs were identified in BmE at 6 d PIE (Fig. 1; Additional Fig. S1), indicating that the responses of miRNAs to 20-HE change with the extending ecdysone incubation. In order to further confirm the time-dependent response, we incubated BmE cells with 5 μM 20-HE and then examined the responding miRNAs in a more detailed time course. The up-regulation of miR-275-3p, miR-305-5p, miR-8-3p and miR-2a-2-3p could be observed at 2 h PIE (Fig. 2B), while the down-regulation of miR-34-5p, miR14-3p and miR-317-3p probably began early before 2 h PIE (Fig. 2C). Clearly, the response intensity of all miRNAs increased with extending exposure to 20-HE.

Juvenile hormone (JH) is a classic antagonist of 20-HE, so we additionally treated BmE cells with the juvenile hormone analogue methoprene (JHA) to investigate whether reverse responses could be observed. Indeed, contrary to their responses to 20-HE, miR-34-5p and miR-14-3p began to be up-regulated by the JHA at 2 h PIE and miR-317-3p was up-regulated by JHA even immediately after the treatment; the response intensity of these miRNAs showed a stronger response in cells with extended exposure to JHA (Fig. 2D). To further address the potential buffering effect between ecdysone and juvenile hormones on the responses of miRNAs, we incubated the BmE cells with the mixture of 20-HE (1 μM) and JHA (10 μM) and examined the expression of miRNAs at 24 h PIE (Fig. 2E). The expression levels of miR-275-3p, miR-305-5p, miR-8-3p and miR-2a-2-3p were not changed in cells when only exposed to JHA but were robustly induced in cells when only exposed to 20-HE. When the cells were exposed to the mixture of these two hormones, the 20-HE-induced expression of miR-275-3p and miR-305-5p was significantly repressed by JHA but the 20-HE-induced expression of miR-8-3p and miR-2a-2-3p was not affected by JHA. Clearly, the expression of miR-275-3p and miR-305-5p was not influenced by JHA alone but JHA could be involved in the accumulation of these two miRNAs when the pathway of ecdysone was activated by 20-HE. Unlike miR-275-3p and miR-305-5p however, the expression of miR-8-3p and miR-2a-2-3p was unrelated to JHA even though in the context of ecdysone. As revealed in the time course test, miR-34-5p and miR-317-3p were significantly upregulated when the cells were only exposed to JHA and were significantly downregulated when the cells were only exposed to 20-HE (Fig. 2C, D). However, their responses to each hormone were obviously buffered when the cells were incubated in the mixture hormones (Fig. 2E). 20-HE and JH regulate the expression of miRNAs either individually or jointly, and display various types of regulation on the production of miRNAs, suggesting that the intricate relationship exists between hormone regulation pathways and miRNA regulation pathways.

As EcR is required by 20-HE to effectively regulate miRNAs [42], we further examined the expression of EcR and USP in the cells incubated with 20-HE, JHA and the mixture hormones at the same time points as those for miRNAs and found that EcR and USP readily responded to 20-HE but did not respond to JHA (Fig. 2E). However, their induction by 20-HE was clearly reduced by JHA in the 20-HE-treated cells. Accordingly, it can be inferred that when cells are treated with a mixture of ecdysone and juvenile hormone, the latter might indirectly regulate the production of miR-275 and miR-305 by affecting ecdysone pathway genes, such as EcR.

miRNAs respond to 20-HE in different cell lines, tissues and individuals

We questioned whether the responses of miRNAs to 20-HE were conserved across cells and tissues in the silkworm. To address this issue, we incubated the ovary cell line (BmNS) with 20-HE, and examined some miRNAs at 6 h PIE and 6 d PIE (Additional Fig. S2). Up-regulation of miR-8-3p and let-7-3p, as well as down-regulation of miR-34-5p, was clearly observed in the ovary cell line at two tested time points; miR-87 did not show up-regulation at 6 h PIE but was strongly up-regulated at 6 d PIE; although in the same cluster miRNAs, let-7-3p responded to 20-HE more strongly and quickly than let-7-5p and miR-100. Like the responses in BmE cell line, miR-281-5p and bantam showed no clear changes in BmNS cell line. However, some miRNAs showed different responses to 20-HE in these two cell lines. For example, miR-14-3p was significantly down-regulated by 20-HE in BmE at 6 h PIE and 6 d PIE but was not changed in BmNS at these two time points; miR-281-3p was significantly up-regulated in BmE but did not show any clear change in BmNS cell line at 6 d PIE; unlike in BmE cell line, no clear up-regulation was observed for let-7-5p and miR-100 in the ovary cell line at 6 d PIE. The same or different responses of miRNAs to 20-HE between silkworm embryonic cell line and ovarian cell line might indicate that common or unique ecdysone-signalling pathways exist for these miRNAs in different cell lines of silkworm.

To reveal whether the response of miRNA is related to the ecdysone amount in the larval body, we separately injected two doses of 20-HE (2 μmol and 5 μmol) into the silkworm at early fourth instar larval stage and collected the whole body of individuals at 24 h post-injection (24 h PIJ) for qRT-PCR assay (Additional Fig. S3A). The results showed that miR-14-3p, miR-100 and miR-2a-2-3p exhibited significant up-regulation in the whole body injected with lower dose (2 μmol) of 20-HE, while miR-275-3p, miR-275-5p and miR-34-5p were significantly down-regulated by both doses of 20-HE; clearly, the responses of miR-100, miR-2a-2-3p and miR-34-5p were consistent with the results revealed in the BmE cell line (Fig. 1 and Additional Fig. S1); we did not find significant responses to 20-HE for let-7-5p, let-7-3p and miRR-8-3p in the whole body. Obviously, however, the responses of miR-14-3p, miR-275-3p and miR-305-5p to 20-HE in the larval body were opposite to their responses in BmE (Fig. 1B, C, Additional Fig. S1).

The response of miRNAs in the larval body should represent the comprehensive accumulation of miRNAs in diverse tissues and cells. To reveal whether and how miRNAs respond to 20-HE in the tissues of silkworm, we performed qRT-PCR assay of some miRNAs in the fat body and haemocyte, which were collected from silkworm larvae at 24 h PIJ of 2 μmol 20-HE (Additional Fig. S3B). Like the responses in the larval body, miR-14-3p was significantly up-regulated and miR-34-5p was significantly down-regulated in the fat body and haemocyte; opposite to the responses in the larval body, but similar to the responses in the BmE, up-regulation of miR-275-3p and miR-305-5p was manifested in the fat body and haemocyte of 20-HE-injected larvae; different from the responses to 20-HE in the BmE (Fig. 1B, C, Additional Fig. S1), no clear expression changes were observed for let-7-5p, let-7-3p, miR-100, miR-8-3p and miR-2a-2-3p in the fat body and haemocyte of 20-HE-injected larvae; similar to the results revealed in BmE and BmNS cell lines, miR-281-3p and bantam did not respond to 20-HE in the fat body and haemocyte. The fat body of insects is the counterpart of mammalian liver and adipose tissue and is the main target organ of ecdysone, in which the transcription factor Myc and miR-8 act as the downstream effectors for gene expression mediating organismal growth [62–64]. Furthermore, the fat body of insect larvae functions as an internal nutrient sensor that coordinates organismal nutritional metabolism and growth in response to nutrient availability [65,66]. To further confirm the responses to 20-HE in the fat body, we checked some miRNAs in the cultured fat body which was isolated from early fourth instar larvae and cultured in Grace medium with 20-HE (5 μM) for 24 hrs (Additional Fig. S3C). No significant expression changes were detected for let-7-5p, miR-100, miR-2a-2-3p and miR-8-3p, significant up-regulation of miR-275-3p and miR-305-5p and significant down-regulation were clearly observed for miR-34-5p in the cultured fat body, which is consistent with their responses to 20-HE in the fat body and haemocyte in vivo. Unlike the up-regulation response in the fat body in vivo however, no clear response was observed for miR-14-3p in the cultured fat body.

Response of miRNAs to ecdysone in the Drosophila S2 cells

In order to explore whether the response of miRNAs to ecdysone was conserved across species, we performed small RNA sequencing to examine Drosophila miRNAs in the 20-HE-treated Drosophila S2 cell line. S2 cell line was originally derived from primary cultures of late-stage D. melanogaster embryos [67] and provides a highly versatile expression system for the study of genes’ functions and interactions in the post-genome applications [68]. The cells were separately incubated with 20-HE for 6 hrs and 48 hrs before collection for RNA extraction (Additional file 1). The sequencing raw data of the 12 libraries have been deposited in the National Center for Biotechnology Information Sequence Read Achieve database (Accession NO. PRJNA632996). The data were filtered and analysed as described in Methods. Finally, 361 known Drosophila miRNAs were identified in S2 cell lines (Additional file 13). Further analysis showed that 9 Drosophila miRNAs exhibited significant responses to 20-HE at 6 h PIE, including 3 up-regulated ones (dme-miR-282-5p, miR-282-3p and miR-278-5p) and 6 down-regulated ones (dme-miR-277-5p, miR-14-5p, miR-980-5p, miR-317-5p, miR-34-5p and miR-34-3p) (Fig. 3A, B; Additional file 14); 78 miRNAs exhibited significant responses to 20-HE at 48 h PIE, including 43 up-regulated ones (dme-let-7-5p, let-7-3p, miR-100-5p, miR-125-5p, etc.) and 35 down-regulated ones (dme-miR-34-5p, miR-34-3p, miR-14-3p, etc.) (Fig. 3A, C; Additional file 15). The differentially expressed miRNAs identified at 6 h PIE remained the same response patterns at 6 d PIE. The three biological replicates of 20-HE incubation were clustered in one clade and were separated from the controls (Fig. 3B, C). Consistent with the silkworm miRNA responses to 20-HE in BmE cell line (Fig. 1B, C; Additional Fig. S1), qRT-PCR results confirmed that some homologous Drosophila miRNAs (dme-let-7-5p, miR-100 and miR-2a-2-3p, etc.) were significantly up-regulated while some others (miR-14-3p and miR-34-5p, miR-34-3p, etc.) were significantly down-regulated in the 20-HE-treated S2 cells (Fig. 4). However, some homologous miRNAs differed in their responses to 20-HE between BmE and S2 cell lines. For example, the cluster miRNAs miR-275-3p and miR-305-5p did not show obvious response in S2 at 6 h and 48 h PIE (Fig. 4A, B) but were significantly up-regulated in BmE cells at all examined time points (Fig. 1B, C; Fig. 2A); miR-8-3p was up-regulated by 20-HE in BmE through 3 h to 6 d PIE (Additional Fig. S1), but displayed no obvious response at 6 h PIE and was significantly down-regulated at 48 h PIE in S2 (Fig. 4). The responses of miRNAs to ecdysone in S2 might also depend on the length of incubation time. For example, let-7-5p showed stronger up-regulation at 48 h PIE than at 6 h PIE; miR-34-5p showed stronger down-regulation at 48 h PIE than at 6 h PIE; miR-14-3p did not respond to 20-HE at 6 h PIE but was significantly down-regulated at 48 h PIE, whereas miR-2a-2-3p was significantly up-regulated at 6 h PIE but was not changed at 48 h PIE (Fig. 4). Collectively, the responses of miRNAs to ecdysone in these two embryonic cells are characterized by similarities and individualities, i.e. some exhibit the same response in both embryonic cell lines, and some exhibit different or even opposite responses in these cell lines. The response differences of miRNAs between BmE and S2 cells indicate that these miRNAs are likely to be mediated by ecdysone under different mechanisms in different species.

Figure 3. Drosophila miRNAs responding to 20-HE in the S2 cell line. (A) The number of responding miRNA at two time points. (B) Hierarchical cluster analysis of differentially expressed miRNAs at 6 h PIE. (C) Hierarchical cluster analysis of differentially expressed miRNAs at 2 d PIE.

Figure 4. Responses of Drosophila miRNAs to ecdysone in the S2 cell line. (A) Responses at 6 h PIE. (B) Responses at 48 h PIE. The data are shown as mean ± SEM. *P <0.05; **P <0.01.

Ecdysone-signalling pathways facilitate the response of miRNAs to 20-HE

To address how the responses of miRNAs to 20-HE were facilitated by ecdysone signalling, we first knocked down BmEcR-B, the receptor of 20-HE in the silkworm (Fig. 5A). Clearly, the transcripts of EcR downstream genes including E75A, E75B, BR-CZ2, BR-CZ4 and E93 were significantly decreased by RNAi of BmEcR-B (Fig. 5B). In response to BmEcR-B RNAi, the ecdysone-upregulated miRNAs, including miR-275-3p, miR-305-5p, miR-2a-2-3p and miR-8-3p, exhibited significant down-regulation, but no obvious expression changes were detected for miR-34, miR-14 and miR-317 (Fig. 5C) although they were significantly down-regulated by ecdysone in BmE cells (Fig. 1B, C; Additional Fig. S1). To clarify whether the response of these miRNAs to 20-HE was dominated by EcR gene alone or by the EcR downstream genes, we then knocked down the EcR downstream genes one by one in BmE cells (Fig. 5D), followed by examining the expression of ecdysone-responsive miRNAs (Fig. 5E). The expression levels of miRNAs were not significantly changed by the knockdown of BmE75A, BmBR-CZ2, BmBR-CZ4 and BmE93 (Fig. 5E), thus showing that these early pathway genes are not responsible for the responses of miRNAs to 20-HE. Unexpectedly, however, miR-275 and miR-305 were significantly induced when BmE75B was knocked down (Fig. 5E). As one downstream gene of BmEcR-B, BmE75B is positively responding to 20-HE and BmEcR-B [28,69]. The effect of BmE75B on these two miRNAs is contrary to the effect of 20-HE, largely due to the existence of a fine E75B-mediated regulatory loop in steroidogenesis [28]. Indeed, the expression of BmEcR-B was significantly increased at 48 h post-RNAi of BmE75B (Fig. 5F). Therefore, ecdysone regulates the cluster miR-275 and miR-305 through the coordination of BmEcR-B and its downstream gene E75B. Our previous work showed that both miRNAs were especially highly expressed in the head, body wall, fat body and whole body of silkworm from cocoon spinning stage to early pupa [53,54]. The similar spatial expression patterns were further confirmed in different tissues of female and male prepupae (Additional Fig. S4A, B). On the contrary, BmEcR-B was expressed at low levels in the head, body wall, fat body and gonads but highly expressed in the midgut and silk gland (Additional Fig. S4C). We checked their temporal expression patterns in the fat body of silkworm and found that both miR-275 and miR-305 were exclusively highly expressed from cocoon spinning stage to early pupal stage (Additional Fig. S4D, E). EcR-B expression in the fat body peaked at early pupa time point (Additional Fig. S4F), just when the developmental changes of ecdysteroid titre occur in haemolymph [70]. These expression data indicated that the miR-275 cluster might play a significant role in the pupa transformation and metamorphosis in the silkworm under the coordinate control of BmEcR.

Figure 5. Relationship between 20-HE signalling pathway genes and miRNAs in the BmE cell line. (A) RNAi knockdown of BmEcR-B. (B) Expression of the 20-HE pathway genes post-knockdown of BmEcR-B. (C) Expression of 20-HE-responsive miRNAs after knockdown of BmEcR-B. (D) RNAi knockdown of BmE75A, BmBR-CZ2, BmBR-CZ4, BmE75B and BmE93. (E) Expression of 20-HE-responsive miRNAs after knockdown of 20-HE pathway genes. (F) Expression of BmEcR-B and BmUSP after knockdown of BmE75B. The cells were collected for qRT-PCR assay after treatment as follows. The cells in each well of 12-well plate were transfected with dsRNA (5 μg/ml) for 36 h and were followed by treatment with 20-HE (5 μM) for 12 h. Three biological replicates were performed for all examined genes and miRNAs. The data are shown as mean ± SEM. *P <0.05; **P <0.01.

EcR element (EcRE) at the upstream of miR-275 cluster

We supposed there might exist some EcR elements (EcREs) at the upstream of miR-275 cluster. In order to identify these potential cis-control elements of miR-275 and miR-305, we cloned the whole fragment of genomic DNA spanning the ~10 Kb bmo-miR-275-C locus into pGL3 basic vector containing the hsp70 minimal promoter (Promega, Madison, WI, USA) (Fig. 6A). The pRL vector (Promega, Madison, WI, USA) containing the vitellogenin promoter for Renilla luciferase was used as the control for normalization (Fig. 6B). The whole fragment was cut into nine fragments and cloned into the hsp70 promoter-contained pGL3 basic vector one by one, followed by dual-luciferase reporter assay in vitro. Fragments 4 and 7 presented significant higher luciferase activities, indicating that they have EcR elements of miRNAs (Fig. 6C). The two fragments were further cut into two shorter fragments and subjected to dual-luciferase assay, respectively. Fragments 4–2 and 7–2 showed higher luciferase activities (Fig. 6D). Similarly, fragments 4–2 and 7–2 were further cut into three and six fragments and subjected to dual-luciferase assay, respectively. The fragments 4-2-4 and 7-2-6 displayed higher luciferase activities (Fig. 6E, F) and were cut into even shorter fragments. Finally, the candidate sites of EcRE were limited to locus between fragments 4-2-4-1 and 4-2-4-2 (−3759_-3158), and locus between fragments 7-2-6-5 and 7-2-6-6 (−591_-538). To identify the potential binding site of EcR on the fragments, we blasted against the loci between fragments 4-2-4-1 and 4-2-4-2, and between fragments 7-2-6-5 and 7-2-6-6 with the EcRE sequence motif of Drosophila, combined with the prediction of transcription factor with online software (http://jaspar.genereg.net) [71]. In all, there existed four EcREs which matched at least 10 bp of the consensus (Fig. 6G). To verify the validity of these EcREs, the fragments containing predicted EcREs and their mutants were submitted to the dual-luciferase reporter assay, respectively. The luciferase activities of each mutated fragment decreased by 80% (Fig. 6H), indicating the authentic existence of these EcREs at these loci.

Figure 6. Ecdysone response elements (EcREs) at the locus of miR-275 cluster. (A) Schematic vector of firefly luciferase reporter. (B) Schematic vector of Renilla luciferase reporter. The promoter of SV40 in the vector pRL-SV40 was replaced by the promoter of vitellogenin (Vg). (C) Luciferase reporter assay of cloned locus fragments. 20-HE-treated BmE cells were transfected with a series of luciferase reporter vectors containing miR-275-C locus fragments 1–9. The empty vector PGL-hsp70 was used as negative control. (D) Fold increase of luciferase activity in 20-HE-treated BmE cells transfected with luciferase reporter vectors containing the fragments (F-4-1, F-4-2, F-7-1, F-7-2). The empty vector PGL-hsp70 was used as negative control. (E) Luciferase activity of truncated fragment F-4-2. (F) Luciferase activity of truncated fragments 7–2. For the convenience of description, the first base of pre-miR-275 cluster was designed as the +1 base in (E) and (F). (G) Predicted EcREs and mutant sequences. The EcRE-1 (wt-1), EcRE-2 (wt-2), EcRE-3 (wt-3) and EcRE-4 (wt-4) sequences were predicted by using the sequences of Drosophila EcRE. (H) Fold increase of luciferase activity in 20-HE-treated BmE cells transfected with luciferase reporter vectors containing mutant miR-275-C EcREs. At least five biological replicates were included for all reporter assays. The cells were transfected with these vectors 6 h before incubation with 20-HE (5 μM) and were then collected for luciferase activity assay at 24 h PIE. The ratio values of firefly luciferase/Renilla luciferase were normalized by setting the ratio value of empty hsp70 vector as 1.

To evaluate whether the heterodimer complex formed by EcR and USP can combine with EcREs on the locus of miR-275-C, we used the sequences containing the predicted binding sites as the probes and performed EMSA assay by incubating the 3ʹ-terminal biotin-labelled probe with the purified proteins of EcR and USP. The result showed that the target band signal gradually increased with the rising concentration of proteins (Fig. 7A–C, lane 2–4; Fig. 7D, lane 2–3). However, the band signal became weaker or even disappeared when the amount of competing probes increased (Fig. 7A–C lane 5, 6; Fig. 7D, lanes 4–6). To verify that the identified site sequence was responsible for the binding of EcR, we generated mutant sites where some nucleotides were replaced by other nucleotides (Fig. 6G). When the mutant probes were used in the incubation, the target bands became weaker or even disappeared (Fig. 7A–D, lane 7). Therefore, these sites of EcRE could specifically bind to EcR and USP proteins. Furthermore, we explored whether the binding of EcREs was limited to EcR by incubating EcRE probes with EcR and USP separately. The bound bands were observed only when the labelled probes were incubated with the EcR protein (Fig. 7A–D, lane 8) but no bands were observed when the probes were incubated with USP protein (Fig. 7A–D, lane 9), thus showing that only EcR protein bound directly to the EcRE on the miRNA locus while USP was likely unrelated to this special binding of EcREs.

Figure 7. Electrophoretic mobility shift assay (EMSA) of EcREs in the pre-miR-275 cluster locus. (A) EMSA of the 3ʹ-biotin-labelled EcRE-1 (wt-1). (B) EMSA of 3-biotin-labelled EcRE-2 (wt-2). (C) EMSA of 3ʹ-biotin-labelled EcRE-3 (wt-3). (D) EMSA of 3ʹend biotin-labelled EcRE-4 (wt-4). All probes were incubated with EcR-B and USP proteins. com-P is the biotin-unlabelled competition probe which is of the same sequence as that of EcRE. mut-P is the 3-biotin-labelled mutant.

Discussion

The moulting hormone ecdysone binds to the heterodimeric nuclear receptor EcR/Usp and activates a cascade of nuclear hormone receptors involved in various developmental processes and as a result provides the temporal signal for various developmental transitions of insects [50,72]. High concentrations of ecdysone are required for its role as a developmental timer, and basal ecdysone levels contribute to the inhibition of tissue growth [63], which is achieved at least in part through EcR activation and Myc repression in the fat body [64]. Functions of hormones might be achieved through the interaction of hormone-signalling genes with miRNAs which manage complex gene networks controlling developmental decisions [73,74]. Drosophila miRNA bantam carries essential and conserved functions in coupling insulin signalling with steroid production to act on tissue growth and body size [51]. This functional similarity between ecdysone and miRNA underscores the involvements of 20-HE in the regulation of miRNA biogenesis. The expression level of silkworm let-7 is higher at the end of each instar larval stage and the beginning of each moulting stage, in accordance with the ecdysone pulse [52]. Drosophila let-7 cluster, which is composed of pri-let-7, pri-miR-100 and pri-miR-125, is a direct target of EcR [75]. Ecdysone signalling reciprocally regulates miR-14 and EcR in Drosophila, forming a positive autoregulatory circuit to ensure robustness in development [76]. The conservative miR-8 is transcriptionally repressed by ecdysone’s early response genes in Drosophila [62]. In mosquitoes, miR-275 and miR-1890 were up-regulated by ecdysone, and may also be regulated by the early genes EcR and BR-C [48,77]. Although the biogenesis and function of an increasing number of miRNAs are regulated by ecdysone hormone, it is far from clear how these hormone signals are transmitted to the downstream miRNAs to coordinate the regulation of life activities in many insects other than Drosophila species.

The biogenesis and function of miRNA are ordinated by various signals inside and outside the cell. It is still an open question as to how miRNAs accumulate and function according to different cells, tissues, organs and individuals, as well as different internal and external environmental factors. In the present study, we systematically studied the responses of silkworm miRNAs to 20-HE in cell lines, tissues and individuals, and compared the responsive miRNAs between BmE and S2 cell lines. It has been documented that the early genes, like E75A, E75A, E75B, E78 and DHR3, begin to respond to 20-HE even at 1 h PIE in Drosophila [78]. Interestingly, we found that some miRNAs responded to 20-HE within 3 h PIE or even early at 0.5 h PIE (Fig. 2B, C; Additional Fig. S1). Reasonably, the quick response of miRNAs to ecdysone might be caused by EcR or the early genes of ecdysone pathway. We knocked down the early genes one by one by RNAi and found that deletion of most early ecdysone-signalling genes did not affect the responses of miRNAs to 20-HE. Therefore, these miRNAs are very unlikely mediated by the early genes in the ecdysone-signalling pathway. We are clear that in this study, we only detected the down-regulation of these pathway genes at the transcription level. Probably they had no significant changes at the protein level, and trace amount of proteins might be enough to maintain the ecdysone-signalling pathway, so it is necessary in the future to confirm the impact of protein level changes of these pathway genes on these miRNAs.

Intriguingly, however, deletion of BmE75B in BmE cells induced the response of miR-275 cluster, showing that BmE75B negatively regulates miR-275 cluster through a feedback loop existing between BmEcR and BmE75B. After BmEcR RNAi, the miRNAs which positively responded to ecdysone, such as miRNA-275 and miRNA-305, were significantly down-regulated, while those which negatively responded to 20-HE, like miR-34-5p and miR-317-5p, were not up-regulated or were unaffected (Fig. 5C). Therefore, those ecdysone-induced miRNAs should be directly regulated by BmEcR or its downstream genes and those ecdysone-deduced miRNAs should be unrelated to BmEcR or its downstream genes. RNAi of BmE75B caused an up-regulation of EcR, further confirming that BmE75B can negatively regulate BmEcR in a feedback loop [28]. Those miRNAs which were down-regulated due to BmEcR RNAi should be up-regulated via RNAi of BmE75B. It is not determined whether EcR could directly interact with the transcription machinery in miRNA biogenesis and control the processivity of RNA polymerase Ⅱ in the silkworm cell lines. We did not find obvious changes in the expression of Drosha, Dicer-1 and AGO-1 in response to 20-HE (Additional Fig. S5). According to our observation that miRNAs responded to 20-HE in diverse ways, conservative across or specific to cells, tissues and organisms, the current model for the relationship between ecdysone and miRNA could be proposed as follows: (a) miRNAs are not up-regulated by 20-HE through the early genes of ecdysone pathway but directly through EcR; (b) miRNAs are down-regulated by 20-HE through the coordination of JH; and (c) E75B negatively regulates some miRNAs, like miR-275 and miR-305, through a feedback effect on EcR (Fig. 8).

Figure 8. Responding model of miRNAs involved in the 20-HE/JH pathway. There existed four EcREs at the upstream of pre-miR-275 cluster.

In our results, some homologous miRNAs showed different response patterns to ecdysone in the BmE and S2 cell lines (Fig. 8). Specifically, miR-275, miR-277, miR-8 and etc. were up-regulated by 20-HE in BmE but were down-regulated in S2, indicating that these miRNAs may have different regulatory pathways. It is not known whether they also respond differently to ecdysone in vivo in the silkworm and fruit fly. Although the response of miRNAs to 20-HE has been extensively studied here in BmE and S2 cell lines, the molecular mechanism regulating the interaction between miRNAs and 20-HE remains less well understood. Although we have revealed the response of miRNAs to ecdysone, their target genes need to be identified in order to understand the regulatory pathways and functions of these responsive miRNAs. Increasing evidences showed that some miRNAs are also involved in regulating hormone pathway genes. Several upstream genes encoding ecdysone biosynthesis enzymes are repressed upon miRNA bantam overexpression, but the molecular mechanism by which bantam represses ecdysone biosynthesis remains elusive [50]. In B. mori, miR-14-3p and miR-14-5p repress BmE75, BmEcR-B and several other genes in ecdysone-signalling pathway to maintain ecdysone homoeostasis for normal development and metamorphosis of this species [60,61]. In D. melanogaster, ecdysone down-regulates miR-14 through EcR, but EcR is a direct target gene of miR-14 [76]. In Blattella germanica, the phenotype from down-regulated Dicer-1 is similar to that of JH treatment, but down-regulation of Dicer-1 does not induce the accumulation of JH [79]. In oriental fruit fly (Bactrocera dorsalis), let-7 controls normal larval-pupal development by regulating the ecdysone-signalling pathway through the exact dose of E75 gene [43]. Another conserved miRNA, miR-34, acts as a node in the intricate interplay between ecdysone signalling and innate immunity in D. melanogaster, and a mutual repression exists between miR-34 expression and ecdysone signalling [80]. 20-HE and JH exhibit opposite effects on the expression of miR-34, miR-100, miR-125 and let-7, and JH inhibits the induction effect of 20-HE on miR-100, miR-125 and let-7 [40,81]. Our results showed that 20-HE and JH individually or jointly regulate the expression of different miRNAs (Fig. 2C–E). Emerging evidences show that miRNAs regulate JH and 20-HE-signalling pathways and these two hormones reciprocally coordinate the expression of miRNAs, thus forming regulatory loops of miRNAs with JH and 20-HE-signalling cascades. However, the fine regulatory mechanism involved in the interplay between hormones and most hormone-responsive miRNAs still remains undetermined, and how their complex relationship contributes to different physiological environments also needs further exploration.

Materials and methods

Insects and cell lines

The larvae of B. mori strain P50 (Dazao) were reared on fresh mulberry leaves at 27°C under a 12 h light/12 h dark photoperiodic regime with a relative humidity of 70–80%. The BmE cells were cultured at 27°C in Grace Insect culture medium (Gibco, Carlsbad, CA, USA) supplemented with 10% heat-inactivated foetal bovine serum (Gibco, Carlsbad, CA, USA) and 1% penicillin/streptomycin. The BmNS cells were cultured at 27°C in TC-100 medium supplemented with 10% heat-inactivated foetal bovine serum (Gibco, Carlsbad, CA, USA) and 1% penicillin/streptomycin. The Drosophila Schneider 2 (S2) cells were cultured at 27°C in Schneider’s Drosophila Medium (Thermo, Waltham, MA, USA) supplemented with 10% heat-inactivated foetal bovine serum (Gibco, Carlsbad, CA, USA) and 1% penicillin/streptomycin.

Incubation of cells with 20-HE and JHA

20-HE and JHA were ordered from Sigma–Aldrich (St. Louis, MO, USA) and were dissolved in Dimethyl sulphoxide (DMSO). The cells were separately cultured in the presence or absence of 20-HE or JHA, and the identical volume of solvent DMSO was used as the control for all treatments. Three biological replicates were set for each treatment and control. For the small RNA sequencing, BmE cells were exposed to 5 μM 20-HE for 6 h and 6 days before collection for RNA extraction, and S2 cells were incubated with 5 μM 20-HE for 6 h and 2 days before collection for RNA extraction. For the time course test, both BmE and S2 cells were exposed to 20-HE (5 μM) or JHA (1 μM) and collected at 0, 0.5, 2, 4, 6, 12, 24, and 48 h of incubation.

In vitro culture of fat body

The fat body was isolated from fourth-instar day-1 larvae of silkworm and was washed in PBS (1%) and Grace medium for 3–5 times, respectively. The fat body was transferred to one 12-well plate and cultured at 27°C in 10% FBS Grace medium for 24 h. The fat body in the wells was then exposed to 20-HE (5 μM) for another 24 h before harvesting for total RNA extraction. The fat body incubated in DMSO was used as the control.

Quantitative Real-Time PCR for miRNAs and protein-coding genes

Total RNA was extracted from tissues and cells by using TRIZOL® reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instruction. Reverse transcriptions for mRNA and miRNA were performed separately with different methods. In detail, 1 μg of total RNA was utilized to generate the oligo (dT)-primed cDNA by using PrimerScript RT Reagent kit with gDNA Eraser (Takara Bio, Shiga, Japan); 10 μg of total RNA was digested for eliminating genome with DNase I (Takara Bio, Shiga, Japan), then used Mir-X miRNA First-Strand Synthesis Kit (Takara Bio, Shiga, Japan) for cDNA synthesis that used to detect miRNAs. The reaction system of qRT-PCR for miRNA: SYBR Advantage Premix (2x) 10 μL, miRNA-specific Primer (10 μM) 0.4 μL, mRQ 3 Primer, cDNA 2 μL and ddH₂O 7.2 μL. The reaction system for mRNAs: SYBR Advantage Premix (2x) 10 μL, gene-specific Forward primer (10 μM) 0.8 μL, gene-specific Reverse primer 0.8 μL, cDNA 2 μL and ddH2O 6.4 μL. As endogenous controls, a house-keeping gene silkworm translation initiation factor 4A (SW22934) and U6 were used to normalize the mRNA and miRNA expression levels, respectively. All primers used for qRT-PCR are listed in Additional Table S1. The qRT-PCR was conducted on a qTOWER 2.2 Real-Time PCR (Analytik Jena AG, Jena, Germany) by using the kit SYBR Premix Ex Taq™ II (Takara Bio, Shiga, Japan) according to the manufacturer’s instructions. Each treatment had at least three biological repeats with three technical repeats. The same qRT-PCR programme was used for the test of miRNAs or mRNAs: initial incubation of 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Melt curve was generated by 60°C for 15 s, 95°C for 1 min and 95°C for 15 s. Statistical analyses were conducted by using Graphpad Prism (7.0) (GraphPad Software, San Diego, CA, USA).

Small RNA sequencing

Total RNA was isolated by TRIZOL® reagent according to the manufacture’s instruction and successively examined by NanoDrop, Qubit 2.0, Agilent 2100. In all, we constructed 12 small RNA libraries with small RNA Sample Pre Kit (NEB Multiplex Small RNA Library Prep Set for Illumina) and purified them with 6% PAGE. All libraries of small RNAs were qualified with Agilent 2100 Bioanalyzer and were then submitted to sequencing by Hiseq 2500 SE 50. After sequencing, the raw data were strictly filtered to generate clean reads which were further analysed with miRDeep 2 [82]. By using the SOAPnuke [83], the raw data were submitted to filtering through removing the adaptors, sequences of low quality and sequences containing N to obtain the clean reads which were then filtered by blasting against several databases, namely Genebank, Repbase and Rfam to remove those sequences of tRNA and rRNA. By using Bowtie2 [84], the filtered sRNA equal to and larger than 18 nt was aligned against the reference genome, and the mapped information and distribution of all reads were evaluated as a whole. To annotate the filtered data, all clean reads were blasted against Gene bank of NCBI to remove rRNA, tRNA, ncRNA, etc., were blasted against the Repeat database to remove the repeated sequence, and were blasted against the Rfam database to remove the sequence related to rRNA, scRNA, snoRNA, snRNA, tRNA, etc. In order to make each unique small RNA have a unique annotation, all clean reads were annotated in the order of priority: miRBase > snoRNA (human/plant) > Rfam > other sRNA. Using Bowtie 2, the filtered small RNA tags were aligned against the miRNA precursor sequences of B. mori or D. melanogaster in the miRBase (the mismatch number was 0), and the known miRNAs were thus identified and were directly used for subsequent analysis. By using miRDeep 2 [82], the novel miRNAs were predicted based on the genome-matched but uncommented reads. We used DESeq2 [85] to analyse the differential expression of miRNAs, and screened the significantly differently expressed miRNA according to the threshold of P value <0.05 or |log₂FC| ≥1 and FDR <0.05.

Synthesis of dsRNA and RNAi

The dsRNAs of ecdysone pathway genes were generated by using T7 RiboMAX^{TM Large} Large-Scale RNA Production System (Promega, Madison, WI, USA). BmE cells cultured in 12-well plates were separately incubated with 5 μM 20-HE or DMSO for 12 h and were then transfected for 36 h with 5 μg of dsRNA by using 6 μL of X-tremeGENE HP DNA Transfection Reagent (Roche, Basel, Switzerland). The dsEGFP was used as the control for RNAi.

Dual-luciferase assay and site-directed mutagenesis of EcREs

The genomic DNA sequence of 10 Kb spanning the upstream of miR-275-C locus was digested by Nhe I and Bgl Ⅱ and cloned into the pGL3-Basic vector which contains the hsp70 minimal promoter. The pRL vector (Promega, Madison, WI, USA) containing the Vg promoter for Renilla luciferase was used as the reference for normalization. BmE cells were cultured in the 96-well cell plates for 24 h when the cell density was about 80%. The same amount of plasmids pLG-hsp70-Luc and pRL-Vg-Rluc (0.1 μg plasmid added to 50 μL serum-free Grace medium) were co-transfected to the cells by using X-tremeGENE HP DNA Transfection Reagent (Roche) with the ratio value of total plasmids to transfection reagent being 1:3. The antibiotic-/serum-free medium was replaced by the complete medium containing antibiotics and 10% serum at 6 h post-transfection when the 20-HE was added to the wells. The cells were incubated for another 24 h before harvesting for dual-luciferase assay. Dual-luciferase reporter assays were performed on GloMax-Multi Detection System Photometer (Promega, Madison, WI, USA) by using Dual-Glo Luciferase Assay System (Promega, Madison, WI, USA). The relative luciferase activity was the ratio between the value of firefly luciferase to the value of reference Renilla luciferase (Rluc). Site-directed mutation of EcREs was acquired by using MutanBEST Kit (Takara Bio, Shiga, Japan).

EMSA

The probes were synthesized by Sangon Biotech (Shanghai, China) and labelled with 3ʹ-terminal biotin. The final concentration of free probes was 1 μM, EcR-B protein concentration is 0.3 μg/μL and USP protein concentration is 0.45 μg/μL. All experiments were performed according to the Chemiluminescent EMSA Kit (Beyotime, Shanghai, China). All probes used for EMSA are listed in Additional Table S1.

Statistics and reproducibility

The mean cycle threshold (Ct) value of RT-qPCR was converted to relative expression level using the 2^−ΔΔCt method [86–88]. All RT-qPCR data from no fewer than three independent experiments were analysed using GraphPad Prism 7.0 (San Diego, CA, USA). Statistical analyses of the expression levels were conducted using a two-tailed unpaired Student t-test. Values with *P <0.05 mean significant difference, **P <0.01 represents very significant difference and ***P <0.001 shows extremely significant difference. All results are presented as a bar graph depicting the mean ± SEM from at least three biological repeats.

Author contributions

X. J., X.W. and L. Z. performed the biological experiments and analysed the data. T. H. and Q. Y. prepared the samples and analysed the data. S.L. conceived this project, performed the small RNA sequencing and data analysis and wrote the manuscript. All authors revised and approved the final version of the manuscript.

Disclosure of Potential Conflicts of Interest

No potential conflict of interest was reported by the authors.

Data availability

The small RNA sequencing data have been deposited into the National Center for Biotechnology Information Sequence Read Achieve database (Accession NO. PRJNA632996). All other data supporting the findings of this study are included in this published article and its additional files.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by grants from National Natural Science Foundation of China [31571334, 31071136]; and Basic and Frontier Research Project of Chongqing [cstc2014jcyjA00025].