Silkworm, Bombyx mori L. is the model insect belonging to Lepidoptera and an important economic insect utilized for silk production (Goldsmith et al., 2005; Xia et al., 2014). Bivoltine silkworm breeds are known for superior productivity and silk quality, while the silk produced by multivoltine breeds or their crossbreeds is characterized by inferior quality. Silkworm hybrid development programmes focus mainly on exploiting bivoltine and multivoltine breeds with an objective to improve silk quality and productivity besides sustainability. The characteristics which determine the silk quality need to be incorporated from bivoltine component to crossbreeds for improvement of overall silk quality in the hybrid. Functional genomics studies has provided several information on diverse animal and plant species for isolating valuable traits or genes employing novel techniques (Bunnik and Roch 2013). These techniques hold potential to improve the quality of the produce offering directions into future gene mining, genetic manipulation or breeding efforts. Several researchers studied genetic loci controlling variations in silk yield to identify the candidate gene related to silk yield in B. mori; however, it still remains to be elucidated (Zhan et al., 2009; Li et al., 2013; Li  et al., 2017; Fang et al., 2020). 

The grading of silk is influenced by several genetic and/or process parameters including size deviation, winding breaks, evenness variation, neatness, cleanness, tenacity, elongation and cohesion. These are qualitative characteristics, also often influenced by process and environmental factors. Genetic factors determine the silk quality, as silkworm breeds of different voltinism/origin produce silks of varying qualities in the similar conditions despite feeding on same quality mulberry leaves. Till date, molecular markers for any of the characteristics that determine silk quality are not identified.

Significant progress has been made in the field of B. mori genetics and genomics (The International Silkworm Genome Consortium, 2008; Xia et al., 2009; Kawamoto et al., 2019) owing to the understandings from transcriptomic studies that facilitated identification of new exons, novel genes, alternative splicing genes, and trans-splicing events (Li et al., 2012; Shao et al., 2012; Xia et al., 2014; Li et al., 2016). Comparative transcriptome analysis of silkworm strains with different silk yields identified the genes associated with silk yield (Xue et al., 2015; Fang et al., 2015, 2020; Li et al., 2016; Wang et al., 2016; Luan et al., 2018). Further, differentially expressed genes (DEGs) accounting vast variations are mainly involved in the processing and biosynthesis of proteins (Li et al., 2016), silk protein synthesis (Wang et al., 2016; Luan et al., 2018). These studies have added information into the insights of complex genetic basis of silk production, but without any substantial clues on gene(s) regulating silk quality and clearly implicating the complexities in the molecular mechanisms of silk synthesis. 

The transcriptome sequencing of silk glands should provide functional and comprehensive information on gene expression in understanding the role of transcripts in silk fibroin synthesis pathway that influence silk quality. The present investigation was carried out by whole transcriptome analysis of silk glands from 3rd day fifth instar larvae to identify differentially expressed genes (DEGs) for genetic characteristics determining silk quality utilizing silkworm breeds of bivoltine and multivoltine, and its crossbreed using Illumina sequencing technology; to also understand molecular mechanism associated with silk synthesis and quality influencing factors. 

Sample Collection: Bombyx mori silkworm breeds contrasting for silk quality  (Bivoltine: CSR2, CSR27; Multivoltine: PM, Nistari; Crossbreed: PM × CSR2) were  reared on fresh mulberry leaves at the Bivoltine Breeding Laboratory, Central Sericultural Research and Training Institute in Mysuru. The rearing condition was  of a temperature of 25°C, 80% relative humidity, and a light-dark cycle of 14 hours of light and 10 hours of darkness, maintained for three generations."Intact silk glands were dissected on 3rd day of the fifth instar larvae and collected in RNA stabilization solution (Invitrogen) following the manufacturer protocol. The anterior  Silk Gland (ASG) and Posterior Silk Gland (PSG) were dissected out for Nistari, CSR27 and Pure Mysore, CSR2, PMxCSR2 (Kolar gold),  respectively. Silk glands from five larvae were pooled and stored for further use at -80°C.

