
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39251726
71503
10.1038/s41598-024-71503-2
Article
Characterization and comparative analysis of sericin protein 150 in Bombyx mori
Wu Bulah Chia-hsiang 12
Zabelina Valeriya 12
Zurovcova Martina 1
Zurovec Michal zurovec@entu.cas.cz

12
1 grid.447761.7 0000 0004 0396 9503 Biology Centre of the Czech Academy of Sciences, Institute of Entomology, 37005 Ceske Budejovice, Czech Republic
2 grid.14509.39 0000 0001 2166 4904 Faculty of Science, University of South Bohemia, 37005 Ceske Budejovice, Czech Republic
9 9 2024
9 9 2024
2024
14 2099031 10 2023
28 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Lepidopteran silk is a complex mixture of proteins, consisting mainly of fibroins and sericins. Sericins are a small family of highly divergent proteins that serve as adhesives and coatings for silk fibers. So far, five genes encoding sericin proteins have been identified in Bombyx mori. Having previously identified sericin protein 150 (SP150) as a major sericin-like protein in the cocoons of the pyralid moths Galleria mellonella and Ephestia kuehniella, we describe the identification of its homolog in B. mori. Our refined gene model shows that it consists of four exons and a long open reading frame with a conserved motif, CXCXCX, at the C-terminus, reminiscent of the structure observed in a class of mucin proteins. Notably, despite a similar expression pattern, both mRNA and protein levels of B. mori SP150 were significantly lower than those of its pyralid counterpart. We also discuss the synteny of homologous genes on corresponding chromosomes in different moth species and the possible phylogenetic relationships between SP150 and certain mucin-like proteins. Our results improve our understanding of silk structure and the evolutionary relationships between adhesion proteins in the silk of different lepidopteran species.

Keywords

Galleria mellonella
Mucin
SP150
Silk glands
CXCXCX
Synteny
Subject terms

Computational biology and bioinformatics
Evolution
Molecular biology
InterregBYCZ01-039 issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Silk is a secretion product of several arthropod groups, including insects. The best-known silk is produced by Bombyx mori using specialized larval salivary glands, known as silk glands (SG). Silk fibers mainly consist of two types of proteins: fibroins and sericins. Fibroins serve as the central structural proteins forming the fibroin core, and imparts strength and durability to the silk. In most moth species, the fibroin complex consists of three protein subunits: fibroin heavy chain (Fib-H), fibroin light chain (Fib-L), and fibrohexamerin/P25 (Fhx/P25)1–3. These subunits are produced in the posterior SG4. Sericins, in contrast, are a small, highly diverse family of adhesive proteins that bind the fibers together, facilitating the construction of intricate structures such as cocoons or feeding tubes5. Sericin is synthesized in the middle of the SG and envelopes the fibroin core in multiple layers6. The use of sensitive proteomic methods and sequencing of the transcriptomes and genomes of various species have revealed that the sericin protein family is larger than previously recognized.

Two major sericin genes in Bombyx mori have been recognized to produce cocoon sericins: sericin 1 (ser1) and sericin 3 (ser3)7,8. Silkworm mutants carrying a truncated ser1 gene are unable to spin or produce coarse cocoons, suggesting that ser1 is crucial for reducing friction during spinning9. Additionally, sericin 2 (ser2) along with two recently identified genes, sericin 4 (ser4) and sericin 5 (ser5), produce proteins found in non-cocoon silk spun by younger larvae (including early last-instar larvae)10–13. The presence of Ser2, Ser4, and Ser5 in these non-cocoon silks suggests that they may have specific functional roles during early developmental stages.

Sericin genes exhibit high sequence divergence, convergent evolution, duplications and deletions, resulting in species-specific sequences that often have no obvious orthologs in other species. Previous studies have shown that several lepidopteran species contain type 1 sericins characterized by a conserved stretch with two cysteine residues at their C-terminus14. In our previous research, we identified sericin protein 150 (SP150) as a major sericin-like protein in the cocoons of the pyralid moths Galleria mellonella and Ephestia kuehniella 14,15. However, a similar protein has not yet been described in B. mori. To investigate the presence of a SP150 ortholog in B. mori, we used advanced omics methods. Our study will help to elucidate the protein structure of silk in B. mori and to understand the complex evolutionary history of the genes encoding these proteins.

Materials and methods

Silkworm strains and datasets used

A non-diapausing B. mori strain, w1-pnd (white egg 1, non-pigmented and non-diapausing egg), was used in the experiments as a wild type (wt) strain. Larvae were reared on mulberry leaves at 25 °C. A list of RNA-seq datasets from NCBI Sequence Read Archive (SRA) is provided in Supplementary Table S1.

RNA extraction, RNA-seq analysis

Total RNA was extracted from the SGs and control tissues (intestine, integument, fat body, head, testes, and ovary) of 3–5-day-old fifth-instar larvae using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. For qPCR: the first cDNA strand was synthesized using 0.5 μg of total RNA as a template with RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA). The cDNA product was used to verify the last exon junction of B. mori SP150 and its expression level in different tissues. Primers were designed using the Geneious Prime software platform (Biomatters, Auckland, New Zealand; version 2021.2.2) and are listed in Supplementary Table S2.

For northern blotting, the 5 µg RNA aliquots were subjected to agarose gel electrophoresis, transferred to a nylon membrane (Hybond N+, Sigma-Aldrich), and hybridized as previously described15. Northern blotting probes were synthesized by RT-PCR using the primers listed in Supplementary Table S2 and labeled with α-32P[dATP] by random priming with a DecaLabel DNA Labeling Kit (Thermo Scientific, Waltham, MA, USA). Hybridization signals were detected by autoradiography using the storage phosphor screen of a STORM 860 phosphor imager (Molecular Dynamics, Ramsay, MN. USA).

The preparation of the (ASG, MSG, and PSG) cDNA libraries was performed as described previously16. A MiSeq (Illumina, San Diego, CA, United States) instrument was used to obtain 150-nt long paired-end reads. Raw RNA-seq reads from B. mori silk glands were trimmed using Trimmomatic (version 0.32)17 and aligned to the reference genome with the software STAR (version 2.7.10b)18. The splice junctions were visualized using the Sashimi plot function in IGV19. Transcript quantification was performed with kallisto (version 0.48.0)20.

