RESEARCH ARTICLE Open Access Expression profile of cuticular genes of silkworm,

Background Insect cuticle plays essential roles in many physiological functions. During molting and metamorphosis tremendous changes occur in silkworm cuticle where multiple proteins exist and genes encoding them constitute about 1.5% of all Bombyx mori genes. Results In an effort to determine their expression profiles, a microarray-based investigation was carried out using mRNA collected from larvae to pupae. The results showed that a total of 6676 genes involved in various functions and physiological pathways were activated. The vast majority (93%) of cuticular protein genes were expressed in selected stages with varying expression patterns. There was no correlation between expression patterns and the presence of conserved motifs. Twenty-six RR genes distributed in chromosome 22 were co-expressed at the larval and wandering stages. The 2 kb upstream regions of these genes were further analyzed and three putative elements were identified. Conclusions Data from the present study provide, for the first time, a comprehensive expression profile of genes in silkworm epidermal tissues and evidence that putative elements exist to allow massive production of mRNAs from specific cuticular protein genes.


Background
Silkworm, a model for Lepidoptera, is a holometabolous insect whose developmental stages include egg, five larval instars, pupa, and adult. During molting and metamorphosis, conspicuous and relatively abrupt changes are seen in its cuticle. Insect cuticle is mainly composed of chitin nanofibres embedded in a matrix of cuticular proteins. In procuticle, a grouping of what has been called the exo-and endocuticle, cuticular proteins bound to chitin and cross-linked with the sclerotizing agents form one of the most infrangible known biological coverings [1]. Generally, cuticle plays essential roles in many physiological functions to protect the insect's body from dehydration, the invasion of pathogens, the penetration of insecticides, and physical injury [2][3][4][5].
As an important component of cuticle, hundreds of cuticular protein sequences have been identified in over 20 species of insects [6]. Many conserved motifs were identified in this data including R&R Consensus [7], CPF&CPFL [8], Tweedle [9], and others. Among them, the cuticular protein sequences containing R&R Consensus (CPR) were extensively studied in Anopheles gambiae, Drosophila melanogaster, Bombyx mori, and Apis mellifera by the annotation of genomic data [10][11][12][13]. Togawa and coworkers subsequently examined the expression profile of 156 CPR genes in A. gambiae by real-time RT-PCR and found that most of them were expressed at single or multiple periods associated with molting [14].
Our bioinformatic analysis and previous work of others have identified more than two hundred cuticular protein genes in the silkworm genome [12], indicating that the silkworm employs more than 1.5% of its estimated protein-coding genes to encode cuticular proteins. These observations led us to focus on the following three questions: 1) How many genes including cuticular protein genes are expressed in silkworm epidermal tissues? 2) Is the expression of a special set of cuticular protein genes metamorphic stage-specific? and 3) Are cuticular protein genes coordinately regulated? The sequencing of the silkworm genome along with microarray technology offered us an opportunity to investigate gene expression profiles on a large scale to answer these questions. Eleven developmental stages were selected, which ranged from day 4 of the fourth instar larva to day 8 of pupa, and microarray-based expression profile analysis of all detectable genes in silkworm epidermal tissues was performed. Our data showed that a total of 6676 genes including the vast majority of silkworm cuticular protein genes were activated in selected stages, with no correlation between expression patterns and the presence of conserved motifs. In addition, twenty-six CPR protein genes distributed on chromosome 22 were co-expressed in larval and wandering stages and three common elements were identified in the 2 kb upstream region of these coexpressed CPR genes.

