Single-Cell Deconvolution Reveals Phenotype-Associated Cellular States in the Silk Glands of Bombyx mori and Its Wild Ancestor

TL;DR

This study shows how domestication changed the cellular states in silkworms to improve silk production, using single-cell analysis to compare domestic and wild silkworms.

Contribution

The study introduces a cellular framework to understand how domestication reshaped silk glands for optimized silk production.

Findings

Domestic silkworms have a 'pro-synthesis' cellular state with activated silk protein genes and metabolism.

Wild silkworms maintain a 'protective–adaptive' state focused on stress response and xenobiotic metabolism.

Pseudotime analysis identified key gene expression changes linked to high silk yield during domestication.

Abstract

This study investigated how domestication altered silk production by comparing the cellular states in silk glands of domestic silkworm and its wild ancestor at single-cell resolution. We identified a “pro-synthesis” cellular state in domestic silkworms, characterized by activated silk protein genes and enhanced metabolism, whereas wild silkworm cells maintained a “protective–adaptive” state geared toward stress response. Trajectory analysis revealed key genetic switches associated with high silk yield. Our results demonstrate that domestication reshaped the silk gland cellular landscape toward optimized production. This work provides a new cellular framework for understanding the evolution of complex traits under selection, with implications for agriculture and functional biology.

Silk production is a classic example of a domestication trait, yet the cell-type-specific driver of its enhancement in the silkworm Bombyx mori remains unresolved. To address this, we integrated extensive bulk RNA-seq data with a single-nucleus RNA-seq atlas of silk glands (SGs) from domestic B. mori and wild B. mandarina for deconvolution analysis. This identified phenotype-associated cell subpopulations (Scissor+ and Scissor− cells) that enrich in B. mori and B. mandarina, respectively. Transcriptomic characterization revealed that B. mori SG cells exhibit a pervasive “pro-synthesis” transcriptional state, with concerted upregulation of silk protein genes and metabolic pathways. Conversely, B. mandarina cells maintained a “protective–adaptive” state, enriched for stress response and xenobiotic metabolism genes. Pseudotime analysis further delineated the cell state transitions, pinpointing key dynamic gene expression linked to high silk yield. Our findings demonstrate that domestication reshaped the silk gland cellular landscape, promoting a systemic shift toward a synthesis-optimized cell state. This study offers a new framework at the cellular level to elucidate the evolution of complex traits under selection.

Figures

Full text

1. Introduction

The domestication of wild species often leads to dramatic improvements in key traits, yet the cell-type-specific molecular mechanisms underlying such changes are frequently obscure [1,2]. The silkworm, Bombyx mori (Lepidoptera), domesticated from its wild ancestor Bombyx mandarina for prolific silk production, presents a powerful model for deciphering the genomic and cellular basis of a selected complex trait [3,4].

Comparative studies between these two species have identified several factors contributing to higher silk yield. Comprehensive transcriptomic analyses of the silk gland (SG), the specialized organ that synthesize silk fibroin and sericin proteins, reveal that B. mori exhibits enhanced expression of genes involved in protein processing, ribosome biogenesis, and energy metabolism [5]. These changes correlate with increased cell numbers in the posterior silk gland (PSG) and elevated translational efficiency [4,5]. Genomic scans further indicate that domestication selected not only core silk protein genes but also a supporting network of genes for nitrogen and amino acid metabolism, potentially fueling the biosynthesis of silk proteins [3,4]. Moreover, the SG undergoes a distinct developmental transition during the final larval instar, characterized by the massive induction of silk protein genes (e.g., fibH, fibL, P25, Ser1, and Ser3) and key transcription factors (TFs) [6,7]. Collectively, this evidence suggests that domestication has reshaped the SG as multiple biological levels.