qPCR analysis

The qPCR was performed using HOT FIREPol EvaGreen qPCR Mix Plus (Solis BioDyne, Tartu, Estonia). The PCR reaction volume of 20 µl contained 5 µl diluted cDNA and 250 nM primer. Amplification was performed using a Rotor-Gene Q MDx 2plex HRM (Qiagen, Hilden, Germany) for 45 cycles (95 °C for 15 s; annealing temperature matched to the primer pair for 30 s; 72 °C for 20 s) after an initial denaturation step (95 °C for 15 min). Each sample was analyzed in triplicate. Results were analyzed using Rotor Gene Q software (version 2.3.5). Elongation factor 1 alpha (EF1a, NM _001044045.1) was used as a reference gene, and the relative expression of target genes was calculated using the Pfaffl method21. Statistical analysis was performed using Student’s t-test in R (version 4.1.1); p values < 0.05 were considered statistically significant. DA detailed statistical analysis is provided in Supplementary Table S3.

Proteomic analysis and data mining

Protein analysis of B. mori cocoons and database searches were performed at the Proteomics Core Facility (BIOCEV, Vestec, Czech Republic) as previously described22. Approximately 10 mg of the silk cocoon was boiled in 8 M urea, and samples were further processed using solid-phase enhanced sample preparation technology (SP3 beads)23. Samples were then digested with trypsin, and the resulting peptides obtained were subjected to liquid chromatography—MS. Four wt cocoons were analyzed in parallel. In addition, raw proteomic data deposited in public databases24 on the composition of individual silk layers in the B. mori cocoon were reanalyzed. The obtained MS/MS spectra were matched against the protein sequences of NCBI B. mori Annotation Release 103, which was enriched for the newly discovered SP150. Quantification was performed using label-free algorithms, and data were analyzed using MaxQuant and Perseus v.1.5.2.425,26.

We also dissolved cocoon silk in LiSCN as follows. The cocoon proteins were dissolved in an approximately 5% protein solution in saturated LiSCN containing 2% 2-mercaptoethanol. The data was analyzed using Data-independent acquisition (DIA) software27.

Chromosomal localization and collinearity analysis

The genome assemblies and annotated information of B. mori (GCF_014905235.1), E. kuehniella14,28, and G. mellonella (GCF_026898425.1) were processed and submitted to the GENESPACE software29 for syntenic analysis. Plots showing the microsyntenic relationships were then generated based on the best mutual hits between the three species and visualized using the R package ggplot230.

Phylogenetic analysis

Codon-based alignment of the SP150 and three-cysteine (CXCXCX) mucin 3′ ends was performed using MEGA7 software following the MUSCLE method31. The phylogram was generated using the IQ-TREE server32,33, which included both the selection of the best substitution model by ModelFinder34 and tree inference using MLE (ultrafast bootstrap, 1000 replicates).

Results

Identification of SP150 gene in B. mori

To identify the homolog of the major sericin gene P150 described previously in G. mellonella and E. kuehniella, we performed a BLAST search against the B. mori genome. We identified two adjacent homologous regions in the genomic sequence that were predicted to be parts of the two B. mori genes. Most homologous sequences belonged to LOC101737213, and the remaining C-terminus was predicted to belong to a separate gene, LOC119629229. Remarkably, the B. mori putative protein sequence shared 47.5% identity with G. mellonella SP150 across 70 C-terminal amino acid residues. We hypothesized that the predicted gene models was incorrect, and that the sequences of both B. mori genes were part of a single large SP150 gene.

To test our hypothesis, we prepared cDNA libraries from the anterior, middle and posterior silk glands (ASG, MSG, and PSG) of last-instar larvae and aligned the silk gland-specific RNA-seq data with the reference genome. Alignment revealed the absence of an intergenic region between LOC101737213 and LOC119629229 (Fig. 1A). To verify that the two putative B. mori genes constituted a single gene, we designed intron-spanning RT-PCR primers to link the last two exons of LOC101737213 to the second exon of LOC119629229 (Fig. 1B). As shown in Fig. 1C, the amplified cDNA fragments supported our hypothesis of a single SP150 gene.Fig. 1 The revised gene model for the B. mori SP150 gene consisting of two putative genes, LOC101737213 and LOC119629229. (A) The RNA-seq reads from B. mori middle silk glands were mapped to the genomic region of Chromosome 12 (NC_051369.1: 1,270,000–1,280,000, strand flipped). The number of reads bridging exon junctions are indicated at the midpoint of the arcs. The absence of bridging reads shown between LOC101737213 and LOC119629229 indicates that the previously identified intergenic region is incorrect. (B) Predicted exon–intron structure of the revised SP150 gene model (green), consisting of four exons and three introns. Gene models from current NCBI annotation are shown for comparison: LOC101737213 (pink); LOC119629229 (blue); MAD (Mothers against dpp; gray); and the KWMTBOMO06993 gene model from SilkBase (brown). The bottom figure shows an approximately 2.2-kb region with chromosomal coordinates and the positions of the primer pairs used in this study. (C) Validation of the 3′ region of the revised SP150 gene model by RT-PCR and agarose gel electrophoresis. The expected product sizes specific for primer pairs: #1F-#1R, #2F-#2R, and #3F-#3R were 1003 bp, 313 bp, and 93 bp, respectively. The electrophoretogram contains a 1 kb ladder (lane 1) and 100 bp ladder (lane 5).

The resulting gene model of B. mori SP150, shown in Fig. 1A, spans approximately 20 kb and comprises four exons and three introns (Fig. 1A). The first two exons encode a signal peptide and are part of a short, non-repetitive N-terminal sequence. The third exon is notably large and comprises 94% of the open reading frame (ORF), featuring two central repetitive regions flanked by unique sequences. The last exon contains a short ORF that ends with a stop codon. Altogether, the gene encodes a protein of 4552 amino acids, including a 19-amino acid signal peptide.

Putative SP150 protein