Note: Total RNA was isolated using TRIzol (Invitrogen) following the manufacturer’s instructions. The quality of RNA preparations was checked on 1% denatured agarose gel and quantified using Nanodrop 8000 Spectrophotometer (Thermo Scientific, USA).  2μg of total RNA was utilized to prepare the libraries using Lexogen SENSE mRNA-Seq Library Prep Kit V2 as per manufacturer’s instructions. The amplified libraries were analysed in Bioanalyzer 2100 (Agilent Technologies) using High Sensitivity (HS) DNA chip as per manufacturer's instructions.

Cluster Generation and Sequencing: After obtaining the Qubit concentration for library and mean peak size, the library was loaded onto Illumina platform for cluster generation and sequencing. The libraries were sequenced using 2 × 150 pair end on Illumina platform.

Mapping and assembly of clean reads: Reference-guided transcript assembly was performed for all five samples; the following steps were taken for mapping and assembling clean reads: Firstly, High-quality (HQ) reads were mapped to the Bombyx mori genome using HISAT2 (Kim et al., 2015) to identify their origin positions. Then, transcript assembly was performed, and transcript abundance was estimated using StringTie (Pertea et al., 2015). The gff compare utility was used to identify potential novel isoforms sharing at least one splice junction with the reference transcript. Lastly, the FPKM (fragments per kilobase of transcript per million mapped reads) values were utilized for the differential expression analysis.

Identification of novel transcripts: To identify novel transcripts not present in the reference GTF file, gff compare utility was run taking the reference GTF and the string-tie merged GTF file. Class code "j" indicates that predicted transcript is a potential novel isoform that shares at least one splice junction with a reference transcript.

Differential Expression of gene analysis: Differential expression of gene analysis was carried out for four different sample sets viz., CSR2 vs PMxCSR2, CSR2 vs PM, CSR27 vs Nistari and PM vs PMxCSR2 and abundances of merged transcripts were estimated using StringTie. The log2 transformed fold change values (log2FC) were used to determine the upregulated and downregulated genes. The criteria used to identify upregulated and down regulated genes include log2FC > 0 for upregulated and log2FC < 0 for downregulated genes. Heatmap for highly significant differential expressed genes (i.e., up and downregulated genes) were generated using the pheatmap package from R software. We analyzed and compared the gene expression levels between the breeds using the method mentioned above. Protein-protein interactions were further analyzed using the STRING database (Szklarczyk et al., 2017).

Pathway analysis: Ortholog assignment and mapping of the transcripts to the biological pathways were performed using KEGG automatic annotation server (KAAS) as described by Moriya et al., (2007).

RT PCR Validation: Among the five key regulators or markers involved in controlling silk quality or silk fibroin synthesis, BGIBMGA008951 - UbiA prenyltransferase domain-containing protein 1 homolog was selected for validation through RT-PCR. Total RNA was isolated from ten silkgland samples using TRIzol reagent from five each of bivoltine (CSR6, CSR16, CSR17, CSR26 and S8) and five multivoltine (BL65, MO6, HB6, FVB1 and L3) breeds on the 3rd day of fifth instar. The quality of isolated RNA was checked on 1% formaldehyde denaturing agarose gel and quantified using Qubit® 2.0 Fluorometer. RT-PCR was carried out on ten samples for the target gene, Ubi-A and Actin-3 gene as a housekeeping gene. Primers used for the reaction include - Ubi-AF: AGACGGGACGAGTTGATTGG and Ubi-AR: ACGGACAGAAGCAACCACAA  and Actin3F: CGGCTACTCGTTCACTACC; Actin3R: CCGTCGGGAAGTTCGTAAG. The mRNA quantity of the UbiA gene was calculated with the CT and ΔΔCq (ΔΔCt) values and normalized to the abundance of the β-Actin gene.

