Transcriptome Analysis Reveals the Gene Expression Changes in the Silkworm (Bombyx mori) in Response to Hydrogen Sulfide Exposure

Simple Summary The fat body is one of the most important tissues in the body of insects due to its number of functions. Nowadays the new physiological function of H2S has gained attention as a novel signaling molecule. H2S performs crucial regulatory functions involving growth, the cardiovascular system, oxidative stress, and inflammation in many organisms. In this study, RNA-seq technology was used to investigate the fat body of the silkworm at the transcriptional level after H2S exposure during the 5th larvae stage. A total of 1200 (DEGs) was identified after 7.5 µM H2S treatment, of which 977 DEGs were up-regulated and 223 DEGs were down-regulated. DEGs were mainly involved in the transport pathway, cellular community, carbohydrate metabolism, and immune-associated signal transduction. Present research provides new insights on the gene expression changes in the fat body of silkworms after H2S exposure. Abstract Hydrogen sulfide (H2S) has been recognized for its beneficial influence on physiological alterations. The development (body weight) and economic characteristics (cocoon weight, cocoon shell ratio, and cocoon shell weight) of silkworms were increased after continuous 7.5 µM H2S treatment. In the present study, gene expression changes in the fat body of silkworms at the 5th instar larvae in response to the H2S were investigated through comparative transcriptome analysis. Moreover, the expression pattern of significant differentially expressed genes (DEGs) at the 5th instar larvae was confirmed by quantitative real-time PCR (qRT-PCR) after H2S exposure. A total of 1200 (DEGs) was identified, of which 977 DEGs were up-regulated and 223 DEGs were down-regulated. Most of the DEGs were involved in the transport pathway, cellular community, carbohydrate metabolism, and immune-associated signal transduction. The up regulated genes under H2S exposure were involved in endocytosis, glycolysis/gluconeogenesis, the citrate cycle (TCA cycle), and the synthesis of fibroin, while genes related to inflammation were down-regulated, indicating that H2S could promote energy metabolism, the transport pathway, silk synthesis, and inhibit inflammation in the silkworm. In addition, the expression levels of these genes were increased or decreased in a time-dependent manner during the 5th instar larvae. These results provided insight into the effects of H2S on silkworms at the transcriptional level and a substantial foundation for understanding H2S function.


Introduction
The silkworm, Bombyx mori, is a pivotal model of Lepidoptera, a holometabolous insect, with a considerable economic importance in the world. The silkworm has been considered a model animal because of its appropriate life cycle, cultivation, and fast reproduction [1,2]. The fat body for insects acts as the dynamic tissue and the metabolic organ, mainly involved in synthesis, nutrient storage, and energy metabolism [3], and it has been widely studied in various studies of silkworm to reflect various biological, physiological, and biochemical processes. The over expression of BmFoxO in the fat body could suppress protein translation in transgenic silkworms [4].
The toxic effects of hydrogen sulfide (H 2 S) as an environmental pollutant have been studied in the past years. However, nowadays, the new physiological function of H 2 S has gained attention as a novel signaling molecule. H 2 S performs crucial regulatory functions involving growth, the cardiovascular system, oxidative stress, and inflammation in many organisms [5,6]. H 2 S has been accepted as the third transmission gas followed by nitric oxide (NO) and carbon monoxide (CO) in past years with applications to various insects. Miller and Roth (2007) found that treatment with high concentrations of H 2 S are toxic to Caenorhabditis elegans; However, exposure to low concentrations of H 2 S increased the thermotolerance and lifespan via increasing the activity of SIR-2.1 [7]. Budde and Roth (2010) found that hif-1 is required when C. elegans were exposed to H 2 S; in turn, H 2 S could increase both HIF-1 protein concentration and nuclear localization [8]. H 2 S is an endogenous regulator of oxidative damage, metabolism, and aging in C. elegans [9]. Wei and Kenyon (2016) found that removing germ cells in C. elegans triggers a level of hydrogen sulfide in nonreproductive somatic tissues. Hydrogen sulfide displayed the protective responses that slow aging [10]. H 2 S treatment significantly increased the survival and lifespan of Drosophila melanogaster under arid and food-free conditions [11,12].
In our previous research, it was found that the development (body weight) and economic characteristics (cocoon weight, cocoon shell ratio, and cocoon shell weight) of silkworm (Bombyx mori) were increased after continuous treatment with 7.5 µM of H 2 S. After exposure to H 2 S, the level of metabolites related to inflammation ((6Z,9Z,12Z)octadecatrienoic acid, hexadecanoic acid, choline phosphate, and malic acid, etc.) in the hemolymph of silkworms were significantly increased compared to control group [13]. However, when the concentrations of H 2 S were higher than 7.5 µM, H 2 S displayed toxic effects to the silkworm, including a decrease in body weight, cocoon weight, cocoon shell weight, etc. RNA-seq technology is widely used in the field of silkworms to investigate physiological and biochemical changes at the transcriptional level under diverse conditions [14]. However, there is no effective report for the effects of H 2 S on silkworms at the transcriptional level.
The present study aimed to investigate the expression level of genes in the fat body of silkworms after exposure to H 2 S. The results based on the transcriptomic and bioinformatic analysis showed significantly different responses to H 2 S and provided useful insights into the role of H 2 S on silkworms.