The predicted protein product of the SP150 gene in B. mori is a large protein weighing 467 kDa and comprising 4552 amino acid residues. It begins with a 19-amino acid signal peptide, followed by a 616-amino acid non-repetitive central segment. This is followed by 45 repeats of a 30-amino acid motif (repeat 1), a 34 amino acid non-repetitive linker, and 73 repeats of a 35-amino acid motif (repeat 2). The protein ends with a non-repetitive C-terminus spanning 388 amino acids. The complete amino acid sequence is listed in Supplementary Text S1. As shown in Fig. 2, SP150 consisted of two types of highly conserved threonine-rich repeat blocks (Supplementary Table S4). The SP150 protein is relatively highly hydrophilic (hydropathy index = − 0.672 compared to − 1.118 of Ser1). The B. mori SP150 contains more than 27% threonine, 14% serine, and 12% alanine residues. The C-terminus (encoded by the last exon) contains a short, conserved three-cysteine motif, CXCXCX (Supplementary Table S5).Fig. 2 Amino acid sequence logos showing conservation pattern in two repeat types of the B. mori SP150 protein. The hydrophobicity of amino acids is indicated by color: hydrophilic (blue; RKDENQ); neutral (dark gray; SGHTAP); hydrophobic (orange; YVMCLFIW). The height of each letter indicates the degree of conservation at that position.

Compared with the SP150 proteins of G. mellonella and E. kuehniella, the SP150 of B. mori was almost three times larger, less hydrophilic, and contained fewer serine residues. Except for the C-terminal amino acids, there were minimal similarities among the SP150 proteins of these species (Supplementary Table S6).

SP150 mRNA is specifically expressed in MSG

To determine the specific expression of B. mori SP150 in the silk glands, we used RNA-seq data from three SG-specific cDNA libraries and quantified transcript abundance using the kallisto software (see “Materials and methods”). As shown in Table 1, the highest transcript abundance of SP150 was found in MSG, but this level remained substantially lower than that of sericins 1–3. In addition, we reanalyzed publicly available RNA-seq data for silk gene expression from previous experiments24,35 using a similar approach to calculate the transcript levels. Consistent with our results, publicly available data also indicated that the highest TPM (transcripts per million) for SP150 was detected in the MSG in its anterior region (Table 1). Table 1 Comparison of RNA-seq data from three studies. Transcript quantification of silk genes in anterior (ASG), middle (MSG), and posterior (PSG) silk glands was computed by the software kallisto in transcript per million (TPM) units. A-MSG denotes anterior MSG, M-MSG denotes middle part of MSG, P-MSG denotes rear part of MSG, "–" denotes "no detection". Sources of the RNA-seq data include this study, BioProject PRJNA559726, and BioProject PRJDB8614 (see Supplementary Table S1 for details).

This study	PRJNA559726	PRJDB8614	
Gene	ASG	MSG	PSG	ASG	MSG	PSG	ASG	A-MSG	M-MSG	P-MSG	PSG	
SP150	0.017	48.024	0.049	0.019	0.599	0.002	0.008	8.632	2.929	0.362	0.002	
MAD	3.263	2.533	2.548	1.764	1.116	1.011	0.936	1.200	1.639	1.482	0.463	
Ser1	11.945	46,294.400	405.991	1.989	5018.353	225.722	4.819	2.249	24,919.433	22,886.967	42.396	
Ser2	167.785	11,006.100	1.232	276.775	3953.410	0.224	0.640	50,552.200	598.646	3.409	1.730	
Ser3	1.155	55,512.700	1.768	0.065	205.281	0.028	0.010	23.712	70.735	0.103	0.029	
Ser4	0.426	2.282	0.347	0.084	0.271	0.122	0.024	1.649	0.392	1.521	0.155	
Ser5	2.258	0.441	–	1.257	0.236	0.004	1.091	1.441	0.265	0.243	0.036	
Muc-12	16.427	1021.620	0.312	0.889	25.841	0.110	0.235	471.462	43.508	2.230	0.124	
FibH	13.409	1161.180	20,647.500	3.484	285.863	6692.817	0.843	0.923	0.658	51.638	2688.293	

To further confirm the tissue specificity of SP150 expression, we isolated mRNAs from different parts of the silk glands and control tissues (intestine, integument, fat body, head, testis, and ovary) from day 3–5 last-instar larvae. We then prepared cDNA and performed qPCR. As shown in Fig. 3, we also analyzed the expression of MAD, a gene adjacent to SP150, and three genes encoding sericins, ser1, ser2 and ser3. The results showed that B. mori SP150 was specifically expressed in the middle silk glands, whereas MAD was ubiquitously expressed.Fig. 3 Quantitative PCR (qPCR) analysis of silk gland-specific gene expression in nine larval tissues of B. mori. The expression levels of the genes ser1, ser2, ser3 and SP150 were measured as well as the MAD gene, which is located adjacent to SP150, for comparison. Statistical differences were evaluated using Student’s t-test (see Supplementary Table S3). The error bars indicate the standard deviation. The results show that B. mori SP150 is specifically expressed in middle silk glands, whereas MAD is ubiquitously expressed.

Finally, we isolated RNA from several tissues of the last-instar larvae (wandering stage) and performed northern blotting. As shown in Fig. 4, the amount of SP150 transcript was very low, close to the detection limit. The SP150 RNA transcript was notably large in length and primarily localized in the middle parts of the SG, similar to the patterns observed for ser1.Fig. 4 Northern blot analysis of the tissue-specific expression of SP150. Total RNA samples from various tissues were analyzed with 32P-labeled cDNA fragments specific to SP150. A Ser1 probe was used as a control. FB fat body, INT integument, GUT intestine, PSG posterior SG, MS2 posterior part of middle SG, MSG1 anterior part of middle SG, ASG anterior SG.

Quantitative proteomic analysis of silk samples

To investigate the presence of SP150 in B. mori cocoons, we performed mass spectrometry (MS) proteomic analyses of wild-type cocoon silk. Data were analyzed using the Andromeda search engine integrated into the MaxQuant software25,26, and the relative protein abundances were determined through label-free quantification. We identified 118 proteins, with a false discovery rate (FDR) of 1% for protein identification. The consistency of protein intensities between biological replicates was robust, as summarized in Supplementary Table 7.