Developmental expression profile of genes in epidermal tissues
In silkworm, oligonucleotide microarrays were employed to examine gene expression profiles of ten tissues on day 3 of the fifth instar larval stage, as reported by Xia et al. [15]. The microarray contained 22,987 70-mer oligonucleotides covering all predicted genes derived from the first draft silkworm genome sequence database [16]. Thereafter, the silkworm genomic database was updated and information refreshed [17,18], and one hundred and forty-seven additional oligonucleotide probes were designed for previously unpredicted genes. The current oligonucleotide microarray contained 23,134 probes, and the complete set of raw and normalized data from this study has been deposited in the Gene Expression Omnibus (GEO) repository (accession number GSE18878). We investigated the developmental expression profile of genes in epidermal tissues using the updated microarrays. A total of 6676 (28.9%) genes were transcribed in at least one selected stage in which clean cuticle tissues were easily isolated. As shown in Figure 1, the clustering analysis was carried out to analyze the expression profile. From the heat map of hierarchical clustering ( Figure 1A), we found that the selected 11 stages were clustered into two groups. Group I contained five stages including V3, V7, W1, Figure 1 Heat map of hierarchical clustering of 6676 genes expressed in epidermal tissues. Genes expressed with signal intensity more than 800 were included in this analysis. Clustering was done using Cluster 3 software (clustering type: hierarchical clustering, Distance metric: Pearson correlation). The colors in the map display the relative values of all tiles within the given 11 developmental stages; green indicates the lowest expression, black indicates the intermediate expression, and red indicates the highest expression. The numerical values give the actual values on a log 2 scale, which were associated with each color. The color scale bar is shown at the top right corner of the figure. A: Genes clustered by the given 11 developmental stages. B: Genes grouped into three clusters on the basis of the similarity of expression.
W2, and P3. The remaining six stages, IV4, IVM, W3, P5, P7, and P8 fell into group II. Stages V7, W1, and W2 of group I were closely clustered into one subgroup, indicating that genes shared a similar expression profile in epidermal tissues in these stages. By contrast, stages V3 and P3 were separated, indicating they had different gene expression profiles from one another which were distinct from the other in group I stages. Unlike group I, group II was further divided into two subgroups, one consisting of stages of IV4 and IVM and the other a tight cluster of stages W3, P5, P7, and P8.
To determine the significance of the expression profiles, K-means clustering was performed ( Figure 1B), resulting in 6676 genes divided into three clusters. Clusters I and III comprised 2450 genes and 2373 genes, respectively. Cluster II, the smaller one, contained 1853 genes. It is worthwhile to mention that genes in the three clusters showed distinctly different expression patterns. Genes in cluster I were highly expressed at stages of W3, P7, and P8, whereas genes in cluster II were highly expressed at IV4, IVM, and P5. Genes in both clusters I and II showed low expression level at stages of V7, W1, and W2. In contrast, genes in cluster III had higher expression levels at those stages.
Diverse expression patterns indicated genes of these three clusters might be involved in distinct functions and physiological pathways. Therefore, Gene Ontology (GO) annotations were determined for genes in the three clusters ( Figure 2). As shown in the GO map, genes in all clusters were mainly involved in four categories including binding, catalytic, structural molecule, and transporter activity. Nucleic acid binding, protein binding, and ion binding functions were predominated in the binding category. In the catalytic category, a large number of genes were expressed to produce three kinds of enzymes including oxidoreductases, transferases, and hydrolases. Furthermore, most of the transporters were transmembrane transporters with substrate-specific transporter activities. Cuticular protein genes and ribosomal genes were predominant in the structural molecule category. Significant differences were observed in the clustering of genes among several GO categories. For instance, the cuticular protein genes showed a different distribution in the three clusters, with the majority found in cluster I instead of clusters II and III. Moreover, molecular transducers that function in signal transduction were mainly contained in clusters I and II. By contrast, more of the nucleic acid binding associated genes were expressed like the pattern shown in cluster III. The GO annotation of 6676 genes is provided in Additional file 1.
Up-regulated and down-regulated genes in epidermal tissues from two stages before ecdysis To understand better what happens in cuticles when the silkworm initiates a molting cycle, we analyzed the gene expression ratios of IVM to IV4 and W3 to W2, trying to find genes with significant expression changes. The cut-off values were set at a ratio of more than 4 or less than 0.25, which represents up-and down-regulated expression, respectively. As shown in Additional file 2, the expression of ninety-four genes appeared to be more than 4-fold higher in the two molting phases. In this up-regulated gene list, we found sixty-eight cuticular protein genes, four juvenile hormone binding protein genes, two putative genes encoding ecdysteroid regulated proteins, six neuropeptide-like binding protein genes, and two genes related to sclerotization. In contrast, only two genes, SP1 and chitinase-related protein 1, were identified with 4-fold less abundant transcripts.
Two hundred and twenty-seven cuticular protein genes were expressed The genome sequencing project for B. mori was completed by Chinese and Japanese groups in 2008 [18]. We identified 255 putative cuticular protein genes in an updated silkworm genome database including 151 CPR genes, 4 Tweedle genes, 1 CPF genes, 4 CPFL genes, 51 CPG genes, and 44 CPH genes. Among them, two hundred and twenty-one genes were reported in a previous study [12], the corresponding accession numbers are listed in Additional file 3. Oligonucleotide probes were designed for all except six of the cuticular protein genes, these were BmorCPR29 (BR000530; EST: BP117575), BmorCPR66 (BR000567; EST: BY914249), BmorCPG1 (BR000422), BmorCPG18 (BR000439; EST: rswdd0005293.y1.abd), BmorCPG29 (BR000450; EST: BP118855), and BmorCPH29 (BR000494). In addition, members of three pairs of genes, namely BmorCPR102 (BR000603) and BmorCPR103 (BR000604), BmorCPR119 (BR000620) and BmorCPR120 (BR000621), and BmorCPR124 (BR000625) and BmorCPR150 (GU070697), shared high similarity in the open reading frame (ORF) region, and no unique probes were to distinguish them. Similarly, a single probe was designed to represent three genes, BmorCPG26 (BR000447), BmorCPG27 (BR000448), and BmorCPG28 (BR000449), which shared high identity. In all we obtained 244 oligonucleotide probes representing 249 cuticular protein genes. According to the microarray data, as shown in Figure 3, the expression signals of 227 cuticular protein genes were detected in at least one stage.

