Transcriptomic analysis of the prothoracic gland from two lepidopteran insects, domesticated silkmoth Bombyx mori and wild silkmoth Antheraea pernyi

The prothoracic gland (PG) is an important endocrine organ of synthesis and secretion of ecdysteroids that play critical roles in insects. Here, we used a comparative transcriptomic approach to characterize some common features of PGs from two lepidopteran species Bombyx mori and Antheraea pernyi. Functional and pathway annotations revealed an overall similarity in gene profile between the two PG transcriptomes. As expected, almost all steroid hormone biosynthesis genes and the prothoracicitropic hormone receptor gene (Torso) were well represented in the two PGs. Impressively, two ecdysone receptor genes, eleven juvenile hormone related genes, more than 10 chemosensory protein genes, and a set of genes involved in circadian clock were also presented in the two PGs. Quantitative real time -PCR (qRT-PCR) validated the expression of 8 juvenile hormone and 12 clock related genes in B. mori PG, and revealed a different expression pattern during development in whole fifth larval instar. This contribution to insect PG transcriptome data will extend our understanding of the function and regulation of this important organ.


the prothoracic gland (pG) is an important endocrine organ of synthesis and secretion of ecdysteroids that play critical roles in insects.
Here, we used a comparative transcriptomic approach to characterize some common features of pGs from two lepidopteran species Bombyx mori and Antheraea pernyi. Functional and pathway annotations revealed an overall similarity in gene profile between the two PG transcriptomes. As expected, almost all steroid hormone biosynthesis genes and the prothoracicitropic hormone receptor gene (Torso) were well represented in the two pGs. Impressively, two ecdysone receptor genes, eleven juvenile hormone related genes, more than 10 chemosensory protein genes, and a set of genes involved in circadian clock were also presented in the two pGs. Quantitative real time -PCR (qRT-PCR) validated the expression of 8 juvenile hormone and 12 clock related genes in B. mori PG, and revealed a different expression pattern during development in whole fifth larval instar. This contribution to insect PG transcriptome data will extend our understanding of the function and regulation of this important organ.
The prothoracic gland (PG) of insect is one of the most important endocrine organs that synthesizes and releases ecdysteroid hormone playing critical roles in regulating growth, moulting and metamorphosis. The insect PG is a CPU-like "decision-making center" that integrates a wide range of systemic cues before permitting the production of an ecdysone pulse 1 . In higher Diptera Drosophila the PG, together with the corpus cardiacum (CC) and corpus allatum (CA), is fused into a complex endocrine structure, known as ring gland. However, in other insects, including the domesticated silkworm Bombyx mori (Lepidoptera: Bombycidae) and wild silkworm Antheraea pernyi (Lepidoptera: Saturniidae), these three endocrine glands form separate structures. In B. mori and A. pernyi, the PGs are a pair of semi-transparent or transparent saccate cell clusters with conjunct theca, respectively, located in the tracheal clusters of the prothorax.
As an important endocrine organ, the insect PG has been considered as a model for steroid hormone biosynthesis and regulation 2 . A recent study has suggested that the local clock is a key driver of steroid hormone production in Drosophila PG 3 . An ultrastructural study in Drosophila has suggested that the PG cells may be performing other roles beyond endocrine synthesis 4 . Previous studies on genes expressed in the insect PG were initially focused on characterization of individual genes, particularly those involved in steroid hormone biosynthesis and regulation as well as circadian clock mechanism 5 . Although the ecdysteroid in numerous insects had been studied for decades, yet advances in understanding this important organ at the molecular level remains largely unknown. At the start of this work, the basic genomic information is lacking, although a proteomic approach had been utilized to investigate the feature of B. mori PG 6 . Very recently, two research teams just released their results for