Variant analysis: Variant analysis/SNP calling was carried out using GATK (https://software.broadinstitute.org/gatk) for analysis of variants like SNPs and indels for the five silkworm breeds/hybrid (CSR2, CSR27, PM, Nistari and PM x CSR2). Mapping reads to the B. mori reference genome, and the corresponding GFF file was downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/assembly/GCF_000151625.1/). The mapped reads were sorted, and mark duplicates of mapped reads were achieved by Picard, followed by the split'N'Trim and reassigned mapping qualities. Further, base recalibration was carried out, followed by variant calling and variant filtering.

Transcriptome sequencing and assembly: To understand the molecular mechanism controlling the silk synthesis and its influence on silk quality among the silkworm breeds (Bivoltine: CSR2, CSR27; Multivoltine: PM, Nistari; Crossbreed: PM × CSR2) for contrasting silk quality parameters, the silkgland Transcriptome at the 3rd day fifth instar were compared (Table 1). Transcriptome sequencing using 2 × 150 pair end chemistry carried out on Illumina Hiseq500 platform resulted in 96,778,970 raw reads in total amounting to a total of 13.77Gb data from CSR2, CSR27, PM, Nistari and PM × CSR2 silkgland samples. The high-quality reads were mapped to B. mori reference using HiSAT2. The RNA sequences raw data was deposited to NCBI Short Read Archive (SRA, http://www.ncbi.nlm.nih.gov.sra/) and obtained the accession numbers SRR8816550 (CSR2), SRR8816551 (CSR27), SRR8816552 (PM), SRR8816553 (Nistari) and SRR8816554 (PM × CSR2).

A total of 85,062,305 reads were mapped from a total of 96778970 HQ reads with a coverage of 82.96% for CSR2, 87.76% for CSR27, 90.13% for PM, 91.28% for Nistari and 91.28% for PM × CSR2. Transcript assembly is carried out using String Tie, which assembles transcripts from RNA reads that align to the genome, a total of 46,183 merged transcripts (CSR2: 14,499; CSR27: 12,867; PM: 21,163; Nistari: 17,526; PMxCSR2: 22,156) were obtained. Further, 4,732 novel transcripts and isoforms were identified employing the gff compare utility considering the reference transcripts and StringTie merged files. The novel isoforms were blast X searched against NCBI database. Functional annotation of novel transcripts and isoforms was carried out, which revealed 55 variants showing similarities with Fibroin H and 210 B. mori Sericin-1 like(LOC101740082).

Analysis of differential gene expression: Abundances of the merged transcripts in the silk gland samples were estimated using String Tie and a Python program to extract fpkm information directly from the files generated by StringTie. Differential expression analysis for these silkworm breeds/hybrid indicated upregulation by 102 genes and 123 genes downregulated for a DEG of CSR2 and Pure Mysore; 81 genes upregulated and 173 downregulated genes in a comparison of CSR27 and Nistari. The fold change annotation results indicated that 86 genes were upregulated and 141 genes downregulated for CSR2 and PMxCSR2. 110 genes and 551 genes have annotated exclusively for CSR2 and PMxCSR2, respectively. For Pure Mysore (PM) and the crossbreed (PM x CSR2), 466 genes were differentially expressed among a total of 1127 genes (Table 2). Mapping of the transcript assembly of PM and the crossbreed (PM x CSR2) revealed a total of 200 upregulated genes and 268 downregulated genes. Further, 345 genes and 317 genes have been annotated exclusively for PM and PMxCSR2, respectively, a total of 82 upregulated genes, along with 174 downregulated genes for CSR27 and Nistari. A total of 109 gene expressions exclusive to CSR27 were observed, wherein two newer transcripts, i.e., fibroin heavy chain-like, were identified. A total of 586 genes were found exclusively expressed in Nistari, including transcription elongation factor, transcription initiation factor, etc.

The transcripts of all five samples (CSR2, CSR27, PM, Nistari and PM × CSR2) were compared to understand percent identity and shared transcripts. The comparison between two bivoltine breeds, CSR2 and CSR27 show 7.1% shared genes; while exclusive genes were 31.6% (CSR2) and 61.3% (CSR27). On the other hand, PM and Nistari (multivoltine breeds) shared 18.1% genes, and the exclusive genes hosted by PM and Nistari were 48.5% and 33.4%, respectively.

A total of 17,754 annotated transcripts were mapped to different reference metabolic pathways; 6091 transcripts represented major biomolecules such as carbohydrates, lipids, nucleotides, amino acids, glycans, cofactors, vitamins, etc., followed by genetic information processing (4317), environmental information processing (3150) and cellular processes (4196); and the distribution. Among the metabolism-related pathways, carbohydrates (1125 transcripts) dominated the other pathways, with 842 transcripts for energy, followed by 810 transcripts for amino acids. Further, translation-related pathways (1735 transcripts) dominated the genetic information process, followed by 1460 transcripts for the folding, sorting and degradation process. The majority of the transcripts were represented by signal transduction (2781) in environmental information processing, whereas transport and catabolism (1914) were found to be dominant in the cellular process pathways category, followed by cell growth and death with 1178 transcripts. Based on the KEGG enrichment analysis, oxidative phosphorylation (ko00190) was overrepresented gene ontologies in both the contrasting silkworm breeds and the crossbreed. RNA biogenesis in eukaryotes (ko03008) and mRNA surveillance pathway (ko03015) were enriched under translation in the genetic information and processing category, followed by enrichment of protein processing in endoplasmic reticulum (ko04141), ubiquitin-mediated proteolysis (ko04120) and RNA degradation (ko03018) under folding, sorting and degradation (Fig. 1). 

The transcripts showing differential expression represented interactions in energy metabolism that are vital for silk fibroin production/synthesis, abundance or biogenesis of levels of ribosomal proteins, proteins related to transportation and folding, factors required for development of the posterior silk gland and most importantly repair system  (Fig. 2). Among these, caspases, Ubiquitin, RNA Pol II, RNA Pol III, Mannosidase and acetylglucosamine are accounted for being the most important regulators involved in silk fibroin synthesis (Fig. 3). The differential expression patterns of these key regulators are represented in Table 3. These regulators, especially Caspases, Ubiquitins and RNA polymerases, play a crucial role in silkworm silk gland development and indirectly influence silk fibroin synthesis. Ubiquitin homologs function as molecular chaperons by activating proteasome synthesis and autophagy, which might influence the degradation of abnormal proteins synthesized in the silk gland in the posterior silk gland. Thus, the UbiA prenyltransferase domain containing protein-1 homolog (BGIBMGA008951) was further validated through RT-PCR. The DEGs were evident here with their upregulation by 2.84 times higher expression of transcript BGIBMGA008951in PSG (CSR2 vs PM), whereas there was insignificant expression in the anterior silk gland (CSR27 vs Nistari). Further, the transcript, UBFD1 like, was found to express 3.521 times more in ASG as compared to 0.273 times in PSG. The CT and ΔΔCq(ΔΔCt) values for all the samples for Ubi-A genes fold-change observations were recorded. Results were analyzed using LightCycler® 480 II instrument (Ver. 1.5.0.39) and indicated that expression of UbiA prenyltransferase domain-containing protein 1 homolog was found highest in bivoltine silkworm breeds, S8 followed by CSR26; further, two multivoltine breeds (BL65 and HB6) also represented similar expression (Fig. 4).

Variant analysis: A total of 40,229 variants were identified among the five samples analyzed, with 39,108 SNPs (97.213%) and 1124 indels (2.79%). The highest number of SNPs were recorded in PMxCSR2 (9864), followed by CSR27 (8314), PM (7771), Nistari (7493) and 6787 in CSR2. The high-quality variants were obtained through filtration using vcf tools. Further analysis identified a total of 192 indels in CSR2, followed by 163, 281, 228 and 260 in CSR27, PM, Nistari and PMxCSR2, respectively. The highest number of multiple allelic sites (17) was identified in PMX CSR2-Crossbreed, whereas the other breeds showed low sites. Transition (Ts) variants were comparatively higher as compared to transversions (Tv); transitions were found highest (68.33%), whereas transversions contributed to about 31.72%, and the ratio of Ts/Tv ranged between 1.97 to 2.3 with the highest recorded in CSR27 and lowest in PM. Interestingly, transitions were found highest in PMxCSR2 (6630), followed by CSR27 (5679), PM (4973), Nistari (4854) and CSR2 (4590). The most common substitution being A>G, followed by G>A, C>T and T>C. The higher complexity in variants indicated that silk fibroin synthesis is under the influence of complex regulatory machinery in B. mori. 

Fig. 1. Heat maps of four comparison sets of silkworm breeds (a-d).

Fig. 2. Protein-protein interaction network of the transcripts and their functional role.

Fig. 3. Interaction network showing roles of key regulators involved in silk fibroin synthesis pathway based on the gene – protein interactions in B. mori

Fig. 4. Expression profile of Ubi-A among the silkworm breeds.

Table 1: Silk parameters of silkworm breeds contrasting for silk quality.

Table 2: Sample wise comparison set and statistics of Differential expression of genes.

Table 3: Differential expression pattern of key regulators in silk fibroin synthesis.

✓: Yes : Present (neither UR nor DR);✕: Absent (neither UR nor DR)

DISCUSSION

Highly accurate and precise genome information of B. mori is a basic requisite for comparative genomics in understanding insect diversity, which is essential for basic studies in silkworm breeding (Kawamoto et al., 2019). Multiple transcription factors are involved in the transcriptional regulation of silk protein genes during the development of silkworms (Ma et al., 2011; Xia et al., 2014). This is the first report on the characterization of transcriptomes of the Indian silkworm breeds (bivoltine and multivoltine) contrasting for silk quality, which could further help to understand the complexities involved in silk fibroin synthesis and silk quality. It has been hypothesized that fibroin production in the posterior silk gland (PSG) of Bombyx is directly proportional to silk yield. This is determined by its gland size and protein synthesis, making it possible to improve silk yield through genetic manipulation.

In the present study, high throughput pair-end RNA sequencing was carried out to understand the transcriptomics of silkworm breeds contrasting for silk quality parameters. A substantial number of DEGs and novel transcript isoforms were identified. A total of 467 upregulated and 703 downregulated transcripts were identified based on DEG analysis. Based on the next-generation sequencing results, studies have found that phenotypic changes play a role in genetic divergence (Cheng et al., 2015). The molecular function and cellular role of a protein can be better understood with knowledge of its interactions with other proteins (Celaj et al., 2017). Protein modifications mediated by ubiquitin is essential for several cellular processes, and ubiquitination is one of the crucial processes which regulates the turnover of proteins required for basic cellular processes like cell cycle and cell death (Lee et al.,  2008; Fang et al., 2015). Thus, the ubiquitin-proteasome pathway (UPP) represents an important role in the programmed cell death of the silk gland (Li et al., 2011; Wang et al., 2016; Cui et al.,  2018). The findings of Cui et al. (2018) indicate that the accumulation of abnormal or mutant proteins in the mutant PSG could potentially result in increased activation of proteasome synthesis and autophagy. This leads to rapid degradation of these abnormal proteins, and the silk gland cells use ubiquitin to tag them.

Among the several transcripts, ubiquitinated and RNA biogenesis transcripts were enriched. DEG analysis indicated 28 ubiquitin-related transcripts among 225 genes, which showed differential expression for CSR2 and Pure Mysore and 254 genes for CSR27 and Nistari. Further, based on the novel isoforms/transcripts obtained through the transcriptomic analysis, 258 ubiquitin-related isoforms were observed. Among the silkworm breeds analysed, two ubiquitin-related genes were exclusively expressed in CSR2, followed by two in CSR27, three in Pure Mysore and one in Nistari.  The gene/transcript E3 ubiquitin ligase RNF126 was identified in both the bivoltine breeds (CSR2 and CSR27), which exhibit superior silk quality in comparison with multivoltines (PM and Nistari). The transcript, UbiA prenyltransferase domain-containing protein 1 homolog, was found exclusively expressed in Nistari, but interestingly, it was upregulated 2.84 times in DEG comparison for CSR2 vs PM.  RT-PCR results indicated that the expression of the UbiA prenyltransferase domain-containing protein 1 homolog was found to be highest in bivoltine silkworm breeds, S8, followed by CSR26.

Ubiquitin homologs function as molecular chaperons that lead to the activation of processes like proteasome synthesis and autophagy (Wang and Wang, 2015), promoting degradation of abnormal proteins synthesized in the silk gland; however, its functional complexity needs to be further studied. Among the interactors of fibroin H, fibroin L and P25, Ubi3 is the only representative interactor, which in turn shows interaction with USP46, Ubiquitin carboxyl-terminal hydrolase 46 protein belonging to peptidase C19 family protein. The exact role of USP46 in fibroin synthesis in the silk gland is yet to be understood. 

Another transcript, Mannosyl-oligosaccharide 1,2-alpha-mannosidase (MAN1) was found to be exclusively expressed and upregulated in multivoltine breed, Nistari, based on DEG and protein-protein interactions determined using STRING analysis. MAN1 is involved in glycoprotein quality control targeting misfolded glycoproteins for degradation (Ogen-Stern et al., 2016), and its main function is protein glycosylation, which is part of protein modification. This upregulation was not found in any other comparison sets and might be indicating an important regulator for fibroin synthesis. MAN1 is hypothesized to capture or accumulate all the misfolded proteins in silk glands and processes by interacting with the ubiquitin ligase complex for final proteolytic activity as chaperons, which in turn leads to the activation of caspases eliciting programmed cell death of unfolded proteins in the cells. This is being predicted for possible involvement in silk protein synthesis and the secretion of functional proteins, thereby leading to quality silk protein synthesis. Thus, based on the observations and identification of upregulation of MAN1 in ASG, it is correlative to improper folding of protein affecting the silk quality in multivoltine silkworm breeds. 

Enrichment of the caspase caspase-dependent pathway was reported as a prerequisite for highly efficient biosynthesis and may suppress the secretion of silk proteins (Li et al., 2016). Ribosome biogenesis underlies the cell's capacity to grow because cell growth or an increase in cell mass requires a large number of ribosomes, which are the molecular factors that carry out protein synthesis (Lempiainen and Shore 2009; Li et al., 2016). Increased expression of ribosome proteins reflects rapid biosynthesis in silk glands and is also crucial for the transport of primary proteins (Li et al., 2016; Zhou et al., 2020). The up-regulation of ribosomal-related proteins, proteins involved in degradation and energy metabolism, is essential for silk fibroin synthesis and also for a highly efficient protein synthesis system (Li et al.,  2011; Wang et al., 2014; Song et al., 2016). For RNA enrichment,11 transcripts were upregulated, and 12 RNA-related transcripts were downregulated for the study groups of DEG.  Further, 51 novel transcripts identified were also further subjected to protein-protein interaction analysis, which revealed that RNA polymerase II is one of the important regulators of silk fibroin synthesis. CSR2 showed the expression of 15 transcripts exclusively, whereas 4 were in CSR27, 8 were in Nistari, and 5 were in PM. The transcript-putative RNA polymerase II subunit B1 CTD phosphatase RPAP2 was upregulated 4.34 folds followed by U3 small nucleolar RNA-associated protein 18 homolog with a fold change of 3.98 in comparison for DEG for CSR2 vs PM. Further, transcript RNA-binding protein, FUS-like, showed a fold change of 3.3 with KEGG annotation of mediator of RNA polymerase II transcription subunit 12. This further indicates that the upregulation of RNA pol. II functions as a crucial factor in silk fibroin synthesis in silk glands of B. mori by activation of processes like proteasome synthesis and autophagy promoting degradation of abnormal proteins synthesized in silk gland. The other transcripts observed include D-aspartate oxidase, Arginine N-methyltransferase, Fibroin heavy chain-like, Carbonyl reductase NADPH1, Tyrosine-Protein phosphate, 5-Hydroxytryptamine Receptor, Carboxylesterase 4A-like, N-acetylgalactosaminyltransferase, N-acetylgalactosamine Kinase, Toll-like receptor 7, Sucrose-6-phosphate hydrolase, ATP-dependent DNA helicase Q1, DNA polymerase epsilon subunit 2, Negative elongation factor B and Actin-related protein 8 represented by functional role in cellular activities line ATP synthesis, splicing, ubiquitinization, mRNA synthesis etc.  The role of these in the silk fibroin synthesis pathway and in controlling silk quality needs to be explored further.