De novo transcriptome sequencing and analysis revealed the molecular basis of rapid fat accumulation by black soldier fly (Hermetia illucens, L.) for development of insectival biodiesel

Background Black soldier fly (BSF, Hermetia illucens L.) can efficiently degrade organic wastes and transform into a high fat containing insect biomass that could be used as feedstock for biodiesel production. Meanwhile, the molecular regulatory basis of fat accumulation by BSF is still unclear; it is necessary to identify vital genes and regulators that are involved in fat accumulation. Results This study analyzed the dynamic state of fat content and fatty-acid composition of BSF larvae in eight different stages. The late prepupa stage exhibited the highest crude fat, with lauric acid being the main component. Therefore, to provide insight into this unexplained phenomenon, the molecular regulation of rapid fat accumulation by BSF larvae was investigated. The twelve developmental stages of BSF were selected for transcriptome analysis, including the eight stages used for investigation of fat content and fatty-acid composition. By Illumina sequencing, 218,295,450,000 nt were generated. Through assembly by Trinity, 70,475 unigenes were obtained with an average length of 1064 nt and an N50 of 1749 nt. The differentially expressed unigenes were identified by DESeq, with 9159 of them being up-regulated and 10,101 of them were down-regulated. The several putative genes that are involved in the formation of pyruvate, acetyl-CoA biosynthesis, acetyl-CoA transcription, fatty-acid biosynthesis, and triacylglycerol biosynthesis were identified. The four vital metabolic genes that are associated with fat accumulation were validated by quantitative real-time PCR (qRT-PCR). The molecular mechanism of fat accumulation in BSF was clarified in this investigation through the construction of a detailed fat accumulation model from our results. Conclusion The study provides an unprecedented level of insight from transcriptome sequencing to reveal the crude fat accumulation mechanism in developing BSF. The finding holds considerable promise for insectival biodiesel production, and the fat content and fatty-acid composition can be altered by genetic engineering approaches in the future for the insect production industry. Electronic supplementary material The online version of this article (10.1186/s13068-019-1531-7) contains supplementary material, which is available to authorized users.


