Transcriptome analysis of the epidermis of the purple quail-like (q-lp) mutant of silkworm, Bombyx mori

A new purple quail-like (q-lp) mutant found from the plain silkworm strain 932VR has pigment dots on the epidermis similar to the pigment mutant quail (q). In addition, q-lp mutant larvae are inactive, consume little and grow slowly, with a high death rate and other developmental abnormalities. Pigmentation of the silkworm epidermis consists of melanin, ommochrome and pteridine. Silkworm development is regulated by ecdysone and juvenile hormone. In this study, we performed RNA-Seq on the epidermis of the q-lp mutant in the 4th instar during molting, with 932VR serving as the control. The results showed 515 differentially expressed genes, of which 234 were upregulated and 281 downregulated in q-lp. BLASTGO analysis indicated that the downregulated genes mainly encode protein-binding proteins, membrane components, oxidation/reduction enzymes, and proteolytic enzymes, whereas the upregulated genes largely encode cuticle structural constituents, membrane components, transport related proteins, and protein-binding proteins. Quantitative reverse transcription PCR was used to verify the accuracy of the RNA-Seq data, focusing on key genes for biosynthesis of the three pigments and chitin as well as genes encoding cuticular proteins and several related nuclear receptors, which are thought to play key roles in the q-lp mutant. We drew three conclusions based on the results: 1) melanin, ommochrome and pteridine pigments are all increased in the q-lp mutant; 2) more cuticle proteins are expressed in q-lp than in 932VR, and the number of upregulated cuticular genes is significantly greater than downregulated genes; 3) the downstream pathway regulated by ecdysone is blocked in the q-lp mutant. Our research findings lay the foundation for further research on the developmental changes responsible for the q-lp mutant.

Introduction nuclear receptors were analyzed between 932VR and q-l p . The results lay a foundation for further research of the q-l p formation mechanism.

The silkworm strains
The silkworm strains 932VR and q-l p were provided by The Sericulture Research Institute, Chinese Academy of Agricultural Sciences (Zhenjiang, China). The silkworms were fed fresh mulberry leaves under standard conditions with alternating 12 hours of illumination and 12 hours of darkness at 25±2˚C. Epidermis samples were collected from every three silkworms after 6 hours of molting in the 4 th instar. All samples were frozen and stored at -80˚C.

RNA extraction and transcriptome sequencing
Total RNA was extracted using RNAiso Plus (TaKaRa, China) and dissolved in RNase-free water. The concentration of total RNA was determined using a NANODROP1000 microspectrophotometer (Thermo, USA) after treatment with DNase. Total mRNA was enriched using Oligo (dT) magnetic beads, after which the mRNAs were sheared into short fragments using hybridization-interruption reagents. These short RNA fragments were then synthesized into double-stranded cDNA using six-base random primers, and terminal modification was performed for the purified double-stranded cDNA; a base (A) tail was added, and the fragments were ligated. The quality of the cDNA library constructed was tested using an Agilent 2100 Bioanalyzer; after quality criteria were met, sequencing was executed using an Illumina HiSeqTM2500 system (Illumina, USA).