Silkworm Rearing and H2S Exposure
P50 strains of silkworm were obtained from the Sericulture Research Institute of the Chinese Academy of Agricultural Sciences, Zhenjiang, China. According to the previous study [13], from the 4th instar larvae to the 5th instar larvae, silkworms were divided into two groups (H 2 S-treated and control) and reared in desiccators using fumigating treatment at 25 ± 1 • C and 12L:12D photoperiod conditions. Silkworms in the H 2 S-treated group were treated by NaHS (Shandong West Asia Chemical Company, Jinan, China), a donor of H 2 S. The experimental concentration of H 2 S was set to 7.5 µM according to the report of Cao et al. (2020) [13]. There were triplicates in each group with 30 silkworms per replicate and the silkworms were fed abundant fresh mulberry leaves twice a day.

RNA Extraction
The total RNA of the silkworm fat body was extracted on ice from the 1st day (5L1D) to the 5th day (5L5D) at the 5th instar larvae using Trizol reagent (Invitrogen, CA, USA). The concentration and purity of total RNA were analyzed by NanoDrop ND-1000 (Nan-oDrop, Wilmington, DE, USA) and the ratios of A260/A280 and 260/230 were used to measure the quality of total RNA.

Sequencing, Data Processing, and Assembly
The total RNA of 5L3D was analyzed in the H 2 S-treated and control groups by RNA-Sequencing. Total RNA (10 µg) of each sample was subjected to Poly(A) mRNA isolation with poly-T oligo-attached magnetic beads (Invitrogen, CA, USA). The poly(A) RNA was fragmented into small pieces using a Magnesium RNA Fragmentation Module (New England Biolabs (Beijing) LTD., cat. e6150, Beijing, China) under 94 • C for 5-7 min. The cleaved RNA fragments were then reverse-transcribed to create the cDNA by SuperScript TM II Reverse Transcriptase (Invitrogen, CA, USA), which were next used to synthesize U-labeled second-stranded DNAs. An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (New England Biolabs (Beijing) LTD., cat.m0280, Beijing, China) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR by the following conditions: initial denaturation at 95 • C for 3 min; 8 cycles of denaturation at 98 • C for 15 s, annealing at 60 • C for 15 s, and extension at 72 • C for 30 s; and final extension at 72 • C for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. Lastly, we performed the 2 × 150 bp paired-end sequencing (PE150) on an Illumina Novaseq TM 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor's recommended protocol.
The transcriptome was sequenced on basis of the Illumina paired-end RNA-seq approach. The average insert size for the paired-end libraries was 300 bp (±50 bp). The low-quality reads with sequencing adaptors, sequencing primer, and nucleotide (q quality score (p value corrected by BH algorithm, also called FDR value, q value, p adj value) <20) were removed using Cutadapt software (https://cutadapt.readthedocs.io/en/stable/, accessed on 1 August 2021) [15] before assembly. We used HISAT2 software (https:// daehwankimlab.github.io/hisat2/,version:hisat2-2.0.4, accessed on 1 August 2021) [16] to map reads to the genome. Three important parameters, such as Q20 (the percentage of bases with mass value ≥20), Q30 (mass value ≥30), and GC content were verified to evaluate all the reads using the FastQC online tool [17].
The Bombyx mori reference genome [18] was used for the alignment of the sample reads using the HISAT package [19]. StringTie (http://ccb.jhu.edu/software/stringtie/, accessed on 1 August 2021) was used for the assembly of the mapped reads of each sample and a comprehensive transcriptome was constructed combining all transcriptomes using Perl scripts [20]. The expression level of mRNAs was used to calculate the fragments per kilobase of exon model per million mapped reads (FPKM) by StringTie (http://www. bioconductor.org/packages/release/bioc/html/ballgown.html, accessed on 1 August 2021) [20]. The FPKM could evaluate the abundance of genes. The DEGs were selected with log 2 (fold change) >1 or log 2 (fold change) <−1 at p-value < 0.05, and q-value > p-value by R package-edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html, accessed on 1 August 2021) [21].

Functional Annotation and Enrichment Analysis of DEGs
The reference sequences for all the silkworm genes were obtained from the NCBI database and BLASTX was used for annotation. All the sequences were aligned against Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using DIAMOND (0.7.12) [22]. GO function or KEGG pathway significant enrichment anal-ysis first mapped all the DEGs to each GO term in the GO database or each KEGG pathway in the KEGG database, and second calculated the number of genes in each GO term or KEGG pathway and finally applied a hypergeometric test to find a GO term or KEGG pathway that were significantly enriched in DEGs compared to the entire genome background.

Quantitative Real-Time PCR (qRT-PCR)
The qRT-PCR was used to validate the accuracy of DEGs in L5D3 at the transcriptome level and the expression pattern of DEGs in the 5th instar larvae was investigated. Reverse transcription was performed by HiScript II QRT SuperMix (Vazyme, Nanjing, China) as follows: total RNA (1 µg) (without genome contamination) and 4 × gDNA Wiper were mixed with 5 × HiScript II QRT SuperMix and then incubated at 25 • C for 10 min, 50 • C for 15 min, and 85 • C for 2 min [14]. The cDNA was used for real-time PCR detection using the UltraSYBR Mixture (Cwbio, Beijing, China) and qRT-PCR was executed using the Light Cycler ® 480 System (Roche, Basel, Switzerland). The primer sequences for target genes were listed in Table S1 using primer 3.0 software (Premier Biosoft International, Palo Alto, CA, USA). The expression level of action A3, a silkworm housekeeping gene, was regarded as an internal reference in standardization [14]. Before Quantitative real-time PCR, the efficiency and specificity of primers were confirmed. The volume of the reaction system was 20 µL, containing 10 µL UltraSYBR Mixture, 1 µL cDNA, 0.5 µL forward and reverse primers and 8 µL ddH 2 O. The qRT-PCR program was as follows: initial denaturation at 95 • C for 10 min; followed by 40 times of denaturation at 95 • C for 10 s, annealing at 60 • C for 30 s, and extension at 72 • C for 32 s; a melting curve program was set at 95 • C for 15 s, 60 • C for 1 min, 95 • C for 15 s, and a60 • C for 15 s with a final cooling step at 4 • C for 30 s [23]. The data were analyzed with the LightCycler ® 96 software (Roche Diagnostics, Indianapolis, IN, USA) using the 2 −∆∆Ct method.

Statistical Analysis
The data were analyzed using one-way ANOVA data analysis. All the data were obtained in triplicate and presented as mean ± standard error (SE) (n = 3). Statistical significance at p < 0.05 was measured using SPSS 20.0 software (IBM, Endicott, NY, USA).

Transcriptome Sequencing and Assembly
The deep sequencing of RNA from the fat body (5L3D) treated with H 2 S using the illumine sequencing and the data processing results were shown in Table 1. The raw data of H 2 S-treated and control groups were approximately 7.75 GB and 7.66 GB, respectively. After filtering and removing the low-quality reads, adaptors and poly-N stretches, the valid data and valid reads ratio of each replicate in H 2 S-treated group were obtained as follows: 7.55 Gb (99.12%), 7.87 Gb (98.91%), and 7.54 Gb (98.39%) and that of the control group were 7.42 Gb (94.00%), 7.37 Gb (98.99%), and 7.50 Gb (98.12%). The quality scores of Q30 exceeded 98% in each sample. Moreover, the mapped ratio in H 2 S-treated and control groups were 97.35%, 97.43%, and 97.72% and 93.19%, 97.35%, and 97.52%, respectively, which highly matched with the silkworm genome (Table S2). The expression interval of genes was mainly focused on the 0.3~3.57 FPKM interval (FI) and the ratio was about 35% gene of the control group and 39% of H 2 S-treated group (Table S3). These results revealed a good quality sequencing analysis and the genes were analyzed in the following section.

DEG Identification and Analysis
DEGs in H 2 S and control groups were calculated according to the FPKM value of the obtained genes and then a large number of genes were differentially expressed after H 2 S exposure. The volcano plot of DEGs was shown in Figure 1 and the detailed information of DEGs was listed in Table S4. A total of 1200 DEGs were obtained between the H 2 S-treated and control groups. The number of up-regulated and down-regulated DEGs was 977 and 223, and the up-regulated genes were higher than down-regulated genes.

DEG Identification and Analysis
DEGs in H2S and control groups were calculated according to the FPKM value of the obtained genes and then a large number of genes were differentially expressed after H2S exposure. The volcano plot of DEGs was shown in Figure 1 and the detailed information of DEGs was listed in Table S4. A total of 1200 DEGs were obtained between the H2Streated and control groups. The number of up-regulated and down-regulated DEGs was 977 and 223, and the up-regulated genes were higher than down-regulated genes. The red points represent the up-regulated genes, the blue points represent the down-regulated genes, and the grey points represent the genes without differential expression.

GO Analysis of DEGs
The GO database can provide standard terms for describing the comprehensive properties of genes in organisms as an international classification system. In this study, DEGs were assigned to different GO categories such as biological processes, cellular components, and molecular functions. The major subcategories of DEGs for biological processes were border follicle cell migration, mitotic spindle elongation, myoblast fusion and microtubule-based processes, as well as others. (Figure 2 and Table S5). Through the results of the GO analysis, it was revealed that H2S could influence these multiple biological processes by regulating the expression of genes in the fat body of silkworms at the molecular level.

GO Analysis of DEGs
The GO database can provide standard terms for describing the comprehensive properties of genes in organisms as an international classification system. In this study, DEGs were assigned to different GO categories such as biological processes, cellular components, and molecular functions. The major subcategories of DEGs for biological processes were border follicle cell migration, mitotic spindle elongation, myoblast fusion and microtubule-based processes, as well as others. (Figure 2 and Table S5). Through the results of the GO analysis, it was revealed that H 2 S could influence these multiple biological processes by regulating the expression of genes in the fat body of silkworms at the molecular level.

KEGG Pathway Analysis
To systematically analyze the biological function and metabolic pathways of D KEGG pathway annotations were obtained. The 561 genes of 1200 DEGs predicted t code enzymes were mapped into 225 KEGG pathways. These annotated genes wer signed to six categories including organismal systems (49 pathways), metabolism pathways), human diseases (39 pathways), genetic information processing (19 pathw environmental information processing (28 pathways), and cellular processes (16 ways) (Figure 3). In the organismal systems, vascular smooth muscle contraction (7 D is the most abundant metabolic pathway followed by the insulin signaling pathw DEGs) and toll-like receptor signaling pathway (5 DEGs). The most abundant meta category in metabolism is inositol phosphate metabolism (12 DEGs) followed by gly sis/gluconeogenesis (11 DEGs), glutathione metabolism (10 DEGs), and the citrate cy DEGs). In human disease, it is mainly related to insulin resistance (17 DEGs). In the gory of genetic information processing, it is mainly enriched in ribosomes (33 D closely related to synthetic genetic material. In environmental information processing most abundant pathways are the Hippo signaling pathway (19 DEGs), MAPK sign pathway (16 DEGs), and phosphatidylinositol signaling system (12 DEGs). In cellular cesses, 36 DEGs were mainly enriched in endocytosis followed by focal adhesio DEGs) and regulation of the actin cytoskeleton (15 DEGs).

KEGG Pathway Analysis
To systematically analyze the biological function and metabolic pathways of DEGs, KEGG pathway annotations were obtained. The 561 genes of 1200 DEGs predicted to encode enzymes were mapped into 225 KEGG pathways. These annotated genes were assigned to six categories including organismal systems (49 pathways), metabolism (74 pathways), human diseases (39 pathways), genetic information processing (19 pathways), environmental information processing (28 pathways), and cellular processes (16 pathways) ( Figure 3). In the organismal systems, vascular smooth muscle contraction (7 DEGs) is the most abundant metabolic pathway followed by the insulin signaling pathway (6 DEGs) and toll-like receptor signaling pathway (5 DEGs). The most abundant metabolic category in metabolism is inositol phosphate metabolism (12 DEGs) followed by glycolysis/gluconeogenesis (11 DEGs), glutathione metabolism (10 DEGs), and the citrate cycle (9 DEGs). In human disease, it is mainly related to insulin resistance (17 DEGs). In the category of genetic information processing, it is mainly enriched in ribosomes (33 DEGs) closely related to synthetic genetic material. In environmental information processing, the most abundant pathways are the Hippo signaling pathway (19 DEGs), MAPK signaling pathway (16 DEGs), and phosphatidylinositol signaling system (12 DEGs). In cellular processes, 36 DEGs were mainly enriched in endocytosis followed by focal adhesion (20 DEGs) and regulation of the actin cytoskeleton (15 DEGs).
A significant enrichment analysis was further performed on the metabolic pathways of DEGs, and a total of 26 significantly enriched metabolic pathways were obtained ( Table 2). The terms were enriched to multiple molecular pathways such as the transport pathway and cellular community, including endocytosis, focal adhesion, regulation of the actin cytoskeleton, tight junction, and adherens junction; amino acid and carbohydrate metabolisms, such as glycolysis/gluconeogenesis, the citrate cycle (TCA cycle), phenylalanine metabolism, and histidine metabolism; immune-associated signal transduction such VEGF signaling pathway, the TNF signaling pathway, NOD-like receptor signaling, the MAPK signaling pathway, the Hippo signaling pathway-fly, the Toll-like receptor signaling pathway, and the NF-κB signaling pathway (Table S6). The top 10 metabolic pathways affected by H 2 S were endocytosis, ribosome, focal adhesion, the Hippo signaling pathway, the MAPK signaling pathway, regulation of the actin cytoskeleton, tight junction, inositol phosphate metabolism, glycolysis/gluconeogenesis, and the TCA cycle (Table 2). These annotations indicated that H 2 S could influence the multiple biological pathways in the fat body of silkworms and provide a new perspective to the effects of H 2 S on silkworms. A significant enrichment analysis was further performed on the metabolic path of DEGs, and a total of 26 significantly enriched metabolic pathways were obtained ( 2). The terms were enriched to multiple molecular pathways such as the transport way and cellular community, including endocytosis, focal adhesion, regulation of t tin cytoskeleton, tight junction, and adherens junction; amino acid and carbohydrat tabolisms, such as glycolysis/gluconeogenesis, the citrate cycle (TCA cycle), phen nine metabolism, and histidine metabolism; immune-associated signal transduction VEGF signaling pathway, the TNF signaling pathway, NOD-like receptor signalin MAPK signaling pathway, the Hippo signaling pathway-fly, the Toll-like receptor s ing pathway, and the NF-κB signaling pathway (Table S6). The top 10 metabolic path affected by H2S were endocytosis, ribosome, focal adhesion, the Hippo signaling path the MAPK signaling pathway, regulation of the actin cytoskeleton, tight junction, in phosphate metabolism, glycolysis/gluconeogenesis, and the TCA cycle (Table 2). annotations indicated that H2S could influence the multiple biological pathways in t body of silkworms and provide a new perspective to the effects of H2S on silkworm

The Validation of DEGs by qRT-PCR
To validate the accuracy of DEGs obtained from the transcriptome between the H 2 Streated and control groups, 28 candidate genes were selected for qRT-PCR analysis. The selected DEGs were mainly annotated to glycolysis/gluconeogenesis, the TCA cycle, focal adhesion, endocytosis, adherens junction, the NF-κB signaling pathway, the MAPK signaling pathway, and the synthesis of fibroin. Good consistency between the results of qRT-PCR and the transcriptome confirmed the accuracy and reliability of the sequencing data and revealed the significant difference of these genes during H 2 S exposure ( Figure 4 and Table S3).   Table S1.

The Expression Pattern Analysis of DEGs
A total of 11 DEGs was representatively selected to analyze the expression leve tern in the 5th instar larvae using qRT-PCR. Genes (relating to the glycolytic pathwa TCA cycle) such as lactate phosphoglycerate kinase (PGK), phosphoglucomutase (P NADP-dependent isocitrate dehydrogenase (IDH), and pyruvate carboxylase (PC) H2S-treated group were up-regulated compared to control group ( Figure 4A, B). M while, the levels of these genes were significantly increased at the 5L3D, 5L4D, and comparing with 5L1D in the H2S group ( Figure 5A-D). In addition, the expression of the members of the Rab family of small GTPases (Rab10) and trehalose transpo (Tret1) in endocytosis were up-regulated after H2S treatment ( Figure 4C). The expr levels of Rab10 and Tret1 were significantly increased in a time-dependent manne  Table S1.

The Expression Pattern Analysis of DEGs
A total of 11 DEGs was representatively selected to analyze the expression level pattern in the 5th instar larvae using qRT-PCR. Genes (relating to the glycolytic pathway and TCA cycle) such as lactate phosphoglycerate kinase (PGK), phosphoglucomutase (PGM), NADP-dependent isocitrate dehydrogenase (IDH), and pyruvate carboxylase (PC) in the H 2 S-treated group were up-regulated compared to control group ( Figure 4A,B). Meanwhile, the levels of these genes were significantly increased at the 5L3D, 5L4D, and 5L5D comparing with 5L1D in the H 2 S group ( Figure 5A-D). In addition, the expression levels of the members of the Rab family of small GTPases (Rab10) and trehalose transporter 1 (Tret1) in endocytosis were up-regulated after H 2 S treatment ( Figure 4C). The expression levels of Rab10 and Tret1 were significantly increased in a time-dependent manner ( Figure 5E,F). On the contrary, transforming growth factor-β activated kinase1 (Tak1) and matrix metalloproteinase-3 (MMP-3) were significantly decreased in the 5th instar larvae after H 2 S exposure ( Figure 5G,H). Moreover, the increased expression of fibroin synthesis genes, such as the fibroin heavy chain (Fib-H), fibroin light chain (Fib-L), and glycoprotein P25 (P25) was observed during the 5th larvae stage after H 2 S exposure ( Figure 4D). At the end of the 5th instar larvae, the expression levels were higher than the early 5th instar larvae ( Figure 5I-K).  Table S1.

Discussion
H2S has been accepted as the third transmission gas followed by nitric oxide (NO) and carbon monoxide (CO) in recent years with applications to various insects, such as C.  Table S1.

Discussion
H 2 S has been accepted as the third transmission gas followed by nitric oxide (NO) and carbon monoxide (CO) in recent years with applications to various insects, such as C. elegans [7][8][9][10], D. melanogaster [11,12], and B. mori [13]. RNA-seq technology is widely used in the field of silkworms to investigate the physiological and biochemical changes at the transcriptional level under diverse conditions [4,14,24]. In our previous research, it was found that the development (body weight) and economic characteristics (cocoon weight, cocoon shell ratio, and cocoon shell weight) were increased after continuous treatment with 7.5 µM of H 2 S. However, there is no effective report for the effects of H 2 S on the silkworm at the transcriptional level [13]. In the present study, the expression levels of genes in the fat body of silkworms after exposure to H 2 S was investigated using transcriptomic and bioinformatic analysis. The results showed significantly different responses to H 2 S and provided useful insights into the role of H 2 S on silkworms.
Using transcriptome, a total of 1200 DEGs were identified in the fat body of silkworms at the 5th instar larvae. After treatment with H 2 S, up-regulated DEGs (977) were higher than the down-regulated DEGs (223). These DEGs were mainly involved in glycolysis/gluconeogenesis, endocytosis, and the TCA cycle. The carbohydrate metabolism is one of the important energy metabolism pathways in organisms, incorporating glycolysis/gluconeogenesis, the TCA cycle, oxidative phosphorylation, and the pentose phosphate pathway [25]. Energy metabolism, especially the metabolic control of carbohydrates, is essential for insect growth and development [26]. Glucose is oxidized to pyruvate in glycolysis/gluconeogenesis [25]. Glycolytic activity is regulated by many enzymes in the glycolytic pathway such as PGK and PGM. PGM can convert glucose 1 phosphate into glucose 6 phosphate, providing the raw materials for the glycolytic pathway [27]. PGK is commonly known as a soluble and membrane-bound protein in different organisms [27]. In one of the glycolytic processes, phosphorylation of 1,3-diphosphoglycerate to 3-phosphoglycerate is catalyzed by PGK [28]. Like glycolysis/gluconeogenesis, the TCA cycle also provides the energy and intermediates for diverse biosynthetic pathways [25]. IDH, one of three rate-limiting enzymes in the TCA cycle, catalyzes the TCA cycle reaction wherein isocitrate is oxidized and decarboxylated to produce α-ketoglutarate, NADH, and CO 2 [29]. PC catalyzes the carboxylation of pyruvate into oxaloacetate for gluconeogenesis, the urea cycle, lipid synthesis, and other pathways which are the primary anaplerotic pathways for the TCA cycle where the intermediates are replenished by alanine and lactate [30]. In this study, the results showed that the levels of PGM, PGK, IDH, and PC were up-regulated after H 2 S exposure. According to the expression pattern analysis of these genes, gene levels had gradually increased in a time-dependent manner. Especially for 5L4D and 5L5D, a significant difference was observed compared to 5L1D. These results indicated that the expression levels of the genes involved in carbohydrate metabolism were increased with H 2 S exposure, suggesting H 2 S can accelerate the carbohydrate metabolism of silkworms and provide much energy for silkworm metamorphosis.
The uptake and sorting of eukaryotes and the following recycling or degradation of fluid, membranes, membrane proteins, and macromolecules are regulated by endocytosis [31]. Endocytosis often depends on clathrin and occurs at clathrin-coated pits [32]. Utilizing endocytosis with clathrin-dependent or clathrin-independent uptake mechanisms, cargo proteins and adaptor molecules are delivered to early endosomes and sorted [33]. The Rab family of small GTPases can mediate membrane the trafficking processes of intracellular compartmentation, generation, and maintenance [34]. Rab10, which is closely related to clathrin-independent endocytosis is largely expressed in many model animals directing the trafficking of proteins. In Caenorhabditis elegans, Rab10 was required in interneuron postsynaptic glutamate receptor recycling and the transport between early endosomes and recycling endosomes [35]. Rab10 in the H 2 S-treated group was up-regulated compared to the control group, which indicated that the H 2 S-treated silkworms were more active in supporting endocytosis for the body. Tret1 participated in the transfer of trehalose, the major hemolymph sugar in most of the insect, which is synthesized in the fat body and released into the hemolymph. The transportation of trehalose with the up-regulated expression level of Tret1 after exposure to TiO 2 nanoparticles was accelerated in the silkworm, showing that TiO 2 nanoparticles could facilitate carbohydrate and the nutrition metabolism [36]. In this study, the level of Tret1 was up-regulated with the H 2 S supplement to accelerate the transport rate of trehalose. Moreover, the expression levels of Rab10 and Tret1 were increased in a time-dependent manner during the 5th larvae. These results showed that H 2 S could accelerate the transmission of mass in the silkworm fat body via up-regulating the expressions of Rab10 and Tret1 to provide material for larva-pupae metamorphosis.
Previous studies reported that H 2 S exerted anti-inflammatory effects on T lymphocytes. Rats injected with NaHS inhibited inflammatory processes such as leukocyte infiltration and adherence to the vascular endothelium [37]. In this study, DEGs were annotated to various pathways related to inflammation by KEGG analysis, including the TNF signaling pathway, the MAPK signaling pathway, and the NF-κB signaling pathway.
Tak1 is a serine/threonine kinase and a member of the mitogen-activated protein kinase kinase (MAP3K) family which is activated by receptors such as TGF-β, TNF-α, and IL-1. MMP-3 is a matrix-soluble protein containing the binding site of NF-κB and can participate in the degradation of ECM and activate the precursors of MMPs. The content of MMP-3 was increased in the mouse arthritis model and the administration of baicalin can reduce the gene expression level of MMPs and alleviate inflammatory processes [38]. After H 2 S exposure, levels of Tak1 and MMP-3 were down-regulated compared to the control group. Moreover, the expression levels of Tak1 and MMP-3 at the end of 5th instar larvae were significantly decreased compared to the early 5th instar larvae. These results indicated that H 2 S might have an anti-inflammatory role in silkworm growth.
In addition, significant changes in the genes related to fibroin synthesis were observed in this study. The domestic silkworm has tremendous economic value and the most important manifestation of its value is silk production [39]. The function of the silk gland is to synthesize silk protein and during the 5th instar larvae it is the biggest organ of B. mori [39]. It was reported that Fib-H, Fib-L, and P25 encoding fibroin were expressed in all the larval instars. However, their highest expression was noticed before the last molting period [40]. The expression levels of Fib-L, P25, and Fib-H were up-regulated in response to H 2 S compared to the control group. In addition, the expression of Fib-H at the end of the 5th instar larvae (5L4D, 5L5D) was rapidly increased and higher than on the other days. Similarly, the expression levels of P25 and Fib-L were also significantly increased. These results indicated that H 2 S could promote the fibroin synthesis by up-regulating the expression levels of Fib-L, P25, and Fib-H, which was beneficial to silk production. These results were similar to the previous study which reported the up-regulated expression of Fib-L, P25, and Fib-H after treatment with titanium dioxide nanoparticles (TiO 2 NPs) [41].

Conclusions
To summarize, this is the first-ever report to investigate the fat body of silkworms at the transcriptional level after H 2 S exposure using the transcriptome analysis during the 5th larvae stage. Transcriptome analysis could detect numerous genes and signaling pathways and most of the DEGs were up-regulated compared to the control group. The expression of key genes related to energy metabolism, transport pathways, and fibroin synthesis in silkworms after H 2 S exposure was up-regulated; conversely, the genes involved in inflammatory processes were down-regulated, which indicated that H 2 S may lead to the promotion of biological processes such as energy and nutrition metabolism and inhibition of inflammatory processes. These results provide a novel approach to comprehending the molecular mechanisms underlying the H 2 S effects in the silkworm fat body.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/insects12121110/s1, Tables S1-S6. Table S1. Primer's sequence of genes used in qRT-PCR analysis. Table S2. Statistical analysis of the RNAs Libraries. Table S3. Distribution of gene expression levels between H 2 S-treated and the control group and the partial DEGs.