Transcriptome analysis of integument differentially expressed genes in the pigment mutant (quail) during molting of silkworm, Bombyx mori.

In the silkworm Bombyx mori, pigment mutants with diverse body colors have been maintained throughout domestication for about 5000 years. The silkworm larval body color is formed through the mutual interaction of melanin, ommochromes, pteridines and uric acid. These pigments/compounds are synthesized by the cooperative action of various genes and enzymes. Previous reports showed that melanin, ommochrome and pteridine are increased in silkworm quail (q) mutants. To understand the pigment increase and alterations in pigment synthesis in q mutant, transcriptome profiles of the silkworm integument were investigated at 16 h after head capsule slippage in the fourth molt in q mutants and wild-type (Dazao). Compared to the wild-type, 1161 genes were differentially expressed in the q mutant. Of these modulated genes, 62.4% (725 genes) were upregulated and 37.6% (436 genes) were downregulated in the q mutant. The molecular function of differently expressed genes was analyzed by Blast2GO. The results showed that upregulated genes were mainly involved in protein binding, small molecule binding, transferase activity, nucleic acid binding, specific DNA-binding transcription factor activity and chromatin binding, while exclusively down-expressed genes functioned in oxidoreductase activity, cofactor binding, tetrapyrrole binding, peroxidase activity and pigment binding. We focused on genes related to melanin, pteridine and ommochrome biosynthesis; transport of uric acid; and juvenile hormone metabolism because of their importance in integument coloration during molting. This study identified differently expressed genes implicated in silkworm integument formation and pigmentation using silkworm q mutant. The results estimated the number and types of genes that drive new integument formation.

Introduction

In insects, widespread pigmentation is involved in mimicry [1], [2], sexual selection [3], [4], thermoregulation [5], [6], cuticle hardening [7], [8], territoriality [9], immunity [10], [11] and preventing ultraviolet damage [12]. Body coloration mainly comes from melanin, ommochromes and pteridines in insects. The molecular mechanisms of pigment pathways were clarified using mutants such as vermilion (v), scarlet (st), brown (bw), white (w), sepia (se) in Drosophila [13]–[15]. The silkworm Bombyx mori was domesticated more than 5000 years and has multiple color mutants. Color mutations are reported throughout the silkworm life cycle in larvae [16], [17], pupae [18], moths [19], [20] and eggs [21], so they are ideal material for elucidating the molecular mechanism of pigment pathways. The larval body color in silkworm is the result of interactions among melanin in the cuticle, xanthommatin, sepialumazine and uric acid in the epidermis [22]. In silkworms, melanism (mln) was used to determine that arylalkylamine-N-acetyl transferase is important in the melanin pathway in Lepidoptera [19]. The mutant gene in the red egg (re) mutant suggests that the mutant product might transport cysteine or methionine into pigment granules, since the mechanisms after 3-hydroxykynurenine is incorporated into pigment granules are still largely unknown. If correct, this would be useful in elucidating the mechanism of the late steps of the ommochrome synthesis pathway [23].

However, these reports mainly focused on single pathway, the relationship of melanin, ommochrome and pteridine pathway is very little. It was reported that melanin, ommochrome and pteridine biosynthesis are increased in q mutant [22]. Thus, it is a good model for studying the interaction of different pigments. The quail (q) mutant has light crescent marking and star spot marking by melanin along the head-to-tail axis, with the star spot marking more dispersive than in wild-type (Dazao). The q mutant results in more small black spots on the dorsal surface of the integument than in Dazao; the larval body color of q mutant is yellowish-brown with light red on the newly molted fifth instar larva that gradually becomes lighter after eating mulberry leaf (Figure 1).

Figure 1

Appearance of quail mutant (q/q) and the wild-type strain (Dazao) at different developmental stages.

(A) Right, q/q mutant; left, wild-type Dazao at fourth molting larval stage, newly molted fifth stage larvae, and day 3 fifth instar insects. Red square, crescent marking; green square, star marking. Bars = 1 cm. (B) Magnified images of larval crescent marking and star spot marking at different developmental stages.