Background
The world's available fuel resources are gradually depleting due to an increase in the energy demand; moreover, the use of fossil fuels is devastating our environment through greenhouse gases (GHG) emission and global warming. Therefore, the search for alternative and sustainable energy sources has attracted widespread attention in recent years to reduce the global environmental problem (GHG emission and global warming) and protect fossil energy resources. Currently, several alternatives are being studied and implement. The biodiesel is the most promising alternative fuel that provides environmental benefits, since there use leads to the decrease in the harmful emission of GHG and their affects [1]. Moreover, it was noted that biodiesel is renewable, nontoxic and biodegradable energy, with performance and characteristic that is very close to fossil diesel fuel [2]. As biodiesel is mainly produced from crop oil, the problems of limited feedstock, expensive production cost and competition with food resources have prevented its largescale applications [3]; therefore, it is important to develop economical feedstock for biodiesel production.
The increase in the global human population resulted in a rise in the level of organic wastes products such as food waste, animal manure, and other agricultural wastes. Moreover, these organic wastes not only contribute to nutrient imbalances of soil results in soil quality deterioration; they also lead to water and air pollution [4]. Currently, the standard practices in the management of these wastes include disposal in landfills, combustion, or agricultural applications such as soil conditioner [5]. These organic wastes managemental practices have caused more critical environmental concerns, such as the creation of leachate and landfill gases, pest attraction, and water and air pollution [6,7] considerably.
A more environmentally friendly and economically feasible approach is to take advantages of insects to convert organic wastes into larval biomass and biofertilizer through a process called biotransformation [6]. The black soldier fly larvae (BSFL) is a widely distributed insect can consume various organic material, from animal manure to vegetables and fruits [7]. The life cycle of BSF includes four phases: egg, larvae, pupa and adult stage (Fig. 1a). Since BSFL can transform various low-value organic wastes into protein and fat [8][9][10], and the fat content have 20-40% [11,12]. So insect fat biomass can be used for biodiesel production that has many valuable properties, such as diesel compatibility, high energy content and cetane number, as was previously demonstrated by Zheng et al. [13,14] and Rehman et al. [6,15]. In our previous studies, fat extracted from BSFL is composed mainly of three fatty acids (FA), including lauric acid (C12:0), palmitic acid (C16:0) and oleic acid (C18:1) [6,14]. Due to a large amount of saturated fatty acid and short carbon chain such as lauric acid (C12:0), the biodiesel produced by BSFL fat has high oxidative stability [6,16]. The performance measurements for biodiesel produced by BSFL fat have shown fuel properties that are consistent with the international standard of biodiesel (EN14214) [15]. Therefore, the high-fat content and waste conversion efficiency of BSFL make it an ideal feedstock for biodiesel production [6].
The advances in sequencing technology have been considerable in recent years [17], such methods, especially those utilizing RNA sequencing (RNA-Seq), allow for examination of mechanisms responsible for metabolic variation [18]. These methods have been utilized on many insects, such as oriental fruit fly, mosquito Culex pipiens, and vinegar flies Drosophila melanogaster, to determine genes that encode enzymes associated with fatty acid (FA) and triacylglycerol (TAG) biosynthesis [19][20][21][22]. Although there are many researches on using BSFL fat for biodiesel production, whereas the molecular regulatory basis of fat accumulation by BSFL is still unclear. Thus, the transcriptomic analysis is essential to clarify the molecular regulatory mechanism of rapid fat accumulation by BSFL for the development of insectival biodiesel.
The present work studied the dynamic changes of fat content and FA composition in developing BSFL at eight stages, including one-day-old larvae (1-d-L), four-dayold larvae (4-d-L), eight-day-old larvae (8-d-L), twelveday-old larvae (12-d-L), early prepupa (E-prepupa), late prepupa (L-prepupa), early pupa (E-pupa), and late pupa (L-pupa). Moreover, the four additional stages are also added for RNA-Seq, including early egg (E-egg), late egg (L-egg), female adult (F-adult), and male adult (M-adult), resulting in a total of twelve stages. The Illumina highthroughput sequencing was performed to obtain unigene annotations, screen differentially expressed unigenes by the DESeq method, and identify temporal expression patterns of putative genes involved in FA and TAG biosynthesis in BSFL. The designated routes of fat accumulation are based on transcriptome data from Illumina sequencing and temporal expression analysis. To reveal the molecular regulatory mechanism of rapid fat accumulation, the FA and TAG biosynthesis pathways in BSF were constructed to form a fat accumulation model that includes detailed information on when and how important genes are expressed and the relationships between its different components. Experimental validation of four vital genes was performed through quantitative realtime PCR (qRT-PCR). The results will provide a useful resource for functional genomics and enrich the public database. These resources will provide a foundation to increase fat accumulation and alter the FA composition of BSFL in the future.

Dynamic patterns of fat content and FA compositions in developing BSF
To investigate the dynamic patterns of fat accumulation in developing BSFL, we analyzed the crude fat (CF) content in eight different stages (1-d-L, 4-d-L, 8-d-L, 12-d-L, E-prepupa, L-prepupa, E-pupa, and L-pupa) (Fig. 1b). It was noted that the CF was lower in BSFL during early stages (1-d-L to 4-d-L) and late stages (E-pupa to L-pupa), whereas higher fat content was recorded at mature larvae stage (12-d-L) and prepupa stages (E-prepupa to L-prepupa). A remarkable increase of CF was recorded at 4-d-L to 8-d-L from 7.5% to 17.8%, with the highest CF (26.1%) at E-prepupa stage. However, a sharp decline (17.2%) in CF was noted at E-pupa to L-pupa after the E-prepupa stage (26.1%).
Moreover, to explore FA composition in developing BSFL, the dynamic spectra of FAs accumulation in various development phases were analyzed (Fig. 2). The relative proportion of lauric acid (C12:0) was found to be large during developmental stages. The decline of lauric acid (C12:0) was recorded from the 1-d-L stage (72.3%) Fig. 2 a Fatty-acid composition of BSF at different stages of development, b fatty acids in small quantity in BSF to the 4-d-L stage (7.7%), at the same time, the proportion of palmitic acid (C16:0), oleic acid (C18:1) and linoleic acid (C18:2) was increased. Whereas the lauric acid (C12:0) showed a rapid increase from 4-d-L to L-prepupa, with a peak value at L-pupa (75.5%), while the proportion of oleic acid (C18:1), palmitic acid (C16:0) and linoleic acid (C18:2) were decreased from 4-d-L to L-prepupa. Moreover, the fluctuation of myristic acid (C14:0) and oleic acid (C18:0) was recorded from 4-d-L to L-pupa. It was noted 90.4% FAs present in developing BSF being short-chain FAs such as C14:0 (myristic acid), C16:0 (palmitic acid) and C12:0 (lauric acid); therefore, it was concluded that BSFL achieves rapid fat accumulation by synthesizing short-chain FAs early in its development, which makes them an ideal feedstock for high-performance biodiesel production.