However, these insights primarily derive from bulk tissue analyses, which are unable to resolve cellular heterogeneity of the SG. Consequently, the specific cell populations that drive the differential silk output and their precise molecular signatures remain unresolved. To dissect the SG ecosystem at cellular resolution, we recently constructed the first single-nucleus RNA sequencing (snRNA-seq) atlas for SGs of both B. mori and B. mandarina [8]. Building on this resource, here, we apply a deconvolution approach (Scissor algorithm [9]) to integrate our single-cell data with a large compendium of bulk SG transcriptomes. This strategy allows us to identify the cell subpopulations most highly associated with B. mori and B. mandarina, compare their cell state trajectories, and delineate the divergent gene expression programs and functional pathways within these phenotype-associated cells. Our work moves beyond a static cell type catalog to define the dynamic, phenotype-associated cellular states that systemically underpin the evolution of high silk productivity during domestication.

2. Materials and Methods

2.1. Data Acquisition

A total of 87 independent public bulk RNA-seq datasets were sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (accessed on 11 November 2024) according to the GSE number. The bulk dataset encompasses samples from different silk gland regions, including anterior silk gland (ASG), middle silk gland (MSG), PSG, and SG, and across multiple developmental stages (see Table S1 for full details). The stages of a subset of these bulk samples overlap with the snRNA-seq stages used in this study (day 3 of the fifth larval instar (L5D3) and day 5 of the fifth larval instar (L5D5) from B. mori, day 6 of the fifth larval instar (L5D6) from B. mandarina). Single-nuclei RNA sequencing (snRNA-seq) data were obtained from our previous study [8] (BioProject ID: PRJCA022374; accession ID: CRA014260).

2.2. Analysis of Bulk RNA-Seq and snRNA-Seq