Although some biochemical and molecular studies of q mutant have been performed [22], [24], [25], the mechanisms that result in body color are not clear. In this study, we used a transcriptome approach to investigate gene networks involved in cuticle formation and pigment synthesis in silkworms. Since cuticular pigmentation occurs 16–18 h after head capsule slippage (HCS) at the fourth molt in P. xuthus and the fourth molting stage is nearly equivalent in B. mori and P. xuthus [26], we performed a comprehensive transcriptome analysis of the integument at 16 h after HCS of the fourth molt. Our analysis targeted two main classes of genes: 1) differently expressed genes involved in melanin, ommochrome and pteridine biosynthesis, and 2) differently expressed genes that might function in cuticle formation. Our results determined differently expressed genes involved in pigment biosynthesis pathways, transport of uric acid, juvenile hormone (JH) metabolism, chitin metabolism, cuticular proteins and nuclear receptors. These data add our understanding of the relationship of different pigment pathways and the mechanisms underlying cuticle formation in B. mori.

Materials and Methods

Silkworm strains and tissue collection

Silkworm q mutant strains and Dazao were obtained from the Silkworm Gene Bank of Southwest University, Chongqing, China. All larvae were reared on fresh mulberry leaves at 25°C under long-day conditions (16 h light: 8 h dark). Staging of molting period was based on head capsule slippage (HCS) timing. The whole integument was collected on ice from 3 larvae of each strain just at 16 h after HCS in the fourth molt and other tissues were cleaned out in PBS solution. The mixture of three integuments were washed quickly with water, blotted on filter paper, and then stored at −80°C until RNA extraction.

RNA preparation and Illumina RNA-seq