Results
Illumina sequencing and transcriptome assembly. Transcriptomic sequence data were generated using two PG cDNA libraries from B. mori and A. pernyi, and Illumina HiSeq 2500 technology. The PGs were collected from a pool of ~30 silkwom larvae of fifth instar. We acquired 28,159,208 and 24,408,498 clean reads from the PG transcriptomes of A. pernyi and B. mori, respectively, after removing adapters, ambiguous nucleotides and low quality sequences. For B. mori, 5.73 Gbp of clean sequence data was generated with a Q30 value of 91.35% and a GC content of 45.75%. The assembly resulted in 49,287 transcripts longer than 200 bp, which were further assembled into 32,302 unigenes, with an N50 of 1,510 and mean length of 798 bp. For A. pernyi, we ultimately obtained 6.60 Gbp of clean sequence data with a Q30 value of 91.73% and a GC content of 44.55%. The trinity assembly of the clean sequence data of A. pernyi resulted in 64,301 transcripts longer than 200 bp, which were further assembled into 44,067 unigenes, with an N50 of 865 and mean length of 549 bp. For each species, at least 6100 unigenes are >1000 bp and 12,500 >500 bp in length (Additional file: Fig. S1). An overview of the sequencing and assembly process is presented in Additional file: Table S1. The sequence data for B. mori and A. pernyi PGs have been deposited in the NCBI Sequence Read Archive (SRA) database under accessions SRX2434881 and SRX2434884, respectively, and the assembled sequences have been deposited in Transcriptome Shotgun Assembly (TSA) database under accessions GFCX00000000 and GFCY00000000 associated with Bioproject PRJNA357974 and PRJNA357975 for A. pernyi and B. mori PGs, respectively.
Functional annotation revealed an overall similarity in gene profile between the two PGs. For functional annotation, we searched all unigene sequences using Blastx tool against NCBI non-redundant protein database (Nr), with a cut-off E-value of 10 −5 . Using this approach, 15,187 (47.02% of all distinct sequences) and 19,035 (43.19%) unigenes for B. mori and A. pernyi returned a Blast hit in the Nr database, respectively, 8,791 and 10,283 unigenes had specific matches in the Swiss-Prot database, and 9,188 and 10,974 unigenes had matches in the Pfam database. Totally, 23,157 (71.69%) and 22,402 (50.84%) unigenes were annotated in at least one database for B. mori and A. pernyi PGs (Additional file: Table S1).
Firstly, we used Blast2GO 11 to perform functional annotation for the PG transcriptome via gene ontology ( Fig. 1). For B. mori, a total of 8,557 unigenes were assigned GO terms, including 6,208 with hits at the Biological Process level, 4,065 at the Cellular Component level and 7,337 at the Molecular Function level. For A. pernyi, 10,621 unigenes were assigned GO terms, including 7,382 at the Biological Process level, 4,662 at the Cellular Component level and 9,170 at the Molecular Function level. Within the Biological Process GO categories, the most abundant transcripts for the two PGs were assigned to "metabolic process" (6,218 in A. pernyi and 4,926 in B. mori), "cellular process" (5,406 in A. pernyi and 4,625 in B. mori), and "single-organism process" (4,006 in A. pernyi and 3,342 in B. mori). "Cell part" (3,337 in A. pernyi and 2,905 in B. mori), "cell" (3,323 in A. pernyi and 2,894 in B. mori), and "organelle" (2,439 in A. pernyi and 2,103 in B. mori) were the most represented GO terms for Cellular Components in both PGs. For Molecular Function, "binding" (5,751 in A. pernyi and 4,587 in B. mori), "catalytic activity" (5,385 in A. pernyi and 4,111 in B. mori) were the most prevalent in the two PGs. Overall, the percentage of Blastx hits distributed among GO categories was highly similar in both PGs.
Secondly, unigenes of the two PGs were characterized by KOG classification to enable conceptualization of its transcripts into potential functional groups. In total, 10,404 and 12,837 unigenes for B. mori and A. pernyi were annotated to 25 KOG categories, respectively (Additional file: Fig. S2). The numbers of each KOG category were similar between A. pernyi and B. mori PG transcriptome. The KOG classification indicated that except 'general function prediction' , genes involved in "signal transduction mechanisms" (1,373 unigenes in B. mori and 1,454 in A. pernyi), "post translational modification, protein turnover, chaperones"(790 unigenes in B. mori and 1,052 in A. pernyi), and "translation, ribosome structure and biogenesis" (638 unigenes in B. mori and 618 in A. pernyi) were most abundant.
Lastly, KEGG orthology (KO) assignments 12 were comparable between the two PGs. The KO assignment analysis resulted in annotation of 179 and 195 pathways, corresponding to 4,665 and 5,820 unigenes in B. mori and A. pernyi, respectively, and the global KO assignments showed similar trends in both PGs ( Fig. 2A). The second-and third-tier pathways (Fig. 2B,C) also indicated a common expressed-gene profile between the two PGs. In the KEGG second-tier pathway hierarchy, "translation", "folding, sorting and degration" and "transport and catabolism" pathways ranked first to third in the two PGs, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ expression of ecdysteroidgenesis genes in two pGs. Insect PG is one of the most important endocrine organs that synthesizes and releases ecdysteroid hormone. As expected, almost all known steroid hormone biosynthesis genes (neverland, spook, phantom, disembodied, shadow, shroud, Cyp6u1) and the prothoracicitropic hormone receptor gene (torso) are well represented in two PGs (Table 1). Local blast search against the transcriptome data of two silkmoths indicated that no expression of the shade gene is detected, which is consistent with ecdysone being activated to 20-Hydroxyecdysone (20-E) in peripheral tissues and not the PG 13 . Note that B. mori neverland was obtained by searching Blastn against a recent released PG transcriptome data (SRX1142589) 8 . The shroud gene in B. mori that encodes a short-chain dehyrogenase/reductase involving in the ecdysteroid biosynthesis pathway is PG-and ovary-specific 14 ; we got the sequences of the shroud gene in both PG transcriptome data. A recent work in D. melanogaster has provided strong evidence that Cyp6u1 may have a role in ecdysteroidogenesis, possibly in the Black Box 15 ; the homologues were also identified in two silkworm PGs. We also identified two genes encoding ecdysone receptor B and ultraspiracle 2 that are expressed in the PG transcriptomes of two silkworms. These genes were confirmed by comparing them with known genes from D. melanogaster using phylogenetic analysis (Additional file: Fig. S3). The RPKM values of these genes were also evaluated, and spook was the most abundant. The high expression of spook in two silkworm PGs was not consistent with the low expression in D. melanogaster 15 . Genes related to juvenile hormone regulation in pGs. Like ecdysteroids, juvenile hormone (JH) is also an important endocrine hormone that determines the nature of molt, and the CA has been considered as the only source of JH in insects 16 . In this study, we identified 11 juvenile hormone related genes that are represented in the PG transcriptome data from two silkworms ( Table 1 and Fig. 3A), including Farnesyl diphosphate phosphatase (FPPP), Aldehyde dehydrogenase (ALDH), juvenile hormone acid methyltransferase (JHAMT), juvenile hormone epoxide hydrolase (JHEH), juvenile hormone esterase (JHE), cytosolic juvenile hormone-binding protein (cJHBP), juvenile hormone binding protein (JHBP), juvenile hormone esterase binding protein (JHEBP), hexamerin, broad and allatostatin receptor. However, we did not get the sequences of the genes such as NADP+ -dependent farnesol dehydrogenase (FOHSDR), methyl farnesoate epoxidase/farnesoate epoxidase (CYP15A1), juvenile hormone diol kinase (JHDK) and sesquiterpenoid omega-hydroxylase (CYP4C7). All these juvenile hormone related genes were confirmed by phylogenetic analysis (Additional file: Fig. S4). www.nature.com/scientificreports www.nature.com/scientificreports/ Genes involved in circadian clock in pGs. Physiological experiments by transplantation have evidenced the presence of a local clock in PG of the saturniid moth Samia cynthia ricini 17 . A local clock has also been suggested to play a key driver of steroid hormone production in Drosophila PG, and all of genes related to circadian clock were well represented in the ring gland transcriptome 15 . In the present study, a set of genes important for circadian clock mechanism were identified in the PGs of two silkmoths (Table 1 and Fig. 3B), including cryptochrome 1, cryptochrome 2, cycle, clock, vrille, timeless, slimb, period, double time, shaggy, PAR-domain protein 1ε, casein kinase 2 alpha, casein kinase 2 beta, which further confirmed the presence of a local clock in the PG, although the expression level of each gene was low. Among these genes, the expression of clock was the least and casein kinase 2 beta was the most abundant. We confirmed these clock related genes by phylogenetic analysis (Additional file: Fig. S5).
Chemosensory protein genes in pGs. The genome of B. mori harbours 16 CSP genes (BmCSP1-16) 18 ; we identified 11 of them and one novel BmCSP17 in the PG transcriptome data ( Table 2; Fig. 4A). BmCSP17 exhibited a highest sequence identity with BmCSP16, with a value of 61%. In A. pernyi PG, we also identified 10 CSPs. These 10 ApCSPs demonstrated 45-77% amino acid sequence identities with corresponding BmCSPs. Sequence comparison between A. pernyi and B. mori showed that there are six pairs of homologous (Fig. 4B).
Changes in expression level of juvenile hormone regulation and circadian clock related genes in PG during development in whole fifth larval instar. We further used quantitative real time -PCR (qRT-PCR) to validate and investigate the changes in expression level of 11 juvenile hormone regulation and 13 circadian clock related genes in B. mori PG during development in whole fifth larval instar. The qRT-PCR results confirmed the expression of 8 genes out of these 11 juvenile hormone regulation related genes (Fig. 5). Among them, 4 genes (ALDH, JHEH, cJHBP and JHEBP) exhibited a similar expression pattern with a trend of rise first Out of 13 clock related genes identified in two PGs, 12 were indeed expressed in the larval PG of B. mori by qRT-PCR detection method (Fig. 6). Among them, 7 genes (cryptochrome 1, cryptochrome 2, timeless, double time, slimb, casein kinase 2 alpha and casein kinase 2 beta,) presented a similar expression pattern with a trend of rise first then fall; 4 genes (cycle, clock, shaggy and PAR-domain protein 1ε) showed a gradual rise tendency with the highest expression at day 10; period gave a distinct expression change with a trend of fall first then rise.