Developmental expression profiles of cuticular protein genes
Hierarchical clustering ( Figure 3A) by both samples and genes was performed to examine the expression profiles of cuticular protein genes. The heat map of hierarchical clustering placed 16 developmental stages into two groups. Eight stages, I3, III3, IV4, V3, V7, W1, W2, and P3 were clustered into group I; the remaining stages were in group II. A marked boundary was observed between the two groups. Genes in group I and group II displayed distinctly different expression profiles. Many more cuticular protein genes were expressed in stages which were clustered in group II than in group I. Kmeans clustering performed to determine the expression patterns of cuticular protein genes yielded five gene clusters with distinct expression patterns, as shown in Figure 3B. Cuticular protein genes contained in cluster A were widely expressed among all 16 developmental stages. Clusters D and E comprised some larval-and pupal-specific expressed cuticular protein genes. Among them, BmorCPR53 was only expressed at four stages including I3, III3, IVM, and V3. The northern blot data for this gene ( Figure 4A) showed a similar expression pattern at stages I3, III3, and IVM, whereas a weak signal was detected at stage V3. Based on the microarray data, BmorCPR125 was identified as a pupalspecific cuticular protein gene expressed at stage P5. The northern blot data ( Figure 4B) indicated that except for P5, it was also expressed at a lower level at stages P3, P7, and P8. In clusters B and C, the cuticular protein genes showed a noticeable coordinate expression pattern. However, genes in cluster C had lower expression levels in pupal stages than those in cluster B. Although well organized expression patterns were shown in clusters B and C, exceptions were still found at several stages.
In silkworm, cuticular protein genes are divided into five families based on conserved motifs. To understand better the expression profiles of cuticular protein genes bearing particular motifs, the number of genes in each cluster and the families to which they belong were summarized. As shown in Table 1, each cluster consisted of more than two families. 81.4% (44/54) RR-1 cuticular protein genes were placed into clusters A, B, and E. Among them, 26% (14/54) RR-1 genes found in cluster A were widely expressed in all tested stages. In contrast, a total of 84 (94.5%) RR-2 genes were found in clusters B, C, D, and E. Only 5 (5.5%) RR-2 genes were identified in cluster A. RR-3 genes were found in clusters A and E. Three of four Tweedle cuticular protein genes shared the same expression pattern as cluster B and only one (BmorCPT1, BR000650) was found in cluster A with a wide expression pattern. A CPF gene was identified in silkworm but its expression level was too low to be detected. Four CPFL genes were found in clusters A, B, and C. For low complexity cuticular protein genes, 17 CPGs and 19 CPHs were found in cluster A. Altogether, cluster A contained mainly RR-1, CPG, and CPH genes, whereas majority of genes in clusters B and C were RR-1 and RR-2 genes.
Twenty-six RR genes distributed in chromosome 22 were strictly co-expressed in larval and wandering molting stages Silkworm cuticular protein genes showed differential patterns of expression, and no obvious correlation was found between expression patterns and the presence of conserved motifs. However, as visualized in the heat map of K-means clustering ( Figure 3B), cuticular protein genes in clusters B and C exhibited distinctive expression pattern. Genes in cluster B had abundant transcripts not only in molting stages but also throughout the four pupal stages. The same expression pattern was found for genes in cluster C from stages I3 to W3. But, the number of cuticular protein genes expressed in pupal stages in cluster C appeared to be lower than in cluster B. Considering their similar patterns of expression in larval and wandering stages, we combined the 54 cuticular protein genes in cluster B with the 47 cuticular protein genes in cluster C for the subsequent analysis. First the chromosomal locations of these 101 cuticular protein genes were examined. As shown in Figure 5, 48.5% (49/101) of these genes were distributed in chromosome 22, including 9 RR-1 and 40 RR-2 genes. Stringent QT (quality threshold) clustering algorithm analysis and manual verification were carried out to analyze the expression profiles of these 49 RR cuticular protein genes. Twenty-six of these 49 genes, as shown in Figure 6, were strictly expressed in larval molting stages and in pharate pupae. No expression signals were detected in stages I3, III3, IV4, and V3 to W2. Such result was also supported by the northern hybridization. As shown in Figure 4C, the microarraybased expression profile of BmorCPR96, one of these 26 RR genes, was completely consistent with the results of the northern blots.
Upstream promoter regions of the twenty-six RR genes shared common putative regulatory elements As noted above, the twenty-six RR cuticular protein genes distributed in chromosome 22 were co-expressed along the larval and wandering stages. To understand the underlying mechanisms of these consistent expression profiles, we analyzed the transcriptional regulatory regions located 2 kb upstream of the putative transcription start sites in the twenty-six genes. The MEME algorithm, which was used to discover similar sequence elements in the promoter regions, identified three putative regulatory elements shared by the 26 RR genes (Figure 7). Elements I, II, and III ( Figure 7A-C) contained 41, 34, and 21 nucleotides, respectively. TOMTOM motif comparison tool was employed to compare these three elements with the known motifs. No motifs were found in element II ( Figure 7B). Only one NKx2-2 transcription factor binding site was found in the middle of element I ( Figure 7A). Nine transcriptional factor  binding sites including Hb, Cdx, Nkx6-1, BR-C Z1, BR-C Z2, PLZF, HNF3, FOXP1, and TCF showed similarity with element III ( Figure 7C). The alignment match logos for elements I and III are illustrated in Additional file 4. In addition, FIMO (Find Individual Motif Occurrences) was applied to determine whether the elements identified in this study were present in the upstream regions of the cuticular protein genes in Anopheles gambiae [19]. Homologues of element II were found upstream of four CPR genes (Additional file 5). Notably, element III was found in the upstream sequences of 59 cuticular protein genes, of which 50 members were CPR genes of A. gambiae (Additional file 6). The TESS (Transcription Element Search System) server was employed to seek known binding sites for transcription factors from the TRANSFAC database in these 2 kb upstream regions [20]. In all the twenty-six upstream sequences, binding sites for at least one of the