To evaluate the effect of solvent on peptide identification, we dissolved the cocoons in 5% protein solution in saturated LiSCN and compared the peptide yields. As shown in Table 2, overall, LiSCN yielded slightly more peptides. However, the number of specific peptides attributed to SP150 remained low. Our proteomic analysis revealed that SP150, similar to Ser2 and Muc-12, was detected in cocoon silk at a low level, close to the detection limit of our instruments. In contrast (Fig. 4A), Ser1 and Ser3 were identified as the most abundant components of cocoon silk, with concentrations at least five orders of magnitude higher than SP150 (Fig. 5A). Table 2 Comparison of the results of proteomic analysis of cocoon samples from B. mori. Counts of detected peptides from silk proteins in cocoon samples dissolved in 8 M Urea and LiSCN. "–" indicates "no detection".

Protein	8 M Urea	LiSCN	
SP150	4	3	
Mad	–	–	
Ser1	88	71	
Ser2	4	14	
Ser3	12	5	
Ser4	33	82	
Ser5	–	–	
Muc-12	4	33	
FibH	10	14	

Fig. 5 Proteomic analysis of B. mori silk proteins (A) Analysis of wt cocoons as previously described22; and (B) Data from individual cocoon layers obtained from a public repository24. Label-free quantification (LFQ) of silk proteins from cocoons was calculated using MaxQuant. LFQ intensities were log2-transformed. Relative protein contents in cocoon silk were analyzed using MaxQuant/Andromeda (eight experiments). Error bars indicate the standard deviation. The proteomic analysis confirmed that Ser1 and Ser3 were the most abundant silk components. The other proteins, including SP150, Ser2, and Muc-12, were detected in cocoon silk at low levels.

To further confirm protein abundance in the cocoons, we reanalyzed existing proteomics data from the public repository24 using our new SP150 annotation. The abundance of SP150, Muc-12, and Ser2 in the cocoons is shown in Fig. 5B. All three proteins were present at very low levels, with SP150 and Ser2 showing highest levels in the innermost cocoon layer (layer 1). Overall, these results confirm that SP150 is present at low levels in cocoons, comparable to the levels observed for Ser2 and Muc-12.

Synteny in regions coding for SP150 genes across Lepidoptera

A previous study on the pyralid moths G. mellonella and E. kuehniella showed that all known sericin genes, except SP150, are located within a cluster of orthologous genes in the same chromosomal region14. In addition, the results revealed some local rearrangements and duplications in this region, including the increase in copy number of several sericin genes in G. mellonella compared to related moth species14.

Our data showed that similar microsynteny was also observed for SP150 between B. mori and G. mellonella or E. kuehniella (Fig. 6). The SP150 gene is located on a different chromosome than other sericin genes in a conserved region between genes encoding metalloprotease 1 and croquemont 1. As shown in Fig. 6, the region on chromosome 12 of B. mori has well-conserved synteny, comprising more than 40 genes, except for an inversion that positions the SP150 region in the reverse orientation relative to the adjacent genes.Fig. 6 Microsynteny maps of SP150 and their flanking genes across species. Horizontal color blocks indicate chromosomal segments in each species. Homologous genes and gene orientation are represented by left- and right-pointing triangles with homologous pairs connected by lines. SP150 homologs are highlighted in red.

SP150 may be related to Muc-12

To investigate the evolution of SP150, we performed a BLAST search for homologous sequences in insect genomes, using the conserved sequence encoding the C-terminal protein end as a query. We found no obvious orthologs in non-lepidopteran insects, suggesting that SP150 is a Lepidoptera-specific gene. Furthermore, there were no SP150 orthologs in members of the superfamily Papilionoidea.

SP150 proteins are highly divergent, making it difficult to align homologous proteins from different lepidopteran families, with the exception of the conserved C-terminus. The most prominent conserved motif is the three-cysteine sequence (CXCXCX), which is located 12–29 amino acids away from the C-terminus.

Interestingly, another silk gland-specific protein also contains three cysteine (CXCXCX) motif, namely mucin-12, which also resembles SP150 due to its size and repetitive structure. To better understand the relationship of both proteins, we constructed a dendrogram using the C-termini of SP150 and Muc-12 from representatives of different lepidopteran families (Fig. 7). The resulting phylogenetic tree was robust and clearly distinguished the SP150 and Muc-12 clades, although the sequences from most primitive species showed less clear separation (Fig. 6). The sequence alignment and consensus sequences are shown in Fig. 6B.Fig. 7 Relationship between lepidopteran SP150 and Muc-12 sequences. (A) Maximum likelihood phylogenetic tree based on the alignment of the C-terminal amino acid sequences of SP150 and Muc-12 homologs from selected lepidopteran species. The Nesw_1 transcript from N. swammerdamellus, (Incurvarioidea), the most primitive in this group, was selected for tree rooting. See Supplementary Table S8 for sequence details. (B) Alignment of the C-terminal regions of SP150 proteins from representative lepidopteran species. Sequences include the characteristic CXCXCX region, which is well conserved between species.

Discussion

Here, we report the discovery of a new sericin-like gene, SP150, in the genome of B. mori. Our results also show that the region on chromosome 12, where this gene is located was misannotated in previous models, which had misidentified the homologous sequence as two separate putative B. mori genes. Here, we present the corrected SP150 gene model and show that its ORF encodes a large protein with a repetitive structure that is specifically expressed in MSG.

Predicting the structures of large genes remains a challenge36. The best results were obtained by comparing genomic and cDNA sequences or by using proteomic data37. Using this approach, SP150 sequences have been successfully identified in G. mellonella and E. kuehniella14,15. However, identification of the SP150 homolog in B. mori is difficult because of its relatively low expression. We identified SP150 based on the conservation of a short conserved motif at the C-terminus. Additionally, the identification of SP150 was supported by microsynteny between the corresponding genomic regions in B. mori and the two pyralid moth species (Fig. 6).