Illumina sequencing and de novo assembly of developmental BSF
To investigate the molecular regulatory mechanism of rapid fat accumulation in developing BSF was investigated; therefore, RNA was extracted from E-egg, L-egg, 1-d-L, 4-d-L, 8-d-L, 12-d-L, E-prepupa, L-prepupa, E-pupa, L-pupa, F-adult, and M-adult. The 24 cDNA libraries were constructed and deeply sequenced by the Illumina HiSeq ™ X Ten system with two replicates in each stage, generating 218,295,450,000 nt from these RNA-Seq samples. After the filtering step, an average of 60,637,625 nt was obtained as clean reads in each sample, with an average Q20 percentage of 97.44% and an average GC percentage of 38.4% (Additional file 1: Table S1). After assembly, 70,475 unigenes were obtained with an N50 of 1749 nt, with the total length of unigenes being 74,988,057 nt and the average length of unigenes being 1064 nt (Additional file 2: Table S2). The length of unigene sequences ranged mainly from 300 to 3000 nt, with 28,121 (39.9%) in the range of 1000 to 3000 nt, 11,335 (16.08%) longer than 2000 nt, and the number of unigenes decreased as the length of sequences increased (Additional file 3: Figure S1). These results indicated that the assembly is of high quality. All clean reads were deposited in the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) database under accession number PRJNA506627.

Functional annotation and classification of BSF unigenes
To investigate the function of assembled unigenes in developing BSF, a total of 70,475 unigenes were matched to public databases, including NR (NCBI nonredundant protein sequences), Gene Ontology (GO), Swiss-Prot (a manually annotated and reviewed protein sequence database), Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups of proteins (COG), and NT (NCBI non-redundant nucleotide sequences). A total of 41,375 (58.7%) unigenes had matches with known genes (Table 1) Similarity analysis between the unigenes of BSF and NR was performed using BLAST (Additional file 4: Figure S2). The results exhibited that 51.9% of the annotated unigenes had strong homology with e value below 1e −45 (Additional file 4: Figure S2. 4A). There were 36.6%, 26.2% and 14.1% of putative proteins showing 40-60%, 60-80% and 80-100% of similarity with known proteins in NR, respectively (Additional file 4: Figure S2. 4B). From the species distribution of NR BLAST matches (Additional file 4: Figure S2. 4C), 24.1% of unigenes had strong homology with Drosophila. When compared to other species within Diptera, 12.3% of unigenes had matches to sequences from Aedes aegypti, followed by Culex pipiens quinquefasciatus (6.7%) and Anopheles gambiae PEST (5.4%) (Additional file 4: Figure S2. 4C), which are all mosquitoes.

Analysis of differentially expressed unigenes of developing BSF
The differential expression patterns of specific unigenes that are associated with BSF fat accumulation were investigated by calculating the 'fragment per kilobase per million' (FPKM) value. A false discovery rate (FDR) ≤ 0.05 was imposed with the absolute value of log 2 (ratio) ≥ 1 to screen differentially expressed unigenes from all assembled BSF unigenes (Additional file 8: Table S3). The differentially expressed unigenes were matched into the GO database (Additional file 9: Table S4) and the KEGG database (Additional file 10: Table S5). Differentially expressed unigenes were concentrated during the early stages and the late stages, with 6911 of them identified during early stages (1-d-L and 4-d-L) and 8793 of them identified during late stages (E-pupa and L-pupa), while very few differentially expressed unigenes were identified at 8-d-L and 12-d-L (Fig. 4a). When the differentially expressed unigenes that are involved in lipid metabolism were screened and analyzed, they were also concentrated during the early stages and the late stages, with 220 of them identified at 1-d-L and 4-d-L, and 262 of them identified at E-pupa and L-pupa (Fig. 4b). As a result, Five unigenes were found to have sustained up-regulation in early stage (1-d-L, 4-d-L and 8-d-L), while four unigenes were found to have sustained up-regulation in late stage (L-prepupa, E-pupa and L-pupa) (Fig. 4c) The distribution of down-regulated unigenes that are involved in lipid metabolism was also analyzed (Fig. 4d). Two unigenes were found to have sustained down-regulation in early stage (1-d-L, 4-d-L and 8-d-L). Among them, one unigene was related to AKR1B, and another unigene was related to diacylglycerol kinase (ATP) (DGK, EC:2.7.1.107) (CL3020.Contig2_All).