Discussion
Despite PG is an important endocrine organ in insects, little information on it at the global level is already known, especially the features beyond endocrine synthesis. The present study, for the first time, use a comparative view to analyze the gene network in the PG of two silkmoths, B. mori and A. pernyi, by taking a transcriptomic approach. Since the complete genomes can be available, their annotation rates reache 71.7% for B. mori PG and 57.8% for D. melanogaster ring glands transcriptomes 15 , respectively. The annotation rate of A. pernyi (50.8%) is lower than that of B. mori, but it is comparable to those of other lepidopteran insects such as Helicoverpa armigera (50.8%), H. assulta (54.0%), Spodoptera frugiperda (51.1%) and Athetis lepigone (41.5%) [19][20][21]  www.nature.com/scientificreports www.nature.com/scientificreports/ and B. mori also indicates that the PG transcriptomes share a common expressed-gene profile. This validates the application of RNA-Seq technology for identification of orthologs in non-model organism. Thus, the PG transcriptomes are useful resources to identify gene networks controlling PG ecdysteroidogenesis and understand the novel features of PG beyond endocrine synthesis.
20-E exerts its effect through binding to its receptor, a heterodimer of ecdysone receptor B (EcR) and ultraspiracle (USP). Two genes encoding EcRB and USP2 are expressed in the PG transcriptomes of two silkworms, providing further evidence that 20-E is involved in feedback loops in the PG 22,23 .
Recent studies have shown that JH plays an important role in regulating PG activity 24,25 ; however, the underlying molecular mechanism is severely limited. Decophering the mechanism underlying cross-talk between JH and ecdysone is key to understanding the control of insect growth and development. JH is degraded predominantly by hydrolytic enzymes, JHE and JHEH 25 . JHEBP might function in JHE transportation and degradation when the JH III titer is high 26 . JHBP serves as a carrier supplying JH to the target tissues. Hexamerin modulates JH availability 27 . Broad is a transcription factor that mediates the effects of ecdysone and juvenile hormone. Our qRT-PCR results confirmed the expression of 8 genes (ALDH, JHEH, cJHBP, JHEBP, JHE, Broad, FPPP and JHBP) in B. mori PG. No detectable expression of 3 genes (JHAMT, hexamerin and allatostatin receptor) by qRT-PCR method was consistent with the RNA-seq results by a very low FPKM value. It has already been shown that in B. mori JHAMT was specificlly expressed in the CA, and trace amounts in the PG 28 . The presence of 8 JH related genes in the PG suggested that through these JH related genes JH cross talk with ecdysteroid hormone. To understand the function of these JH related genes in the PG, it is necessary to investigate whether they are functionally important by performing qRT-PCR with CA and other tissues (for example, fat body) in the future work.
Circadian clock is an important regulator of behavior and physiology in insects. In addition to the central clock in brain, peripheral clocks reside in various organs and tissues 29 . For examples, the peripheral clocks in the Malpighian tubules, antenna, fat body of Drosophila have been well studied. These peripheral clocks are independent of the central clock, and the oscillatory machinery and entrainment mechanism of peripheral clocks vary between different tissues and organs. In Drosophila, the eclosion rhythm is set by a local clock residing in the PG that is a key driver of steroid hormone production 3 . Compared with Drosophila, very little attention has been paid to the peripheral clocks in lepidopteran insects 17 . In the PG transcriptomes of two silkmoths, we also identified 13 genes important for circadian clock mechanism and qRT-PCR results confirmed the expression of 12 genes in the B. mori PG, thus extending our understanding of the local peripheral clock residing in the PG of lepidopteran insects 17 that may also involve in steroid hormone production. To address this issue, we will investigate the fluctuations in expression levels of circadian clock related genes of night and day PG samples in the future.
Insect CSPs are small soluble acidic proteins that are believed to be involved in chemical communication, including perception, identification, transport and transduction of semiochemicals from environment (olfaction, taste and others) and may be associated with regulation of circadian rhythms and maturation of tissue or appendage 30 . CSPs are expressed not only in insect sensory organs, but also in other tissues that lack gustatory and www.nature.com/scientificreports www.nature.com/scientificreports/ olfactory neurons 31 . These non-chemosensory tissues included cuticle, legs, labial palp, pheromone gland, tarsi, proboscis, wings, testes, ovaries, compound eyes, hemolymph and ejaculatory 30 . The PG transcriptome data of two silkworms offered us an opportunity to investigate gene expression profiles of CSP genes on a large-scale in  www.nature.com/scientificreports www.nature.com/scientificreports/ PGs. Our study evidenced the presence of 12 and 10 CSPs in B. mori and A. pernyi PG, suggesting that there have a link between CSPs and ecdysteroids 32 .
Earlier study in B. mori demonstrated that KK-42 can reduce the incidence of embryonic diapause when administered to the mother during her final larval instar 33 . By contrast, our recent results showed that KK-42 can delay termination of the pupal diapauses in A. pernyi and Helicoverpa zea, and boost pupal diapause incidence when administered to larvae of H. zea 34 . The mechanism is that KK-42 appears to act by inhibiting ecdysteroid biosynthesis within the PG, without killing the PG cells 35 . Previous studies suggested that a KK-42 binding protein might be a receptor of an endogenous signaling compound 36 ; the expression of this gene could be detectable in many organs of A. pernyi 37  In conclusion, this present work provided a comparative analysis of the PG transcriptomes of two silkmoths, whose associated expressed-gene profile were highly similar. Our results uncovered the presence of at least 8 juvenile hormone related genes, 12 circadian clock genes, and 10 chemosensory protein genes in both PGs. This contribution to insect PG transcriptome data will extend our understanding of the function and regulation of this important organ.