The P150 protein of B. mori, with a mass exceeding 467 kDa, stands out from its putative homologs in other lepidopteran species because of its considerable size. For instance, the size of other SP150 homologs ranges from 110.4 kDa in Pectinophora gossypiella (family Gelechiidae, LOC126368327) to 196 kDa in Manduca sexta (family Sphingidae, XP_037298329), and 307.8 kDa in Mythimna loreyi (family Noctuidae, KAJ8704851). Such size differences are not so unusual, as it has been shown that genes coding for silk proteins with repetitive sequences often expand or shrink12,38. For example, the size of fibroin genes varies from 216.1 kDa in Antheraea pernyi (Lepidoptera) to 998.8 kDa in Limnephilus lunatus (Trichoptera)38,39.

Unlike spiders, which have several fibroin genes, moths possess only one fibroin gene that functions as the primary structural component of silk. However, silk composition varies among moth species, with sericins being particularly important components for building three-dimensional silk structures, contributing to its strength and toughness40. For example, G. mellonella builds dense feeding tubes and cocoons required to protect the larvae in the hive. In G. mellonella, sericins and other soluble silk components constitute approximately 48% of the cocoon mass, whereas in B. mori cocoons, they make up only about 26% of the mass. Additionally, some moths, including those of Samia ricini, contain as little as 16% soluble proteins12,41.

The large divergence among sericin proteins suggests that they evolved to perform slightly different functions depending on the specific needs of the silk-producing organism. The three-cysteine motif (CXCXCX) encoded by the C-termini of some SG-specific genes is reminiscent of the “cysteine knot” motif described in some mammalian growth factors, including the VEGF family42. Such motifs contribute to the structural integrity of proteins.

The similarity between SP150 and Muc-12 suggests that they share a common origin. Both are large, highly divergent proteins with repetitive sequences encoding ORFs composed of simple amino acids, and both likely function as silk adhesives15. Our phylogram (Fig. 7) clearly separates the two clades SP150 and Muc-12 with good support, although the placement of the sequences of the most primitive species, which contained only one such gene, was difficult. This suggests that the possible duplication and divergence between the P150 and mucin lineages may have occurred later in the evolutionary history of Lepidoptera. The question of whether SP150 and Muc-12 share a common ancestor and diverge widely, or whether they are the result of convergent evolution, remains an important question for future research.

Sericin P150 is one of the most abundant silk proteins in G. mellonella and E. kuehniella15. In contrast, our results showed that it is only present in trace amounts in the silk of B. mori. The SP150 protein in B. mori is mainly located in the inner cocoon layer at the end of the last larval instar and is also present in the non-cocoon silk produced by larvae from earlier instars, as suggested by the reexamined data of Dong et al.35. Non-cocoon silk contains Ser2, Ser4, and Ser5 proteins10–12, which play a role in the initial stages of silk spinning, and are needed for securing the molting larva and anchoring the cocoon to a suitable substrate12. Thus, SP150 is a minor component of B. mori silk that is present in both cocoon and non-cocoon silk.

The low abundance of SP150 protein in B. mori contrasts with its high abundance in pyralid moths. This difference may be partly attributed to the unusually high levels of sericin proteins found in pyralid moths compared to those in B. mori15. It has also been hypothesized that gene expression inversely correlates with gene length, as the process of transcription is time-dependent43,44. Thus, the production of extremely large proteins such as SP150, may be repressed in B. mori. In addition, the sericin sequences in B. mori cocoons are relatively simple in structure with a relatively high serine content (more than 38%, Supplementary Table S5), facilitating their digestion by cocoonase. Sericins with lower serine contents (less than 15%) are therefore more likely limited to non-cocoon silk or present in low amounts45. This is particularly important in B. mori as failure to open the cocoon is lethal to adults, as evident in mutants lacking cocoonase45. Detailed functional studies in different lepidopteran species are required to determine the cause of the observed differences in the expression of SP150.

Overall, our results contribute to a deeper understanding of silk structure and the molecular mechanisms regulating silk production and provide a basis for future studies investigating the evolution of silk proteins in Lepidoptera.

Supplementary Information

Supplementary Information.

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-71503-2.

Acknowledgements

This research was supported by the European Community’s Program Interreg Bayern—Tschechien BYCZ01-039. This publication is also supported by the project “BIOCEV—Biotechnology and Biomedicine Centre of the Academy of Sciences and Charles University” (CZ.1.05/1.1.00/02.0109), from the European Regional Development Fund. Computational resources were provided by the e-INFRA CZ project (ID:90254), which is supported by the Ministry of Education, Youth and Sports of the Czech Republic and by the ELIXIR-CZ project (ID:90255), part of the international ELIXIR infrastructure. We thank Dr. Miluše Hradilova for her help with RNA sequencing.

Author contributions

Mi.Z.: supervision writing; B.C.W.: investigation, data analysis; writing; Ma.Z.: phylogenetic analysis; V.Z.: B. mori rearing and staging.

Data availability