Expression patterns of enzymes involved in pyruvate and acetyl-CoA formation in developing BSF
To investigate the expression patterns of genes associated with pyruvate formation, putative genes that are related to enzymes required for glycolysis were obtained from Illumina sequencing analysis. Among the 122 putative genes that are associated with glycolysis in BSF, 12 of them were related to hexokinase (HK, EC:2.7.1.1, e value: 6e −9 to 0), 5 of them were related to ADP-dependent glucokinase (ADPGK, EC:2.7.1.147, e value: 2e −18 to 5e −140 ), 12 of them were related to glucose-6-phosphate isomerase (GPI, EC:5. Temporal transcript analysis was performed to examine the dynamic expression patterns of putative genes that are involved in glycolysis. The putative genes from glycolysis pathway were highly expressed during early and late stages of BSFL development (Fig. 5a). When compared the temporal profile of FPKM for putative genes encoding isozymes HK, ADPGK, gpmA and gpmI that are involved in glycolysis, HK had higher expression in early stage, while ADPGK was up-regulated at L-prepupa (Fig. 6a), meanwhile, the expression patterns of gpmI showed more consistent with CF accumulation patterns (Fig. 6b).
Pyruvate dehydrogenase complex (PDC) is an important enzyme for acetyl-CoA formation. There were 8, 7, 6 and 11 putative genes that are related to PDC subunits E1-α, E1-β, E2 and E3, respectively, with e value ranging from 6e −16 to 3e −170 , 5e −25 to 4e −172 , 5e −6 to 0 and 2e −6 to 0, respectively. The expression levels of these genes gradually decreased before 8-d-L, but increased from 8-d-L, with the highest expression level at L-prepupa (Fig. 5b). Since the putative genes that are associated with acetyl-CoA formation are also putative genes involved in glycolysis (Fig. 5a), the PDC subunits can respond to BSFL development, with acetyl-CoA production for CF accumulation occurring mainly in early stage. Although the PDC subunits had high expression level in late stages such as E-pupa and L-pupa, the CF content did not increase significantly. As a result, acetyl-CoA was used to provide energy rather than CF accumulation in late stage.

Expression patterns of enzymes involved in acetyl-CoA transportation and FA biosynthesis in developing BSF
By Illumina sequencing analysis, 8 and 5 putative genes were identified relation to the citrate synthase (CS, EC:2.3.3.1, e value: 3e −6 to 0) and ATP-citrate lyase (ACLY, EC:2.3.3.8, e value: 0). Acetyl-CoA can be catalyzed by CS to produce citric acid by condensation of oxaloacetic acid, and citrate is preferentially exported to the cytosol via the tricarboxylate transporter. Similarly, we analyzed the dynamic expression patterns of putative genes that are associated with acetyl-CoA transportation, the putative genes of CS and ACLY had similar expression patterns as pyruvate and acetyl-CoA biosynthesis (Fig. 5c).

Expression patterns of enzymes involved in triacylglycerol synthesis in developing BSF
Triacylglycerol (TAG) biosynthesis begins with acyl-CoA formation. Two isozymes were identified by Illumina sequencing analysis in the acyl-CoA biosynthesis pathway, with 10 putative genes that were related to long-chain-fatty-acid-CoA ligase (ACSBG, EC:6.2.1.3, e value: 7e −11 to 0) and 6 putative genes that were related to long-chain acyl-CoA synthetase (ACSL, EC:6.2.1.3, e value: 1e −17 to 0). ACSBG had high expression level in developing BSFL from the temporal transcript analysis of ACSBG and ACSL (Fig. 6c). By Illumina sequencing analysis, 54 putative genes for TAG biosynthesis were identified. Temporal transcript analysis showed that the putative genes for TAG biosynthesis were highly expressed during early stages (1-d-L to 4-d-L) and late stages (E-prepupa to E-pupa) (Fig. 5e). Since differential expression analysis indicated that TAG degradation occurs mainly during the early stages, these results indicated that rapid TAG accumulation occurs mainly during the late stages.
In the third step of TAG biosynthesis, phosphatidate is dephosphorylated to provide 1,2-diacylglycerol (DAG) for the biosynthesis of TAG. Two isozymes of phosphatidate phosphatase were identified by Illumina sequencing analysis, with 11 putative genes that were related to phosphatidate phosphatase (PLPP1_2_3, EC:3.1.3.4, e value: 2e −27 to 8e −105 ), and 8 putative genes that were related to phosphatidate phosphatase LPIN (LPIN, EC:3.1.3.4, e value: 3e −12 to 0). Temporal transcript analysis of PLPP1_2_3 and LPIN showed that PLPP1_2_3 had high expression level, and the expression patterns are consistent with the ones from putative genes involved in TAG biosynthesis (Fig. 6f ).

Experimental validation and analysis of key genes involved in BSF fat accumulation
To assess the accuracy of sequencing and assembly of the BSF transcriptome, the relative expression levels and temporal transcript patterns of the putative genes which involved in fat accumulation were analyzed. Four putative genes of vital enzymes, including FAS, ACC , ACSBG and DGAT1, were selected for qRT-PCR (Additional file 12: Table S6). The results from qRT-PCR showed that the relative expression levels of these selected genes were mostly consistent with the FPKM comparative ratios (with 1-d-L as the control) (Fig. 7). These results indicated that the unigene assembly is accurate and reliable, and it is feasible to use the DESeq method to select differentially expressed gene. Both enzymes ACC and FAS that are involved in FA biosynthesis had high expression level at the 4-d-L stage.

Discussion
BSFL develops on organic wastes, reducing environmental pollution and converting waste into biodiesel from the larval fat body. Therefore, BSF fat as feedstock for biodiesel production is studied widely; moreover, the investigation of the molecular regulatory basis of rapid fat accumulation by BSF is pivotal for biodiesel development. Furthermore, the identification of vital genes and regulators that are involved in fat accumulation will optimise the conversion to biodiesel. The fat content and FA composition of BSFL in different stages were compared (Figs. 1, 2). Larvae early in the development process have not accumulated fat, while those late in the development process (i.e., E-pupa and L-pupa) have utilized fat reserves to complete metamorphosis, which is in agreement with previous results [11]; meanwhile, the sudden decline of lauric acid from the 1-d-L stage (72.3%) to the 4-d-L stage (7.7%) may be due to its utilization as stored energy at 1-d-L. The dynamic variation of lauric acid (C12:0) during BSFL development indicated that it could be important in energy and metabolism. Similar results showed that other flies (Diptera) were characterized to have a very high proportion of lauric acid (C12:0), palmitoleic acid (C16:1), and decanoic acid (C10:0) [23,24].
The fat accumulation and regulation in BSF at the molecular level, sequencing of BSF transcriptome were performed to comprehensively analysis the transcriptome profile in developing BSF. Since the reference genome of BSF is not available, to obtain the reliable results of assembly and annotation, exception the eight stages (1-d-L, 4-d-L, 8-d-L, 12-d-L, E-prepupa, L-prepupa, E-pupa and L-pupa) that we investigated the dynamic patterns of fat content and FA composition, four additional stages (E-egg, L-egg, F-adult and M-adult) are also added, resulting in a total of twelve stages was performed deep sequencing, the greater sequencing depth and accurate sequence reads provide reliable and effective data of BSF transcriptome profile.
The analysis of differentially expressed gene is a powerful approach to identify the key regulator underlying important biological processes. In this study, we identified five enzymes (lip, LPL, CES1, UGT and GLB1) that have sustained up-regulation in early stage, three enzymes (AKR1B, ELOVL4 and HSD17B4) that have sustained up-regulation in late phase (Fig. 4c). It is noteworthy that the enzymes up-regulated expression in early stage primary involvement lipid and carbohydrate metabolism, these results indicated that lip, LPL, CES1, UGT and GLB1 play vital role in energy supply. The previous researches showed that UGT is a very important detoxifying enzyme that participates in host defense against endogenous toxins and xenobiotic chemicals, and is discovered in various living organisms from bacteria to humans [25]. In humans, CES1 was found to have fattyacid ethyl ester (FAEE) synthesizing activity, in which long-chain fatty acids are transesterified with ethanol to generate FAEEs [26]; moreover, it was formerly noted that CES1 appears to be participated in some processes important to cell biology, like cholesterol trafficking [27]. In addition to its catalytic actions, it was recorded that CES1 appears to play critical role for protein release and retention from the endoplasmic reticulum (ER) [28], and the UGT phase II drug metabolism enzymes can form a complex with CES1 to hold UGTs in the ER lumen [29]. These results indicated that UGT and CES1 play vital roles in BSFL lipid metabolism and defense against toxins. AKR1B was proved to play a significant role in the biosynthesis of ascorbic acid and develop scurvy [30]. Whereas in this investigation AKR1B was discovered sustained down-regulation in early stage and supported up-regulation in late stage. Since it is well known that ascorbic acid is an antioxidant that suppresses osteoclast activity, the ALR1B may play a vital role in BSF eclosion. Although ELOVL4 and HSD17B4 were proved involvement in fatty-acid elongation and biosynthesis of unsaturated fatty acids, and both enzymes exhibited upregulation in late stage. Whereas increase in the content of long-chain FA and unsaturated fatty acids (Fig. 2) was not observed, because there is not enough substrate for the biosynthesis of long-chain fatty acids and unsaturated fatty acids, as BSF stops eating from the prepupa stage [31].
The down-regulated unigenes that are associated with lipid metabolism was also analyzed (Fig. 4d). Two enzymes (ydfG and ALDH) were found to have sustained down-regulation in late stage (E-prepupa, L-prepupa, E-pupa and L-pupa). Although the major function of ALDH is NAD(P)+-dependent aldehyde oxidation, previous studies indicated that ALDHs exhibit multiple functions, such as ALDH1A1, ALDH2, ALDH3A1 and ALDH4A1 are known to catalyze ester hydrolysis, suggesting that ALDHs may have more than one catalytic role [32]. Since both ydfG and ALDH convert aldehydes into acids in developing BSFL, the down-regulation of ydfG and ALDH in late-stage indicated that the conversion of aldehydes to acids was important for lipid metabolism in early and middle stage.
With pyruvate and acetyl-CoA as the substrate for FA biosynthesis, the putative genes from pyruvate and acetyl-CoA formation pathway and acetyl-CoA transportation pathway showed similar expression patterns (Fig. 5a-c). The substrate formation and transportation for FA biosynthesis in the early stage may play an essential role in fat accumulation in developing BSFL. The putative genes encoding isozymes from pyruvate formation pathway were identified, the first step in metabolism of glucose is phosphorylation that is catalyzed by HK and ADPGK [33,34], the temporal profile of HK and ADPGK suggested that HK may play a vital role in pyruvate formation and fat accumulation, while ADPGK plays an essential role in energy supply, because BSF needs energy for tissue differentiation and eclosion in late stages (Fig. 6a), whereas the expression patterns of gpmI are consistent with CF accumulation, and gpmI may play an important role in glycolysis (Fig. 6b).
Several putative genes from FA biosynthesis in developing BSF were also identified. The analysis of canonical pathways of fatty-acid biosynthesis showed that FASN catalyzes a series of reactions in this pathway (Additional file 11: Figure S6), this phenomenon indicated that FASN plays a significant role for FA biosynthesis in developing BSF. Although the putative genes from FA biosynthesis pathway had higher expression level in late stage, CF content did not increase. Since BSF does not eat at the prepupa stage, there is not enough substrate to accumulate fat. At the pupa stage, energy is needed for tissue differentiation, and fat is consumed to provide energy. Although the putative genes had high expression level in E-pupa to L-pupa stages, CF content is low.
The putative genes from pyruvate and acetyl-CoA formation pathway, acetyl-CoA transportation pathway, FA and TAG biosynthesis pathway showed similar expression patterns (Fig. 5), as above-mentioned the temporal profile of vital genes involved in fat accumulation, these results indicated BSF possesses the regular fat metabolism pattern, and the fat accumulation in early stages play an important role in BSF development.
The expression profiles of metabolic enzymes that are involved in the biosynthesis of acetyl-CoA, FA and TAG were analyzed to investigate the fat accumulation mechanism in BSF. A detailed fat accumulation model was constructed from our results (Fig. 8). It starts with glycolysis and results in the formation of pyruvate and acetyl-CoA, followed by the transportation of acetyl-CoA to the cytosol. The finding of this investigation will help to understand the molecular regulatory basis of rapid fat accumulation in BSF, and provide a foundation to increase content of fat without compromising growth by metabolic engineering, it is an essential aspect of advancing the economic feasibility of BSF waste conversion technology. It was previously noted the knockdown of a multifunctional lipase/phospholipase/acyltransferase in the diatom Thalassiosira pseudonana increased lipid content without affecting growth [35]. Moreover, to redirect carbon flux toward fat biosynthesis, the deletion of glycerol-3-phosphate dehydrogenase isomer in Yarrowia lipolytica exhibited a threefold remarkably increase in lipid content [36]. Thus, the gene data sets were developed in this investigation lay a foundation for metabolic engineering to increase the content of fat in BSF and also provide a foundation to develop molecular markers for fat breeding studies.

Conclusion
The dynamic patterns of fat accumulation and FA compositions in different development stages of BSFL were analyzed. It was noted that the E-prepupa showed the highest fat content during BSF metamorphosis. Moreover, FA analysis indicated that lauric acid (C12:0) is the major component in BSF fat body. Furthermore, 24 samples selected for in-depth transcriptome sequencing analysis, a total of 70,475 unigenes were assembled with an average length of 1064 nt from clean reads. By applying the DESeq method, 14,597 differentially expressed unigenes were identified.
The expression profiles of metabolic enzymes that are involved in the biosynthesis of acetyl-CoA, FA and TAG were analyzed to investigate the fat accumulation mechanism in BSF. The regulatory enzymes that play important roles in fat accumulation in BSF, including HK, ADPGK, PFK, CS, ACLY, ACC, FASN, ACSBG, AGPAT3_4, PLPP1_2 and DGAT1 were characterized (Additional file 13: Table S7). The FA and TAG biosynthesis pathways were reconstructed in BSFL, and the temporal expression levels of four vital genes were validated through quantitative real-time PCR (qRT-PCR).
A detailed fat accumulation model was constructed from our results (Fig. 8). It starts with glycolysis and results in the formation of pyruvate and acetyl-CoA, followed by the transportation of acetyl-CoA to the cytosol. The rapid accumulation of predominantly short-chain FAs are produced from acetyl-CoA and ligated to form long-chain acyl-CoA, lipogenesis process ends with TAG formation. Whereas many genes that are involved in these processes are both highly expressed during early stages (1-d-L to 4-d-L) and late stages (E-pupa to L-pupa), it was noted that the CF content does not increase significantly during late stages. Therefore, most of these processes up to FA biosynthesis occur during the early stages. Moreover, it was observed that the TAG degradation occurs mainly during the early stages by differentially expressed analysis; this evidence proved that the rapid TAG accumulation occurs mainly during the late stages.
The finding of this investigation will help to understand the molecular regulatory basis of rapid fat accumulation in BSF. It will provide a useful resource for functional genomics studies, and enrich the public database. Meanwhile, the gene data sets were developed in this investigation lay a foundation for metabolic engineering to increase the content of fat in BSF and also provide a foundation to develop molecular markers for fat breeding studies.

Insect materials and extraction of BSF fat
The colony of BSF was maintained in a greenhouse at the National Engineering Research Center of Microbial Pesticides, Huazhong Agricultural University (HZAU), Wuhan, Hubei, China, for about 10 years. The commercial chicken feed used to rear BSF was produced by Charoen Pokphand Group, Wuhan, China. The nutritional composition of this chicken feed was listed in (Additional file 14: Table S8).
Fresh eggs of BSF were collected from the colony. The eggs were immediately frozen in liquid nitrogen for 90 s, marked as E-egg (< 12 h), and stored at − 80 °C until they were used. A batch of eggs was placed in incubators at 26 °C with 65-70% relative humidity (RH) for hatching, with neonate larvae developed on the chicken feed. The eggs marked as L-egg (< 72 h) were collected after 2 days from hatching eggs and frozen in liquid nitrogen for 90 s and stored at − 80 °C. The egg hatchings were monitored, and the larvae that hatched within 3 h were collected and frozen in liquid nitrogen for 90 s, marked as one-day-old (1-d-L), and stored at − 80 °C. The remaining neonate larvae were separated into five groups and reared with 150 g of moist chicken feed having 70% water content at 26 °C with 60% RH. Larvae samples at various developmental stages (4-d-L, 8-d-L and 12-d-L) were collected using the previously described method. The newly emerged prepupae were removed promptly every day. Prepupae were collected from the 2nd day after preparation started and frozen in liquid nitrogen for 90 s, marked as E-prepupae, and stored at − 80 °C. The remaining prepupae were left for pupation at 26 °C with 60% RH. The pupae were observed, immediately removed and collected. The freshly emerged pupae were selected and marked as L-prepupae and E-pupa, frozen in liquid nitrogen for the 90 s and stored at − 80 °C. The remaining pupae were monitored for adult emergence. Once adult emerged, late pupae were collected and marked as L-pupa, frozen in liquid nitrogen for 90 s and stored at − 80 °C. The adult female and male flies that emerged within 12 h were selected and marked as F-adult and M-adult, frozen in liquid nitrogen for the 90 s stored at − 80 °C. The fat content and fatty-acid composition of BSF were analyzed concerning China National Standards GB/T 5009.6-2003 and GB/T 17376-2008. In brief, to measure the fat content, the 2 g sample was put into the 50 ml test tube with 10 ml of hydrochloric acid and 8 ml of water. The sample was digested completely in 70-80 °C water bath, the 10 ml ethanol was added into the digested sample, following the digested sample was transferred to a 100 ml mix cylinder with stopper, the cylinder was added of 5 ml petroleum ether-ether mixture and allowed to stand. The supernatant was taken out from the cylinder, the residue was collected and transferred to the flask on a water bath, the flask was evaporated repeatedly to constant weight, finally, the total fat content was measured (GB/T 5009. . To analysis the fatty-acid composition, 0.5 g sample was put into 50 ml flask with 6 ml of sodium hydroxide methanol solution, the sample was saponified. To obtain the methyl ester solution, seven ml of boron trifluoride was added in the flask and extracted for 3 min. Finally, the FAs were determined by gas chromatography analysis using 0.2 μl of methyl ester solution (GB/T 17376-2008).

cDNA library preparation and sequence data analysis and assembly
Twelve developmental stages (E-egg, L-egg, 1-d-L, 4-d-L, 8-d-L, 12-d-L, E-prepua, L-prepua, E-pupa, L-pupa, F-adult, and M-adult) in BSF were selected for transcriptome sequencing and analysis. With two biological replicates per stage, 24 samples were constructed. RNA was extracted with Trizol reagent (Vazyme BioTech Co. Ltd, Nanjing, China) and treated with DNase I. The mRNA was isolated by magnetic beads with Oligo (dT) and mixed with the fragmentation buffer to obtain the short fragments. The cDNA synthesis was performed using the mRNA fragments as templates. The short fragments purified from the cDNA synthesis were used for end reparation, single nucleotide A (adenine) addition and adaptors connection. Following, the fragments with suitable length were selected for the PCR amplification as templates; meanwhile, the quantification and qualification of the sample libraries were measured by Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System. Finally, the libraries were sequenced using Illumina HiSeq ™ X Ten.
Since the raw reads produced from sequencing machines contain dull reads with adaptors, unknown or low-quality bases, they were filtered before bioinformatics analysis: (1) remove adaptor sequences; (2) remove reads with more than 5% unknown nucleotides; and (3) remove low-quality reads with more than 20% of nucleotides having quality value < 10. Clean reads were assembled by the Trinity software (version: release-20130026) [37].