Discussion
Cuticle is generally considered as a protective cover for insects [2][3][4][5]. A sclerotized and tanned integument layer consisting of chitin and particular cuticular proteins synthesized and secreted by the epidermal cells play vital roles in insects' lives [1]. In the present study, we monitored gene expression profiles of silkworm epidermal tissues isolated from eleven developmental stages ranging from day 4 of fourth instar larvae to day 8 pupae. In these stages, the expression of 6676 genes was detected. This represents the first global gene expression analysis of Lepidoptera insect epidermal tissues, which provides important functional insight. Chitin and cuticular proteins are the major components of insect cuticle [21,22]. It is noteworthy to mention that the gene coding for chitin synthase, a critical enzyme for chitin synthesis [23][24][25], had peak expressions at the same stages in which the majority of cuticular protein genes were also highly expressed, indicating that the processes of cuticle formation were activated concurrently. Our analysis detected the expressions of seven yellow protein genes, several of which are involved in cuticle pigmentation [26][27][28]. Furthermore, two wellstudied melanin syntheses related genes tyrosine hydroxylase (TH) and dopa decarboxylase (DDC) [29,30], were strongly induced when molting was initiated. These results suggested that many genes involved in cuticle tanning were activated in epidermal tissues before ecdysis. In addition, four juvenile hormone binding protein genes (JHBP) and two ecdysteroid regulated proteins (ERP) were highly expressed before molting. Previous studies showed that JHBPs protect juvenile hormone from degradation [31,32], and the mRNA of Manduca sexta ERP20, a homolog of ERP, is abundant in epidermal tissues during molting [33]. Although the functions of these genes have not yet been clarified, it can be assumed that JHBPs and ERPs are involved in the molting process. The expression of six neuropeptide-like protein genes (NPLP) was also stimulated prior to molting. In silkworm, EST evidences showed that neuropeptide-like protein genes were present in epidermis, but their functions were unclear. We speculate that these NPLPs might participate in the regulation of molting. In contrast to numerous up-regulated genes, only two down-regulated genes were detected. SP1 is considered as a storage protein in hemolymph and ovary associated with the development of Lepidoptera insects [34,35]. Another down-regulated gene, chitinase-related protein 1 showed high identity to chitinase. However, without a glutamate residue in the catalytic sites chitinase-related protein 1 did not hydrolyze chitin [36].
The results of the GO analysis presented in Figure 2 showed that genes expressed in epidermal tissues were involved in different pathways, which indicated that cuticle functions not only as a protective cover against external threats but also as a place for active metabolism. In our data, the expressions of silkworm chitinase and genes related to the degradation of proteins were detected. Moreover, numerous transporter genes were expressed in epidermal tissues. One of the conclusions made from these observations was that the synthesis of new cuticle was concomitant with the degradation of the old cuticle. Interestingly, the majority of heat shock protein genes of silkworm were expressed in epidermal tissues throughout the developmental stages, suggesting their role as molecular chaperone was necessary for synthesis and degradation of epidermal proteins [37][38][39].
Silkworm is a holometabolous insect that develops from larva to pupa, and then pupa to adult. In order to grow and change the appearance, it must molt and cast its old cuticles. So, during molting and metamorphosis dramatic changes occur in cuticles. As the major component of cuticle, cuticular proteins are obvious choices for studying development and metamorphosis. Here, the expressions of cuticular protein genes at 16 stages ranging from day 3 of the first instar larvae to day 8 of pupae were investigated. Of the 244 available probes for cuticular protein gene, 227 genes (93%) had expression signals and were distributed in diverse gene expression patterns. One of goals of this research was to learn which families of cuticular protein genes were expressed in each metamorphic stage. Our data clearly showed no correlation between the expression profiles of cuticular protein genes and the presence of conserved motifs. This result is consistent with what Togawa et al found in the expression profile of putative CPR cuticular protein genes of A. gambiae by qRT-PCR [14]. On the other hand, our data explicitly revealed massive expression of many cuticular protein genes at molting stages, when the silkworm was building its cuticle. A reasonable explanation for this observation is that these cuticular protein genes were transcribed and immediately translated to proteins participating in the construction of cuticles. Okamoto and colleagues reported that massive numbers of ESTs for BmorCPR32, BmorCPR39 and BmorCPG3 were only isolated during the fourth larval molt, whereas more ESTs of BmorCPR41 and BmorCPR46 genes were identified in the intermolt stage [40]. In our microarray, BmorCPR32, BmorCPR39, BmorCPR41, and BmorCPR46 showed similar patterns to the Okamoto et al.'s [40]. However, BmorCPG3 showed a high level at both molting and intermolt stages in our data, which was somewhat different from the patterns of Okamoto et al. This can be explained by the differences in sensitivity between EST sequencing and microarrays. Our data also revealed that although the cuticles of silkworm larva, pupa, and adult were distinct, a number of cuticular protein genes were commonly expressed at all stages, indicating that the properties of cuticle depend on the amount of cuticular proteins, and their spatial distribution, the degree of sclerotization and tanning [1], rather than simply the types of cuticular proteins present.
The silkworm used more than 1.5% of its total estimated genes to encode the cuticular proteins, and at each molting stage massive cuticular protein genes were activated. Another goal of this research was to learn how the silkworm regulates the transcription of cuticular protein genes. Coordinate regulation using the same conserved motifs or localizing genes on the same chromosomes might be the simplest and most efficient mechanism. Our data revealed no evidence for coordinate regulation in the majority of cuticular protein genes. However, twenty-six RR cuticular protein genes coexpressed in larval stages were found to be distributed in chromosome 22. Togawa and colleagues also found that a portion of A. gambiae cuticular protein genes distributed in a narrow chromosome region showed highly similar expression patterns [14]. Three similar elements were further identified in the 2 kb upstream promoter regions of the 26 silkworm cuticular genes on chromosome 22. Notably, the similar sequences of element II and III were found in the upstream regions of A. gambiae CPR genes, and a binding site of transcriptional factor NKx2-2 was identified in the middle of element. Nkx2-2 acts cooperatively with Pax6, whose function was conserved from invertebrate to vertebrate for dorsal and ventral patterning [41,42]. A recent study showed that Broad-Complex and βFTZ-F1 positively regulated the transcription of wing cuticular protein gene in silkworm [43]. Interestingly, element III found in the present study was much longer than the binding sequence of any Broad Complex isoform. The three elements might account for the coordinate regulation of the 26 cuticular protein genes. Future analysis of these elements will increase the understanding of transcriptional regulation of cuticular protein genes.
Besides, many binding sites for known transcription factors were identified in the upstream regions of these cuticular protein genes. Binding sites for EcR and E74A were discovered in the upstream regions of 25, and 23 of the 26 cuticular protein genes, respectively. EcR formed heterodimer with USP (ultraspiracle protein) to function as the receptor of ecdysone and E74A was known as the early genes induced by ecdysone and functioned as one of ecdysone signal transducers [44][45][46]. Binding sites for bZIP, bHLH, and C/EBP were discovered in upstream sequences of all 26 cuticular protein genes. Previous study showed that transcription factors with bZIP and bHLH domain played roles in transcriptional regulation of neuropeptides and peptide hormone [47], which were engaged in the regulation of insect molting. C/EBP was the factor bound to the promoters of silkworm chorion genes and regulated their precise spatial and temporal expression [48,49]. In silkworm, POU factors played a critical role in transcriptional regulation of neuropeptides and silk genes [50][51][52].
Generally, the expression of cuticular protein genes is regulated by two hormones, ecdysone and juvenile hormone. Ecdysone induces the transcription of the primary-response genes [46,53], including Broad Complex genes, E74 isoforms, βFTZ-F1, and orphan nuclear receptors. Primary-response genes activate and regulate the transcription of the secondary-response target genes. Although a number of juvenile hormone binding proteins have been identified [31,32], little was known about the juvenile hormone receptor. Cuticular protein genes are usually considered to be located downstream of the hormones action hierarchy [44]. The transcription of cuticular protein genes could be activated by various signaling pathways, which provides a possible explanation for the diverse expression patterns demonstrated in this study.