RNA-Seq data analysis
Quality analysis was conducted for the original data obtained using the Illumina HiSeqTM2500 system; after quality criteria were met, clean reads were screened by filtering out low-quality reads. The sequences were aligned to the silkworm genome database SilkDB (http://silkworm. swu.edu.cn/silkdb/); after a second quality analysis for alignment, analysis of the distribution and coverage of the clean reads on the reference sequence was conducted. RPKM (Reads Per Kb per Million reads) [24] was used to calculate the expression level of genes, with RPKM = mapped reads of gene/(the total mapped reads of all genes Ã the length of this gene) Ã 10^9. The RPKM of a gene ranged up to 5, and difference in expression was considered at P < 0.05; the fold change of q-l p and 932VR RPKMs was greater than 2. The function of the differentially expressed genes and pathways related to pigment were analyzed using BLASTGO and KEGG, respectively.

Quantitative reverse transcription PCR
Epidermis samples were collected from every three silkworms after 6 hours of molting in the 4 th instar, and total RNA was prepared using RNAiso Plus and reverse transcribed using the Reverse Transcriptase M-MLV (RNase H-) kit (TaKaRa, China) after treatment with DNase. The cDNA was diluted to 20 ng/μL and used as the template for qRT-PCR. The 20-μL reaction included 1 μL primer (10 μmol/L, S1 Table), 1 μL cDNA, 10 μL 2×SYBR 1 Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, China) and 8 μL ddH 2 O. After a rapid centrifugation step, quantitative reverse transcription PCR (qRT-PCR) was performed using a LightCycle 96 real time PCR system (Roche, Switzerland) with the reaction program below: a three-step reaction protocol of 45 cycles of 95˚C for 10 s, 58˚C for 10 s and 72˚C for 10 s after a 10-min step of pre-degeneration, followed by melting. Relative expression was calculated using the 2 -ΔΔCt method [25] with glyceraldehyde-3-phosphate dehydrogenase (GAPDH, BGIBMGA007490) as the reference gene. Melting analysis was performed to guarantee the qRT-PCR quality. The relative expression of 932VR and q-l p were compared, and RNA-Seq and qRT-PCR data were compared to analyze several key pathways.

General information of RNA-Seq
Prior to bioinformatic analysis, we assessed the quality of the transcriptome sequencing data, including quality evaluation of bases, base distribution, pretreatment of the data quality and pollution detection of reads (Fig 2). As expected, all of test results met the requirements. For 932VR, we obtained 19,992,852 paired-end reads, of which 13,792,700 reads were aligned to the silkworm genome database (68.99%) and 11,357,633 to silkworm genes (56.80%). For q-l p , we obtained 16,048,696 paired-end reads, of which 10,987,505 reads were aligned to the silkworm genome database (68.46%) and 8,639,063 to silkworm genes (53.83%). Expression of 14,623 silkworm genes was calculated based on the RPKM method (S2 Table), where an RPKM for one gene less than 5 was regarded as scarcely expressed. At a P value < 0.05, a difference in expression was considered to exist when the ratio of the RPKM for q-l p and 932VR was up to 2 or less than 0.5. The results indicated 515 differentially expressed genes between 932VR and q-l p , of which 234 were upregulated and 281 downregulated (S3 Table).
Based on BLASTGO analysis, we examined the top 10 upregulated and downregulated genes in q-l p compared to 932VR. The upregulated genes mainly encode structural constituents of the cuticle, membrane proteins, transport related proteins, protein binding proteins, enzymes catalyzing oxidation-reduction process, hydrolases, and ATP-binding proteins. The downregulated genes mainly encode protein-binding proteins, structural constituents of the membrane, transport related proteins, enzymes catalyzing oxidation-reduction process, hydrolases, ATP-binding proteins, nucleic acid-binding proteins, and transcription factors. (Fig 3).
Verifying the accuracy of the RNA-Seq data by qRT-PCR First, we analyzed several housekeeping genes. Although no significant difference between the wild-type 932VR and the mutant q-l p was observed for glyceraldehyde-3-phosphate dehydrogenase (BGIBMGA007490, fold = 1.34) or ribosomal protein L3 (BGIBMGA013567, fold = 1.05), a significant difference for BmActin3 (BGIBMGA005576, fold = 0.26) was found. We then determined that the accuracy of the RNA-Seq data reached the required level. Based on manifestations of the mutants and the differentially expressed genes according to RNA-Seq, we selected certain genes for qRT-PCR testing, including differentially expressed genes, key biosynthesis genes for the three types of pigments and chitin and genes coding cuticular proteins and several related nuclear receptors. The results are listed in Table 1 and S1 Table. Differentially expressed genes After analysis of differentially expressed genes revealed by RNA-Seq, we selected 12 strongly differentially expressed genes for qRT-PCR testing; the results are listed in Table 1. Of the 12 genes, 6 were upregulated by RPKM and 6 downregulated. qRT-PCR results showed that 10 genes were consistent with the RPKM values. The remaining two genes were upregulated according to RPKM but downregulated according to qRT-PCR. These results indicated an effective rate of RNA-Seq of 83.33% and that the RNA-Seq data were reliable. Among the 10 genes strongly differentially expressed between 932VR and q-l p , four appear to encode cuticle proteins, one encodes the heavy chain of myosin, one encodes trypsin 1, one encodes nesprin 1, one encodes mucin, one encodes a member of the cytochrome P450 family and one is unknown.
The main phenotype of q-l p is an epidermis that is similar in appearance to the pigment mutant quail; the pigment is one of the important characteristics of the epidermis [26], and there are many cuticle proteins in the epidermis [27]. Indeed, almost half of the 10 strongly differentially expressed genes encode cuticle proteins, and they may play key roles in formation of the epidermis and pigmentation. Members of the cytochrome P450 family comprise monooxygenases belonging to the CYP gene family [28]. In insects, cytochrome P450 is involved in many anabolism processes, including metabolism of exogenous substances such as plant secondary metabolites and pigments, as well as endogenous substances such as juvenile hormone, ecdysone, fatty acid, and pheromones [29]. Downregulation of cytochrome P450 in q-l p may influence the synthesis and signaling of a series of hormones. Trypsin is a protein hydrolase, a large class of differentially expressed genes in q-l p . Mucin, mainly distributed in epithelial cells, has a protective role inside cells and can also contribute to a protective extracellular mucin gel [30]. Downregulation of mucin in q-l p may be related to its high death rate, as q-l p is more easily damaged by adverse environmental factors. Nesprin 1, a component of the karyotheca,