Total RNA was prepared using TRIzol (Invitrogen) from the mixture of three whole integument of each strain. RNA integrity was analyzed using an Agilent bioanalyzer (Agilent Technologies 2000). mRNAs were purified with the oligo(dT) magnetic beads, fragmented and used to synthesize cDNA following the TruSeq RNA Sample Preparation v2 Guide(Illumina). Sequencing adaptors were ligated to cDNA fragments by PCR amplification. Sequencing raw data were generated using an Illumina HiSeq 2000 system (Illumina, USA). The raw data presented in this publication have been deposited in NCBI Short Read Archive (http://www.ncbi.nlm.nih.gov/sra/) and are accessible through SRA accession number: SRR1177794 and SRR1177795.

RNA-seq data analysis

The quality analysis of RNA-seq data was determined using a CLC Genomics Workbench 5.5 following the manufacturer's instructions. The silkworm genome sequence was obtained from the SilkDB (http://silkworm.swu.edu.cn/silkdb/). Reads were mapped to the silkworm genome sequence using CLC Genomics Workbench 5.5. Gene expression levels were calculated using the reads per kb per million reads (RPKM) method [27]. P value<0.05, false discovery rate (FDR)≤0.001 and RPKM>5 were thresholds for gene expression. The RPKM of q was divided by the RPKM of Dazao at 16 h after HCS, and the ratio was used to determine differentially expressed genes. Ratio higher than 2.0 or lower than 0.5 was the thresholds for prominent changes. Genes with ratios higher than 2.0 were regarded as upregulated; genes with ratios lower than 0.5 were regarded as downregulated. Differently expressed genes were aligned by BLASTp and results were analyzed for mapping, annotation, enzyme code and combined graph using BLAST2GO software (version 2.6.6, http://www.blast2go.org) with default settings.

Quantitative reverse transcriptase PCR

RNA was isolated using TRIzol (Invitrogen) from the mixture of three whole integument of each strain, purified with absolute alcohol and treated with with DNase. RNA was reverse transcribed using PrimeScript RT reagent kits (RR037A, Takara) following manufacturer's instructions. Quantitative reverse transcriptase PCR (qRT-PCR) was performed using an ABI7500 real-time PCR system (Applied Biosystems) with a two-step reaction protocol of 40 cycles of 94°C for 3 s and 60°C for 30 s, followed by dissociation for quality control. qRT-PCR mixtures was performed in 15 μL with 1.5 μL of 2 μM specific primers, 7.5 μL SYBR Premix Ex Tag II (Tli RNaseH Plus, Takara) (2×), 0.3 μL ROX Reference Dye II (50×) and 4 μL cDNA. mRNA quantity of each gene was calculated with the 2−△△Ct method [28] and normalized to the silkworm housekeeping gene, ribosomal protein L3 (RPL3, BGIBMGA013567). Ratios after normalization were expressed as fold-change compared with samples from the Dazao integument. The primers are listed in Table S1.

Results and Discussion

General transcription patterns

Whole genome mRNA sequencing was used to monitor global changes in gene expression in q mutant and Dazao. RNA-seq was performed on integument samples at 16 h after HCS. Quality analysis indicated good RNA-seq data (Figure S1). A total of 33,815,986 paired-end reads (100 bp) were obtained for Dazao, of which 19,724,505 (58.3%) mapped uniquely to the silkworm genome; In q mutant, 33,585,319 (61.3%) of 54,754,452 total paired-end reads matched to the silkworm genome. Based on RPKM method [27], the level of 14,623 annotated silkworm genes were determined (Table S2). Applying a FDR cut-off of 0.1% and P value<0.05 resulted in 2905 genes for analysis. To ensure that genes were expressed in q and Dazao, RPKM>5 in q or Dazao was applied as a threshold. Genes were considered not expressed when the RPKM value was less than 5 in both q and Dazao. Differentially expressed genes between q and Dazao were identified by expression ratios higher than 2.0 or lower than 0.5. This resulted in 1161 genes that were differentially expressed: 725 upregulated and 436 downregulated genes in q (Table S3).

We analyzed whether transcripts in q mutant were enriched in specific molecular functions using Blast2GO (Table S4). Most upregulated gene transcripts in the q strain were in protein binding, small molecule binding, transferase activity, nucleic acid binding, specific DNA-binding transcription factor activity, chromatin binding; specific down-regulation of genes in q were in oxidoreductase activity, cofactor binding, tetrapyrrole binding, peroxidase activity and pigment binding (Figure 2).

Figure 2

Functional categories of differentially expressed genes in quail mutant.

Differentially expressed genes were categorized by molecular function. Bars, number of genes. Upregulated or downregulated genes based on Blast2GO analysis were tested for enrichment in specific functional categories.

Validation of data reliability by qRT-PCR

To verify the RNA-seq results, the expression of two silkworm housekeeping genes, RPL3 and glyceraldehyde-3-phosphate dehydrogenase (GAPDH, BGIBMGA007490), were examined in the RNA-seq data. Almost no difference in expression between q and Dazao was observed. Previous studies found that the expression of GTP-CH I was higher in q mutant [22]. Silkworms have two isforms of GTP-CH I, a and b [18]. The RNA-seq results showed significant upregulation of GTP-CHI a (17.4-fold) and b (12.9-fold) in q mutant relative to Dazao, similar to previous studies. To corroborate the RNA-seq results, we selected all differentially expressed genes listed in Table 1 and top ten in up-/down-regulated differentially expressed cuticular protein genes in q mutant relative to Dazao for secondary confirmation using qRT-PCR. The similarity of changing trend was 92.6% in Table S5. Figure 3 show that qRT-PCR confirmed the reliability of RNA-seq data. Fold-changes for all genes as determined by qRT-PCR were similar to fold-changes observed by RNA-seq.

Table 1

Differentially expressed genes involved pigment pathway, transporting of uric acid, JH metabolism, chitin metabolism and differentially expressed nuclear receptor genes in q mutant.

		Fold (q/Dazao)		Concordant
Genes	Name	Fold a	Fold b	up	down
BGIBMGA003866	PAH	2.04	1.34	Yes
BGIBMGA000031	ebony	9.36	1.83	Yes
BGIBMGA001149	yellow	0.41	0.31		Yes
BGIBMGA014032	yellow-f2	0.09	0.03		Yes
BGIBMGA010917	yellow-f4	0.09	0.25		Yes
BGIBMGA007253	yellow-e	2.61	1.48	Yes
BGIBMGA007255	yellow-h2	2.73	5.74	Yes
BGIBMGA007254	yellow-d	2.04	3.16	Yes
BGIBMGA003918	yellow-f4-2	3.32	3.27	Yes
BGIBMGA014224	yellow x	3.15	0.17	No
BGIBMGA001235	GTP-CH I a	17.40	1.62	Yes
BGIBMGA008134	GTP-CH I b	12.94	1.67	Yes
BGIBMGA007856	KFM	5.26	9.79	Yes
BGIBMGA006740	PHS	2.70	1.49	Yes
BGIBMGA007285	OBP1	4.95	1.40	Yes
BGIBMGA007286	OBP2	2.45	2.98	Yes
BGIBMGA002922	Bm-w-3	2.86	3.53	Yes
BGIBMGA003864	Bm-os	2.24	1.50	Yes
BGIBMGA002581	Bm-ok	4.12	2.90	Yes
BGIBMGA010392	JHAMT	64.51	819.88	Yes
BGIBMGA000772	JHE	13.26	4.22	Yes
BGIBMGA013930	JHEH	2.31	2.37	Yes
BGIBMGA008815	JHDK	2.27	2.37	Yes
BGIBMGA013971	JH-inducible protein	2.48	1.32	Yes
BGIBMGA011646	β-N-acetylglucosaminidase	2.27	0.80	No
BGIBMGA004221	GPI	3.08	3.70	Yes
BGIBMGA007517	GFAT	4.91	1.86	Yes
BGIBMGA001609	UDPAP	3.66	1.74	Yes
BGIBMGA006767	EcR	4.46	1.10	Yes
BGIBMGA002964	HR38	6.50	1.16	Yes
BGIBMGA000716	βFTZ-F1	3.02	1.05	Yes
BGIBMGA007914	HR39	2.99	3.04	Yes
BGIBMGA007970	E74A	2.37	0.73	No

Fold: Fold change (q/Dazao) of gene expression in RPKM.

Fold: Fold change (q/Dazao) of gene expression in qRT-PCR.

Abbreviations: PAH, phenylalanine hydroxylase; GTPCH I, GTP cyclohydrolase I; KFM, kynurenine formamidase; PHS, phenoxazinone synthetase; OBP, Ommochrome-binding protein; JHAMT, JH acid methyltransferase; JHE, JH esterase; JHEH, JH epoxide hydrolase; JHDK, JH diol kinase; GPI, Glucose-6-phosphate isomerase; GFAT, Glutamine:fructose-6-phosphate aminotransferase; UDPAP, UDP-N-acetylglucosamine pyrophosphorylase.

Figure 3

Relative gene expression.

(A) Sequencing results. (B) Real-time quantitative PCR results. RPKM, reads per kb per million reads; PAH, phenylalanine hydroxylase; GTPCH I, GTP cyclohydrolase I; KFM, kynurenine formamidase; PHS, phenoxazinone synthetase.

Differentially expressed genes in melanin synthesis

Based on enzymes reported for D. melanogaster [29], homologous genes involved in melanin synthesis pathway genes in silkworms were examined in RNA-seq data. Phenylalanine hydroxylase (PAH) was upregulated (2.04-fold, Table 1) and ebony was significantly upregulated (9.36-fold, Table 1), while the expression of yellow, yellow-f2and yellow-f4 decreased (Table 1 and Figure 4A). In addition, yellow relate gene were upregulated in q mutant: yellow-e 2.61-fold, yellow-h2 2.73-fold (Table 1). Tyrosine can be obtained from the diet and L-phenylalanine (phenylalanine) generated by PAH for melanin synthesis in insects [30]. With new cuticle synthesized, sclerotization and pigmentation occur during molting, so conversion of phenylalanine to tyrosine by PAH is critical for melanin formation.

Figure 4

Differentially expressed genes in pigment pathways.

(A) Simplified and synthetic melanin; (B) pteridine; (C) ommochrome pathways. Red, upregulated genes; green, downregulated genes in quail mutant. PAH, phenylalanine hydroxylase; TH, tyrosine hydroxylase; DAT, dopamine acetyltransferase; DOPA, dihydroxyphenylalanine; NADA, N-acetyl dopamine; NBAD, N-β-alanyl dopamine; PO, phenoloxidases; GTP, guanosine triphosphate; GTPCH I, GTP cyclohydrolase I; NBH2, dihydroneopterin triphosphate; BH4, tetrahydrobiopterin; TRPO, tryptophan pyrrolase; KFM, kynurenine formamidase; KMO, kynurenine 3-monooxygenase; PHS, phenoxazinone synthetase.

Our results suggested that upregulation of PAH provided more tyrosine for melanin formation in q mutant. It was reported that the yellow is required to produce black melanin in Drosophila and Bombyx mori [18], [31]. Previous studies found that yellow-f and yellow-f2 catalyze the conversion of dopachrome into 5,6-dihydroxyindole during dopa melanization [32] and yellow-f4 is similar to Drosophila yellow-f [33]. The downregulation of these genes might depress dopa melanization in q mutant. During molting, dopamine is involved in the synthesis of light-colored pigments and sclerotization. The ebony gene encodes N-β-alanyl-dopamine (NBAD) synthase, which converts dopamine to NBAD. NBAD is a precursor for the yellow and reddish-brown cuticular pigments of Drosophila and P. xuthus larvae [31], [34] and is oxidized to NBAD-ortho-quinone which crosslinks cuticle proteins for sclerotization [35]. Drosophila ebony mutants have more darkly pigmented thorax and wings [31]; in silkworm ebony mutants, the body color is smoky in the larval and adult stage [18]. It was reported that yellow-e and yellow-h2 may inhibit melanin pigmentation [36]; in silkworm, yellow-e disruption promoted melanin pigmentation in the larval head and tail [37]. Our results suggested that dopa melanization was prevented and dopamine was converted to NBAD or NADA in q mutant, resulting in light black crescent markings, light star spot markings and a brown cuticle. In addition, we found that other yellow-related genes (yellow-d, yellow-f4-2 and yellow-x) were upregulated in q mutant. The yellow gene family has been identified in many insects, but the function of the gene products is still largely unknown. Three yellow-related genes were upregulated in q mutant, so these genes might have other functions in this mutant.

Differentially expressed genes in pteridine biosynthesis

The genes for enzymes in pteridine biosynthesis were investigated. We found that GTP-CH I was notably upregulated in q mutant: GTP-CH I a by 17.40 fold and GTP-CH I b by 12.94 fold (Figure 4B). GTP-CH I is the first enzyme in the pteridine biosynthesis, but both isoforms were upregulated in q mutant, suggesting increased conversion of GTP to dihydroneopterin triphosphate compared to wild-type resulting in more BH4. GTP-CH I levels are higher in q mutant compared to a background strain during the fourth molt [22], consistent with our results. The increase of GTP-CH I in q mutant might increase BH4 content. Also this was conformed in previous study, which said that the level of BH4 was higher in q mutant than wild type strain [22]. As BH4 is essential as a TH and PAH cofactor during the cuticle screlotization and melanization, the increase of BH4 content in q mutant might increase PAH and TH catalysis of phenylalanine and tyrosine in the melanin pathway(Figure 4). In P. xuthus, GTP-CH I a form is expressed in black regions and is involved in cuticular pigmentation during the fourth molt [2]. Therefore, GTP-CH I a might reinforce melanin synthesis by TH in q mutant. Meanwhile, since Dopa melanization might be prevented, it resulted that more dopamine were convert to NBAD or NADA, which increased cuticle sclerotization and brown color.

Differentially expressed genes in ommochrome biosynthesis

The ommochrome pathway is an important route for eliminating excessive tryptophan metabolites [22] and ommochromes are widely distributed in integument, eyes, testes, wing and extreta in insects [38]. Many mutants with mutations that affect pathway are identified in Drosophila and B. mori, but the molecular mechanisms of the pathway after 3-hydroxykynurenine is largely unknown. The q mutant accumulate high concentrations of xanthommatin in the epidermis [25], so q mutant might be affected in later ommochrome pathway steps. Xanthommatin is a yellow-brown pigment derived from condensation of two molecules of 3-hydroxy-kynurenine via phenoxazinone synthetase (PHS) [39]. However, no reports have described the sequence of phenoxazinone synthetase in insects. To identify genes homologous to PHS in B. mori, we used the Streptomyces antibioticus PHS amino acid sequence (U04283) to query genes predicted in the B. mori genome using BLASTP. We found 6 genes homologous to PHS and belonging to the Cu-oxidase family in SilkDB. Among these genes, one was upregulated (2.70-fold, Table 1) in q mutant, which might result in increased conversion of 3-hydroxykynurenine to xanthommatin (Figure 4C). In addition, kynurenine formamidase, which converts formylkynurenine into kynurenine, was also upregulated in q mutant. Two ommochrome-binding proteins with affinity for phenoxazinone pigments and localized in pigment granules [24] were upregulated in q mutant (Table 1). These results suggested that kynurenine formamidase provided more kynurenine, an ommochrome precursor, and upregulated PHS catalyzed 3-hydroxykynurenine to increase xanthommatin resulting in accumulation in q mutant.

Differentially expressed genes in uric acid transport

In addition to melanin, pteridine and ommochrome, uric acid is excreted as a nitrogen metabolite in insects and is indispensable role for white larval epidermis in B. mori. Uric acid is synthesized in fat bodies, the midgut and Malpighian tubules [40], transported to the epidermis and accumulates as urate granules, making the silkworm larval skin white and opaque [41]. We found three genes (BGIBMGA002922, BGIBMGA003864 and BGIBMGA002581) involved in transport of uric acid were upregulated (2.86-fold, 2.24-fold and 4.12-fold, Table 1) in q mutant. Bm-w-3 (BGIBMGA002922) is a w-3 mutant silkworm gene. Mutants have translucent larval skin from a deficiency in uric acid transport [42]. BGIBMGA003864 is an amino acid transporter of solute carrier proteins involved in uric acid transport in insects and is responsible for the sex-linked translucent (os) phenotype [43]. Mutation of an ABC transporter gene (BGIBMGA002581) is responsible for the failure to incorporate uric acid in the epidermis of ok mutants of B. mori [44]. In addition, Bm-w-3 is involved in pteridine formation and might transport pteridines in the silkworm [44]. Bm-w-2 and Bm-w-3 form a heterodimer responsible for transport of ommochrome precursors. These data suggested that the three genes are involved in transport of uric acid and Bm-w-3 has important role for transporting pteridine and ommochrome. The expression of these genes was significantly increased in q mutant, suggesting that increased uric acid, pteridine and ommochrome were transported to the epidermis and in combination with melanin, resulted in the characteristic body coloration of q mutant.

Differentially expressed genes in juvenile hormone metabolism

Regulation of body coloration by JH is reported in several insects. In swallowtail butterfly, P. xuthus, JH regulates larval pattern switches: young caterpillars (from the first to the fourth instars) are mimics of bird droppings, whereas fifth instar larvae have a pattern that camouflages the larvae among host plant leaves [45]. In the tobacco hornworm Manduca sexta, the normally black mutant larvae do not turn black when treated with JH during molting [46] because JH prevents ommochrome synthesis [47]. Conversely, JH induces xanthommatin in the silkworm larval epidermis [48]. JH acid methyltransferase (JHAMT) converts JH acids or inactive JH precursors to active JHs in the final step of JH biosynthesis in insects. JHMAT levels correspond to JH biosynthesis rates and JHAMT is detected only in the corpora allata of the silkworm fourth instar larvae [49].

We found that JHAMT was highly expressed in integument of q mutant but not in the Dazao strain, which was confirmed by q-PCR (Table S5). Key enzymes (JH esterase, JH epoxide hydrolase and JH diol kinase) (Table 1 and Figure S2) [50], in the JH degradation pathway were also upregulated in q mutant. Xanthommatin content is markedly increased after JH injection in silkworms [48]. Therefore, More JH was likely synthesized and degraded in q mutant during molting, and might be involved in pigment formation in q mutant. A juvenile hormone-inducible protein (BGIBMGA013971) was also upregulated (2.48-fold, Table 1) in q mutant.

Differentially expressed genes in cuticular protein and chitin metabolism pathway

Degradation of old cuticle and regeneration of new cuticle occur during insect molting. The insect cuticle is a thin outer epicuticle with a thicker procuticle of proteins and chitin [51]. In the silkworm, 220 putative cuticular protein genes, expressed mainly in the epidermis and wing disc, have been identified using genome sequences [52]. Of these genes, 62 were differentially expressed in q mutant: 15 upregulated and 47 downregulated (Table S6). The expression profiles of related enzymes in chitin metabolism were examined. The results showed that β-N-acetylglucosaminidase, glucose-6-phosphate isomerase, glutamine: fructose-6-phosphate aminotransferase and UDP-N-acetylglucosamine pyrophosphorylase were significantly upregulated in the q strain (Table 1). Chitin, a linear polymer of β-(1,4)-N-acetyl-D-glucosamine (GlcNAc), is a structural component of the insect trachea, cuticle, cuticular lining of the foregut, hindgut, and peritrophic membrane that lines the lumen of the midgut [53]. Chitin is digested in the cuticle to GlcNAc in a binary enzyme system of a chitinase and a β-N-acetylglucosaminidase in molting fluid [54]. The chitinase hydrolyzes chitin into oligosaccharides and the β-N-acetylglucosaminidase degrades the oligomers to monomers. β-N-acetylglucosaminidase was upregulated in the q mutant, suggesting increased chitin digestion. This result was consistent with degradation of old cuticle. Glucosamine: fructose-6-phosphate aminotransferase (GFAT) is a critical enzyme in chitin synthesis. GFAT is a rate-limiting enzyme and provides the GlcN precursor in chitin biosynthesis [55]. The key enzyme was upregulated in q mutant, indicating more chitin in q strains. The mechanical properties of cuticle are regulated by factors such as the relative amounts of chitin and proteins, protein composition, and degree of sclerotization [56]. We found that many cuticular proteins were downregulated and some enzymes in chitin synthesis were upregulated in q mutant, suggesting higher amounts of chitin and a relatively small amount of cuticular protein than in the Dazao strain. The q mutant has more black spots on the dorsal surface of the integument than Dazao, indicating that melanin in the cuticle increased in q mutant, facilitating sclerotization.

Differentially expressed nuclear receptor genes

Nuclear receptors are ligand-regulated transcription factors involved in a variety of biological processes [57]. We surveyed 19 nuclear receptors in the B. mori genome [58] in our RNA-seq data. Among them, EcR, βFTZ-F1, HR38 and HR39 were upregulated in q mutant (Table 1). Key enzymes in chitin biosynthesis are significantly downregulated and chitin contents in the cuticle are significantly decreased in Spodoptera exigua after injection of EcR dsRNA [59], suggesting that EcR is an important regulator in chitin biosynthesis. EcR and key enzymes were upregulated in q mutant, which might result in increased chitin. Mutants of βFTZ-F1 show that βFTZ-F1 is necessary for larval molting and is crucial for cuticle formation [60]. Cuticular proteins are expressed in different regions of the B. mori epidermis [61]. Three cuticle genes are significantly downregulated in HR38 mutant pupae, disrupting cuticular integrity [62]. These findings indicated HR38 has an important function in cuticle formation. DHR39 is a nuclear receptor with high sequence similarity to FTZ-F1 [57], so DHR39 might have a function similar to FTZ-F1. E74A, an ecdysone-inducible transcriptional factor, was upregulated in q mutant. βFTZ-F1 and E74A regulate the promoter of target cuticular protein to determine the time and location of expression [62].

Quality analysis of RNA-seq data. (A) Distribution of sequence lengths of q and (B) Dazao at 16 h after HCS of fourth molt. Vertical axis, number of sequences. Length normalized to total number of sequences. (C) Distribution of average sequence quality scores for q and (D) Dazao at 16 h after HCS of fourth molt. Sequence quality was calculated as the arithmetic mean of its base qualities. Vertical axis, number of sequences at a quality score normalized to total number of sequences. (E) Coverage for the four DNA nucleotides and ambiguous bases of q and (F) Dazao at 16 h after HCS of fourth molt. Vertical axis, number of nucleotides per type normalized to total number of nucleotides at that position.

(TIF)

Click here for additional data file.

Differentially expressed genes in JH degradation pathway. Red box represents upregulated genes in quail mutant.

(TIF)

Click here for additional data file.

Primer pairs for q-PCR.

(XLS)

Click here for additional data file.

The RPKM of 14,623 annotated silkworm genes in q and Dazao at 16 h after HCS of fourth molt.

(XLS)

Click here for additional data file.

Differentially expressed genes by RNA-seq at 16 h after HCS in q mutant.

(XLS)

Click here for additional data file.

Functional categories of upregulated and downregulated genes in q mutant.

(XLS)

Click here for additional data file.

Differentially expressed genes were verified by qRT-PCR.

(XLS)

Click here for additional data file.

Differentially expressed cuticular protein genes.

(XLS)

Click here for additional data file.