Materials and Methods
Insect materials and samples collection. Samples were derived from strain Shenhuang no. 2 of A. pernyi and Dazao of B. mori maintained at the Department of Sericulture, Shenyang Agricultural University in Shenyang. Larvae of A. pernyi strain were reared on oak trees in the field until the fifth larval stage, and then reared at room (25 °C, natural humidity) using oak branches with leaves. Larvae of B. mori were reared using mulberry leaves during the whole larval stage at 25 °C and 50-70% relative humidity in natural light. To generate samples, ~30 worms of B. mori and A. pernyi were collected in the third day and the fifth day of the fifth instar, www.nature.com/scientificreports www.nature.com/scientificreports/ respectively. The PGs were carefully removed from the worms in insect Ringer physiological saline buffer under dissecting microscope, then immediately placed in 500 μl TRIzol reagent (Beijing Sinogene, China) and stored at −80 °C. Frozen tissues in TRIzol were shipped to Biomarker Technologies in Beijing, China for RNA extraction, library preparation and DNA sequencing.
RNA extraction and sequencing. Total RNAs were extracted from frozen PGs with TRIzol. RNA integrity and concentration were assessed with a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries were prepared using NEBNext Ultra RNA Library Prep Kit for Illimina (NEB, E7530) and NEBNext Multiplex Oligos for Illumina (NEB, E7500). Two RNA samples were sequenced on the Illumina HiSeq ™ 2500 sequencing platform (paired-end, 125 bp reads).
Transcriptome assembly. High-quality clean reads were obtained by removing the adaptor sequences, duplicated sequences, ambiguous reads ('N'), and low-quality reads. For transcriptome data of A. pernyi that has no complete genomic data available, the clean reads were pooled for assembly using Trinity (http://trinityrnaseq. sourceforge.net/) 38 , and the related contigs were then clustered using the TGICL software 39 to yield unigenes (without N) that cannot be extended on either end, and redundancies were removed to acquire non-redundant unigenes. For transcriptome data of B. mori that has complete genomic data available, the clean reads were mapped to genome using Tophat2 software 40 . Transcript expression levels were estimated with FPKM values (fragments per kilobase of exon per million fragments mapped) by the Cufflinks software 41 . N50 and mean lengths of the transcripts associated with each sample were calculated. The values for N50 length and mean length indicated high quality samples, sequences and assemblies for the PGs of two silkmoths.
Annotation. The unigenes of the two silkmoths were compared against public databases, including NR (non-redundant), GO (gene ontology), KOG (eukaryotic ortholog groups), KEGG (Kyoto Encyclopedia of Genes and Genomes), Swiss-Prot and TrEMBL databases using BLASTx with an E-value cutoff at 10 −5 to retrieve protein functional annotations with the highest sequence similarity. High-priority databases (followed by Nr, Swiss-Prot, and KEGG) were selected to determine the direction of the unigene sequences. The best aligning results were used to predict the coding region sequences from unigenes, and the coding sequences were translated into amino sequences using the standard codon table. phylogenetic analysis. The accession numbers of sequences used for phylogenetic analysis are listed in Additional file: Table S2. Amino acid sequences were aligned with ClustalX 1.83 42 and unrooted trees were constructed with MEGA6.0 43 using the neighbour-joining method, with Poisson correction of distances and bootstrap replications set at 1000.
Quantitative real time -pCR (qRt-pCR). Two PGs per larva were used as one sample to extract the total RNA from 1-day-old fifth instar larva to 10-day-old fifth instar larva (matured silkworm) in the present study. All the extracted RNA (at least 2.6 μg/per fifth instar larva) using TRIpure (Beijing Aidlab biotechnologies Co. Ltd.) was converted into cDNA using the oligo(dT) 15 primer with the PrimeScript RT reagent Kit with gDNA Eraser (Takara Biotechnology Dalian Co. Ltd.). The total volume of qRT-PCR reactions was 10 µl, containing 3.6 µl of TB Green Premix Ex Taq (TaKaRa), 0.4 µl of specific primers (10 µM), 1 µl of cDNA and 5 µl of ddH 2 O. qRT-PCR was performed with a BIO-RAD CFX Connect Real-Time System, and the conditions were as follows: 95 °C for 30 s followed by 40 cycles in 95 °C for 5 s and 60 °C for 30 s. Gene-specific primers used for qRT-PCR analysis are listed in Additional file: Table S3. The mRNA expression levels of the genes of interest were calculated with the 2 −ΔΔCt method and nomalized to the abundance of a house-keeping gene, ribosome protein 49 (rp49). The relative mRNA levels of each gene were represented as folds over the expression levels of rp49.