Differentially expressed pigment biosynthesis genes
Melanin is one of the most important pigments, and it has been studied in detail in Drosophila melanogaster [4, [9][10][11][12][13][14]. We analyzed key genes of melanin biosynthesis based on the synthesis pathway of D. melanogaster (Fig 4A). Melanin is produced from tyrosine by a series of enzymatic reactions. Tyrosine is produced from phenylalanine via phenylalanine hydroxylase (PAH) [33] and is then transformed to dopa via tyrosine hydroxylase (TH). The results for PAH (fold = 0.77) and TH (fold = 0.63) expression in q-l p were not obviously different from the results for 932VR. However, yellow (fold = 0.02) and yellow-f2 (fold = 0.21) were downregulated in q-l p ; these two genes play key roles in the processing of dopa to dopa melanin [12,15]. Dopa decarboxylase (Ddc, fold = 19.38) catalyzes the step from dopa to dopamine, and dopamine acetyltransferase (dat, fold = 10.46) catalyzes the step from dopamine to N-acetyl dopamine (NADA); both were upregulated in q-l p , leading to accumulation of dopamine and NADA. Yellow-e catalyzes the reaction from NADA to NADA sclerotin, and downregulation of yellow-e (fold = 0.05) leads to further accumulation of dopamine and NADA. Silencing ebony, which encodes the enzyme catalyzing the reaction from dopamine to N-β-acetyl dopamine (NBAD), also leads to further accumulation of dopamine and NADA. A lack of ebony might be a reason for the downregulation of tan (fold = 0), which catalyzes the step from NBAD to dopamine and the upregulation of blank (fold = 47.4), which catalyzes the step from uracil to β-alanine [34]. Upregulation of Ddc leads to most available dopa being transformed into dopamine, and a decrease in dopa and downregulation of yellow and yellow-f2 lead to a lack of dopa melanin. However, accumulation of dopamine causes production of more dopamine melanin, resulting in brownness, and might be responsible for the mutant q-l p . GTP cyclohydrolase (GTP-CH) is the first enzyme in the biosynthesis of pteridine [20,35]. Two types of GTP-CHs are expressed differentially; no difference in GTP-CH I a (fold = 0.78) was observed between 932VR and q-l p , though GTP-CH I b (fold = 39.31) was markedly upregulated in q-l p , which is consistent with quail [20]. Upregulation of GTP-CH I would cause an increase in BH2, followed by an increase in BH4, and eventually cause an increase in biological pteridine, which might lead to pigment deposition in the epidermis (Fig 4B).  Transcriptome analysis of the purple quail-like mutant (q-l p ) in silkworm, Bombyx mori Ommochrome is a metabolite of tryptophan and can be transformed into different colors varying from red, purple, black to yellow by participating in redox reactions [36]. As an excess of tryptophan is poisonous to an organism, transforming tryptophan to ommochrome avoids tryptophan accumulation [7]. Tryptophan 2,3-dioxygenase plays a role in converting tryptophan to formylkynurenine. Kynurenine is then produced from formylkynurenine via kynurenine formamidase (KFM) and converted to 3-hydroxy kynurenine. Xanthommatin is a ommochrome derived from two molecules of 3-hydroxy kynurenine through the action of phenoxazinone synthetase (PHS) [37,38]. In q-l p , these two key genes, KFM (fold = 1.68) and PHS (fold = 2.85), were upregulated, which might result in more xanthommatin production. Ommochrome-binding protein (OBP) is a key pigmentation-related protein in the epidermis [39]. According to our RNA-Seq results, two OBPs are downregulated in q-l p . However, after confirmation using qRT-PCT, OBP1 (fold = 1.57) was shown to be upregulated in q-l p , with no difference in OBP2 (fold = 0.96) between 932VR and q-l p . These results also correspond to the quail mutant [20]. Upregulation of KFM and PHS leads to accumulation of ommochrome, and upregulation of OBP enhances this process, which might facilitate the observed pigmentation in q-l p (Fig 4C).

Differentially expressed genes encoding cuticle proteins
The epidermis is a very important organ of the silkworm, providing protection for external adverse environmental factors. The silkworm must go through several molting processes in larval development stages, during which, the old epidermis is replaced by a new epidermis. To maintain epidermal markings, pigment is deposited in every molt. The epidermis is composed of a large number of cuticle proteins and chitin. There are more than 200 cuticle protein genes in B. mori [27] (S4 Table), of which 107 were not expressed (RPKM<5 in both 932VR and q-l p ) and 101 were expressed. Of the 101 expressed genes, 78 were upregulated and 23 downregulated in q-l p . In addition, 66 and 13 genes were upregulated and downregulated !2-fold in q-l p ; for fold changes !5 in q-l p , we found 42 upregulated genes and 6 downregulated genes (Fig 5).
The results above reveal more upregulated genes (77.2%) than downregulated genes (22.8%), and for fold changes !5, the ratio of upregulated genes reached 87.5%. In q, the inverse is true: of 62 differentially expressed genes, only 15 genes are upregulated and 47 downregulated, with a ratio of upregulated genes of 24.2% [20]. By comparing samples, we found that the control q strain Dazao exhibits normal markings but that the control q-l p 932VR strain is a plain silkworm. Therefore, some of these upregulated genes might be associated with normal marking formation.
By comparing differentially expressed cuticular genes between q/Dazao and q-l p /932VR, we found 5 upregulated genes and 3 downregulated genes in q and q-l p as fold 0.5 or fold!2 (Table 2). These genes might play key roles in the pigmentation of q and q-l p .
Increasing chitin biosynthesis and decreasing chitin metabolism might lead to accumulation of chitin. Correspondingly, many chitin-binding proteins were increased, consistent with the fact that more than 83.5% of differentially expressed cuticle protein genes were upregulated in q-l p . There were also more cuticle protein genes expressed in q-lp (96 genes expressed at RPKM!5) than in 932VR (73 genes expressed at RPKM!5). These results suggest that q-l p requires more cuticle proteins than 932VR to produce its epidermis.

Differentially expressed nuclear receptor genes
As transcription factors and signal transduction-related genes were abundant among the differentially expressed genes, we examined several nuclear receptors. Nuclear receptors, a gene superfamily, interact with transcription factors to regulate the growth and development of eukaryotic cells [45]. There are two types of nuclear receptors, one of which has a known endogenous ligand; the other is termed an orphan receptor, as the ligand has yet to be identified [46]. In this study, we identified several common nuclear receptors, including FTZ-F1, E74a, HR38, HR39 and EcR [47], and the signal transduction pathway controlled by 20hydroxyecdysone (20E) was exhibited in Fig 6. All of these receptors were upregulated, except for the ecdysone receptor (EcR, fold = 0.52).
HR38 belongs to the NR4A gene family, which can block the binding of heterologous dimers of EcR (ecdysone receptor) and USP (Ultra spiracle) [48,49]. HR38 is also associated with epidermis formation [50,51]. Upregulation of HR38 (fold = 2.51) might block the structure of EcR-USP heterologous dimers and affect the binding of EcR-USP and ecdysone [52], thus altering the downstream pathway regulated by ecdysone with correspondingly downregulation of the ecdysone receptor. Based on the above results, we hypothesized that upregulation of HR38 causes a reduction in EcR-USP and then blocks ecdysone signaling pathways, which might be responsible for developmental abnormalities phenotype, such as slow growth, low food consumption, and high mortality. Moreover, upregulation of HR38 also strengthens the synthesis of cuticle proteins, corresponding to the upregulated cuticular genes compared to downregulated cuticular genes in q-l p and consistent with the finding that key enzymes of the chitin synthesis pathway were upregulated.
HR39 is necessary in female D. melanogaster flies due to its roles in copulation and fertilization. Additionally, HR39 also controls specific expression of cytochrome P450 [53,54].
FTZ-F1 participates in larvae during molting and is associated with epidermization [55]. FTZ-F1 is highly homologous to HR39 [56], with the same binding sites in alcohol dehydrogenase (Adh), fushi tarazu (ftz), and new glue (ng) [57], indicating that FTZ-F1 and HR39 have other similar functions. E74a contains an E26 transformation-specific (ETS) domain and is induced by steroids to regulate transcription of genes such as L71-6 and L71-6, which are related to epidermis formation in pupae [58]. This steroid receptor was also upregulated in ql p . Upregulation of these nuclear receptors is associated with epidermis formation, indicating that more cuticle proteins are synthesized in q-l p than in 932VR. This was also consistent with the genes differentially expressed between 932VR and q-l p and with the quantity and content of cuticle proteins in these strains.

Conclusion
In this study, we analyzed the transcriptome of the epidermis silkworm (B. mori) purple quaillike mutant (q-l p ) 4 th -instar larvae during molting using RNA-Seq and validated the reliability by qRT-PCR. Our results allowed us to reach several conclusions. First, genes involved in biosynthesis of the pigments melanin, ommochrome and pteridine were all increased in the q-l p mutant, which might be responsible for the pigmentation of the epidermis of this strain. In the melanin synthesis pathway, a large amount of dopamine melanin would be produced but less dopa melanin, NBAD sclerotin and NADA sclerotin. Second, a greater abundance of cuticle proteins were expressed in q-l p than in 932VR, with significantly more upregulated than downregulated cuticular genes, indicating that the q-l p mutant requires more cuticle proteins to produce its epidermis than 932VR. Third, the downstream pathway regulated by ecdysone appeared to be blocked by HR38 upregulation, and correspondingly we found that the ecdysone receptor was downregulated, which might be responsible for the developmental abnormalities phenotype of q-l p .
Supporting information S1