Raw bulk RNA-seq reads were first subjected to quality control using FastQC (v0.12.1) [10] (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) (accessed on 22 January 2025) to assess the sequencing quality. Low-quality bases and adaptor sequences were subsequently removed with Trimmomatic (v0.39, RWTH Aachen, Aachen, Germany) [11]. The cleaned reads were subsequently aligned to the reference genome (https://silkbase.ab.a.u-tokyo.ac.jp/cgi-bin/download.cgi) (accessed on 29 January 2025) using STAR (v2.7.11b, CSHL, Cold Spring Harbor, NY, USA) [12] with the default parameters. Gene- and transcript-level expression was quantified with RSEM (v1.3.3, UW-Madison, Madison, WI, USA) [13] using the rsem-calculate-expression module, generating both count and TPM matrices. The ‘ComBat_deq’ function of Sva (v3.50.0, JHSPH, Baltimore, MD, USA) [14] was used to remove batch effects from the samples. The input for ComBat_deq was the log2 (TPM+1)-transformed expression matrix of the bulk RNA-seq data. The model included only batch information as a known covariate; no other covariates were included. The differentially expressed genes (DEGs) identified in the bulk RNA-seq data were analyzed using the limma-voom pipeline (v3.58.1, UoM, Parkville, Victoria, Australia) [15] package with the criteria of logFC > 0.8 and adj.p.adj < 0.05. The snRNA-seq datasets were processed and analyzed using the Seurat (v5.3.0, NYU, New York, NY, USA) [16] R package. We used Seurat’s FindMarkers function with the Wilcoxon rank-sum test. p-values were corrected for multiple testing using the Bonferroni method within each cell group comparison. FindIntegrationAnchors and IntegrateData with default parameters. LogNormalize with a scale factor of 10,000.

2.3. Scissor Analysis

Single-cell identification of subpopulations with bulk sample phenotype correlation analysis was performed with Scissor (v2.0.0, OHSU, Portland, OR, USA) [9]. The expression matrix, phenotype information (0 for B. mandarina and 1 for B. mori) of the bulk RNA-seq data, and snRNA-seq data were analyzed, with the parameter alpha set to 0.3. Alpha = 0.3 was selected after preliminary testing as it provided a balanced trade-off between identifying a sufficient number of phenotype-associated cells and maintaining a stringent correlation threshold. The bulk RNA-seq TPM matrix was log2 (TPM+1) transformed and then Z-score normalized per gene. The snRNA-seq data were normalized using Seurat’s LogNormalize method. The intersection of highly variable genes (HVGs) identified separately in the bulk and single-nucleus datasets was used as the common feature set for Scissor analysis. The bulk RNA-seq data consisted of 87 SG samples, of which 16 were from B. mandarina and 71 were from B. mori. The snRNA-seq data consisted of 3 SG samples, including L5D3 and L5D5 from B. mori and L5D6 from B. mandarina (Table S2). L5D3 and L5D5 for B. mori represent the beginning and sustained high-synthesis phases of silk production, respectively. L5D6 for B. mandarina was chosen based on its physiological similarity to B. mori L5D5.

We subsequently calculated the fractions of Scissor+ cells (positively associated with B. mori), Scissor− cells (positively associated with B. mandarina), and background cells in different samples and cell types. We further identified DEGs in the Scissor+ cells and Scissor− cells using the Seurat package ‘FindMarkers’ with the criteria of |avg_log2FC| > 0.5 and p.adj < 0.05.

2.4. Pseudotime Analysis

Pseudotime analysis was performed on the scissor cells using the Monocle 2 (v2.32.0, UW, Seattle, WA, USA) [17] R package. The structure of the trajectory was plotted in two-dimensional space using the DDRTree dimensionality reduction algorithm, and the Scissor cells were ordered according to pseudotime, with pseudotime rooted at the L5D3 developmental stage.

Pseudotime analysis was performed on cocoon sericin-synthesizing and secreting cells (FBSs) using the Monocle 3 (v1.3.1, UW, Seattle, WA, USA) [18] R package. A trajectory graph was built on the FBSs UMAP using the “learn_graph” function with the default parameters. Pseudotime was calculated with the “order_cells” function, after which the start of pseudotime was oriented on the basis of the 4th molting (L4M) development stage. Afterward, the cells were grouped into 50 evenly sized bins throughout each trajectory, and loess regression was applied to the average module eigengene (ME) of Modules 4 and 7 in these bins to inspect the dynamics of each module throughout the FBSs trajectory using the “PlotModuleTrajectory” function.

2.5. Gene Set Enrichment Analysis (GSEA)

Gene set enrichment analysis was performed on the gene matrix provided in Table S3 using the desktop tool obtained from the GSEA website (http://software.broadinstitute.org/gsea/index.jsp) (accessed on 8 March 2025).

2.6. Quantitative Real-Time PCR (qRT-PCR) Analysis

The SGs of L5D5 B. mori (Dazao, under standard experimental conditions and using fresh mulberry leaves, was raised at a temperature of 25 °C) and L5D6 B. mandarina (Collected from Anhui Province, China) were dissected, and total RNA was extracted using a tissue RNA extraction kit (Omega, Irving, TX, USA)). cDNA was synthesized from the extracted RNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Osaka, Japan) according to the manufacturer’s instructions. Three independent biological replicates were performed for each group, with each replicate consisting of a pooled sample. Each 20 μL reaction contained 2 μL of cDNA template, 10 μL of 1 × SYBR Green I master mix (TB Green Premix Ex Taq II, Takara), 0.4 μL of ROX Reference Dye (Takara), and 0.8 μL each of forward and reverse primer (sequences are listed in Table S4). BmeIF-4a (SilkDB Probe: sw22934) was used as an internal control for normalization. Amplification was carried out on a 7500 FAST Real-Time PCR System (ABI, Los Angeles, CA, USA). Relative gene expression levels were calculated using the 2^−ΔΔCt^ method. Statistical comparisons between groups were performed using unpaired two-tailed Student’s t-tests, and data are presented as mean ± standard error of the mean (SEM). All analyses and graphical representations were conducted with GraphPad Prism software (v9.0).

3. Results

3.1. Dynamic Reshuffling of Cellular Composition During SG Development

The SG undergoes coordinated changes in cellular composition to support stage-specific functions. Using the snRNA-seq atlas of B. mori and B. mandarina from our previous study [8], we first quantified the proportional dynamics of major SG cell types across key developmental stages (Figure 1a). As anticipated, cell types constituting the ASG—primarily involved in silk protein processing rather than synthesis—maintained stable proportions after tissue remodeling at L4M, indicating functional stabilization prior to spinning (Figure 1b, Table S5). In contrast, the proportions of the primary silk protein-synthesizing cell types exhibited marked stage-specific heterogeneity. FBSs progressively increased in abundance throughout the fifth instar. Conversely, cocoon sericin-synthesizing and secreting cells (CSSs) peaked at L5D3 and subsequently declined at L5D5 (Figure 1b). A critical divergence between species was evident: B. mori possessed a significantly higher proportion of FBSs but a lower proportion of CSSs compared to B. mandarina at matched time points (Figure 1b, Table S5). This differential cellular architecture aligns with and provides a cellular-resolution explanation for their distinct physiological priorities: the domesticated silkworm allocates more cellular resources to fibroin production for cocoon bulk, whereas its wild ancestor commits a greater fraction of cells to sericin synthesis, likely crucial for constructing a more protective cocoon structure.