Conclusions
This study describes the expression profile of genes in silkworm epidermal tissues for the first time. Microarray data showed activation of a total of 6676 genes involved in various functions and physiological pathways. More than 93% cuticular protein genes were expressed in selected developmental stages, displaying diverse expression patterns.
The majority of cuticle proteins showed no evidence of coordinate regulation as a function of common cuticle protein motifs. However, 26 RR genes distributed in chromosome 22 were co-expressed at larval and wandering stages, and three putative elements were identified in the 2 kb upstream region of these 26 RR genes. Extensive expression data and the analysis of transcriptional factor binding sites provided novel insights into the functional coordination of these genes.

Microarray design and construction
Based on a previously designed silkworm oligonucleotide microarray [15], we added 147 novel oligonucleotide probes for the cuticular protein genes that were not present in the original array. All the probes were designed by CapitalBio Corporation (Beijing, China) and were synthesized by MWG Biotech (Ebersberg, Germany). The microarray slide contained 48 blocks, each with 22 rows and 23 columns. Five housekeeping genes and eight yeast intergenic sequences were dotted in one block as positive and external controls, respectively. Dual channel microarray hybridization was performed with a Cy3-labeled control sample and Cy5-labeled test sample. Total RNAs extracted from the whole body of silkworm at day 3 of the fifth instar larvae served as a normalization control for data analysis.