Data will be made available on request at zurovec@entu.cas.cz or bulah@entu.cas.cz.

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Tanaka K Mori K Mizuno S Immunological identification of the major disulfide-linked light component of silk fibroin J. Biochem.-Tokyo 1993 114 1 4 10.1093/oxfordjournals.jbchem.a124122 8407860
Tanaka, K., Mori, K. & Mizuno, S. Immunological identification of the major disulfide-linked light component of silk fibroin. J. Biochem.-Tokyo 114, 1–4 (1993).8407860 10.1093/oxfordjournals.jbchem.a124122
2. Zurovec M Kodrik D Yang C Sehnal F Scheller K The P25 component of Galleria silk Mol. Gen. Genet. 1998 257 264 270 10.1007/s004380050647 9520260
Zurovec, M., Kodrik, D., Yang, C., Sehnal, F. & Scheller, K. The P25 component of Galleria silk. Mol. Gen. Genet. 257, 264–270 (1998).9520260 10.1007/s004380050647
3. Zurovec M Vaskova M Kodrik D Sehnal F Kumaran AK Light-chain fibroin of Galleria mellonella L Mol. Gen. Genet. 1995 247 1 6 10.1007/BF00425815 7715595
Zurovec, M., Vaskova, M., Kodrik, D., Sehnal, F. & Kumaran, A. K. Light-chain fibroin of Galleria mellonella L. Mol. Gen. Genet. 247, 1–6 (1995).7715595 10.1007/BF00425815
4. Lucas, F. & Rudall, K. M. in Comprehensive Biochemistry, vol. 26 (eds Florkin, M. & Stotz, E. H.)s 475–558 (Elsevier, 1968).
5. Sehnal F Sutherland T Silks produced by insect labial glands Prion 2008 2 145 153 10.4161/pri.2.4.7489 19221523
Sehnal, F. & Sutherland, T. Silks produced by insect labial glands. Prion 2, 145–153. 10.4161/pri.2.4.7489 (2008).19221523 10.4161/pri.2.4.7489
6. Shibukawa A Studies on the silk substance within the silk gland in the silkworm, Bombyx mori L Bull. Sericult. Exp. Sta. 1959 15 401
Shibukawa, A. Studies on the silk substance within the silk gland in the silkworm, Bombyx mori L. Bull. Sericult. Exp. Sta. 15, 401 (1959).
7. Okamoto H Ishikawa E Suzuki Y Structural analysis of sericin genes. Homologies with fibroin gene in the 5' flanking nucleotide sequences J. Biol. Chem. 1982 257 15192 15199 10.1016/S0021-9258(18)33412-4 6294094
Okamoto, H., Ishikawa, E. & Suzuki, Y. Structural analysis of sericin genes. Homologies with fibroin gene in the 5’ flanking nucleotide sequences. J. Biol. Chem. 257, 15192–15199 (1982).6294094 10.1016/S0021-9258(18)33412-4
8. Takasu Y Identification and characterization of a novel sericin gene expressed in the anterior middle silk gland of the silkworm Bombyx mori Insect. Biochem. Mol. Biol. 2007 37 1234 1240 10.1016/j.ibmb.2007.07.009 17916509
Takasu, Y. et al. Identification and characterization of a novel sericin gene expressed in the anterior middle silk gland of the silkworm Bombyx mori. Insect. Biochem. Mol. Biol. 37, 1234–1240. 10.1016/j.ibmb.2007.07.009 (2007).17916509 10.1016/j.ibmb.2007.07.009
9. Takasu Y Iizuka T Zhang Q Sezutsu H Modified cocoon sericin proteins produced by truncated Bombyx Ser1 gene J. Silk Sci. Technol. Jpn. 2017 25 35 47 10.11417/silk.25.35
Takasu, Y., Iizuka, T., Zhang, Q. & Sezutsu, H. Modified cocoon sericin proteins produced by truncated Bombyx Ser1 gene. J. Silk Sci. Technol. Jpn. 25, 35–47. 10.11417/silk.25.35 (2017).10.11417/silk.25.35
10. Dong ZM Identification of Bombyx mori sericin 4 protein as a new biological adhesive Int. J. Biol. Macromol. 2019 132 1121 1130 10.1016/j.ijbiomac.2019.03.166 30928374
Dong, Z. M. et al. Identification of Bombyx mori sericin 4 protein as a new biological adhesive. Int. J. Biol. Macromol. 132, 1121–1130. 10.1016/j.ijbiomac.2019.03.166 (2019).30928374 10.1016/j.ijbiomac.2019.03.166
11. Guo K Identification and characterization of sericin5 reveals non-cocoon silk sericin components with high beta-sheet content and adhesive strength Acta Biomater. 2022 150 96 110 10.1016/j.actbio.2022.07.021 35902035
Guo, K. et al. Identification and characterization of sericin5 reveals non-cocoon silk sericin components with high beta-sheet content and adhesive strength. Acta Biomater. 150, 96–110. 10.1016/j.actbio.2022.07.021 (2022).35902035 10.1016/j.actbio.2022.07.021
12. Kludkiewicz B Structure and expression of the silk adhesive protein Ser2 in Bombyx mori Insect Biochem. Mol. 2009 39 938 946 10.1016/j.ibmb.2009.11.005
Kludkiewicz, B. et al. Structure and expression of the silk adhesive protein Ser2 in Bombyx mori. Insect Biochem. Mol. 39, 938–946. 10.1016/j.ibmb.2009.11.005 (2009).10.1016/j.ibmb.2009.11.005
13. Michaille JJ Garel A Prudhomme JC Cloning and characterization of the highly polymorphic Ser2 gene of Bombyx-Mori Gene 1990 86 177 184 10.1016/0378-1119(90)90277-X 2323571
Michaille, J. J., Garel, A. & Prudhomme, J. C. Cloning and characterization of the highly polymorphic Ser2 gene of Bombyx-Mori. Gene 86, 177–184. 10.1016/0378-1119(90)90277-X (1990).2323571 10.1016/0378-1119(90)90277-X
14. Wu BC Characterization of silk genes in Ephestia kuehniella and Galleria mellonella revealed duplication of sericin genes and highly divergent sequences encoding fibroin heavy chains Front. Mol. Biosci. 2022 9 1023381 10.3389/fmolb.2022.1023381 36523651
Wu, B. C. et al. Characterization of silk genes in Ephestia kuehniella and Galleria mellonella revealed duplication of sericin genes and highly divergent sequences encoding fibroin heavy chains. Front. Mol. Biosci. 9, 1023381. 10.3389/fmolb.2022.1023381 (2022).36523651 10.3389/fmolb.2022.1023381
15. Kludkiewicz B The expansion of genes encoding soluble silk components in the greater wax moth, Galleria mellonella Insect Biochem. Mol. Biol. 2019 106 28 38 10.1016/j.ibmb.2018.11.003 30448349
Kludkiewicz, B. et al. The expansion of genes encoding soluble silk components in the greater wax moth, Galleria mellonella. Insect Biochem. Mol. Biol. 106, 28–38. 10.1016/j.ibmb.2018.11.003 (2019).30448349 10.1016/j.ibmb.2018.11.003
16. Rouhova L Silk of the common clothes moth, Tineola bisselliella, a cosmopolitan pest belonging to the basal ditrysian moth line Insect Biochem. Mol. 2021 10.1016/j.ibmb.2021.103527
Rouhova, L. et al. Silk of the common clothes moth, Tineola bisselliella, a cosmopolitan pest belonging to the basal ditrysian moth line. Insect Biochem. Mol.10.1016/j.ibmb.2021.103527 (2021).10.1016/j.ibmb.2021.103527
17. Bolger AM Lohse M Usadel B Trimmomatic: A flexible trimmer for Illumina sequence data Bioinformatics 2014 30 2114 2120 10.1093/bioinformatics/btu170 24695404
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. 10.1093/bioinformatics/btu170 (2014).24695404 10.1093/bioinformatics/btu170
18. Dobin A STAR: Ultrafast universal RNA-seq aligner Bioinformatics 2013 29 15 21 10.1093/bioinformatics/bts635 23104886
Dobin, A. et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. 10.1093/bioinformatics/bts635 (2013).23104886 10.1093/bioinformatics/bts635
19. Robinson JT Integrative genomics viewer Nat. Biotechnol. 2011 29 24 26 10.1038/nbt.1754 21221095
Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24–26. 10.1038/nbt.1754 (2011).21221095 10.1038/nbt.1754
20. Bray NL Pimentel H Melsted P Pachter L Near-optimal probabilistic RNA-seq quantification Nat. Biotechnol. 2016 34 525 527 10.1038/nbt.3519 27043002
Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34, 525–527. 10.1038/nbt.3519 (2016).27043002 10.1038/nbt.3519
21. Pfaffl MW A new mathematical model for relative quantification in real-time RT-PCR Nucleic Acids Res. 2001 29 e45 10.1093/nar/29.9.e45 11328886
Pfaffl, M. W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 29, e45. 10.1093/nar/29.9.e45 (2001).11328886 10.1093/nar/29.9.e45
22. Zabelina V Mutation in Bombyx mori fibrohexamerin (P25) gene causes reorganization of rough endoplasmic reticulum in posterior silk gland cells and alters morphology of fibroin secretory globules in the silk gland lumen Insect. Biochem. Mol. Biol. 2021 135 103607 10.1016/j.ibmb.2021.103607 34102294
Zabelina, V. et al. Mutation in Bombyx mori fibrohexamerin (P25) gene causes reorganization of rough endoplasmic reticulum in posterior silk gland cells and alters morphology of fibroin secretory globules in the silk gland lumen. Insect. Biochem. Mol. Biol. 135, 103607. 10.1016/j.ibmb.2021.103607 (2021).34102294 10.1016/j.ibmb.2021.103607
23. Hughes CS Single-pot, solid-phase-enhanced sample preparation for proteomics experiments Nat. Protoc. 2019 14 68 85 10.1038/s41596-018-0082-x 30464214
Hughes, C. S. et al. Single-pot, solid-phase-enhanced sample preparation for proteomics experiments. Nat. Protoc. 14, 68–85. 10.1038/s41596-018-0082-x (2019).30464214 10.1038/s41596-018-0082-x
24. Zhang Y Comparative proteome analysis of multi-layer cocoon of the silkworm, Bombyx mori PLoS One 2015 10 e0123403 10.1371/journal.pone.0123403 25860555
Zhang, Y. et al. Comparative proteome analysis of multi-layer cocoon of the silkworm, Bombyx mori. PLoS One 10, e0123403. 10.1371/journal.pone.0123403 (2015).25860555 10.1371/journal.pone.0123403
25. Cox J Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction, termed MaxLFQ Mol. Cell Proteom. 2014 13 2513 2526 10.1074/mcp.M113.031591
Cox, J. et al. Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction, termed MaxLFQ. Mol. Cell Proteom. 13, 2513–2526. 10.1074/mcp.M113.031591 (2014).10.1074/mcp.M113.031591
26. Tyanova S The Perseus computational platform for comprehensive analysis of (prote)omics data Nat. Methods 2016 13 731 740 10.1038/nmeth.3901 27348712
Tyanova, S. et al. The Perseus computational platform for comprehensive analysis of (prote)omics data. Nat. Methods 13, 731–740. 10.1038/nmeth.3901 (2016).27348712 10.1038/nmeth.3901
27. Lou R Benchmarking commonly used software suites and analysis workflows for DIA proteomics and phosphoproteomics Nat. Commun. 2023 14 94 10.1038/s41467-022-35740-1 36609502
Lou, R. et al. Benchmarking commonly used software suites and analysis workflows for DIA proteomics and phosphoproteomics. Nat. Commun. 14, 94. 10.1038/s41467-022-35740-1 (2023).36609502 10.1038/s41467-022-35740-1
28. Visser S Volenikova A Nguyen P Verhulst EC Marec F A conserved role of the duplicated Masculinizer gene in sex determination of the Mediterranean flour moth, Ephestia kuehniella PLoS Genet. 2021 17 e1009420 10.1371/journal.pgen.1009420 34339412
Visser, S., Volenikova, A., Nguyen, P., Verhulst, E. C. & Marec, F. A conserved role of the duplicated Masculinizer gene in sex determination of the Mediterranean flour moth, Ephestia kuehniella. PLoS Genet. 17, e1009420. 10.1371/journal.pgen.1009420 (2021).34339412 10.1371/journal.pgen.1009420
29. Lovell JT GENESPACE tracks regions of interest and gene copy number variation across multiple genomes Elife. 2022 10.7554/eLife.78526 36083267
Lovell, J. T. et al. GENESPACE tracks regions of interest and gene copy number variation across multiple genomes. Elife.10.7554/eLife.78526 (2022).36083267 10.7554/eLife.78526
30. Wickham, H. ggplot2: Elegant graphics for data analysis. Use R, 1–212. 10.1007/978-0-387-98141-3 (2009).
31. Kumar S Stecher G Tamura K MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets Mol. Biol. Evol. 2016 33 1870 1874 10.1093/molbev/msw054 27004904
Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874. 10.1093/molbev/msw054 (2016).27004904 10.1093/molbev/msw054
32. Hoang DT Chernomor O von Haeseler A Minh BQ Vinh LS UFBoot2: Improving the ultrafast bootstrap approximation Mol. Biol. Evol. 2018 35 518 522 10.1093/molbev/msx281 29077904
Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: Improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, 518–522. 10.1093/molbev/msx281 (2018).29077904 10.1093/molbev/msx281
33. Nguyen LT Schmidt HA von Haeseler A Minh BQ IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies Mol. Biol. Evol. 2015 32 268 274 10.1093/molbev/msu300 25371430
Nguyen, L. T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274. 10.1093/molbev/msu300 (2015).25371430 10.1093/molbev/msu300
34. Kalyaanamoorthy S Minh BQ Wong TKF von Haeseler A Jermiin LS ModelFinder: Fast model selection for accurate phylogenetic estimates Nat. Methods 2017 14 587 10.1038/Nmeth.4285 28481363
Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: Fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587. 10.1038/Nmeth.4285 (2017).28481363 10.1038/Nmeth.4285
35. Dong Z Comparative proteomics reveal diverse functions and dynamic changes of Bombyx mori silk proteins spun from different development stages J. Proteome Res. 2013 12 5213 5222 10.1021/pr4005772 24093152
Dong, Z. et al. Comparative proteomics reveal diverse functions and dynamic changes of Bombyx mori silk proteins spun from different development stages. J. Proteome Res. 12, 5213–5222. 10.1021/pr4005772 (2013).24093152 10.1021/pr4005772
36. Wang J Vertebrate gene predictions and the problem of large genes Nat. Rev. Genet. 2003 4 741 749 10.1038/nrg1160 12951575
Wang, J. et al. Vertebrate gene predictions and the problem of large genes. Nat. Rev. Genet. 4, 741–749. 10.1038/nrg1160 (2003).12951575 10.1038/nrg1160
37. Levin M Butter F Proteotranscriptomics—A facilitator in omics research Comput. Struct. Biotechnol. 2022 20 3667 3675 10.1016/j.csbj.2022.07.007
Levin, M. & Butter, F. Proteotranscriptomics—A facilitator in omics research. Comput. Struct. Biotechnol. 20, 3667–3675. 10.1016/j.csbj.2022.07.007 (2022).10.1016/j.csbj.2022.07.007
38. Sezutsu H Yukuhiro K Dynamic rearrangement within the silk fibroin gene is associated with four types of repetitive units J. Mol. Evol. 2000 51 329 338 10.1007/s002390010095 11040284
Sezutsu, H. & Yukuhiro, K. Dynamic rearrangement within the silk fibroin gene is associated with four types of repetitive units. J. Mol. Evol. 51, 329–338. 10.1007/s002390010095 (2000).11040284 10.1007/s002390010095
39. Rouhova L Comprehensive analysis of silk proteins and gland compartments in Limnephilus lunatus, a case-making trichopteran BMC Genom. 2024 25 472 10.1186/s12864-024-10381-4
Rouhova, L. et al. Comprehensive analysis of silk proteins and gland compartments in Limnephilus lunatus, a case-making trichopteran. BMC Genom. 25, 472. 10.1186/s12864-024-10381-4 (2024).10.1186/s12864-024-10381-4
40. Li YR Wei YK Zhang GZ Zhang YS Sericin from fibroin-deficient silkworms served as a promising resource for biomedicine Polymers-Basel. 2023 10.3390/polym15132941 38232044
Li, Y. R., Wei, Y. K., Zhang, G. Z. & Zhang, Y. S. Sericin from fibroin-deficient silkworms served as a promising resource for biomedicine. Polymers-Basel.10.3390/polym15132941 (2023).38232044 10.3390/polym15132941
41. Prasong S Yaowalak S Wilaiwan S Characteristics of silk fiber with and without sericin component: A comparison between Bombyx mori and Philosamia ricini silks Pak. J. Biol. Sci. 2009 12 872 876 10.3923/pjbs.2009.872.876 19803122
Prasong, S., Yaowalak, S. & Wilaiwan, S. Characteristics of silk fiber with and without sericin component: A comparison between Bombyx mori and Philosamia ricini silks. Pak. J. Biol. Sci. 12, 872–876 (2009).19803122 10.3923/pjbs.2009.872.876
42. Muller YA Heiring C Misselwitz R Welfle K Welfle H The cystine knot promotes folding and not thermodynamic stability in vascular endothelial growth factor J. Biol. Chem. 2002 277 43410 43416 10.1074/jbc.M206438200 12207021
Muller, Y. A., Heiring, C., Misselwitz, R., Welfle, K. & Welfle, H. The cystine knot promotes folding and not thermodynamic stability in vascular endothelial growth factor. J. Biol. Chem. 277, 43410–43416. 10.1074/jbc.M206438200 (2002).12207021 10.1074/jbc.M206438200
43. Brown JC Role of gene length in control of human gene expression: Chromosome-specific and tissue-specific effects Int. J. Genom. 2021 10.1155/2021/8902428
Brown, J. C. Role of gene length in control of human gene expression: Chromosome-specific and tissue-specific effects. Int. J. Genom.10.1155/2021/8902428 (2021).10.1155/2021/8902428
44. Castillo-Davis CI Mekhedov SL Hartl DL Koonin EV Kondrashov FA Selection for short introns in highly expressed genes Nat. Genet. 2002 31 415 418 10.1038/ng940 12134150
Castillo-Davis, C. I., Mekhedov, S. L., Hartl, D. L., Koonin, E. V. & Kondrashov, F. A. Selection for short introns in highly expressed genes. Nat. Genet. 31, 415–418. 10.1038/ng940 (2002).12134150 10.1038/ng940
45. Gai TT Cocoonase is indispensable for Lepidoptera insects breaking the sealed cocoon PLoS Genet. 2020 10.1371/journal.pgen.1009004 32986696
Gai, T. T. et al. Cocoonase is indispensable for Lepidoptera insects breaking the sealed cocoon. PLoS Genet.10.1371/journal.pgen.1009004 (2020).32986696 10.1371/journal.pgen.1009004
