Transcriptome sequencing and assembly results
Illumina transcriptomes of total RNA were sequenced and assembled in five different developmental stages of M. separata. The assembly results totaled 140562 contigs, with an average length of 1201.07bp, N50 of 2159bp, and GC of 41.94%. A detailed description of the assembly results is shown in Table 1.
Through Blast alignment at NCBI (e value < e-3), more than half of annotated genes (74,059)were obtained from the target alignment sequences. Among them, 6120, 351851, 47632, 48973, 27817, and 31569 unigenes were effectively annotated in NR, NT, Swiss-Prot, KEGG, COG, and GO, respectively (Figure S1). Finally, a total of 19733 unigenes were collectively annotated in five different databases (Figure S2).
Species homology analysis results show that the most paired species are Amyelois transitella 12958 (21.17%), B.mori 10635 (17.38%), Papilio xuthus 5457 (8.92%), Papilio machaon 4205 (6.87%), Danaus plexippus 3085 (5.04%), Operophtera brumata 2920 (4.77%), Papilio polytes 2348 (3.84%), Plutella xylostella 2335 (3.82%), H. armigera 1148 (1.88%). Other insect species showed the rate of unigenes hit less than 1% (Figure S3).
Global analysis of differentially expressed unigenes
There were 22884 (10990 up-regulated and 11894 down-regulated), 23534 (11534 up-regulated and 12000 down-regulated), 26643 (12718 up-regulated and 13925 down-regulated), and 33238 (18001 up-regulated and 15237 down-regulated) DEGs identified between ML and W, W andP1, P1 and P5, and P5 and P10, respectively (Figure 2A). The number of DEGs between other groups is shown in Figure S4. A total of 3135 (0.76%) commonly DEGs were found to be common to five developmental stages (Figure 2B).
GO distribution analysis of DEGs potentially involved in the prepupal-pupal transition in M. separata
GO distribution analysis results showed that these DEGs were categorized into 59 sub-categories belonging to three main GO categories: biological processes, cellular components, and molecular functions. These categories consisted of 17, 22, and 21 sub-categories respectively (Figure 3).
Under biological process, the oxidation-reduction process was the most significantly prepupal-pupal transition-related GO term in ML vs W and P5 vs P10, with 700 (6.06%, 193 up-regulated and 507 down-regulated genes) and 1023 (4.77%, 585 up-regulated and 438 down-regulated genes) DEGs, respectively, while the most significantly prepupal-pupal transition-related GO term was found in the proteolysis in W vs P1, P1 vs P5, with 629 (5.95%, 233 up-regulated and 396 down-regulated genes) and 622 (5.83%, 215 up-regulated and 407 down-regulated genes) DEGs, respectively (Figure 3).
GO enrichment analyses of DEGs showed that the membrane and integral component of membrane in the cellular component category were significantly enriched. A total of 2388 (8.11%, 1082 up-regulated and 1306 down-regulated genes), 2245 (9.49%, 901 up-regulated and 1344 down-regulated genes), 2404 (9.41%, 1115 up-regulated and 1289 down-regulated genes) and 3615 (8.03%, 2024 up-regulated and 1591 down-regulated genes) DEGs were enriched in the membrane sub-category in ML vs W, W vs P1, P1 vs P5, and P5 vs P10, respectively. Integral component of membrane was enriched with 2346 (7.97%, 1064 up-regulated and 1282 down-regulated genes), 2208 (9.33%, 892 up-regulated and 1316 down-regulated genes), 2372 (9.28%, 1093 up-regulated and 1279 down-regulated genes) and 3514 (7.80%, 1975 up-regulated and 1539 down-regulated genes) DEGs in ML vs W, W vs P1, P1 vs P5, and P5 vs P10, respectively (Figure 3). These DEGs were mainly focused on solute carrier family, MFS transporter, low-density lipoprotein receptor, glucuronosyltransferase, ATP-binding cassette (ABC), P450, alcohol-forming fatty acyl-CoA reductase, ATPase, NADH dehydrogenase, stearoyl-CoA desaturase, chitin, fatty acids, E3 ubiquitin-protein(Table S1). Also, the CPs, 20E and JH related genes that are closely related to the metamorphic development of Lepidoptera insects (Willis et al., 2010; Liu et al., 2015), were also identified in the membrane and integral component of membrane sub-categories (Table S2), suggesting these genes are likely to play a role in the prepupal-pupal transition of M. separata.
The hydrolase activity was the most significantly affected GO term under molecular function in ML vs W, W vs P1, P1 vs P5, with 784 (6.79%, 286 up-regulated and 498 down-regulated genes) and 727 (9.49%, 201 up-regulated and 526 down-regulated genes), 689 (6.46%, 271 up-regulated and 418 down-regulated genes) DEGs, respectively. The metal ion binding was significantly enriched in P5-vs-P10, with 1262 (5.88%, 596 up-regulated, and 666 down-regulated genes) (Figure 3).
KEGG pathway analysis of DEGs potentially involved in the prepupal-pupal transition in M. separata
KEGG pathway analysis indicated that the DEGs were enriched in the metabolic pathways during the prepupal-pupal transition of M. separata. These DEGs were mainly focused on starch and sucrose metabolism, pyrimidine metabolism, protein digestion and absorption (Figure 4). In the metabolic pathways, there were 2218 (16.85%, 745 up-regulated and 1473 down-regulated genes), 2053 (16.26%, 890 up-regulated and 1163 down-regulated genes), 2046 (14.64%, 914 up-regulated and 1132 down-regulated genes), and 3037 (16.14%, 1890 up-regulated and 1147 down-regulated genes) DEGs in ML vs W, W vs P1, P1 vs P5 and P5 vs P10, respectively.
Besides, there were only 35 (0.27%), 26 (0.21%), 28 (0.20%), and 37 (0.2%) DEGs enriched in insect hormone biosynthesis, which was closely associated with JH and 20E, and 78 (0.59%), 73 (0.58%), 79 (0.57%), 143 (0.76%) and 60 (0.46%), 50 (0.40%), 70 (0.50%), and 73 (0.39%) DEGs enriched in CPs belongs to the GnRH signaling pathway and the ErbB signaling pathway in ML vs W, W vs P1, P1 vs P5 and P5-vs-P10, respectively (Figure S5).
We further elaborated on the gene families of particular interest as potential targets for the prepupal-pupal transition of M. separata below.
The number of DEGs detected in the CPs, 20E and JH biosynthesis
A total of 33 genes encoding CPs, and 18 and 7 genes involved in 20E and JH biosynthesis, respectively, were detected during the prepupal-pupal transition of M. separata (Table 2). The number of genes in CPs, 20E, and JH in the five different development stages of M. separata is similar(Table S3). To further confirm whether those genes were well clustered into their relevant phylogenetic branches based on the phylogenetic analysis of the CP (Figure 5), 20E (Figure 6) and JH (Figure 7) from M. separata and eight lepidopteran insects including S. litura, A. dissimilis, T.ni, S. frugiperda, G. mellonella, B.mori, H.armigera, and S.exigua.
Expression levels of CP genes
There were 11 up-regulated and 2 down-regulated CP genes with a maximum up-regulated and down-regulated |log2Ratio| value of 3.3 (CL1279.Contig108_All) and 4.9 (CL7899.Contig3_All) in ML vs W(Table 2). From P1 to P5, one gene (CL9103.Contig5_All) was significantly increased (|log2Ratio|=2), and 9 genes were decreased (the maximum |log2Ratio|=3.8). In the comparison between W-vs-P1 and P5-vs-P10, 5 genes were up-regulated whereas 6 were down-regulated (Table 2).
Five genes (CL12980.Contig3_All, CL1585.Contig16_All, CL2596.Contig51_All, CL14397.Contig3_All, and CL16856.Contig2_All) were lowest expressed in ML(Figure 8A). The expression levels of 8 genes (unigene12087_All, CL16856.Contig2_All, unigene8645_All, CL18975.Contig2_All, CL6991.Contig6_All, CL10379.Contig10_All, CL9076.Contig3_All, and CL1337.Contig49_All) in W and P1 were higher than in other stages, while 3 genes (CL1585.Contig16_All, CL2596.Contig51_All, and CL1279.Contig108_All) were highly expressed in W(Figure 8A). The expression levels of 3 genes (Unigene15668_All, CL10689.Contig4_All, and CL1337.Contig49_All) in P5 were high than in other stages. Moreover, 3 genes, including CL12980.Contig3_All, CL4397.Contig3_All and CL9739.Contig5_All were highly expressed in P10 while the expression level of unigene15668_All was lower in P10 than in other stages (Figure 8A).
Expression levels of genes involved in JH biosynthesis
There were 11 and 6 down-regulated genes in ML vs W and P1 vs P5, respectively, and both of them have two up-regulated genes. Compared to W, 4 genes were significantly increased (the maximum |log2Ratio|=3.8) and 5 genes were decreased (the maximum |log2Ratio|=4.3) in P1. There were 7 up-regulated and 3 down-regulated genes in P5 vs P10; of these genes, CL696.Contig45_All and CL9391.Contig2_All showed the maximum up-regulated and down-regulated |log2Ratio| value of 5.8 and 2.3, respectively (Table 2).
Six genes (CL17942.Contig2_All, CL14984.Contig2_All, CL15275.Contig2_All, CL817.Contig4_All,CL4540.Contig4_All, and CL9391.Contig2_All) in the JH pathway were mainly expressed in ML(Figure 8B). The expression levels ofCL18477.Contig4_All and CL11190.Contig1_All in W was higher than in other stages. Among five stages, CL5146.Contig76_All was highest expressed in P1(Figure 7B). While CL17942.Contig2_All and CL817.Contig4_All were lowly expressed in P5, the expression level of CL16665.Contig4_All in P5 was higher than in other stages(Figure 8B). The expression levels of CL4540.Contig4_All and CL9391.Contig2_All in P10 was lower than in other stages. Meanwhile, unigene19996_All was mainly expressed in P10 (Figure 8B).
Expression levels of genes involved in 20E biosynthesis
Except for P1-vs-P5 which has one up-regulated gene (CL17221.Contig1_All, |log2Ratio|=2.1), there are no DEGs in the other three groups. A total of 2 (Unigene4431_All and CL17221.Contig1_All), 1 (CL14010.Contig2_All) and 3 (CL14010.Contig2_All, CL17221.Contig1_All and CL19241.Contig2_All) genes were down-regulated in W-vs-P1, P1-vs-P5 and P5-vs-P10, respectively (Table 2). The expression of DEGs in other groups is shown in Table S4. The expression levels of CL17221.Contig1_All in ML, W, and P5 were higher than in P1 and P10 (Figure 8C). Besides, CL19241.Contig2_All was mainly expressed in the P5 stage and the expression level of CL14010.Contig2_All was lowest in P10 than in other stages (Figure 8C).
Validation of the expression of DEGs using qRT-PCR
Correlation analysis uncovered that the relative expression levels of 10 selected CP, 20E, and JH-related genes obtained by real-time PCR quantitative (RT-qPCR) matched well with their FPKM values derived from RNA-seq (Figure 9).