Silkworm strain and reagents
Silkworm larvae (p50 strain) maintained at the Institute of Sericulture and System Biology (Southwest University, China) were reared on mulberry leaves at 25°C~26°C. Silkworms grow through five instars until cocoon spinning which begins at the end of the fifth instar larva day 7. After spinning for three days, silkworms develop into the pupal stage, which takes about 10 days, followed by emergence from the cocoon, mating and egg lay. We selected sixteen time points around the molting phases, ranged from the first instar larva day 3 to pupa day 8. Considering the small body size, we used whole larval bodies from the first to the third instar larvae to isolate total RNA. From the fourth instar larva day 4 to pupa day 8, we collected epidermal tissues to isolate total RNA. TRIzol reagent was obtained from Invitrogen (Carlsbad, CA, USA). Reverse transcriptase was made in Promega (Madison, WI, USA). ECL direct nucleic acid labeling and detection system was from GE Healthcare (Buckinghamshire, UK).

RNA isolation, amplification, labeling and array hybridizations
Total RNAs were isolated using TRIzol reagent and further purified using a NucleoSpin RNA clean-up kit (Macherey-Nagel, Germany). The amplification and labeling of mRNA were performed as described in previous studies [15,54]. Five micrograms of total RNA were primed with 1 μl of 100 μM primer containing T7 RNA polymerase promoter sequence at 70°C for 10 min, then reversed transcribed at 42°C for 2 h in the presence of 200 U CbcScript (CapitalBio Corp, China). The second strand of cDNA was synthesized at 16°C for 2 h in the presence of RNaseH and DNA polymerase. cRNA was synthesized by T7 Enzyme Mix (CapitalBio Corp, China) using the cDNA template. 2 μl of cRNA were primed with 1 μl random primer at 65°C for 10 min, then reverse transcribed at 25°C for 10 min and 37°C for 1.5 h in the presence of CbcScript II (CapitalBio Corp, China). The Cy3-and Cy5-dCTP double-stranded cDNA was labeled using a CapitalBio cRNA Amplification and Labeling Kit (CapitalBio, Beijing, China). Cy5-dCTP or Cy3-dCTP were added at a final concentration of 120 μM of each dATP, dGTP, and dTTP and 60 μM dCTP and 40 μM Cy5-dCTP for test samples. For reference samples, Cy3-dCTP was used. The Cy3and Cy5-dCTP double-stranded cDNA was dissolved in 80 μl hybridization solution containing 3 × SSC, 0.2% SDS, 5 × Denhart's, and 25% formamide. The slides were covered with a LifterSlip™ coverslip (Erie Company, Portsmouth, NH, USA) and hybridized in a closed chamber at 42°C over-night. After hybridization, slides were washed three times in 0.2% SDS, 2 × SSC at 42°C for 5 minutes and three times in 0.2 × SSC at room temperature for 5 minutes before signal scanning.