3.2. Identification of Phenotype-Associated Cell Subpopulations via Deconvolution

To move beyond static cell type classifications and identify specific cell subpopulations whose states are highly linked to the domestic or wild phenotype, we employed the Scissor algorithm. This approach integrates single-cell data with bulk transcriptomic signatures, here using 87 external SG RNA-seq datasets as a reference to deconvolute phenotype-associated cells. From a total of 39,477 high-quality cells, Scissor selected 2549 “Scissor+” cells and 1617 “Scissor−” cells as being most highly associated with the B. mori and B. mandarina phenotypes, respectively (Figure 2a, Table S6).

Examining their distribution across samples revealed a clear biological pattern. Scissor+ cells were overwhelmingly dominant in B. mori samples at L5D3, coinciding with the previously documented peak period of cell–cell communication and the onset of high-efficiency silk protein synthesize [8]. Conversely, Scissor− cells reached their highest absolute proportion in B. mandarina at L5D6, underscoring their association with the wild phenotype (Figure 2b, Table S6). Analysis across cell types showed that the distribution of Scissor+ and Scissor− cells was not uniform but exhibited significant bias in three key types: ESMs, CSSs, and NCSs (Figure 2c, Table S6). This biased distribution suggests divergent functional priorities at the cellular level. The significant enrichment of Scissor+ cells within ESMs implies an elevated metabolic or support role in B. mori, potentially required to fuel its high-rate biosynthesis of silk proteins. In contrast, the skew of Scissor+/− cells toward CSSs and NCSs, respectively, indicates enhanced cellular commitment to the cocoon sericin and noncocoon sericin production pathways in B. mori and B. mandarina. This observation reveals, at single-cell resolution, differences in sericin synthesis between the domesticated silkworm and its wild ancestor.

3.3. Convergent Cell State Transitions and Terminal State Activation in Phenotype-Associated Cells

We next sought to understand the cell state of the Scissor+/− cells. Pseudotime trajectory analysis of these cells revealed a branching path that culminated in three distinct transcriptional states (Figure 3a). Based on the distribution of cells from early and late developmental stages, we designated state 2 as an initial state and state 1 as the terminal state (Figure 3b and Figure S1a). Intriguingly, both Scissor+ and Scissor− cells showed increased abundance along the pseudotime trajectory and were predominantly localized to the terminal State 1 (Figure 3c and Figure S1b, Table S7). This terminal state (state 1) was also precisely where the core silk-producing cell types (FBSs and CSSs) were most highly enriched (Figure 3d). This convergence indicates that despite their molecular differences, cells from both species progress towards a shared, highly specialized “synthesis-competent” state at the endpoint of SG maturation.

Gene expression dynamics along this trajectory were categorized into three major patterns (cluster 1–3), each peaking in a corresponding state (Figure S1c, Table S8). The gene program activated in the terminal State 1 was particularly revealing. It included not only core silk protein genes (e.g., P25, Ser1, Ser3) but also genes with diverse auxiliary functions: Cbp (involved in pigment formation) [19], serpin-15 (implicated in antimicrobial activity) [20], and a myosin heavy chain gene (KWMTBOMO04757) [21] that potentially related to cellular contractility or vesicle transport (Figure 3e, Table S8). The coordinated upregulation of this functionally diverse suite of genes underscores a comprehensive, conserved developmental program that activates the multifaceted machinery necessary for final silk production.

3.4. Molecular Signatures Defining the Pro-Synthesis and Protective–Adaptive Cellular States

