98 ISJ 18: 98-107, 2021 ISSN 1824-307X RESEARCH REPORT Transcriptomic analysis of strain-specific and gender-specific response of silkworm to BmNPV infection S He, J Xu, Y Fan, F Zhu, K Chen* Institute of Life Sciences; Jiangsu University, Zhenjiang 212013, China This is an open access article published under the CC BY license Accepted June 18, 2021 Abstract Bombyx mori nucleopolyhedrovirus (BmNPV) is one of the main pathogens causing serious economic losses in sericulture. However, the molecular mechanism of silkworm resistance to BmNPV is still largely unclear, and the differences in the anti-BmNPV response between silkworms of different genders have been rarely studied. In this study, BmNPV resistant strain NB and BmNPV sensitive strain 306 of different genders were used as experimental materials to inoculate BmNPV, and their transcriptomes were sequenced to analyze their response to BmNPV. Eighteen genes specifically differentially expressed in NB after BmNPV inoculation were finally obtained through transcriptomic analysis, fourteen of which were up-regulated and four were down-regulated, suggesting that they might be related to BmNPV resistance. Among them, the expression abundance of eight genes were higher in males than in females, and one gene was in the contrary. These genes suggested that there were certain differences in the anti-BmNPV response between silkworms of different genders. This study provided a new understanding of the molecular mechanism of silkworm resistance to BmNPV and the differences in the anti-BmNPV response between silkworms of different genders, and laid a foundation for future prevention and control of BmNPV. Key Words: Bombyx mori; BmNPV resistance; transcriptome sequence; strain-specific response; gender-specific response Introduction The silkworm, Bombyx mori has been domesticated for more than 5000 years and is an important economic insect in many developing countries (Goldsmith et al., 2005). It is also a good model Lepidoptera insect that is often used in the study of insect genetics and immunology (Tanaka et al., 2009; Guo et al., 2015). Bombyx mori nucleopolyhedrovirus (BmNPV) is a baculovirus that infects silkworms. Although some silkworm strains with BmNPV resistance like NB and CVDAR have been cultivated (Chen et al., 1991; Qing et al., 2019)，their cocoon yield and quality still do not meet the production requirements. So far, the main silkworm strains used for production like Yuke9 and Yuncan10 are BmNPV-sensitive strains (Liu et al., 2013; Yang et al., 2019). There is no effective method or drug to control this virus disease, which brings huge losses to silkworm industry every year. Therefore it is very urgent and important to identify _________________________________________ Corresponding author: Keping Chen Institute of Life Sciences Jiangsu University Zhenjiang, Jiangsu 212013, China E-mail: kpchen@ujs.edu.cn the BmNPV resistant mechanism of silkworm. There have been many researches on silkworm resistance to BmNPV, and some resistance genes have been successively identified. For example, the increased expression of suppressor of cytokine signaling 2 (BmSOCS2) has been reported to be correlated with suppression of BmNPV replication in silkworm (Yuan et al., 2020). Arginine kinase (BmAK) and trypsin (Bmtryp) have also been reported to be involved in the antiviral process (Kang et al., 2011; Ponnuvel et al., 2012). However, the molecular mechanism of silkworm resistance to BmNPV has not been clearly interpreted yet. In nature, male and female individuals of many species have great differences in many aspects (Matthews et al., 2019). For example, there is sexual dimorphism in genetic loci for anthropometric traits (Randall et al., 2013). And the fungal infection could induce sex-specific transcriptional changes in Silene latifolia (Zemp et al., 2015). Similarly, there are certain differences in the development, cocoon characters and disease resistance between silkworms of different genders (Qin et al., 2014; Chen et al., 2016). However, the differences in the anti-BmNPV response between silkworms of different genders have been rarely studied. 99 Table 1 Primers used in RT-qPCR GENE Forward Primer Reverse Primer Product Length (bp) BGIBMGA004955 5' ATCGAGATTCGGTCACGAGC 3' 5' TTCAGCTGATCGCGCCAATA 3' 169 BGIBMGA009799 5' ACCTGGCTAGGACTGAACGA 3' 5' GGTGGTGCCACACCTTTGTA 3' 220 BGIBMGA006251 5' GGGAGGTTGCAACGGTGAT 3' 5' CGGGGCTGATACATTTGCCT 3' 173 BGIBMGA004809 5' ACTCGAAGTGGCTCTCAACG 3' 5' AGCTTCAATAGCTGCCGTGT 3' 205 BGIBMGA012486 5' AGCGGTTGCAGTTTCTACGA 3' 5' CAGTTCGTGCATACCCCAGT 3' 168 BGIBMGA010514 5' TTGCAAGCCGTTTTTGCTGT 3' 5' TCTCCGGCCTCCAGTAGTAG 3' 188 BGIBMGA010863 5' TGCTGCATAAGATGCGGAGA 3' 5' TGTAGCACCTGGCTACTTGG 3' 219 BGIBMGA011868 5' TGCCAGGCTAACACAGACAG 3' 5' CTCTGAGCCTGCACTTCCAT 3' 150 BGIBMGA009091 5' CAGCCAGGGTTCGGTTGAAA 3' 5' TCAACGGGTACGCATTTTCC 3' 180 BmRPS3 5' CGATTCAACATTCCAGAGCA 3' 5' GAACACCATAGCAAGCACGAC 3' 142 The primers used in RT-qPCR ,and their product length were listed above. BmRPS3 was used as reference gene in this study In this study, we conducted transcriptome sequencing of NB (BmNPV-resistant strain) and 306 (BmNPV-sensitive strain) in a gender specific manner. Through comparative transcriptomic analysis, differentially expressed genes in silkworms after BmNPV inoculation were obtained. And some candidate genes that may be related to BmNPV resistance were finally obtained according to Venn analysis. This study provides a new understanding of the molecular mechanism of silkworm resistance to BmNPV and the differences in the anti-BmNPV response between silkworms of different genders. Materials and methods Virus and silkworm The BmNPV and silkworms were all preserved in our laboratory. The virus was obtained from the hemolymph of infected larvae and purified by repeated and differential centrifugation according to the method of Rahman et al. (2004), and the virus concentration (OB/mL) was measured by the hemocytometer. The half lethal concentration (LC50) of NB reached 6.39×108 OB/mL, which was nearly 1000 times higher than that of 306. The first four instar larvae were reared with fresh mulberry leaves. On the first day of fifth instar, the PBS solution containing BmNPV polyhedrons (1×10^7 OB/mL) was evenly smeared on the surface of fresh mulberry leaves and fed to the experimental group. After 24 hours feeding with virus, the experimental group was refed with normal mulberry leaves. And the control group was set simultaneously. On the third day of the fifth instar, the whole silkworms NB♂-V, NB♀-V, 306♂-V, 306♀-V of the experimental group and NB♂-C, NB♀-C, 306♂-C, 306♀-C of the control group were stored in the -80 °C refrigerator for later use. RNA extraction and transcriptome sequencing Trizol method was used to extract total RNA from the eight samples, with three biological replicates for each sample. The purity, concentration and integrity of RNA were determined by NanoDrop ONE/ONEC (Thermo, USA), Invitrogen Qubit 4 (Thermo, USA), Agilent 2100 BioAnalyzer (Agilent Technologies, USA) and agarose gel electrophoresis. After that, mRNA was enriched with Oligo(dT) magnetic beads to construct a cDNA library. After the library was qualified with Agilent 2100 BioAnalyzer (Agilent Technologies, USA), transcriptome sequencing was performed using Illumina Novaseq 6000 (Illumina, USA). Transcriptome assembly, annotation and gene expression analysis The raw reads were filtered by removing reads with adaptor, reads with a N ratio of more than 10 % (N denotes bases that cannot be identified), and reads containing more than 50 % bases with a Qphred ≤ 20 to obtain the clean reads. Then the clean reads were aligned with the reference genome (http://metazoa.ensembl.org/Bombyx_mori/Info/Inde x) by HISAT2 to obtain the mapped reads for subsequent transcriptome assembly and gene expression analysis. Stringtie software was used to assemble the mapped reads based on the reference genome. GO, KEGG, COG, NR, Swiss-Prot and Pfam databases were used to annotate the assembled transcriptome. The expression abundance of the annotated genes was analyzed using RSEM software. Read counts data of genes was obtained from the alignment result and annotation file. Standardized gene expression abundance was obtained by FPKM (Fragments per kilobase of transcript per million Fragments) transformation (Mortazavi et al., 2008). http://metazoa.ensembl.org/Bombyx_mori/Info/Index http://metazoa.ensembl.org/Bombyx_mori/Info/Index 100 Table 2 Gene annotation results statistics Database Expressed Gene number percent All Gene number percent GO 4388 32.54 % 4531 28.78 % KEGG 7162 53.12 % 7501 47.64 % COG 12652 93.84 % 13372 84.93 % NR 13125 97.34 % 14097 89.53 % Swiss-Prot 9403 69.74 % 9609 61.03 % Pfam 10083 74.78 % 10311 65.49 % Total-anno 13187 97.80 % 14219 90.31 % Total 13483 100.00 % 15745 100.00 % Expressed Gene number: Numbers of expressed genes in this study in each database All Gene number: Numbers of all expressed genes of silkworm in each database Total-anno: Total expressed genes that have been annotated in all databases Total: Total expressed genes in all databases Differentially expressed genes (DEGs) analysis Gene expression abundance was compared between the experimental group and the control group using DEseq2 software to obtain up-regulated DEGs and down-regulated DEGs. The comparison schemes were as follows: NB♂-V vs NB♂-C, NB♀-V vs NB♀-C, 306♂-V vs 306♂-C, 306♀-V vs 306♀-C. P-adjust < 0.05 and |log2FC| ≥ 1 were set as standard to screen DEGs. Venn diagrams were constructed to obtain unique DEGs of NB♂-V and NB♀-V, as well as DEGs that were common to NB♂-V and NB♀-V but not to 306♂-V or 306♀-V. Further analysis was made on the DEGs that were common to NB♂-V and NB♀-V but not to 306♂-V or 306♀-V to select genes that differentially expressed between NB♂-V and 306♂-V, as well as between NB♀-V and 306♀-V. Genes with low expression abundance (FPKM<10) in all samples were not considered, and genes with a ratio ≥ 1.35 were considered to be up-regulated, genes with a ratio ≤0.74 were considered to be down-regulated. The genes differentially expressed in NB♂-V and NB♀-V were screened from these genes to analyze the differences in the anti-BmNPV response between NB of different genders. The functions of DEGs and candidate genes were classified by GO analysis, including molecular functions, cellular components and biological processes. Real-time quantitative PCR (RT-qPCR) analysis Nine genes that might be related to BmNPV resistance were selected and their relative expression levels in NB-V, NB-C, 306-V and 306-C (mixed samples including male and female silkworms) were detected by RT-qPCR. And relative expression levels of 4 genes were detected in NB♂-V and NB♀-V as well. All the Primers were listed in Table 1. RT-qPCR reactions were prepared with the HiScript Q RT SuperMix Kit (Vazyme, Nanjing, China), following the manufacturer’s instruction. Reactions were carried out in Applied Biosystems QuantStudio 3 Real-Time PCR System (Thermo Fisher, Nanjing, China). All samples were performed in triplicate. B.mori ribosomal protein s3 (BmRPS3) gene was used as a reference gene. And the relative expression levels were calculated using the 2-∆∆Ct method (Livak et al., 2001). Results Overview of transcriptome sequencing results A total of 172.4 GB of clean reads were obtained for transcriptomic analysis. For all libraries, the GC content was about 46 %, and the Q30 was all higher than 92.12 %, indicating that the sequencing data were of sufficient quality and accuracy for further analysis. Most reads of all samples were mapped successfully with the reference genome, with a uniquely mapped ratio of about 80 % (S1 Table). There was no significant difference in the sequencing results among the biological replicates, indicating that the library construction and sequencing results were qualified. Table 3 DEGs statistics Experiment Control Up-regulated Down-regulated DEGs 306♀-V 306♀-C 741 1410 2151 306♂-V 306♂-C 251 363 614 NB♀-V NB♀-C 328 222 550 NB♂-V NB♂-C 644 279 923 Total 3205 Standard: P-adjust < 0.05 & |log2FC| ≥ 1 101 Fig. 1 Venn diagrams showing the DEGs of different samples. (A) Venn diagram of up-regulated DEGs; (B) Venn diagram of down-regulated DEGs Gene annotation and expression analysis Reference sequence alignment was performed in GO, KEGG, COG, NR, Swiss-Prot and Pfam databases, and 13,187 of the 13,483 expressed genes were successfully annotated (Table 2). And gene expression abundance was then analyzed. Differentially expressed genes (DEGs) analysis A total of 3205 DEGs were obtained by comparing the gene expression of the experimental group and the control group (Table 3). In both male and female, NB had more up-regulated DEGs than down-regulated DEGs after BmNPV inoculation, while 306 were on the contrary. And the number of DEGs was higher in NB male than in female, while that in 306 was higher in female than in male. In order to further analyze which DEGs were related to BmNPV resistance, Venn diagrams were performed on the up-regulated and down-regulated DEGs respectively (Fig 1). According to the Venn diagrams, 588 specific DEGs were found in NB♂-V (S2 Table), of which 474 were up-regulated and 114 were down-regulated. While NB♀-V had 305 specific DEGs (S3 Table), of which 244 were up- r egulated 102 Fig. 2 Gene ontology (GO) functional annotation of specific DEGs. (A) GO functional annotation of specific DEGs of NB♂-V; (B) GO functional annotation of specific DEGs of NB♀-V. The X axis shows the number of genes, and the Y axis shows the categories of gene functions, including molecular functions, cellular components, and biological processes and 61 were down-regulated. GO analysis showed that specific DEGs of NB♂-V and NB♀-V had similar functional categories, both of which were mainly involved in binding, metabolic process, membrane, catalytic activity and others (Fig 2). However, the number of specific DEGs involved in all functions of NB♂-V was about twice that of NB♀-V, which was consistent with the ratio of their total DEGs (NB♂-V : NB♀-V = 923 : 550). Meanwhile, 68 DEGs (53 up-regulated and 15 down-regulated) that were common to NB♂-V and NB♀-V but not to 306♂-V and 306♀-V were also obtained according to the Venn diagrams (S4 Table). Go analysis revealed that 8 genes were related to metabolic process, 5 genes were related to membrane part, 8 genes were related to binding, and 12 genes were related to catalytic activity (Fig 3). By analyzing the expression abundance of these 68 genes, 18 genes that differentially expressed between NB♂-V and 306♂-V as well as between NB♀-V and 306♀-V were obtained. Among them, 4 genes were down-regulated and 14 were up-regulated in NB males and females after BmNPV inoculation (Table 4). These genes might be involved 103 Fig. 3 Gene ontology (GO) functional annotation of 68 specific genes of NB males and females Table 4 Expression abundance of 18 genes Gene ID NR annotation 306♂-V FPKM 306♀-V FPKM NB♂-V FPKM NB♀-V FPKM NB♂-V vs 306♂-V ratio NB♀-V vs 306♀-V ratio BGIBMGA010563 juvenile hormone acid O-methyltransferase-like 26.31 51.06 2.29 3.82 0.09 0.07 Down-regulated BGIBMGA004955 serine protease inhibitor 13 precursor 17.48 13.3 4.84 5.25 0.28 0.39 BGIBMGA002176 None 26.38 311.93 7.22 5.26 0.27 0.02 BGIBMGA001147 facilitated trehalose transporter Tret1-like 20.3 17.44 1.19 1.2 0.06 0.07 BGIBMGA014348 uncharacterized protein LOC106069352 3699.11 108.64 7868.94 3157.47 2.13 29.06 Up-regulated BGIBMGA012002 sericin 3 precursor 5463.96 179.71 16584.36 6671.33 3.04 37.12 BGIBMGA001213 None 108.52 2.26 1161.74 112.25 10.71 49.67 BGIBMGA009261 uncharacterized protein LOC105842476 0.1 0 94.66 9.37 916.06 BGIBMGA009799 aldo-keto reductase AKR2E4-like 112.76 93.66 288.16 182.82 2.56 1.95 BGIBMGA006251 zonadhesin-like 8.47 3.61 16.05 8.16 1.89 2.26 BGIBMGA004219 uncharacterized protein LOC101743399 7.81 8.29 48.88 47.61 6.26 5.74 BGIBMGA004809 uncharacterized protein LOC105842986 70.93 66.23 121.04 90.48 1.71 1.37 BGIBMGA012486 cytochrome b5-related protein 20.06 42.64 44.91 95.09 2.24 2.23 BGIBMGA010285 uncharacterized protein LOC101738995 4.13 7.83 18.49 24.06 4.48 3.07 BGIBMGA010514 glucose dehydrogenase [FAD, quinone] isoform X1 1.36 0.81 10.95 1.72 8.05 2.14 BGIBMGA010863 aspartate--tRNA ligase, mitochondrial 6.5 6.33 53.93 40.12 8.29 6.34 BGIBMGA011868 uncharacterized protein LOC101742906 1.53 1.67 17.83 17.47 11.63 10.48 BGIBMGA009091 fungal protease inhibitor F-like 925.76 186.12 1814.67 1295.84 1.96 6.96 The ratio represents the fold change of FPKM values: genes with a ratio ≥ 1.35 were considered to be up-regulated, genes with a ratio ≤ 0.74 were considered to be down-regulated 104 Table 5 Expression abundance of 9 genes Gene ID NR annotation NB♂-V NB♀-V NB♂-V vs FPKM FPKM NB♀-V ratio BGIBMGA012486 cytochrome b5-related protein 44.91 95.09 0.47 Down-regulated BGIBMGA014348 uncharacterized protein LOC106069352 7868.94 3157.47 2.49 Up-regulated BGIBMGA012002 sericin 3 precursor 16584.36 6671.33 2.49 BGIBMGA001213 None 1161.74 112.25 10.35 BGIBMGA009261 uncharacterized protein 94.66 9.37 10.1 BGIBMGA009799 aldo-keto reductase AKR2E4-like 288.16 182.82 1.58 BGIBMGA006251 zonadhesin-like 16.05 8.16 1.97 BGIBMGA010514 glucose dehydrogenase [FAD, quinone] isoform X1 10.95 1.72 4.61 BGIBMGA009091 fungal protease inhibitor F-like 1814.67 1295.84 1.4 in the anti-BmNPV response of silkworm. By comparing the expression abundance of these 18 genes in NB♀-V and NB♂-V, nine genes were found to have different expression abundance, of which 8 genes had higher expression abundance in NB♂-V and 1 gene had higher expression abundance in NB♀-V (Table 5). These genes denote that NB of different genders have certain differences in the anti-BmNPV response. RT-qPCR validation of DEGS The relative expression levels of 9 genes in NB-V, NB-C, 306-V, and 306-C were detected by RT-qPCR (Fig 4A), as well as the relative expression levels of 4 genes in NB♀-V and NB♂-V (Fig 4B). The RT-qPCR results were consistent with the transcriptome data, confirming the reliability of transcriptome sequencing result in this study. Discussion In this study, we obtained some genes involved in the silkworm response to BmNPV infection through transcriptome sequencing, so as to study the mechanism of silkworm resistance to BmNPV and the differences in the anti-BmNPV response between silkworms of different genders. By comparing the gene expression between the experimental group and the control group, some differentially expressed genes that might be related to the BmNPV resistance were successfully identified. Among them, BmTret1 has already been reported to be involved in the anti-BmNPV response (Yang et al., 2016), but its resistance mechanism has not been clarified. In many non-mammals, Tret1 can transport exogenous trehalose into cells and induce autophagy (Sarkar et al., 2007), while the autophagy of silkworm cells is conducive to BmNPV infection (Wang et al., 2017). Therefore, we speculated that the down-regulation of BmTret1-like may play a role in the anti-BmNPV response by regulating autophagy. A variety of serine proteases have also been proved to have anti-BmNPV activity (Nakazawa et al., 2004; Li et al., 2017). Among the candidate genes obtained in this study, both serine protease inhibitor 13 and fungal protease inhibitor F-like were the inhibitors of serine proteases (Pham et al., 1996). And they might be involved in the anti-BmNPV process by interacting with serine proteases. Viruses must rely on cellular proteins to complete replication in cells, so the protein metabolism of the host plays an important role in the game between host and virus (Emmett et al., 2005). In this study, some candidate genes related to protein metabolism were also identified. For example, aspartate-tRNA ligase could regulate protein translation (Ibba et al., 2000; Ribas et al., 2000). Cytochrome B5 could significantly regulate the function of Cytochrome P450, thus extensively affecting the metabolism of exogenous and endogenous compounds (Zhang et al., 2007; Im et al., 2011). GO analysis showed that LOC105842986 might have lyase activity and participate in the transport and metabolism of amino acids, and LOC101738995 might be involved in coenzyme binding and catalysis activity. The up-regulated expression of these genes suggested that they might influence the replication of BmNPV by regulating protein metabolism of silkworm. 105 Fig. 4 A) RT-qPCR results of 9 genes. B) RT-qPCR results of 4 genes. The X axis represents different samples, and the Y axis shows relative gene expression levels. The significance of changes in different gene expression levels was also indicated. Samples used for qRT-PCR of nine genes were mixed samples of female and male silkworms, and samples used for qRT-PCR of four genes were female and male silkworm samples respectively Ecdysone and juvenile hormone are the two most important hormones in insects, which can jointly regulate the growth, development and metamorphosis process of insects (Herboso et al., 2015; Riddiford et al., 2020). In this study, Aldo-keto reductase-like (AKR2E4-like) was up-regulated, and juvenile hormone acid O-methyltransferase-like (JHAMT-like) was down-regulated in NB after BmNPV inoculated. AKR2E4 could catalyze the production of ecdysone (Yamamoto et al., 2017), while JHAMT was the enzyme that catalyzes the last step of juvenile hormone biosynthesis in lepidopteran insects (Shinoda et al., 2003). Therefore, changes in their expression abundance may influence the biosynthesis of ecdysone and juvenile hormone. Interestingly, in this study, the expression abundance of Sericin 3 precursor and BGIBMGA009261 (contained Fibroin P25 domain) in NB were up-regulated after BmNPV inoculation, which might be caused by the changed titer of ecdysone and juvenile hormone. In drosophila, steroid hormone signaling played a key role in regulating innate immunity and fighting bacterial infection (Regan et al., 2013). Ecdysone could closely regulate innate immunity to recognize and defend against bacterial infection by controlling the expression of pattern recognition receptor （PGRP-LC）in drosophila (Rus et al., 2013). While juvenile hormone was an immunosuppressive factor that could strongly interfere with this ecdysone-dependent immune enhancement (Flatt et al., 2008). However, these two hormones have not 106 been found to be involved in the immune system of silkworm. In this study, the expression abundance of AKR2E4-like and JHAMT-like were changed in NB after BmNPV inoculation, suggesting that juvenile hormone and ecdysone might participate in the regulation of immune response and the anti-BmNPV process in silkworm. The pupation time of NB males was about 24 hours earlier than that of females under normal rearing conditions, and this phenomenon still maintained in NB of different genders after BmNPV inoculation. Therefore, a number of genes must be differentially expressed between NB of different genders after BmNPV inoculation. However, it is not clear whether there are differences in the anti-BmNPV response between silkworms of different genders. In this study, NB of different genders produced different numbers of DEGs and different specific DEGs after BmNPV inoculation, indicating that they did have different responses to BmNPV. Although GO analysis showed that the functional classifications of their specific DEGs were similar, there were differences in the genes involved in the anti-BmNPV response and their expression abundance. We screened nine genes that differentially expressed in NB of different genders after BmNPV inoculation. Among them, The expression abundance of cytochrome b5-related protein was lower and fungal protease inhibitor F-like and AKR2E4-like was higher in NB♂-V than those in NB♀-V. These genes have been discussed above.,and they might participate in the anti-BmNPV response by regulating protein metabolism, protein interaction or other ways. While the higher expression abundance of Sericin 3 precursor, BGIBMGA009261 (contained Fibroin P25 domain), glucose dehydrogenase [FAD, quinone] isoform X1 and zonadhesin-like in NB♂-V might be related to the faster growth and development of NB males. These differences in genetic background between silkworms of different genders were still maintained after BmNPV inoculation, but whether they were involved in the anti-BmNPV response remains to be further studied. In conclusion, this study provides some insights into the mechanism of silkworm resistance to BmNPV and the gender-specific differences in anti-BmNPV response of silkworm. Acknowledgments This work was supported by National Natural Science Foundation of China (31861143051 and 31572467). Reference Chen KP, Lin CQ, Wu XD, Yao Q, Fang QQ. Resistance of preservative Bombyx mori strains to nuclear polyhedrosis virus (in Chinese). Sci Sericul, 17: 45-46, 1991. Chen SL, Yang RG, Gao JH, Liao PF, Li YM, Gao X, et al. Difference analysis of different varieties and genders of silkworm infected by Microsporidium isolated from Pieris rapae L(in Chinese). Hubei Agricul Sci, 55(19): 5108-5110, 2016. Emmett SR, Dove B, Mahoney L, Wurm T, Hiscox JA. The cell cycle and virus infection. Methods Mol Biol. 296: 197-218, 2005. Flatt T, Heyland A, Rus F, Porpiglia E, Sherlock C, Yamamoto R, et al. Hormonal regulation of the humoral innate immune response in Drosophila melanogaster. J Exp Biol. 211(Pt 16): 2712-2724, 2008. Goldsmith MR, Shimada T, Abe H. The genetics and genomics of the silkworm, Bombyx mori. Annu Rev Entomol. 50: 71-100, 2005. Guo R, Wang S, Xue R, Cao G, Hu X, Huang M, et al. The gene expression profile of resistant and susceptible Bombyx mori strains reveals cypovirus-associated variations in host gene transcript levels. Appl Microbiol Biotechnol. 99(12): 5175-5187, 2015. Herboso L, Oliveira MM, Talamillo A, Pérez C, González M, Martín D, et al. Ecdysone promotes growth of imaginal discs through the regulation of Thor in D. melanogaster. Sci Rep. 5: 12383, 2015. Ibba M, Soll D. Aminoacyl-tRNA synthesis. Annu Rev Biochem. 69: 617-650, 2000. Im SC, Waskell L. The interaction of microsomal cytochrome P450 2B4 with its redox partners, cytochrome P450 reductase and cytochrome b(5). Arch Biochem Biophys. 507(1): 144-153, 2011. Kang L, Shi H, Liu X, Zhang C, Yao Q, Wang Y, et al. Arginine kinase is highly expressed in a resistant strain of silkworm (Bombyx mori, Lepidoptera): Implication of its role in resistance to Bombyx mori nucleopolyhedrovirus. Comp Biochem Physiol B Biochem Mol Biol. 158(3): 230-234, 2011. Li G, Zhou Q, Qiu L, Yao Q, Chen K, Tang Q, et al. Serine protease Bm-SP142 was differentially expressed in resistant and susceptible Bombyx mori strains, involving in the defence response to viral infection. PLoS One. 12(4): e0175518, 2017. Liu M, Wu KJ, Tian MH, Huang P, Shen ZL, Dong ZP，et al. Breeding of silkworm variety ‘Yuncan 10’ with double sex-limited markings for spring rearing. Southwest China J.Agric.Sci. 26(05): 2147-2152, 2013. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method. Methods. 25(4): 402-408, 2001. Matthews G, Hangartner S, Chapple DG, Connallon T. Quantifying maladaptation during the evolution of sexual dimorphism. Proc Biol Sci. 286(1908): 20191372, 2019. Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 5(7): 621-628, 2008. Nakazawa H, Tsuneishi E, Ponnuvel KM, Furukawa S, Asaoka A, Tanaka H, et al. Antiviral activity of a serine protease from the digestive juice of Bombyx mori larvae against nucleopolyhedrovirus. Virology. 321(1): 154-162, 2004. 107 Pham TN, Hayashi K, Takano R, Nakazawa H, Mori H, Ichida M, et al. Expression of Bombyx family fungal protease inhibitor F from Bombyx mori by baculovirus vector. J Biochem. 119(6): 1080-1085, 1996. Ponnuvel KM, Nithya K, Sirigineedi S, Awasthi AK, Yamakawa M. In vitro antiviral activity of an alkaline trypsin from the digestive juice of Bombyx mori larvae against nucleopolyhedrovirus. Arch Insect Biochem Physiol. 81(2): 90-104, 2012. Qin L, Shi H, Xia H, Chen L, Yao Q, Chen K. Comparative proteomic analysis of midgut proteins from male and female Bombyx mori (Lepidoptera: Bombycidae). J Insect Sci. 14: 226, 2014. Qin Q, Dong ZQ, Lei XJ, Cao MY, Tang L, Shi MN, et al. Analysis of resistance characteristics to BmNPV of silkworm strain CVDAR (in Chinese). Acta Microbiology. 59(12): 2390-2400, 2019. Randall JC, Winkler TW, Kutalik Z, Berndt SI, Jackson AU, Monda KL, et al. Sex-stratified genome-wide association studies including 270,000 individuals show sexual dimorphism in genetic loci for anthropometric traits. PLoS Genet. 9(6): e1003500, 2013. Rahman MM, Gopinathan KP. Systemic and in vitro infection process of Bombyx mori nucleopolyhedrovirus. Virus Res. 101(2): 109-118, 2004. Regan JC, Brandão AS, Leitão AB, Mantas Dias AR, Sucena E, Jacinto A, et al. Steroid hormone signaling is essential to regulate innate immune cells and fight bacterial infection in Drosophila. PLoS Pathog. 9(10): e1003720, 2013. Ribas de Pouplana L, Schimmel P. A view into the origin of life: aminoacyl-tRNA synthetases. Cell Mol Life Sci. 57(6): 865-870, 2000. Riddiford LM. Rhodnius, Golden Oil, and Met: A History of Juvenile Hormone Research. Front Cell Dev Biol. 8: 679, 2020. Rus F, Flatt T, Tong M, Aggarwal K, Okuda K, Kleino A, et al. Ecdysone triggered PGRP-LC expression controls Drosophila innate immunity. EMBO J. 32(11): 1626-1638, 2013. Sarkar S, Davies JE, Huang Z, Tunnacliffe A, Rubinsztein DC. Trehalose, a novel mTOR-independent autophagy enhancer, accelerates the clearance of mutant huntingtin and alpha-synuclein. J Biol Chem. 282(8): 5641-5652, 2007. Shinoda T, Itoyama K. Juvenile hormone acid methyltransferase: a key regulatory enzyme for insect metamorphosis. Proc Natl Acad Sci U S A. 100(21): 11986-11991, 2003. Tanaka H, Sagisaka A, Fujita K, Kaneko Y, Imanishi S, Yamakawa M. Lipopolysaccharide elicits expression of immune-related genes in the silkworm, Bombyx mori. Insect Mol Biol. 18(1): 71-75, 2009. Wang L, Xiao Q, Zhou XL, Zhu Y, Dong ZQ, Chen P, et al. Bombyx mori Nuclear Polyhedrosis Virus (BmNPV) Induces Host Cell Autophagy to Benefit Infection. Viruses. 10(1): 14, 2017. Yamamoto K, Ozakiya Y, Uno T. Localization of an Aldo-Keto Reductase (AKR2E4) in the Silkworm Bombyx mori (Lepidoptera: Bombycidae). J Insect Sci. 17(5): 94, 2017. Yang JH, Li G, Qian HY, Zhao GD, Xu AY. Identification and Expression Analysis of Tret1-like in Bombyx mori (in Chinese). Abstracts of Jiangsu Genetic Society 2016 Annual Conference papers, China Conference. 129, 2016. Yang RK, Lei T, Wang Y, Yang SY, Huang CS, Zhang SH, et al. Breeding of new double sex-limited marking silkworm variety Yuke 9 for summer and autumn rearing. Sci Sericul. 45(3): 451-455, 2019. Yuan Y, Zhu F, Xiao R, Ge Q, Tang H, Kong M, et al. Increased expression of Suppressor of cytokine signaling 2 (BmSOCS2) is correlated with suppression of Bombyx mori nucleopolyhedrovirus replication in silkworm larval tissues and cells. J Invertebr Pathol. 174: 107419, 2020. Zhang H, Im SC, Waskell L. Cytochrome b5 increases the rate of product formation by cytochrome P450 2B4 and competes with cytochrome P450 reductase for a binding site on cytochrome P450 2B4. J Biol Chem. 282(41): 29766-29776, 2007. Zemp N, Tavares R, Widmer A. Fungal Infection Induces Sex-Specific Transcriptional Changes and Alters Sexual Dimorphism in the Dioecious Plant Silene latifolia. PLoS Genet. 11(10): e1005536, 2015.