Microarray data processing and analysis
The slides were scanned with a confocal LuxScan scanner (CapitalBio Corp.) and the raw data were extracted using LuxScan™ 3.0 software (CapitalBio Corp.). For dual-channels microarray data, the scanning setting for Cy3 and Cy5 channels were balanced by visual inspection of the external control spots. The LOWESS (Locally Weighted Scatterplot Smoothing) method was used to normalize the dual channel data using all the signals from the Cy3-labeled sample. The ratios of signal intensity of test and control samples were used to perform clustering analysis. The one with a fluorescence intensity higher than 800 after subtracting the background was considered as an expressed gene since the signal greater than that detection level was more reliable. The expression of a cuticular protein gene was defined by the ratio of the original signal intensity divided by 800. The Xfold values were used in the subsequent clustering analysis to display the expression of cuticular protein genes at different developmental stages. HCL (Hierarchical Clustering) analysis was carried out using both Cluster 3.0 software and Mev software (version 4.2.01) [55,56]. In addition, Cluster 3.0 software was used for K-mean clustering analysis. Mev software was also used for QT (quality threshold) clustering. The parameter setting for clustering analysis was based on the distance metric of the Pearson correlation and the average linkage method. TreeView software was used to display heat map of clustering results. Gene ontology analysis was performed at the BGI WEGO website [57].