To dissect the transcriptional basis of the phenotype-associated states, we identified DEGs across both Scissor+ and Scissor− cells, as well as between different cell types (Table S9). The expression patterns of top DEGs were corroborated in B. mori and B. mandarina SGs at the tissue level, confirming the robustness of our deconvolution approach (Figure 4a,b). These DEGs were enriched for several key functional categories.

3.4.1. Genes Linked to Domestication, Growth, and Environmental Response

We interrogated the DEG list for genes with known roles in silkworm biology. Several genes previously associated with domestication or SG development showed differential expression (Figure 4c). Scissor+ cells upregulated KWMTBOMO14482 (linked to cell volume regulation [22,23]), KWMTBOMO09538 (essential for influencing cell size [24]) and KWMTBOMO08097 (associated with posterior SG expansion [25,26]). In contrast, Scissor− cells upregulated KWMTBOMO15198, an SG-specific gene that enhances the UV-protective capacity of the cocoon, a trait likely critical for wild survival [27,28]. Furthermore, Scissor+ cells exhibited higher expression of positive regulators of cell cycle and growth (e.g., Cdt1, Ras1), while Scissor− cells upregulated potential growth inhibitors such as Geminin [29] and fibroblast gene (KWMTBOMO13997) [30].

3.4.2. Metabolic Pathway Divergence

Gene set enrichment analysis (GSEA) uncovered profound metabolic differences. Despite comparable abundance in CBCs (Figure 2c), Scissor− cells were significantly enriched for pathways like “metabolism of xenobiotics by cytochrome P450” (Figure 5a, Table S3). Despite B. mandarina possessing a smaller genomic repertoire of cytochrome P450 genes [31], the significant enrichment of this pathway in Scissor− cells suggests a heightened functional emphasis on detoxification-related transcription in the wild species’ SG, potentially reflecting a key adaptive difference. Conversely, within ESMs, Scissor+ cells specifically upregulated BmDJ-1β (Figure 5b,c, Table S9), a gene implicated in maintaining mitochondrial homeostasis under oxidative stress [32,33], which may support the energy balance required for sustained high-rate biosynthesis.

3.4.3. Regulation of Silk Protein Synthesis

KWMTBOMO14478 and KWMTBOMO12520, which are highly positively correlated with silk protein synthesis [4], were upregulated in Scissor+ cells (Figure 4c). Focusing on sericin-synthesizing cells (CSSs and NCSs), we found distinct expression patterns for sericin genes and their regulators (Figure 5b,c, Table S9). Ser2 and Ser3 were upregulated in Scissor− cells of NCSs and CSSs, respectively. Notably, in CSSs, Scissor+ cells co-upregulated Ser1 and its positive transcriptional regulator, BmSuc1. Given that BmSuc1 is known to promote Ser1 expression, silk fiber mechanics, and overall yield [34,35], this coordinated activation provides a precise molecular correlate for the superior Ser1 production and cocoon weight observed in B. mori.

In summary, the molecular profiles of Scissor+ and Scissor− cells encapsulate the fundamental trade-off shaped by domestication. Scissor+ cells in B. mori are molecularly tuned for enhanced growth, anabolic output, and efficient production of key silk proteins, defining a “pro-synthesis” state. Scissor− cells in B. mandarina retain signatures of environmental resilience, including stress response and detoxification pathways, constituting a “protective–adaptive” state. Our deconvolution framework thus successfully maps the domestication-driven divergence in silk production to specific cellular subpopulations and their defining transcriptional programs.

4. Discussion

Single-cell spatiotemporal omics provides a high-resolution research perspective for insect science [36]. This study leverages a high-resolution snRNA-seq dataset, combined with a deconvolution approach, to dissect the cellular and molecular basis of silk protein synthesis divergence between B. mori and B. mandarina. Our analysis reveals that the evolution of high silk yield in B. mori involves a system-level reorganization, characterized by recalibrated cellular proportions, distinct transcriptional states, and rewired functional genes within the SG ecosystem.