Computational identification of putative regulatory elements
The MEME algorithm (Multiple Expectation maximization for Motif Elicitation) was used to identify common elements present in the 2 kb promoter regions upstream of the transcription start sites of cuticular protein genes [58]. TOMTOM motif comparison tool was used to compare the elements identified in this study to known motifs [59]. In TOMTOM analysis, the TRANSFAC database was selected and the Pearson correlation coefficient was employed to survey the Motif Column Comparison Function. FIMO (Find Individual Motif Occurrences) was applied to search for whether the identified regulatory elements existed upstream of other genes [19]. In FIMO analysis, the Anopheles_gambiae_EnsEMBL_upstream database was selected as the reference and the p-value output threshold was set at l × e-5. TESS (Transcription Element Search System) was applied to search the binding sites for known insect transcription factors from the TRANSFAC database [20].

Northern hybridization
Northern hybridization was performed to confirm the microarray data. The sequences of cuticular protein genes used to design the hybridization probes were obtained from the Silkworm Genome Database [17]. DEPC water was employed to prepare the related solutions and to clean the associated equipments. Five micrograms total RNA per sample was loaded to perform denaturing formaldehyde gel electrophoresis. The transfer of RNA from gel to Hybond+ (GE) membrane was completed in 2 hr by using Transfer Equipment (Amersham Biosciences). All reagents used in prehybridization, probes labeling, hybridization and signal detection were provided by the Amersham ECL Direct Nucleic Acid Labeling and Detection Systems (GE Healthcare, Cat No: PRN 3001), which is based on enhanced chemiluminescence. The optimized temperature of hybridization mixture 42°C was adopted to protect the activity of the horseradish peroxidase. The cDNA probes, which were labeled with the enzyme horseradish peroxidase, were completely denatured to single-strand form to hybridize the target RNA on Hybond+ membrane. Membranes were washed to a stringency of 0.1 × SSC, and labeling and detection were carried out according to the manufacturer's instructions.
Additional file 1: The GO (Gene Ontology) annotation of 6676 probes. This file contains the GO number of 6676 probes expressed in our microarray. It has multiple columns-the first of which is the probe ID and the rest with GO numbers. Some probes representing new genes have no GO number. Click here for file [ http://www.biomedcentral.com/content/supplementary/1471-2164-11-173-S1.XLS ] Additional file 2: The list of ninety-four up-regulated and two down-regulated genes from two stages before ecdysis. This file contains the ninety-four up-regulated and two down-regulated genes from two stages before ecdysis. It has a table which contains five columns. The first column is the function terms which these genes belong to; the second to the fifth with names of these genes, symbol or homolog, accession number and E-value, respectively. Abbreviations I3: day 3 of the first instar larva; IM: larval molting stage from 1 st to 2 nd instar; IIM: larval molting stage from 2 nd to 3 rd instar; III3: day 3 of the third instar larva; IIIM: larval molting stage from 3 rd to 4 th instar; IV4: day 4 of the fourth instar larva; IVM: larval molting stage from 4 th to 5 th instar; V3: day 3 of the fifth instar larva; V7: day 7 of the fifth instar larva; W1: day 1 of the wandering phase; W2: day 2 of the wandering phase; W3: day 3 of the wandering phase; P3: day 3 of the pupa; P5: day 5 of the pupa; P7: day 7 of the pupa; P8: day 8 of the pupa; CPF: cuticular protein with 44-amino acid motif; CPFL: CPF-like protein; CPG: glycine-rich cuticular protein; CPH: hypothetical cuticular protein; CPR: cuticular protein with the R&R Consensus; CPT: cuticular protein with a Tweedle motif; GO: gene ontology; QT clustering: quality threshold clustering; MEME: multiple expectation maximization for motif elicitation (motif discovery tool); TOMTOM: motif comparison tool for searching a database of motifs with a given query motif; FIMO: find individual motif occurrences; TESS: transcription element search system.