A key finding is the altered allocation of cellular resources. The significant expansion of fibroin-synthesizing FBSs in B. mori, alongside a relative contraction of sericin-producing CSSs compared to B. mandarina, provides a cellular basis for the domesticated silkworm’s fibroin-rich, high-yield cocoons. This aligns with the breeding goal of maximizing silk output [37]. Conversely, the wild silkworm maintains a larger proportion of CSSs, consistent with an ecological need for enhanced cocoon protection. This fundamental trade-off between “productivity” (fibroin) and “protection” (sericin) is reflected in the population dynamics of SG cell types (Figure 1). This notion is further supported by metabolite studies showing that wild silkworm cocoons contain higher levels of protective metabolites, such as citric acid and long-chain hydrocarbons, which may act as antimicrobials and water repellents [38].

Applying the Scissor algorithm, we identified discrete, phenotype-associated subpopulations (Scissor+ and Scissor− cells) linked to B. mori and B. mandarina, respectively (Figure 2). This indicates that molecular signatures of domestication are concentrated within specific cellular states. The preferential association of Scissor+ cells with ESMs in B. mori suggests that heightened silk output is supported by an augmented metabolic support system. This observation aligns with transcriptomic studies indicating that improved silkworm strains exhibit repressed basic nitrogen biosynthesis but enhanced dynamics of protein post-translational modification and resource reallocation to the SG, likely for efficient silk production [37].

Pseudotime trajectory analysis revealed that both Scissor+ and Scissor− cells progress toward a similar terminal cellular state (State 1). This state is enriched for key silk-producing cells (FBSs and CSSs) and exhibits peak expression of core silk protein genes (e.g., P25, Ser1, Ser3) and secretion-related factors (e.g., Cbp, Myosin). This points to a conserved functional configuration for efficient synthesis (Figure 3). However, the paths to this state and the molecular composition of the phenotype-associated cells exhibit species-specific differences. Differential expression of genes governing cell cycle progression (Geminin/Cdt1), growth signaling (Ras1), and cellular expansion (TBC1 domain genes) in Scissor+ versus Scissor− cells suggests divergent developmental tuning, which may contribute to the larger SG size in B. mori (Figure 4). These findings are consistent with genetic studies that have identified quantitative trait loci associated with cocoon yield traits [39]. Furthermore, the significant enrichment of xenobiotic metabolism pathways, particularly those involving cytochrome P450 enzymes, in B. mandarina Scissor− cells likely reflects evolutionary pressures on wild populations to detoxify environmental compounds. Interestingly, although B. mandarina possesses a smaller cytochrome P450 gene repository than B. mori, our analysis suggests that specific P450 genes may be crucial for this function, a selective constraint that appears relaxed during domestication [31]. This aligns with the idea that wild silkworms require metabolic adaptations to cope with natural environments.

Overall, this study advances our understanding of insect secretory organs from a static catalog of cell types to a dynamic model of a cellular ecosystem shaped by selection. Methodologically, we demonstrate the power of deconvolution frameworks like Scissor in entomology, where they can unlock cellular insights from the rich legacy of bulk transcriptomic studies, especially when comprehensive single-cell atlases are still being built. Our findings establish a cell-state perspective for understanding the evolution of complex economic traits, with implications beyond sericulture for the study of animal domestication and specialized tissue function. From a practical standpoint, the cell-state signatures we identified, such as the ‘pro-synthesis’ transcriptional program in ESMs and FBSs, provide a new set of candidate markers for molecular breeding programs aimed at enhancing silk yield or quality. The specific regulatory genes (e.g., BmSuc1) and pathways highlighted offer precise targets for functional validation via genetic perturbation to directly test their roles in silk production.

We acknowledge several technical limitations. First, incomplete developmental stage matching between species may confound species-specific interpretations, and pseudotime inferences lack robustness checks such as alternative root selection or dimensionality reduction. Second, tissue-level qRT-PCR does not provide cell-state-resolved validation, which would benefit from future in situ single-molecule imaging. Third, despite standard normalization and the use of pooled samples, technical differences in snRNA-seq (e.g., capture efficiency) and limited data availability for B. mandarina remain considerations. Therefore, while our deconvolution framework identifies coordinated cellular states associated with domestication, claims of ecosystem-level reconfiguration should be interpreted with caution and regarded as hypothesis-generating pending validation with stage-matched, replication-rich datasets and orthogonal spatial approaches.