Single-cell and spatial transcriptomics define 20E-driven developmental reprogramming in silkworm wing disc

Domain: chemistry
Published: 2026-02-23
Editorial outcome: minor_revision

Abstract
Insect wing development involves tissue patterning, cell fate transitions, and hormone signaling, yet its spatiotemporal logic remains unclear. The silkworm, with large wing discs and defined stages, provides an ideal model for high resolution analysis. Here, we construct a spatiotemporal single-cell atlas of the silkworm wing disc across 10 time points, identifying 12 major cell types and their developmental transitions. Wing morphogenesis (Wm) cells act as central progenitors, differentiating into epithelial and cuticle lineages under lineage-specific transcription factors. Time‑resolved snRNA‑seq reveal hierarchical transcriptional reprogramming, with Wm cells functioning as early signaling hubs. Functional modules and signaling pathways were activated in spatiotemporal controlled manner. 20‑hydroxyecdysone treatment rapidly accelerates fate transitions and gene expression, recapitulating natural development within hours. Integration of morphology, hormone levels, and gene expression supports a five-stage Gene Transition Model describing progressive fate resolution. This work reveals wing development in silkworm and provides insights into hormone-driven organogenesis and potential manipulation of insect development in agriculture.

Introduction
Insects are the only invertebrates that can fly, a key evolutionary innovation that allows them to travel long distances, escape predators, locate mates, and expand their foraging range. This ability has enabled insects to occupy diverse ecological niches and become the most species-rich group of animals on Earth1. Understanding how insect wings develop is thus fundamental to developmental biology and has practical value for pest control, evolutionary studies, and bio-inspired engineering. In holometabolous insects, their wings develop from specialized larval structures known as imaginal discs, particularly wing discs2. Wing discs are composed of diverse cell types that undergo continuous morphogenetic changes, ultimately forming the fully developed insect wing3,4. Knowing the spatial distribution and functional roles of these cells is important for understanding how wings form. Previous studies have shown that hormones and several signaling pathways play key roles in wing disc development5,6. Among them, 20-hydroxyecdysone (20E) is the main hormone, driving the metamorphosis from larvae to adult stages by guiding cell proliferation, differentiation, and apoptosis5,7. Meanwhile, classic pathways such as Hedgehog (Hh), Wnt, and Notch also play roles in cell fate, tissue organization, and organ formation8–10.

Although single-cell analyses have been conducted for the wing disc of Drosophila melanogaster (D. melanogaster) and Bicyclus anynana (B. anynana)11,12, they have mainly focused on cell type composition and gene expression patterns at isolated stages. These studies fail to capture the full continuum of cellular transitions across developmental time and often lack spatial resolution. Consequently, it remains unclear how different cell types interact, how their spatial arrangement evolves, and how they collectively respond to dynamic 20E signaling. A comprehensive, spatiotemporal view of wing disc development is imperative to fill these gaps.

The silkworm Bombyx mori (B. mori), is an important lepidopteran model of economic and developmental significance and serves as an ideal system for investigating developmental processes13,14. Bulk transcriptomic studies have identified temporal gene expression dynamics involving regulators such as BmBlimp-115, BmIGFLP16, BmChiNAG17, and BmSd18. However, these studies often focus on either temporal or spatial dimensions alone. Furthermore, bulk transcriptomic approaches average signals across heterogeneous tissues, lacking the resolution to uncover cell-type-specific dynamics or spatial regulatory cues. This makes it hard to understand how different cells change and interact within the tissue. Recent advances in single-cell RNA sequencing (scRNA-seq) and spatial transcriptomic technologies provide powerful tools to resolve cellular heterogeneity and gene expression dynamics at high resolution19–21. By integrating these approaches, it is possible to see how different cell types are arranged, how they develop over time, and how they communicate22–26.

In this work, we combined scRNA-seq, single-nucleus RNA sequencing (snRNA-seq), and spatial transcriptomics (Stereo-seq) to construct a high-resolution spatiotemporal transcriptomic atlas of B. mori wing discs from larva to pupa. We systematically mapped the dynamics of 12 functional cell populations in wing disc and revealed molecular mechanisms underlying their differentiation trajectories, which are regulated by key transcription factors such as Rfx, Blimp-1, Dll, and Pur-alpha. Performing snRNA-seq on wing disc subjected to time-constrained 20E treatment, we discovered a “time-axis compression” effect. In this process, 20E caused rapid reprogramming of cell fates. Finally, we built a Gene Transition Model (GTM) explaining how cellular signal cascades respond to 20E levels changes to guide wing development. In summary, our work reveals the dynamic processes that build insect wings and provides new insights in studying wing development and lepidopteran pest insect’s control.

Results
Temporalandspatialatlasofsilkwormwingdisccells
Thesilkwormwingdiscundergoesprogressivemorphologicaltransitions,shiftingfromtissueremodelingduringthelarvalstage(fifthinstar)tospinningandultimatelyformingmaturewingsduringthepupalstage(Fig.1a).Tosystematicallyresolvedevelopmentaltrajectoriesofwingdisccellpopulations,weexploitedthescRNA-seqtechnology,generating11librariesfromsamplesacross10timepointsspanning3developmentalstages:larval(L5D1toL5D7),spinning(W1,W2),andpupalstages(P6,thesixthdayofthepupa)(Fig.1b andSupplementaryData1).Afterqualitycontrol,atotalof126,041high-qualitysingle-celltranscriptomeswereretained,covering16,122geneswithamediangenecountpercellof1989(SupplementaryFig.1a,b).UsingdimensionreductionandUniformManifoldApproximationandProjection(UMAP),thecellsweregroupedinto12clusters.

Based on known marker genes of B. mori22, and other lepidopterans11,12, a Drosophila single-cell RNA-seq dataset27, and enriched Gene Ontology(GO)processes ofthe top 200ranked differential expressed genes from each cell cluster(Fig.1c,SupplementaryFig.1andSupplementaryData2),12majorcelltypesweredefined(Fig.1d).Amongthese,threeclustersshowedspecificexpressionofcuticle-relatedgenes.Theseclusterswerecharacterizedas:Cuticle1(C1)withhighserpingenes(serpin-5andserpin-1028)expressionanditsproportiongraduallydecreasedfromlarvalstage(L5D4)topupalstage(P6),Cuticle2(C2)withhighfibroblastgrowthfactorreceptor(Btl11)expression,andCuticle3(C3)withhighcuticleprotein(CPH28,CPG3 and CPG929) expression and a high proportion of pupalstage (SupplementaryFig.1c,d).ThreeepithelialcelltypeswereidentifiedbasedonmarkerssuchasMf-cpa30(SupplementaryFig.1e).Epithelial1(Ep1)showedastableproportionfromL5D1toW2andalmostdisappearedatP6.Epithelial2(Ep2)wasprominentlyabundantatW2andalmostdisappearedatP6too.Epithelial3(Ep 3) very rare in whole stages exception(SupplementaryFig.1f)11,12,22,27.Wealsoidentifiedsixothermajorcelltypes,includingWingmorphogenesis(Wm),Immune(Im),Matrix (Ma), Metabolic process (Mp), Axon development (Ad), and Ciliated(Ci)cells(SupplementaryFig.1andSupplementaryData2),givingatotalof12majorcelltypes(Fig.1d).Wmcells,markedwithkeygenessuchasFat31andNID32,weremainlyenrichedinpathwaysrelatedtocellmorphogenesisandwingpatterning,includingtubemorphogenesisandimaginaldisc-derivedwingmarginformation(Fig.1candSupplementaryFig.1g).Imcellswereenrichedforhematopoiesisandimmuneorgandevelopment,suggestingtheirinvolvementinimmunecelllineageformation.Ma cells were associated with extracellular matrix organization, including basement membrane and respiratory tubedevelopment.Mpcellsshowedenrichmentinpathwaysregulatingcuticledevelopmentandmetabolicprocesses.Adcellswerelinkedto axon development, neuron projection, and cell–cell recognition, reflectingrolesinneuralwiring.Cicellsexhibitedsignificantenrichmentinciliumassemblyandorganization,microtubule-basedtransport,andcellprojectionassembly.Theseenrichmentprofilesindicatethateachcelltypecontributesspecializeddevelopmentalorfunctionalrolesduringwingdiscdevelopment.

To characterize the spatial distribution of these cell types, we generated Stereo-seq spatial transcriptome data for wing discs at L5D2, L5D3,and W2timepoints (SupplementaryData 3).Stereo-seq data were annotated through deconvolution mapping of gene signaturesderivedfrom12celltypes,atbin20(10µm×10µm)resolution(Fig. 1e and Supplementary Fig. 2a). Combing spatial information (Supplementary Fig. 2b) and histological features (Supplementary Fig.3a–d), werecapitulatedthespatialdistributionofeachcelltypeatL5D2,includinglarvacuticle(Cuticle1–3),epidermis(Epithelial1–3),wingbud(Wmcells),andhematopoieticorgan(ImmuneandMatrixcells) (Fig. 1f). The spatial distribution pattern of cells shifts during wing disc development. At the L5D2 and L5D3 stages, cuticle cells (Cuticle1andCuticle2)weremainlydistributedattheoutermostlayerofthewingdisc,followedbyepithelialcells(Epithelial1andEpithelial2). In contrast, Wm cells were primarily localized in the wing bud region(Fig.1eandSupplementaryFig.3e).BytheW2stage,cuticlecellsstillmadeuptheoutermostlayer,buttheproportionofWmcellsinthecentralregionofthewingdischadsignificantlydecreased,andepithelialcellsbegantooccupyabroaderarea,whichmayberelatedtocellulardifferentiationandmigration(Fig. 1e and SupplementaryFig.3e).

Toexemplifytheutilityoftheannotateddataset,weexamined53genesknowntoplayimportantrolesduringwingdiscdevelopment(Supplementary Data 4). These genes were mapped onto both our single-cellandspatialtranscriptomicdatasetsandexhibitedcell-typeandregion-specificexpressionpatterns,providingrefinedinsightsintohowtheyfunctionatthecellularandspatialdistributionlevelsinwingdiscdevelopment. For instance, STAT33,34 was highly expressed in Cuticle1 cells located in the larval cuticle, Dpp35,36 showed peak expressioninWmcellswithinthewingbudatL5D5(SupplementaryFig.4),andFng37–39wasenrichedinEpithelial1cellswithinthewingbudandepidermisatL5D7(SupplementaryFig.4).Overall,bycombining Stereo-seq and scRNA-seq-based spatiotemporal transcriptomics, we revealed the heterogeneity and dynamic spatiotemporalchangesincelltypecompositioninthesilkwormwingdisc,alongwiththeirdistinctspatialdistributionpatternsacrossdevelopmentalstages.

Wmcellsasadifferentiationhubinwingdevelopment
Basedontheclassificationofmajorcelltypesinthewingdisc(Fig.1d),weanalyzedhowtheirproportionschangedacrossdevelopmentalstages.Fromthelarvalstage(L5)tothepupalstage(P6),severalcelltypes—including Wm, epithelial, cuticle, immune, and matrix cells—showed notable shifts in abundance (Supplementary Fig. 5a). For instance, cuticle cells increased rapidly during the pupal stage, suggesting active morphological transformation (Supplementary Fig.5b,c).Somecelltypes,suchasCuticle3andEpithelial2,appearedat later stages and becamemore prominent during wing formation (SupplementaryFig.5c).Thisisconsistentwiththemajorstructuralchanges that occur near pupation. Notably, Wm cells gradually decreasedfromL5toP6andwerealmostabsentbythepupalstage(SupplementaryFig.5d,e).

To investigate how cell fate transitions contribute to wing disc development,weanalyzedthedifferentiationtrajectoriesofvariouscelltypesusingTOME-tradj40.Mostcelltypesremainedinstablestateswithout clear differentiation trajectories (Supplementary Fig. 6). However, Wm cells showed complex and continuous transitions, suggesting a dynamic differentiation process primarily toward epithelial and cuticle lineages (Supplementary Fig.6g). Additionally, Wm, epithelial,andcuticlecellswerefoundincloseproximityintheUMAPembeddingspace(Fig.1d),whichindicatessimilarityingeneexpressionprofilesamongthesecells.

TofurtherinvestigatethefateofWmcells,wereconstructedtheirdifferentiationtrajectorieswithMonocle2andidentifiedthreedistinctcellularstates(Fig.2a).Wmcellsweremainlyenrichedinstate1,whileepithelial and cuticle cells were primarily found in states 2 and 3, respectively(Fig.2b).Basedontheirtemporaldistribution,state1wasidentifiedasthestartingpoint(Fig.2c).Pseudotimeanalysisfurtherconfirmedabifurcatedtrajectory,withsequentialemergenceofEpithelial1–3andCuticle1–3subtypes,showinganintrinsicmaturationgradientwithineachlineage(Fig.2dandSupplementaryFig.7a,b).Earlysubtypes(Epithelial1,Cuticle1)appearedinlarvalstages,whilelatersubtypes(Epithelial3,Cuticle3)peakedduringthe wandering and pupal stages (Fig. 2e). These findings indicate that Wm cells differentiate into two main cell lineages: epithelial and cuticle cells.

Togainmechanisticinsightsintotheregulationofthesecellfatetrajectories,weanalyzeddifferentiallyexpressedgenes(DEGs)alongpseudotime (Supplementary Data 5). In fate 1, DEGs were mainly enriched in Toll/Imd signaling, longevity-regulating pathways, and endoplasmic reticulum (ER)-related processes (Fig. 2f and SupplementaryFig.7c).Thesepathwaysareassociatedwithtissueremodeling during development. Toll/Imd signaling showed low activity in earlylarvalstagesbutbecamemoreactiveinthepupalstages,suggestingaroleinimmuneregulationandtissuemorphogenesisduringlaterdevelopment(Fig.2g).

Incontrast,fate2cells,whichbecamedominantduringtheW1 andW2period,showed enrichmentinribosome,Hippo,andMAPKpathways (Fig. 2f and Supplementary Fig. 7c). Notably, the Hippo pathway, which is known to regulate somatic growth and cell proliferation41,showedamarkedincreaseinactivityfromL5D5toW2,followedbyarapiddeclineatP6(Fig.2h).Consistentwiththis,Hippo-relatedgenessuchasA3andLix1showeddistinctexpressionpatternsacrossthetwofatebranches(SupplementaryFig.8a).ThispatternsuggeststhatHipposignalingmayregulateWmcellproliferationandtissuemorphogenesisbeforefinaldifferentiation.

To identify upstream regulatorsof Wm differentiation, weperformedSCENICanalysistoidentifykeytranscriptionfactors(TFs)withhighactivityacrossdevelopmentaltimepointsinWmcells(Supplementary Fig. 8b). Based on the Regulon Specificity Score (RSS), we selectedthetopfiveTFsattwokeystagesduringdifferentiation:early(Dll,sqz,Pur-alpha,Rfx)andWperiod(Blimp-1)(Fig.2i).EnrichmentanalysisoftheirdownstreamtargetgenesindicatedthatgenesregulatedbyDllandsqzweremainlyassociatedwithwingmorphogenesisand displayed elevated RSS values at L5D1 and P6, suggesting their functionalroleinguidingWmdifferentiationtowardthecuticlefate(SupplementaryFig.8c,d).Incontrast,targetgenesofPur-alphawereprimarilyassociatedwithepithelial differentiation, whereas those regulated by Blimp-1 showed enriched in Wm-related genes and exhibitedhighRSSscoresfromL5D7toW1,suggestingapotentialroleinpromotingthetransitionfromWmtoEpithelialcells(SupplementaryFig.8e,f).Additionally,Rfxwasidentifiedasacell-type-specificcentral regulator in Wm cells with the highest RSS value at L5D1 (SupplementaryFig.8g).

To functionally validate the role of transcriptional regulators implicatedinWmdifferentiation,wefocusedonRfx,atranscriptionfactor specifically expressed in Wm cells. RNA interference (RNAi)-mediated knockdown of Rfx significantly disrupted its downstream regulatorynetwork,includingseveralkeytargetgenesinvolvedinthedifferentiationofWmcells(Fig.2jandSupplementaryFig.9).Amongthese,severalgeneshavewell-characterizedrolesintissueorganizationandmorphogenesis,suchasbsh42,dpy43,futsch44andpnr45.Within48hoursafterRNAitreatment,wingdiscsshowedclearmorphologicaldefects,includinglobulationinthecentralregion,irregularepithelialinvagination,anddisorganizedtrachealbranching(Fig.2k).Inadults,Rfx knockdown caused severe wing deformities, such as shortened veins,malformedmarginalbristles,andjaggedposteriorwingedges(Fig.2landSupplementaryData6).Inaddition,Rfxishighlyconservedacrosslepidopteratermsofphylogeneticrelationship,motiforganization,anddomainarchitecture(SupplementaryFig.10).Knockdownof its ortholog in the fall armyworm Spodoptera frugiperda (SfRfx, XM_050701463.1)resultedinsimilarmorphologicaldefects,indicatingthat Rfx plays a conserved and essential role in wing disc morphogenesisacrosslepidoptera(Fig.2m–o).Furthermore,RNAi-mediatedknockdownofothercandidateTFs,includingDll,andPur-alpha,also caused structural abnormalities in the wing disc (Supplementary Fig.11).Blimp-1,whichhaspreviously been implicated in regulating wingdiscdevelopment,showedconsistentfunctionaleffectsonwingpatterning15.Together,thesefindingssupportthatRfx,togetherwithBlimp-1,Dll,andPur-alpha,formsatranscriptionalregulatorymoduledrivesthedifferentiationofWmcellsduringwingmorphogenesis.

20Eregulateswingdiscdevelopmentaldynamics
20Eisacentralhormonalregulatorininsectdevelopment,knowntocontrolcelldifferentiation and tissue morphogenesis2,46.To explore how 20Einfluencessilkwormwingdiscdevelopment,wefirstexaminedtheexpressionpatternsof20E-relatedgenesinnaturallydevelopingtissues.Thesegenesdisplayedstrongspatiotemporalspecificityacrosscelltypes(SupplementaryData7).Forexample,20Ereceptor genes such as USP47 were predominantly expressed in Wm cells, whereascuticlecellshighlyexpresseddownstreamsignalingcascadeslike E7448 (Supplementary Fig. 12a). Integration with spatial transcriptomics further confirmed cell-type-specific expression of 20E-responsive genes (Supplementary Fig. 12b, c), suggesting that 20E playsakeyroleinguidingcellfatetransitionswithinthewingdisc.

To capture the rapid and transient transcriptional responses inducedby20E—beyondtheresolutionofnaturaldevelopment—weestablishedaninvitro20Estimulationsystemandperformedtime-resolved snRNA-seq across nine time points, ranging from 0 to 360minutes (mins) after treatment (Fig. 3a and Supplementary Data8).WingdiscsweredissectedattheL5D2stage,whichischaracterizedbylowendogenous20Elevels49,providingacleanbaselineforhormonalresponse.A20Econcentrationof2µg/mLwasselectedbasedondose–responseoptimization,asitrobustlyinducedcanonicaltargetgenessuchasEcRandUSP(SupplementaryFig.13a).Afterqualitycontrol,atotalof125,717high-qualitycellswereretainedwithamediangenenumberof1304.UsingscRNA-seqdatafromnaturallydevelopingwingdiscsasareference,weidentifiedsixmajorcelltypes in the 20E-treated samples: Wm, Cuticle 1, Cuticle 2, Epithelial 1, Immune,andMatrixcells(Fig.3bandSupplementaryFig.13b,c).Gene expressionexhibitedaclearcascaderesponsefollowing20Eexposure. Withinthefirst20min,genesinvolvedinlarvalcuticledevelopment were significantly upregulated. At 40min, genes related to lipid responsebecameactivated.Between50and60min,genesinvolvedinepithelialcellandhemocytedifferentiationwerehighlyexpressed.By 4–6hours (h), genes related to cytoskeleton organization were stronglyupregulated(Fig.3candSupplementaryData9).Allcelltypesrespondedto20Eatdifferentlevelsduringthisperiod.Intheearly phase (0–30min), Wm, Cuticle, and Epithelial cells showed the strongest transcriptional responses. At later stages, the response shifted toward Ma and Im cells. This pattern indicates a time-dependentregulatorysequence(SupplementaryFig.13d).

To assess how 20E influences intercellular communication, we analyzedthedynamicsofsignalingpathwayactivityamongcelltypes. Notably,theTENASCINpathwayshowedasharpincreaseinsignaling strength at 40min. In contrast, the BMP, COLLAGEN, Notch, PARs, FGF,andSEMA4pathwaysbecamestronglyactivatedatlaterstages, especiallyat4and6haftertreatment(Fig.3d).Amongthese,theFGF pathwaywasactivatedearliest,withstrongactivityobservedinWm cells. BMP and Notch signaling peaked around 30min and stayed active,althoughtheirlevelsgraduallydiminishovertime.By6h,most pathwaysreachedtheirhighestactivitylevels.Theseresultssuggest that 20E triggers a stepwise and coordinated activation of multiple signalingpathways(SupplementaryFig.14a).Before40min,Wmcells were the main source of outgoing signals, particularly through the BMP,Notch,PARs,COLLAGEN,andFGFpathways.After40min,Im cellsbecamethedominantsignalingsource.Thisindicatesashiftin intercellular communication patterns between different cell types duringthe20Eresponse(SupplementaryFig.14b).

Tofurtherinvestigatethegeneexpressionprogramsregulatedby 20E, we identified two gene expression modules: Module 1 (M1) and Module 2 (M2) (Supplementary Fig. 15a). Genes in M1 was highly expressed in Wm, epithelial, and cuticle cells during the early phase (0–30min). It was enriched in genes related to morphogenesis and cuticle development. This indicates that Wm differentiation begins before 30min and is driven by morphogenetic signaling pathways (Fig.3e,fandSupplementaryFig.15b).Incontrast,genesinM2were mainlyexpressedinMaandImcellsduringthelaterstages(40–360min). (Fig. 3e, f and Supplementary Fig. 15c). This shift in module activity occurred alongside a reduction in morphogenetic cell types and an increaseinMaandImcell(SupplementaryFig.13c).Theseresultsindicate that20Econtrolsastage-specificswitchingeneexpressionprograms.

Toidentifykeyregulatoryfactors,weusedSCENICtoanalyzed TFswhoseactivitychangedovertimeinresponseto20E.BasedonRSS values,weselectedthetopsixkeyTFs,whichshowspecificenrichment in different cell types (Supplementary Fig. 16a). In the early stage (before30min),Hr38(LOC693064),Nubbin(LOC101744796),andSix2 (LOC101735990)showedstronglyactivity.TheseTFsregulatedgenes involvedinWm,Epithelial1,andCuticle1cells(Fig.3gandSupplementaryFig.16b–d).After30min,theactivityofCREB(LOC692871), FoxN3 (LOC101739133), and GATA-A (LOC10174013) were involved in ImmuneandMatrixcells.Thesefactorstargetedgenesrelatedtocell migration and immune responses (Fig. 3g and Supplementary Fig.16e–g).Theseresultsrevealatime-dependentshiftintranscriptionalregulationcontrolledby20E.

Compared with gene expression profiles of wing discs under natural development (scRNA-seq dataset), 20E-treated wing discs (snRNA-seqdataset)exhibitedadualtranscriptionalsignature,showing both L5D2-like baseline features and characteristics of later developmentaltimepoints(Fig.3h).Withincreasingdurationsof20E exposure, the induced transcriptional profile progressively transitioned from L5D4 (before 30min) to L5D5 (30min), and further to W1stage(50minto6h),suggestingthatthe30-mintimepointmay represent a critical tipping point for 20E-induced transcriptional reprogramming(Fig.3h).Pseudotimetrajectoryanalysisrevealedthat Wm cells from the treated group closely followed a trajectory like natural differentiation (Supplementary Fig. 17a, b), gradually transitioned into epithelial and cuticle lineages. Moreover, the key Wm regulators showed a distinct shift in expression before and after 30min in both datasets (Supplementary Fig. 17c). Together, these findings suggest that 20E induces natural developmental-like transcriptionalprogramsbutinafaster,time-compressedmanner.

Althoughthemajorcelltypesremainedthesameinbothconditions,theirrelativeabundancesandactivationtimingsweredifferent. Forinstance,theproportionsofWm,Cuticle,andEpitheliallineages graduallychangedfromL5D1toW1innaturaldevelopment.Incontrast,20Etreatmentcausedarapidandshort-termexpansionofthese lineageswithinthefirst30min.After40min,MatrixandImmunecells became dominant. Matrix cells showed a continuous increase over timeinbothgroups(SupplementaryFig.18a).GOandKEGGenrichment analyses revealed that genes enriched in endocytosis, mTOR signaling pathway, structural constituent of cuticle, and DNA-binding transcriptionfactorbindingwerespecificallyactivatedduringnatural development, reflecting sustained biosynthetic and developmental activityduringnaturalwingdiscmaturation.Incontrast,genesenriched inribosomebinding, mRNA binding, and oxidative phosphorylation were activated upon 20E stimulation, indicating a stronger emphasisonproteinsynthesis,energymetabolism,andtranscriptional reprogramming (Supplementary Fig. 18b, c). A set of 273 shared upregulatedgenesbetween20E-treatedandnaturallydevelopingwing disc were identified as candidate direct targets of 20E stimulation (SupplementaryFig.19aandSupplementaryData10).GOandKEGG enrichmentanalysisofthesegenesrevealedsignificantinvolvementin developmentalprocessessuchascellularandorganmorphogenesis, andgrowth,aswellasinstressandimmuneresponses.Notably,severalpathwaysassociatedwithproteinhomeostasisandsignaltransduction—such as autophagy, MAPK signaling, and ER protein processing—werealsoenriched,suggestingthat20Emaycoordinate bothstructuralremodelingandcellularstressresponsesduringwing discdevelopment(SupplementaryFig.19b).Wefurtheranalyzedthe dynamicexpressionpatternsofthesegenesacrossnaturaldevelopmentand20E-treated(SupplementaryFig.19c).Thesegenesshoweda progressive,stage-specificactivationfromL5D1toP6undernatural development,transitioningfromclustersC1toC10.Upon20Etreatment,genesinclustersC3–C6respondedrapidly.Incontrast,asubset ofgenesinclustersC8–C10wasactivatedwithinthefirst30min,while the remaining genes in these clusters showed gradual activation thereafter.Together,theseresultsprovideahigh-resolutionviewof 20E-inducedtranscriptionalreprogrammingandunderscoreitsrolein acceleratingdevelopmentaltransitionsinthewingdisc.

Genetransitionmodeldefineswingdiscdevelopmentalstages
To systematically characterize the transcriptional dynamics and developmentaldrivenregulatedby20Einthewingdisc,bycombing geneexpressionpatternssimilarity(Fig.4a),local20Econcentrations (Fig. 4b, Supplementary Fig. 20a and Supplementary Data 11), and morphological changes (Supplementary Fig. 20b–e), we propose a five-stage gene transition model of wing disc development: Stage 1 (L5D1–L5D2), Stage 2 (L5D3–L5D6), Stage 3 (L5D7), Stage 4 (W1–W2), and Stage 5 (P6). DEGs analysis showed distinct gene expression peaks ineachstage(Fig.4c,SupplementaryFig.21a–c,andSupplementary Data 12). This model reflected the molecular shifts underlying each developmentalphase.

Specifically, in Stage 1, 20E concentrations were low (SupplementaryFig.20a),andwingdiscmorphologyremainedunexpanded (SupplementaryFig.20b).L5D1andL5D2exhibitedstrongtranscriptional similarity, indicating closely coordinated gene expression in early stage. This period featured co-activation of the Toll/Imd and Hippopathways,initiatingearlymorphogenesis,alongsidechromatin remodeling that reflected enhanced transcriptional plasticity (SupplementaryFig.21d).As20ElevelsbegantoriseinStage2,wingdisc cells became more structure, and the tissue displayed increasing organization, with visible formation of distinct layers and regions (SupplementaryFig.20a,c).Transcriptionalsimilarityremainedhigh from L5D3 to L5D6, suggesting sustained gene activity. Pathways involved in homologous recombination and DNA replication were stronglyenriched(SupplementaryFig.21e),likelytomaintaingenome stabilityduringrapidcellproliferation.Stage3representedamajor transition point. L5D7 displayed strong transcriptional separation, coincidingwithpeak20Elevelsandwingdiscenteredaremodeling phase(SupplementaryFig.20a,d).MorphogenesisandMAPKsignalingpathwayswerelikelyinvolvedtotissueboundaryreorganization (Supplementary Fig. 21f). In Stage 4, W1 and W2 formed separate transcriptional clusters, showing stage-specific transcriptome shifts. 20Elevelsremainedelevated,andpathwayssuchasoxidativephosphorylation and ribosome biogenesis were activated, providing the energyandproteinsynthesisrequiredfortissuematuration(Supplementary Fig. 21g). Morphologically, the wing disc showed defined regionalstructure(SupplementaryFig.20e).ByStage5,P6displayed strongtranscriptionalseparation,suggestingamajorgeneexpression shift.20Elevelswerehighestatthisstage,andcoordinatedactivation of Wnt and MAPK signaling helped refine tissue architecture, while longevity-associated pathways contributed to maintaining tissue homeostasis (Supplementary Fig. 21h). Cuticle cells became the dominant population (Fig. 4d), consistent with high expression of cuticle-related genes such as CPH32, CPH42, CPG46, and CPG39 (Fig.1c).Thesegeneslikelysupportthefinalformationoftheadult wingcuticle.Thewingdiscmorphologyreacheditsfinalform,withall major structures and cell types fully differentiated and functionally mature. GO and KEGG enrichment analyses of stage-specific genes revealedacascadeofsignalingpathwaysthatweresequentiallyactivatedduringdevelopmentaltransitions.

Spatialmappingofstage-specificDEGsrevealedtheirenrichment in distinct regions, indicating functional compartmentalization (Fig.4e).Forexample,geneexpressionwasconcentratedinthewing budatStage1,shiftedtotheepidermisatStage2,andreturnedtothe wingbudatStage4.Thesespatialdynamicsofgeneexpressionclosely mirrored tissue morphology, highlighting the tight coordination betweengeneactivityandstructuraldevelopment(Fig.4c,e).Meanwhile,theexpressionof20E-responsivegenes,suchasUSPandE75, also showed dynamic and stage-specific expression peaks (Fig. 4f). Notably,the20EreceptorgeneUSPwassubstantiallyupregulatedin Stage3,when20Elevelsreachedtheirpeak(Fig.4b),suggestingitskey roleintriggering20E-drivenmetamorphicremodeling.

To dissect how cell–cell interactions (CCIs) regulate wing disc development,weperformedCCIanalysisacrossstages.Overall,the interactionnetworkgrewmorecomplexfromStage1to4andthen simplified at Stage 5, suggesting a shift from active remodeling to stabilization(Fig.4l).Sevensignalingpathways—Wnt,Notch,Hh,BMP, FGF, PARs, and SEMA4—showed stage-specific activation (Fig. 4g–k andSupplementaryFig.22).DuringStage1,intercellularsignalingwas mainly concentrated in Wm, Im, and Ma cells. Notch signaling was prominently initiated byMacells and received by Wm, Ep1, and C1 cells, with key ligand–receptor pairs, including DLK1–NOTCH3 (Ma-secreted)andJAG1–NOTCH3(Wm-secreted)(Fig.4g,mandSupplementaryFig.22a).TheseinteractionsindicatethatWmcellsservedas primaryrecipientsofmorphogeneticcuesduringearlypatterning.The preferential Notch reception by structural cell types (Wm, Ep1, C1) highlightedNotch’scentralroleinestablishingspatialorganizationat thisstage.

InStage2,Imcellsemergedasmajorsignalinghubsbyintegrating multiplepathwayssuchasBMP,Hh,Notch,SEMA4,andFGF.Meanwhile, PARsignalingwaspredominantlyinitiatedbyWm,Ep1,andImcells,with Macellsasprincipaltargets.Notably,ligand–receptorinteractionslike PRSS1–PARD3(fromWmorImtoMa)suggestedaconservedsignaling axispotentiallyinvolvedinimmune–metaboliccoordination(Fig.4h,n andSupplementaryFig.22b).InStage3,Notchsignalingbecamemore diverse,withWm,Ep1/Ep2,C1/2,andImcellsactingassources.Wmcells retainingdualrolesasbothprimarysendersandreceivers.Wntsignals mainlytargetedC1,Wm,andEp1,whileFGFligandsfromMa,Im,and WmcellsactedonEp1andC1/2.AkeyFGF1–FGFR2axis(WmtoEp1)was identified, potentially driving epithelial morphogenesis at this stage (Fig.4i,oandSupplementaryFig.22c).

AtStage4,WmandEp1/Ep2cellsdominatedthereceptionand transmission of BMP, Notch, SEMA4, and Wnt signals. Notably, Ep2 cellsstoodoutbysendingWNT2andSEMA4GligandstoWmandEp1 through WNT2–(FZD2+LRP5) and SEMA4G–PLXNB2 pairs. These interactions likely promote epithelial layering and boundary formation.(Fig.4j,pandSupplementaryFig.22d).ByStage5,C1cellsmainly receivedWntsignaling,whileMacellsreceivedHhsignals.BMPsignals from Ad cells targeted C2, and C3 showed strong intra-type BMP activity via GDF5–(BMPR1B+ACVR2B). PAR signaling was largely reducedatthisstage.(Fig.4k,qandSupplementaryFig.22e).Notably,MacellscontinuedtoreceivePARsignalsacrossallstages.BothMa and Im cells remained active participants in Hh signaling and had extensivecontactswithothercells.ThesefindingsindicatethatMaand Imcellsserveaspersistentsignalinghubsthroughoutdevelopment (SupplementaryFig.22).Together,thesefindingsrevealadynamicand lineage-dependent signaling framework that ensures stage-specific coordinationofsilkwormwingdiscmorphogenesis.

Tosummarizethesecoordinatedchanges,weproposetheGene Transition Model (GTM),whichdivides wing discdevelopmentinto fivestages:(1)Developmentalblueprint(L5D1–L5D2),characterizedby earlymorphogeneticsignalingandhightranscriptionalplasticity;(2) Cellularfoundation(L5D3–L5D6),markedbycellgrowthandgenomic stabilization;(3) Remodeling and sculpting (L5D7), defined by peak 20Eandstructuralboundaryreorganization;(4)Structuralformation (W1–W2),characterizedbytissuelayeringandincreasedmetabolism; and(5)Maturationandstability(P6),wherefinaldifferentiationand homeostasis shape the adult wing (Fig. 5). This model outlines a hormone-gradient-driventemporalaxis,spatialdecodingviaregion-specificactivationof20E-responsivegenes,andmodularcoordination through dynamic signaling networks such as Wnt, Notch, Hh, and MAPK.Together,thesefindingsprovideaunifiedframeworklinking hormonal regulation, gene expression, and tissue morphogenesis duringinsectorgandevelopment.

Methods
Experimentalanimalpreparationandsamplecollection
TheB.moristrainDazaowasrearedundercontrolled conditionsat 27±1°C with 60% relative humidity under a natural photoperiod. Larvaewerefedfreshmulberryleavesuntilthespinningstage.Wing discs were collected from both larval and pupal stages. Larvae and pupae ofbothsexeswererandomlyselectedforallexperiments.S. frugiperdalarvaewererearedinplasticcontainersat26±2°C,66±5% relativehumidityundera16-hlight/8-hdarkcycle.Theywerefedan artificialdietuntilpupation.Afterpupation,individualsweremovedto separate containers and maintained until adult emergence. All silk-worm and S. frugiperda were provided by the Biological Science ResearchCenteratSouthwestUniversity,China.

Fortissuecollection,gentlepressurewasappliedtotheregion belowthespiraclesofthesecondandthirdthoracicsegmentstomake thewingdiscareaslightlyprotrude.Then,theepidermiswaspierced using the tipof forceps (Dumont, 0209-5-PO,Switzerland), and the wing disc was carefully extruded. Each disc was immediately transferredintophosphate-bufferedsaline(PBS)at4°C.Underaste-reo microscope (Leica, dm2700 p, Germany), visible fat body, tracheae,andotherunwantedtissueswerecarefullyremoved.Theentire procedurewascompletedwithin5mintopreservetissueintegrity.

Microscopicobservationofwingdisc
Wingdiscswereobservedusingastereomicroscope(Olympus,Bx53, Japan).Eachdiscwascarefullymountedonaglassslidetomaintain stabilityandensureclearimaging.High-resolutionimageswerecapturedandanalyzedusingPhotoshop(v22.4.3,Adobe,USA).Themajor (long)andminor(short)axesofeachdiscweremeasuredtoassesssize andmorphology.

WingdiscsamplecollectionandscRNA-seqsequencing
WingdiscsampleswerecollectedfromB.moriacrossdevelopmental timepointsfromL5D1toL5D7,aswellasW1,W2,andP6,toensure comprehensivedataforscRNA-seqanalysis.Duetothesmallersizeof the wing discs, 200, 100, and 50 independent wing discs were collectedfortheL5D1–L5D7,W-1/W-2,andP6stages,respectively.The wingdiscsamplescollectedfromeachtimepointweredisassociated into single-cell suspensions using PBS (BBI, Cat# E607008-0500, China).Beforesequencing,cellcountingandviabilityassessmentwere performed, ensuring a cellviability rate of ≥80%.The prepared cell suspensions were processed using a DNBelab C Series High-Throughput Single Cell RNA Library Prep Kit V2.0 (MGI, Cat# 940-000519-00) according. Briefly, magnetic beads containing Oligonucleotidesequenceswereencapsulatedtogetherwithwingdisccellsin droplets,formingsingle-cellgelbeadsinemulsions.Withinthedroplets,cellsunderwentlysistoreleasemRNA,whichwassubsequently capturedbythecellbarcodesequencesonthebeads.Afterdroplet breaking, reverse transcription was performed to synthesize cDNA. SequencinglibrairieswereconstructedforthebarcodedcDNAmoleculesandsequencingplatform(MGI,MGISEQ2000,China).

scRNA-seqdatapre-processing
TherawFASTQsequencingdatawerealignedtotheNCBIB.mori reference genome (GCF_014905235.1_Bmori_2016v1.0) using the DNB-C4 pipeline (https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_HT_scRNA-analysis-software, version3.1.5). The alignmentandqualitycontrolstepsresultedinthegenerationofa gene expression matrix. Subsequent filtering was conducted for each sample using the Seurat package (v4.3.0)51 based on the followingcriteria:
1) Potential doublets identified using the DoubletFinder (version2.0.3)52algorithmwereremoved;2)cellswithreasonableUMI counts (600≤nUMI<max nUMI×0.95), expressed genes (nGene≥ 500wereretained.3)Cellswithlog10(genes)perlog10(nUMI)≤0.8 wereexcluded.Inaddition,mitochondrialgeneswereexcludedfrom downtream analysis. Dimensionality reduction and clustering were performedusingtheScanpypackage(v1.9.3)53.Highlyvariablegenes (HVGs) were identified, retaining the top 2500 genes. Batch effects werecorrectedusingtheHarmony54algorithm.Principalcomponent analysis (PCA) was applied, retaining the top 30 PCs, followed by neighbor graph construction with n_neighborssetto 30. Clustering wasconductedusingtheLeidenalgorithmwitharesolutionof0.4.
DEGswereidentifiedusingthet-testapproach,andthetop200DEGs wereselectedforGOenrichmentanalysis.Celltypeannotationwas performed based on the functional properties of top markers andDEGs.

WingdisctissuecollectionandStereo-seqsequencing
Spatial transcriptomic profiling of wing disc tissue was performed using Stereo-seq (Spatial Enhanced Resolution Omics Sequencing) technology55. Briefly, multiple wing disc tissues from B. mori at differentdevelopmentalstages(L5D2–L5D3andW2)wereembedded in OCT compound (Sakura, Cat# 4583, Japan), sectioned, and subjectedtomorphologicalqualitycontrolandRNAintegrityassessment. The sections were then mounted onto spatiotemporal FC chips for permeabilization. Fluorescent imaging at various time intervals was performedtodeterminetheoptimalpermeabilizationtime.
After identifying the optimal permeabilization time, sections with larger exposed areas of the wing disc were selected and mounted onto spatialSCchips. Themounted sectionswerefixed and stained with Qubit™ ssDNA Assay Kit (Thermo Fisher, Cat# Q10212, USA) for fluorescence imaging, which recorded tissue morphology and position. The tissues were permeabilized in an incubator at 37°C for the previously determined optimal time to release RNA, followed by reverse transcription at 42°C for 3h to generatecDNA.AftertissueremovalandRNArelease,thesamples underwent magnetic bead purification and cDNA amplification. Next, the amplified cDNA products were fragmented and ligated with PCR adapters for further amplification and size selection. Subsequently,DNAnanoballswerepreparedandsequencedusing the DNBSEQ-T7RS platform. The spatial chips are equipped with captureprobesthatbindtomRNAmoleculesreleasedfromtissue cells.Thecapturedgeneinformationwasspatiallyresolvedbased ontheprobe’suniqueCID(CaptureID),whichcorrespondstothe originallocationwithinthewingdisc.

Stereo-seqdatapre-processing
RawStereo-seqFASTQdatawerealignedtotheNCBIB.morireference genome(GCF_014905235.1_Bmori_2016v1.0)usingtheSAWpipeline56 (https://github.com/STOmics/SAWversion6.0.0)forqualitycontrol. Downstream analysis was performed at a bin20 resolution. Quality controlofthedatawasconductedusingtheScanpypackage(v1.9.3)53. Cells were filtered using the criteria of min_cells=20 and min_-genes=max (genes)×0.05. Cell types annotated in the scRNA-seq datasetwereusedasreferences.Spatialannotationwasachievedusing theCell2location(v0.1.3)57algorithm.

20-hydroxyecdysonegeneset
To identify the 20E hormone-related genes in silkworm wing discs (WDs),weconducted a comprehensive reviewofrelevant literature availableontheWebofScience.Basedontheirfunctionsandroles within signaling pathways, these genes were classified into several categories: 20E Receptor, 20E Signaling Cascade Proteins, Signal Cascade-Related Enzymes, Immune-Related Factors, Functional Proteins, Protein Receptors, Transcription Factors (TFs), and Cuticular Proteins.AdetailedlistofthesegenesisprovidedinSupplementary Data 4. This classification offers a systematic framework for understanding the roles of these genes in WD development and their involvementinthe20Esignalingnetwork.

Trajectoryanalysisandpseudotimeanalyses
Cellular trajectories were constructed using the Tome pipeline (https://github.com/ChengxiangQiu/tome_code)40.PseudotimeanalysiswasperformedinMonocle2(v2.28.0)58usingtheDDRTreemethod. Tofocusonspecificcellpopulations,expressionmatrixesforepithelial,cuticle,andWm-relatedcellswereextracted.DEGswereidentified usingtheFindAllMarkersfunction,withsignificantgenesfilteredusing anadjustedp-valuethresholdofp_val_adj<0.05.PseudotimevisualizationwasconductedusingCellRank2(v2.0.0)59.

Pathwayenrichmentanalysisofgeneset
UsingtheRpackage“ClusterProfiler”(v4.6.2,Bioconductor,China)60 toperformGOenrichmentandKEGGpathwayenrichmentanalyses, with a p-value cutoff set at 0.05. Data visualization was primarily conductedusingtheggplot2(v3.5.1)package.

RNAinterferenceexperiment
Double-stranded RNA (dsRNA) targeting specific genes was synthesizedusingtheRibomaxT7RNAProductionSystem(Promega,Cat# P1300, USA) according to established protocols61,62 and stored at –80°Cforsubsequentexperiments.Double-strandedprimersequenceswereprovidedinSupplementaryData6.Size-matchedfifth-instar B.morilarvae,attheearlystageofday2,wererandomlydividedinto experimentalandcontrolgroups(threebiologicalreplicates,with30 independentwingdiscperreplicate).Eachlarvaintheexperimental groupwas injectedwith10μgdsRNAintothe third abdominalsegmentusingafinecapillaryneedle,whilecontrollarvaewereinjected withanequivalentamountofdsRNAtargetingEGFP.After48h,five larvaefromeachgroupweredissectedtoassesswingdiscmorphology andevaluateRNAiefficiencyviaRT-qPCR.Theremaininglarvaewere rearedtothepupalstage,duringwhichwingdevelopmentwasmonitored. The incidence of wing disc segmentation defects and adult wing abnormalities was statisticallycompared between groups.The sameprocedurewasappliedtoassesstheroleofRfxinS.frugiperda. Fifth-instar larvae with uniform growth conditions (three biological replicates,with20independentwingdiscperreplicate)wereinjected with dsRNA targeting SfRfx. Wing discs were dissected 48h post-injection to examine morphological changes and determine knockdown efficiency using RT-qPCR. Surviving larvae were allowed to pupate, and adult wing morphology was analyzed and compared statistically.

Rfxconservationanalysis
ThegeneandproteinsequencesofRfxwereobtainedfromSilkDB3.0 (https://silkdb.bioinfotoolkits.net)andNCBI.Phylogenetictreeswere builtinMEGA8.0.Conserveddomainanalysiswasperformedusing the CDD website (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/ bwrpsb.cgi), and motif analysis was conducted online using the MEMEwebsite(https://meme-suite.org/meme/tools/meme).
Using blastp to perform homology alignment between human (GCF_000001405.40)andB.moritoobtainhomologousgenecorrespondences.B.morigenesweremappedtotheirhumanhomologsto enable downstream analysis. CCIs were inferred using the CellChat package(v1.6.1)65.Foreachdevelopmentalstage,aseparateCellChat objectwascreated,andCCInetworkswereindependentlycalculated. These stage-specific CellChat objects were subsequently merged to facilitatecomparativeanalysisacrossstages.

Invitrocultureofwingdiscsand20Etreatment
Toexaminethetime-dependenttranscriptionalresponseofwingdiscs to20E,weestablishedaninvitroculturesystemandapplied20Eat different time points. Wing discs were dissected using forceps (Dumont, 0209-5-PO,Switzerland)andwashedinPBStoremovesurrounding tissues. Cleaned wing discs were transferred into Grace’s medium(Gibco,Cat#11300-043,USA)supplementedwith10%fetal bovineserum(FBS)(Gibco,Cat#A5669701,USA)and20E(Yuanye, Cat#S25531,China).
Todeterminetheoptimal20Econcentration,wefirstperformeda dose–response test (see Supplementary Fig. S12). Wing discs were culturedinmediumcontaining0,0.1,0.2,0.5,1,2,5,or10μg/mLof 20E,eachwith10%FBS.Wemeasuredtheexpressionlevelsofknown 20E-responsivegenes,includingEcRandUSP.Theresultsshowedthat 2μg/mLof20E-inducedsignificantlyincreasedtheexpressionofthese genes. Previous studies also reported that this concentration effectivelyactivated20Eresponsesinculturedwingdiscs.
Silkwormshaveanopencirculatorysystem,soendogenous20Ein hemolymphcaninfluencemultipletissues.However,theefficiencyof 20E uptake into isolated wing discs remains uncertain. To ensure strongandconsistentstimulation,wechose2μg/mLof20Etoinducea robusttranscriptionalresponse.Thisallowedustocapturetheactivationtimingof20Etargetgenesandcompareitwithnaturaldevelopment.Therefore,weusedthisconcentrationintheformalexperiment.
Intheformalexperiment,wingdiscswereculturedincomplete mediumcontaining2μg/mLof20Efor10,20,30,40,50,or60min, andfor4or6h.Ablankcontrolgroup(immediatelyfrozenwithout culture) was also included. After incubation, we transferred each sampleintoasterilemicrocentrifugetube.Weremovedthemedium, washed the discs twice with PBS, blotted off excess liquid, and immediatelyfrozethesamplesinliquidnitrogen.Thesesampleswere usedforsnRNAsequencing.

WingdisctissuecollectionandsnRNA-seqsequencing
Thefresh-frozenwingdisctissuesweregroundandlysedtoprepare single-nucleussuspensions.Nucleiwereextractedandisolatedfrom fresh-frozen wing disc tissue using PBS (BBI, Cat# E607008-0500, China)andresuspended.Nucleisuspensionsweresubjectedtoquality control,onlythosewithcellviability>85%andcellcount>100thousand were retained for downstream experiments. Nuclei were processedusingtheDNBelabCSeriesSingle-CellLibraryPrepSet(MGI, 1000021082,China)63.Single-nucleusRNAlibrarieswereconstructed usingtheDNBSEQ-T7RSHigh-ThroughputSequencingKitV3.0(MGI, 940-000269-00, China) and sequenced on the MGI DNBSEQ-T7RS sequencingplatform.

snRNA-seqdatapre-processing
TherawsequencingdatawerealignedtotheNCBIB.morireference genome (GCF_014905235.1_Bmori_2016v1.0) using the DNB-C4 pipeline (https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_HT_scRNA-analysis-software) for quality control and alignment, resulting in the generation of a gene expression matrix. The gene expressionmatrixwasfilteredusingtheSeuratpackage(v4.3.0)using theaforementionedcriteria.Dimensionalityreductionandclustering wereperformedusingScanpy(v1.9.3)HVGswereidentified,retaining thetop3000genes.BatcheffectswerecorrectedusingtheHarmony algorithm, and PCA was performed, retaining 30 PCs. A neighbor graphwasconstructedwithneighborssetto10,andclusteringwas performedusingtheLeidenalgorithmwitharesolutionof0.4.Cell type annotation was conducted using SingleR64, with scRNA-seq annotatedcelltypesservingasreferences.

Inferringcell–cellcommunication
Usingblastptoperformhomologyalignmentbetweenhuman (GCF_000001405.40)andB.moritoobtainhomologousgenecorrespondences.B.morigenesweremappedtotheirhumanhomologsto enable downstream analysis. CCIs were inferred using the CellChat package(v1.6.1)65.Foreachdevelopmentalstage,aseparateCellChat objectwascreated,andCCInetworkswereindependentlycalculated. These stage-specific CellChat objects were subsequently merged to facilitatecomparativeanalysisacrossstages.

Transcriptionfactoranalysis
UsingblastptoperformhomologyalignmentbetweenD.melanogaster (GCF_000001215.4) and bombyx to obtain homologous gene correspondences. B. mori genes were mapped to their D.melanogasterhomologstoenableTFregulonanalysis.TheSCENIC package(v1.3.1)66withdefaultsettingsusedtoinferstatisticallyactive TFsandtheirtargets.
UsedwiththeD.melanogasterTFdatabasetoidentifyTFregulons andtheirassociateddownstreamtargetgenes.TopTFswerefiltered basedonanRSS>0.01andazThresholdof2.TheresultingTF–target regulatorynetworkswerevisualizedusingCytoscape(v3.9.1)67.

Genenetworkconstruction
Single-cell weighted gene co-expression network analysis (WGCNA) was performed using the hdWGCNA package (v0.2.26)68. Stage-specifichubgeneco-expressionnetworkswereidentifiedandvisualizedusing Cytoscape (v3.9.1) to explore the dynamic regulatory relationshipsacrossdevelopmentalstages.

Similarityanalysis
TheHarmonyalgorithmwasappliedtointegratescRNA-andsnRNA-seqdatasets.Toevaluatethecorrelationbetweencellsfromdifferent time points in the scRNA- and snRNA-seq datasets, the Lorentzian distancemetricfromtheRpackage“philentropy”(v0.8.0)69wasused tocalculatepairwisedistancesbetweencells.

Quantificationof20Econtentineachwingdisc
Toaccuratelymeasure20Elevelsduringwingdiscdevelopment,we directlyquantifiedthehormoneconcentrationinthewingdisctissueinsteadofthehemolymph.Thistissue-specificmethodavoided thedilutioneffectcausedbythesilkworm’sopencirculatorysystem. Wingdiscswerecollectedatdifferenttimepoints.Threebiological replicates were used for each timepoint, with 30, 10, and 5 independentwingdiscsperreplicatefortheL5D1–L5D7,W-1/W-2,and P1/P3/P5/P7/M stages, respectively. The concentration of 20E was measured using an ELISA kit (Boshen Biotechnology, Cat# BS-E19221O2,China).Foreachsample,wingdiscswerehomogenizedin 100μL of PBS. The homogenate was centrifuged at 2000×g for 10min, and the supernatant was collected for analysis. The ELISA platewasequilibratedat25°Cfor20minbeforeuse.Standardswere prepared, and 50μL of each standard solution was added to the correspondingwells.Forsamplewells,10μLofthesamplesupernatantand40μLofsamplediluentwereadded.Blankwellswereleft empty.Then,100μLofHRP-conjugatedantibodywasaddedtoeach well,excepttheblanks.Theplatewasincubatedat37°Cfor60min. Afterincubation,thewellswerewashedfivetimes.Next,50μLeach ofsubstratesolutionsAandBwereaddedtoallwells.Theplatewas incubatedinthedarkat37°Cfor15min.Thereactionwasstopped byadding50μLofstopsolutiontoeachwell.Opticaldensity(OD) values were measured within 15min at 450nm with a reference wavelength of 630nm. A standard curve was generated in Excel (MicrosoftCorporation,Redmond,WA,USA)byplottingODvalues against the standard concentration. The curve followed the equation:Y=−0.001086×2+0.1162X+0.07753,wasusedtocalculatethe 20Econcentrationineachsample.Thetotal20Econtentperwing discwascalculatedbydividingthetotal20Econtentbythenumber ofdiscsineachsample.

Drosophilamelanogasterwingdiscsingle-cellRNAsequencing (scRNA-seq)dataanalysis
SequencingdataandalignedmatrixesweredownloadedfromtheGEO database(accessioncodeGSE155543).Thestandardanalysispipeline wasthesameasthatusedinthearticlemethods11.Clustersweresplit intoAMPandepithelialcells.Finally,weobtained12,872cellsand6552 homologousfeatureswiththesilkworm.Correlationswerecalculated usingSpearman’smethod.

Statisticsandreproducibility
Each statistical test used to analyze data is indicated in the figure legendswherepvalueswerecalculated.Theanimalexperimentdatain thefigurearepresentedasmean±SD.Two-tailedunpairedStudent’st testwasusedforstatisticalanalysisofthedifferencesbetweenthetwo groups. Data in Fig. 2g, h were analyzed by one-way ANOVA. ***p<0.001,**p<0.01,and*p<0.05wasacceptedasstatisticallysignificant. All analysis assumed both populations were normally distributedwhenparametrictestswereused.Aconfidenceintervalof95% wasused.Specificquantificationmethodsandpvaluesareindicatedin thefigure’slegends.
To identify differentially expressed genes (DEGs) between cell typeswithineachdevelopmentalstage,weperformedcluster-based differential expression analysis. Specifically, for the subset of cells fromeachdevelopmentalstage,andweemployedtheFindAllMarkers function from the Seurat package (v4.3.0). This function performs pairwise comparisons between each cell cluster and all other cells withinthatstagetoidentifyitsspecificallyupregulatedgenemarkers. Thekeyparametersusedintheanalysiswereasfollows:test.use=“ MAST”,only.pos=TRUE.Thissettingretainsonlygenesthataresig-nificantly upregulated in the target cluster. min.pct=0.01 and min.cells.feature=3.Agenewasincludedintheanalysisonlyifitwas detectedinatleast1%ofcellswithintheclusterofinterestandinat least3cellsofthatcluster.logfc.threshold=0.1.Onlygeneswithan average log fold change (base 2) greater than 0.1 were retained. return.thresh=0.01. Only genes with an adjusted p-value (p_val_adj) less than 0.01 were considered statistically significant DEGs. The functionautomaticallycorrectsformultiplehypothesistestingusing theBonferronimethod,adjustingallrawp-values(p_val)basedonthe total number of genes tested across the entire dataset, yielding adjustedp-values(p_val_adj).Foreachdevelopmentalstage,theanalysisgeneratedadataframelistingthepositivemarkergenesforeach cellclusteroreachtimepoint,alongwithmetricssuchastheaverage logexpression,thepercentageofcellsexpressingthegene,p-values, andadjustedp-values.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Liu etal. present a comprehensive multi-omics study that integrates scRNA-seq, snRNA-seq, and spatial transcriptomics to\nconstruct a high-resolution spatiotemporal atlas of the Bombyx mori wing disc across multiple developmental stages. Based\non this resource, the authors propose a five-stage Gene Transition Model (GTM) regulated by 20-hydroxyecdysone (20E) to\ndescribe the temporal dynamics of cell fate reprogramming. This is an interesting and timely study that offers novel\ninsights into hormonal regulation during insect organogenesis and provides a theoretical framework for understanding genetic\ncontrol and potential molecular intervention strategies in Lepidoptera. In addition to its biological relevance, the dataset\nrepresents a valuable resource for future research into wing disc development and\nhormonal regulation. Overall the study is well performed and reads clearly. Prior to publications, I’d like to see some clarifications and revisions.\nMain concerns:\n1. In Fig. 3h, the comparison between naturally developing and 20E-induced samples could be expanded. It would be\nvaluable to clarify the similarities and differences in cell composition and pathway activity between the two conditions,\nparticularly considering how closely 20E treatment mimics natural developmental processes.\n2. For Fig. S17, please provide functional enrichment results for the shared upregulated gene set between the natural and\n20E-induced groups. This would help clarify the core biological processes directly regulated by 20E.\n3. In Fig. 4h–i, the analysis of cell–cell interactions (CCI) does not clearly reflect stage-specific dynamics. The authors are\nencouraged to highlight temporal shifts in interaction patterns across developmental stages.\n4. The adult wing morphology following Rfx RNAi in Spodoptera frugiperda is not shown. Including representative images of\nadult wing phenotypes would strengthen the phenotypic validation. Given its importance in supporting the evolutionary\nconservation of Rfx, the S. frugiperda RNAi results should be moved to the main figures. Additionally, a comparative\nanalysis of Rfx conservation across insects would further support this claim.\n5. Beyond Rfx, the authors are encouraged to functionally validate additional key transcription factors identified through\nSCENIC analysis (e.g., Dll or Blimp-1) via RNAi or similar approaches.\n6. In Fig. 3g, it would be helpful to indicate which cell types express the highlighted transcription factors and how their\nactivity contributes to cell-type-specific transcriptional regulation and heterogeneity.\n7. The evolutionary comparison with Drosophila and B. anynana is relatively brief. A more detailed discussion of conserved\nand lineage-specific developmental strategies across these species would enhance the impact of the comparative analysis.\n8. While the methods section is generally well-detailed, additional clarification on the threshold settings and statistical\nparameters used in transcription factor enrichment, module identification, and pathway analysis would improve transparency\nand reproducibility.\nBelow are some minor comments for the authors to consider:\n9. In Fig. 2j, the figure legend should specify “dsRfx” to match the experimental description. The labels in Fig. 3a–f should be\nadjusted for better alignment, and the RSS values in Fig. 3j are too small to read clearly and should be enlarged.\n10. Several supplementary figures appear to lack sufficient resolution. Please ensure all images are high-quality and include\nappropriate legends and annotations for ease of interpretation.\n11. Reference formatting should be carefully reviewed. Inconsistencies in journal name abbreviations and missing page\nnumbers should be corrected for uniformity.\n12. In the author list, a comma should follow “Ying Lin” to maintain consistent punctuation.\n13. Line 27: The full name of 20E should be written as “20-hydroxyecdysone”.\n14. Line 90: “Fig. S1” should include a period to follow journal formatting guidelines.\n15. Lines 105–106: The font style of reference citations appears inconsistent (non-italic). Please standardize formatting\nacross the manuscript.\n16. Lines 180 and 185: The correct formatting should be “Fig. 2f and Fig. S7c”.", "reviewer_2": "The current manuscript outlines a spatial transcriptomic approach to map the development of wing discs in Bombyx mori.\nThe authors have compiled an impressive dataset, and undertaken extensive analysis, to identify the cell types within this\ntissue, their spatial location, the changes in cell fates over time, and changes in gene expression that accompany these\nchanges. In addition, they explore how the steroid hormone ecdysone affects wing disc development, demonstrating that it\nregulates the pace of wing disc development. We commend the authors on the scale of this work, and feel that it makes a\nvaluable contribution to our understanding of organ development. We have no experience with the transcriptomics\napproaches used and thus will limit our comments to other aspects of the work, which we think would help strengthen the\nmanuscript.\nThe acronyms for the different cell types are really hard to keep track of. Can you use the whole word to make it easier to\nfollow?\nThe experiments knocking down Rfx in the wing disc were very interesting and highlighted the role of this gene in Wing\nMorphogenesis cell type development. However, it wasn’t clear why this particular cell type was the focus?\nThis might sound a bit picky, but the order of the methods did not follow the order of experiments. This isn’t necessary, but\ndoes help when trying to find the relevant method to match to the experiment.\nFor the wing disc culture experiments, more detail is required in the methods to understand how these were performed. What\nmedium did you use for culture and did you add any insulin or bovine serum? Preferably, this would have it’s own section.\nFor the experimental design, why did you choose to change the time of exposure instead of 20E dose? Lastly, 2µg / ml of\n20E is a high dose for wing disc cultures. It would be useful to explain why you chose this dose.\nWhy were 20E titres collected from the wing discs and not from the hemolymph?\nSome of the paragraphs were in really big blocks of text. It would help with readability to break up these paragraphs into\ndiscrete ideas and to provide a bit more context throughout with regards to the aims of your different analyses.\nMinor comments:\n• Line 38-39: not all insects have wing discs (ametabolous and hemimetabolous insects have wing buds which are\nanalogous)\n• Line 48-50: Rephrase this sentence as it sounds like you are saying that Bicyclus anynana is not a lepidopteran.\n• Figure 2K legend: what do the arrowheads point to? Also, there are two sets of dsEGFP and two sets of dsRfx panels. How\nare these different?\n• Line 242: You mention Usp is expressed in Wm cells, but EcR isn’t in the list of ecdysone response genes?\n• In the text, you refer to 12 supplementary tables. These tables are found in the excel file. However, in the pdf there are only\n3 supplementary tables and the formatting of these tables appears to be off.\n• For Figure 3d (and the accompanying supplementary figure), it wasn’t clear what was meant by relative information flow.\n• Figure S18: What do BF and Trans refer to?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts."}, "author_response": "Response to reviewers\nWe are grateful for the time and thought the reviewers have given to\nassessment of our manuscript. We have made extensive changes including new\nexperiment results, figure parts, additional clarifications and methodology details\nto the text. In addition to tracked changes, we have highlighted modified text in\nyellow in the revised manuscript to help identification of them in the merged file.\nOur point-by-point responses are listed below:\nReviewer #1\nComment: 1. Liu et al. present a comprehensive multi-omics study that integrates\nscRNA-seq, snRNA-seq, and spatial transcriptomics to construct a high-resolution\nspatiotemporal atlas of the Bombyx mori wing disc across multiple developmental\nstages. Based on this resource, the authors propose a five-stage Gene Transition Model\n(GTM) regulated by 20-hydroxyecdysone (20E) to describe the temporal dynamics of\ncell fate reprogramming. This is an interesting and timely study that offers novel\ninsights into hormonal regulation during insect organogenesis and provides a\ntheoretical framework for understanding genetic control and potential molecular\nintervention strategies in Lepidoptera. In addition to its biological relevance, the dataset\nrepresents a valuable resource for future research into wing disc development and\nhormonal regulation. Overall the study is well performed and reads clearly. Prior to\npublications, I’d like to see some clarifications and revisions. Main concerns:\nResponse: We sincerely thank the reviewer for the positive evaluation and thoughtful\nsummary of our work. We have carefully addressed each of the main concerns raised\nby the reviewer and have made corresponding revisions throughout the manuscript.\nPage references and detailed point-by-point responses are provided below.\nComment: 2. In Fig. 3h, the comparison between naturally developing and 20E-\ninduced samples could be expanded. It would be valuable to clarify the similarities and\ndifferences in cell composition and pathway activity between the two conditions,\nparticularly considering how closely 20E treatment mimics natural developmental\nprocesses.\nResponse: We sincerely thank the reviewer for this constructive suggestion. We agree\nthat a detailed comparison between natural development and 20E-induced conditions is\nessential for evaluating the fidelity of the hormone induction system. To address this\npoint, we expanded our analyses of both cellular composition and pathway activity in\nthe revised manuscript.\nAt the cellular level, the major cell types captured in both conditions were similar.\nHowever, their relative abundances and activation timings were different. During\nnatural development, the proportions of Wm, Cuticle, and Epithelial lineages changed\ngradually from L5D1 to W1. In contrast, 20E treatment triggered a rapid and transient\nexpansion of these lineages within the first 30 min, followed by a pronounced shift\ntoward Matrix and Immune cell enrichment after 40 min. Notably, Matrix cells\nexhibited a steady increase in both conditions (Fig. S18a).\nAt the transcriptional level, GO and KEGG enrichment analyses revealed distinct\nfunctional states. Natural development (sc group) showed enrichment in Endocytosis,\nmTOR signaling pathway, structural constituent of cuticle, and DNA-binding\ntranscription factor binding, reflecting sustained biosynthetic and morphogenetic\nactivity. Conversely, 20E-induced samples (sn group) were enriched in Ribosome\nbinding, mRNA binding, and Oxidative phosphorylation, indicating enhanced protein\nsynthesis, energy metabolism, and rapid transcriptional reprogramming (Fig. S18b–c).\nThese findings demonstrate that 20E treatment recapitulates the natural developmental\ntrajectory but in a time-compressed manner, accelerating lineage activation and\ntranscriptional remodeling. We have included these results in the revised manuscript\n(Page 13, Lines 326–337) and added Supplementary Fig. S18 to clarify these\ncomparisons.\nComment: 3. For Fig. S17, please provide functional enrichment results for the shared\nupregulated gene set between the natural and 20E-induced groups. This would help\nclarify the core biological processes directly regulated by 20E.\nResponse: We sincerely thank the reviewer for this valuable suggestion. According to\nyour suggestion, we performed Gene Ontology (GO) and KEGG pathway enrichment\nanalyses on the 273 shared upregulated genes identified in both the natural development\nand 20E-induced datasets. The enrichment results revealed that these shared genes are\nsignificantly involved in developmental processes such as cellular morphogenesis,\norgan growth, and tissue remodeling, as well as in stress and immune responses. In\naddition, several pathways related to protein homeostasis and signal transduction—\nincluding autophagy, MAPK signaling, and endoplasmic reticulum (ER) protein\nprocessing—were also enriched. These findings suggest that 20E regulates a core\ntranscriptional program that coordinates both structural remodeling and cellular stress\nadaptation during wing disc development (Fig. S19b). We also analyzed the dynamic\nexpression patterns of these shared genes across natural development and 20E-treated\nsamples (Fig. S19c). The results showed a progressive, stage-specific activation from\nL5D1 to P6 under natural development, transitioning sequentially from clusters C1 to\nC10. Upon 20E treatment, many genes in clusters C3–C6 and C8–C10 responded\nrapidly within 20 minutes, while genes in cluster C1, which were activated early during\nnatural development, were induced later under 20E exposure. These observations\nprovide a high-resolution view of 20E-induced transcriptional reprogramming and\nunderscore its role in accelerating developmental transitions within the wing disc. We\nhave incorporated these analyses into the revised Supplementary Figure S19 and\nupdated the corresponding description in the main text (Page 14, Lines 337–353). The\nlegend of Fig. S19 has also been revised accordingly.\nComment: 4. In Fig. 4h–i, the analysis of cell–cell interactions (CCI) does not clearly\nreflect stage-specific dynamics. The authors are encouraged to highlight temporal shifts\nin interaction patterns across developmental stages.\nResponse: We appreciate the reviewer’s careful work. We fully agree that highlighting\nthe stage-specific dynamics of cell–cell interactions (CCIs) would make the results\nmore intuitive. Accordingly, we revised and expanded the CCI analysis to better\nillustrate temporal changes in intercellular communication across the five\ndevelopmental stages. In the revised version, we visualized the overall interaction\nnetworks for each stage and emphasized their dynamic complexity (Fig. 4l). The cell-\ncell interaction network gradually became more complex from Stage 1 to Stage 4 and\nthen simplified at Stage 5, indicating a transition from active remodeling to stabilization\n(Fig. 4l). To improve clarity, we restructured Figs. 4m-n by highlighting the main\nsignaling sources with bold text and labeling each stage-specifically enriched\ninteractions with the corresponding stage label. This presentation makes the temporal\npattern of CCIs and the stage-specific signaling hubs easier to follow and interpret.\nPlease see Figure 4 for the revised visualization. The corresponding explanations have\nbeen added to the main text of the revised manuscript (Page 17-19, Lines 414–448).\nComment: 5. The adult wing morphology following Rfx RNAi in Spodoptera\nfrugiperda is not shown. Including representative images of adult wing phenotypes\nwould strengthen the phenotypic validation. Given its importance in supporting the\nevolutionary conservation of Rfx, the S. frugiperda RNAi results should be moved to\nthe main figures. Additionally, a comparative analysis of Rfx conservation across\ninsects would further support this claim.\nResponse: We sincerely thank the reviewer for this important and constructive\nsuggestion. In the revised manuscript, we have added representative images of the adult\nwing phenotypes in S. frugiperda following Rfx RNAi to Figure 2 (Fig. 2o and Page\n10, Line 231).\nThese images clearly display the morphological defects, such as shortened veins\nand irregular wing margins, which support the functional disruption of Rfx. To highlight\nthe evolutionary conservation of Rfx, we have moved the Rfx RNAi knock down results\nof S. frugiperda from the Supplementary Materials to the main figures (Fig. 2m–n).\nFurthermore, to strengthen the evidence for evolutionary conservation, we\nperformed a comparative analysis of Rfx orthologs across representative insect species.\nWe constructed a phylogenetic tree using Rfx protein sequences from 23 representative\ninsect species and examined their motif organization and conserved domain\narchitecture. The analysis revealed that Rfx is highly conserved across Lepidoptera and\nother insect lineages, suggesting a conserved regulatory function in wing development\n(Fig. S10). Please see Supplementary Figure S10 below for the phylogenetic and motif\nanalyses. These new results and corresponding figure have been incorporated into the\nrevised manuscript (Supplementary Fig. S10 and Page 9, Lines 220–222). The\nlegends for Fig. 2 and Fig. S10 have also been updated accordingly.\nSupplementary Fig. 10 Conservation and evolutionary analysis of Rfx. a Phylogenetic analysis\nof Rfx orthologs in representative insects. b Conserved motif architecture of Rfx proteins. c Domain\ncomposition of Rfx proteins based on CDD annotation.\nComment: 6. Beyond Rfx, the authors are encouraged to functionally validate\nadditional key transcription factors identified through SCENIC analysis (e.g., Dll or\nBlimp-1) via RNAi or similar approaches.\nResponse: We thank the reviewer for the valuable suggestion. In the revised manuscript,\nwe have included additional functional validation experiments for key transcription\nfactors identified through the SCENIC analysis. Specifically, we focused on Distal-less\n(Dll) and Pur-alpha, both predicted to play central roles in Wm cell differentiation and\nwing disc development. RNAi-mediated knockdown of these genes in Bombyx mori\nresulted in noticeable morphological defects in the wing disc, confirming their essential\nregulatory roles in morphogenesis (Fig. S11). Please see Supplementary Figure S11\nbelow. As for Blimp-1, previous studies have reported that RNAi-mediated knockdown\nof Blimp-1 in Bombyx mori leads to disrupted wing development (Wu SY et al., 2019)1.\nThis literature evidence supports the functional importance of Blimp-1 identified in our\nanalysis and further exemplify the reliability of our transcription factor (TF) screening\nstrategy. Together, these results provide additional support for the functional relevance\nof the key TFs identified by SCENIC. The new experimental results and related\nreferences have been incorporated into the revised manuscript (Fig. S11 and Page 9,\nLines 225–228).\nSupplementary Fig. 11 Effects of dll and pur-α gene knockdown on larval wing discs and adult\nwing development in silkworm. a qPCR analysis of gene changes after Dll and Pur-α genes\ninterference, Data represented by the mean ± SD, were analyzed for intergroup comparison using\nunpaired t-tests, with significance at **, p < 0.01. b Morphological changes in wing discs of\nsilkworn after Dll and Pur-α genes knockdown, showing developmental fissures and defects.\nBright-field imaging revealed the overall tissue architecture. Phalloidin staining labeled actin\nfilaments to visualize cytoskeletal organization. DAPI staining marked nuclear DNA for cellular\nlocalization. Merged images integrated Phalloidin and DAPI signals to delineate spatial\nrelationships between cytoskeletal elements and nuclei. c Knockdown of Dll and Pur-α genes\nresulted in atrophy and malformation of wing development in adult silkworm. Dll knockdown\ninduced wing atrophy with reduced size, Pur-α knockdown led to morphological distortions in adult\nwings.\nComment: 7. In Fig. 3g, it would be helpful to indicate which cell types express the\nhighlighted transcription factors and how their activity contributes to cell-type-specific\ntranscriptional regulation and heterogeneity.\nResponse: We thank the reviewers for pointing out this issue. According to your\nsuggestion, we have provided the cell-type–specific expression patterns and temporal\nactivities of the transcription factors in Supplementary Figure S16a. Based on regulon\nspecificity scores (RSS) and expression dynamics across developmental stages, we\nidentified six key TFs that exhibit distinct temporal activation patterns and cell-type\npreferences. In the early stage (before 30 min), Hr38 (LOC693064), Nubbin\n(LOC101744796), and Six2 (LOC101735990) showed strong activity, regulating gene\nnetworks enriched in Wm, Epithelial 1, and Cuticle 1 cells. After 30 min, CREB\n(LOC692871), FoxN3 (LOC101739133), and GATA-A (LOC10174013) became\npredominant, with activities mainly associated with Immune and Matrix cells. These\nfindings suggest a temporal transition of TF activity from early patterning factors to\nregulators of morphogenesis and tissue remodeling. We have incorporated these\nanalyses and explanatory descriptions into the revised manuscript (Figs. 3g and S16,\nPage 12, Lines 304–312), and the legends of Fig. 3g and Fig. S16 have been updated\naccordingly.\nComment: 8. The evolutionary comparison with Drosophila and B. anynana is\nrelatively brief. A more detailed discussion of conserved and lineage-specific\ndevelopmental strategies across these species would enhance the impact of the\ncomparative analysis.\nResponse: We sincerely thank the reviewer for this valuable suggestion. We fully agree\nthat a deeper evolutionary perspective would enhance the impact of our comparative\nanalysis. Comparing our findings in B. mori with that of D. melanogaster and B.\nanynana revealed both conserved regulatory modules and lineage-specific\ndevelopmental strategies (Fig. S23a and b). Core pathways such as Wnt, Notch, and\nMAPK were consistently involved in cell fate determination and tissue patterning\nacross species. However, our findings highlight several specific adaptations in B. mori.\nSpecifically, the Wm cell in B. mori followed a bifurcated developmental trajectory\ntoward epithelial and cuticle fates. This distinct progenitor identity has not been clearly\nreported in Drosophila, suggesting a silkworm-specific mechanism for lineage\nallocation. In B. anynana, a more diversified set of scale precursor cells supports\nelaborate wing coloration and patterning, highlighting another lineage-specific\nadaptation (Fig. S23c). Furthermore, although 20E functions as a conserved hormone\nregulating metamorphosis, our data point to a silkworm-specific mode of spatial and\ntemporal hormone responsiveness. In B. mori, Wm cells responded rapidly and locally\nto 20E, initiating early morphogenetic programs. This contrasts with the more systemic\nand less spatially resolved 20E response observed in Drosophila. Taken together, these\ncomparisons demonstrate that while fundamental signaling frameworks are preserved,\neach lineage has evolved tailored regulatory strategies—shaped by ecological pressures\nand morphogenetic needs. Our single-cell and spatial atlas of B. mori not only confirms\nthese shared principles but also provides a valuable reference for exploring how\nevolutionary flexibility in hormone integration and progenitor dynamics contributes to\nwing diversity in insects. We have expanded the Discussion section to include this\ndetailed comparison in the revised manuscript (Pages 23-24, Lines 554–572).\nComment: 9. While the methods section is generally well-detailed, additional\nclarification on the threshold settings and statistical parameters used in transcription\nfactor enrichment, module identification, and pathway analysis would improve\ntransparency and reproducibility.\nResponse: We thank the reviewer's careful work. We have thoroughly gone through the\nMethods section and revise accordingly. To be specific, we have added detailed\ndescriptions of the threshold settings and statistical parameters used in transcription\nfactor enrichment, gene module identification, and pathway enrichment analysis. In\naddition, we have further included important methodological details to enhance\ntransparency and reproducibility, such as catalog numbers of key reagents, names and\nversion of software, database sources, and parameter configurations used throughout\nthe analyses. These additions have been incorporated into the revised Methods section\n(Pages 25–31).\nComment: 10. Below are some minor comments for the authors to consider: In Fig. 2j,\nthe figure legend should specify “dsRfx” to match the experimental description. The\nlabels in Fig. 3a–f should be adjusted for better alignment, and the RSS values in Fig.\n3j are too small to read clearly and should be enlarged.\nResponse: We thank the reviewers for pointing out this issue. All suggested corrections\nhave been implemented in the revised manuscript and figure, including but not limited\nto: “dsRfx” has been added to the legend of Fig. 2j to accurately reflect the experimental\ncondition; the labels in Fig. 3a–f have been realigned for improved clarity and visual\nconsistency; the RSS values in Fig. 3j have been enlarged to enhance readability. These\nmodifications have been incorporated into the revised figures and corresponding\nlegends (Page 10, line 231 and Page 15, line 354).\nComment: 11. Several supplementary figures appear to lack sufficient resolution.\nPlease ensure all images are high-quality and include appropriate legends and\nannotations for ease of interpretation.\nResponse: Thank you for your careful review. We have addressed this concern by\nadjusting the resolution of all supplementary figures to ensure high-quality\nvisualization. In addition, we have thoroughly checked the entire manuscript and\nsupplementary materials to confirm that all images are of sufficient resolution (>300\ndpi) and include clear legends and annotations for ease of interpretation. Please see the\nrevised Supplementary data with changes marked.\nComment: 12. Reference formatting should be carefully reviewed. Inconsistencies in\njournal name abbreviations and missing page numbers should be corrected for\nuniformity.\nResponse: We sincerely thank the reviewer for this helpful comment. We have carefully\nreviewed and standardized the reference formatting throughout the manuscript to\nensure full consistency. Please see the revised manuscript with changes marked.\nComment: 13. In the author list, a comma should follow “Ying Lin” to maintain\nconsistent punctuation.\nResponse: We thank the reviewer's careful work. We have corrected the punctuation\nby adding a comma after “Ying Lin” to ensure consistency in the author list (Page 1,\nline 6).\nComment: 14. Line 27: The full name of 20E should be written as “20-\nhydroxyecdysone”.\nResponse: Thank you for your careful review. We have revised the text to spell out the\nfull name of 20E as “20-hydroxyecdysone” (Page 1, line 26).\nComment: 15. Line 90: “Fig. S1” should include a period to follow journal formatting\nguidelines.\nResponse: We thank the reviewer's careful work. We have corrected the formatting\nissue by adding a period to “Fig. S1” (Page 3, line 93).\nComment: 16. Lines 105–106: The font style of reference citations appears\ninconsistent (non-italic). Please standardize formatting across the manuscript.\nResponse: Thank you for your careful work. We have corrected the font style of the\nreference citations (Page 4, line 106). Additionally, we have reviewed the entire\nmanuscript to standardize reference formatting throughly.\nComment: 17. Lines 180 and 185: The correct formatting should be “Fig. 2f and Fig.\nS7c”.\nResponse: We thank the reviewers for pointing out this issue. We have corrected the\nformatting in Lines 180 and 185 to “Fig. 2f and Fig. S7c” (Page 8, line 192). We have\nalso carefully reviewed the entire manuscript to ensure that all figure references follow\nconsistent and proper formatting.\nReference\n1. Wu SY, et al. BmBlimp-1 gene encoding a C2H2 zinc finger protein is required\nfor wing development in the silkworm Bombyx mori. Int J Biol Sci 15, 2664-\n2675 (2019).\nReviewer #2\nComment: 1. The current manuscript outlines a spatial transcriptomic approach to map\nthe development of wing discs in Bombyx mori. The authors have compiled an\nimpressive dataset, and undertaken extensive analysis, to identify the cell types within\nthis tissue, their spatial location, the changes in cell fates over time, and changes in gene\nexpression that accompany these changes. In addition, they explore how the steroid\nhormone ecdysone affects wing disc development, demonstrating that it regulates the\npace of wing disc development. We commend the authors on the scale of this work, and\nfeel that it makes a valuable contribution to our understanding of organ development.\nWe have no experience with the transcriptomics approaches used and thus will limit our\ncomments to other aspects of the work, which we think would help strengthen the\nmanuscript.\nResponse: We deeply appreciate the respected reviewer for their meticulous\nexamination and valuable suggestions on our manuscript. According to your\nsuggestions, we have made extensive modifications to our manuscript and\nsupplemented extra data to make our results convincing. Thank you again for your\npositive comments and valuable suggestions to improve the quality of our manuscript.\nComment: 2. The acronyms for the different cell types are really hard to keep track of.\nCan you use the whole word to make it easier to follow?\nResponse: We thank the reviewers for pointing out this issue. We fully understand that\nthe numerous acronyms may be difficult to follow. To improve clarity and readability,\nwe have carefully revised the text and figures to ensure that all abbreviations are clearly\ndefined and consistently explained.\nIn our study, due to the large number of cell types identified through single-cell\nanalysis, we adopted standard practice by using abbreviations (e.g., Wm for wing\nmorphogenesis, Im for immune, etc.) to enhance figure readability and minimize label\novercrowding. Nonetheless, we fully recognize the importance of clarity. Therefore, we\nhave adopted the following strategy:\nFor cell types with longer or more complex names (e.g., wing morphogenesis), we\nuse abbreviations such as “Wm” to simplify figure labels while maintaining readability.\nFor short or self-explanatory names (e.g., Cuticle 1–3), we now use the full terms\ninstead of abbreviations.\nFor figure clarity and layout efficiency, some figures are abbreviated as needed;\nall abbreviations are defined in the corresponding figure legends.\nIn addition, we have carefully revised all figures, figure legends, and relevant parts\nof the main text to ensure that all abbreviations are clearly defined and explained at first\nmention. We have also added explicit acronym-to-full-name mappings. All related\nfigures, legends, and relevant text sections have been carefully updated to maintain\nconsistency. Please see the revised manuscript for these modifications with changes\nmarked.\nComment: 3. The experiments knocking down Rfx in the wing disc were very\ninteresting and highlighted the role of this gene in Wing Morphogenesis cell type\ndevelopment. However, it wasn’t clear why this particular cell type was the focus?\nResponse: We sincerely thank the reviewer for this insightful question. Our focus on\nthe Wing Morphogenesis (Wm) cell type was based on several layers of observations\nduring wing disc (WD) development.\nFirst, when examining the dynamic changes in the proportion and differentiation\ntrajectory of cell types across developmental stages, we observed that most cell types\nin the WD maintained a relatively stable state and did not exhibit significant lineage\nprogression (as shown in Fig. S6a–f). In contrast, Wm cells demonstrated a markedly\ndynamic and continuous differentiation trajectory, with clear transitions toward cuticle\nand epithelial (EP) lineages (Fig. S6g). This suggests that Wm cells are a pivotal cell\ntype in the developmental progression of the wing disc. Second, spatial transcriptomics\ndata further corroborated this finding. We identified Wm cells as being specifically\nenriched in the wing bud region, the future wing-forming domain of the WD (Fig. 1g–\nh). This spatial localization aligns with their predicted role in driving the formation of\nwing-specific structures. Motivated by the above observations, we sought to uncover\npotential transcriptional regulators specifically involved in wing disc cell fate decisions.\nThrough the reconstruction of a transcription factor regulatory network and RSS\nscoring, we found that Rfx emerged as the top-ranked specific transcription factor at\nthe L5D1 stage. Which leading us to hypothesize that Rfx may serve as a core regulator\nin Wm cell lineage commitment.\nTo validate this hypothesis, we performed RNAi-mediated knockdown of Rfx and\nobserved striking disruptions in wing disc morphology and gene expression associated\nwith Wm-derived lineages (Fig. 2j–l and Fig. S9), confirming its critical role in\nmodulating Wm cell fate and wing development.\nComment: 4. This might sound a bit picky, but the order of the methods did not follow\nthe order of experiments. This isn’t necessary, but does help when trying to find the\nrelevant method to match to the experiment.\nResponse: We thank the reviewer for this valuable suggestion. We fully agree that\norganizing the Methods section to follow the experimental sequence described in the\nResults enhances readability and makes it easier for readers to locate relevant\nmethodological details. To address this, we have reordered the Methods subsections in\nthe revised manuscript to match the chronological flow of experiments presented in the\nmain text—from sample collection and sequencing, through data processing and\nanalyses, to functional validation and hormone treatment assays. Please see the revised\nMethods section (pages 25–31) in the revised manuscript.\nComment: 5. For the wing disc culture experiments, more detail is required in the\nmethods to understand how these were performed. What medium did you use for\nculture and did you add any insulin or bovine serum? Preferably, this would have it’s\nown section. For the experimental design, why did you choose to change the time of\nexposure instead of 20E dose? Lastly, 2µg / ml of 20E is a high dose for wing disc\ncultures. It would be useful to explain why you chose this dose.\nResponse: Thank you for your valuable suggestions. We have added a dedicated\nsection in the Methods to comprehensively describe the in vitro wing disc culture and\n20E treatment experiments. Specifically, the wing discs were cultured in Grace’s\nmedium supplemented with 10% fetal bovine serum.\nRegarding the experimental design, we chose to vary the duration of 20E exposure\nrather than its concentration. This decision was based on two main considerations. First,\ndue to the open circulatory system of Bombyx mori, endogenous 20E levels in the\nhemolymph can influence multiple tissues simultaneously. Second, the permeability\nand uptake efficiency of exogenous 20E by the wing disc in vitro is unclear. This\napproach allowed us to capture the temporal sequence of gene activation in response to\n20E and to compare it with the natural developmental timeline.\nTo ensure a robust and consistent hormonal response will be induced in cultured\nwing discs, we performed a preliminary dose-response experiment to select optimized\nconcentration for 20E. Wing discs were cultured with various concentrations of 20E (0,\n0.1, 0.2, 0.5, 1, 2, 5, and 10 μg/mL) and relative expression of key 20E responsive genes\nwere measured using quantitative PCR. We found that at 2 µg/mL, the transcriptional\nlevels of EcR and USP, were markedly induced (Fig. S13a). This concentration is also\nsupported by previous studies reporting effective hormonal stimulation at similar levels\nin cultured wing discs (Ali et al., 2019)1. Hence, we used 2 µg/mL 20E to stimulate\nwing disc in vitro to investigate gene expression dynamics. All these details have been\nincorporated into the revised Methods section (Page 11 lines 266–269 and pages 28-\n29 lines 707–730).\nSupplementary Fig. 13 Single-nucleus RNA sequencing (snRNA-seq) analysis of 20E treated\nwing disc. a Changes in relative expression levels of EcR and USP genes in wing discs after 4 h of\nin vitro culture under different concentrations of 20E hormone. b UMAP plots showing the\ndistribution of different cell types at various time points (0 min–6 h). c Temporal changes in the\nproportions of different cell types. d The proportion of time point-specific upregulated genes in\ndifferent cell types.\nComment: 6. Why were 20E titres collected from the wing discs and not from the\nhemolymph?\nResponse: Thank you for your thoughtful question. We chose to analyze 20E titers\ndirectly from the wing discs (WDs) rather than the hemolymph for the following\nreasons. First, the dynamic changes of 20E concentration in the hemolymph during\nsilkworm development have already been well documented in previous studies. Second,\nthe silkworm possesses an open circulatory system, and 20E in the hemolymph is\ndistributed to various tissues and organs, not exclusively acting on the wing discs.\nTherefore, measuring 20E in the hemolymph would not provide precise information\nabout the local hormonal milieu experienced by the developing WDs. Our goal was to\ninvestigate the tissue-specific 20E dynamics during wing morphogenesis, and thus we\ndirectly measured the hormone levels in the WDs. We have clarified this rationale in\nthe Method section of revised manuscript to improve transparency (Page 30, line 774–\n778).\nComment: 7. Some of the paragraphs were in really big blocks of text. It would help\nwith readability to break up these paragraphs into discrete ideas and to provide a bit\nmore context throughout with regards to the aims of your different analyses.\nResponse: We thank the reviewer for the valuable suggestion. We have broken up the\nlong paragraphs into shorter, focused units and added brief contextual sentences to\nclarify the purpose of each analysis, improving overall readability and flow. Please see\nthe version of revised manuscript with changes marked.\nComment: 8. Minor comments: Line 38-39: not all insects have wing discs\n(ametabolous and hemimetabolous insects have wing buds which are analogous)\nResponse: We thank the reviewer's careful work. To address this, we have revised the\nsentence to clarify that wing discs are specific to holometabolous insects (Page 2,\nlines 38-39).\nComment: 9. Line 48-50: Rephrase this sentence as it sounds like you are saying that\nBicyclus anynana is not a lepidopteran.\nResponse: We agree that the original sentence may unintentionally imply that Bicyclus\nanynana is not a lepidopteran. We have revised the sentence as “Although single-cell\nanalyses have been conducted for the wing disc of Drosophila melanogaster (D. melanogaster)\nand Bicyclus anynana (B. anynana), they have mainly focused on cell type composition and\ngene expression patterns at isolated stages.”, to improve clarity and accuracy (Page 2, lines\n48–50).\nComment: 10. Figure 2K legend: what do the arrowheads point to? Also, there are two\nsets of dsEGFP and two sets of dsRfx panels. How are these different?\nResponse: We greatly appreciate your attention to this detail. We apologize for the lack\nof clarity in the original figure legend. The arrowheads in Figure 2K highlight regions\nof morphological changes in the wing disc before and after RNAi treatment (Page 11,\nLines 245–252). As for the two sets of dsEGFP and dsRfx panels: they represent the\nsame experimental groups but illustrate different levels of observation. The upper\npanels display the overall morphology of the wing discs. The lower panels show\ntransverse histological sections of the same samples, allowing visualization of\nstructural differences in tissue organization and cellular morphology after RNAi\ntreatment. We have revised both Figure 2 and its legend to clarify this distinction (Fig.\n2 and Page 11, Lines 246–248).\nComment: 11. Line 242: You mention Usp is expressed in Wm cells, but EcR isn’t in\nthe list of ecdysone response genes?\nResponse: We thank the reviewer’s careful work. We found that EcR gene expression\nwas barely detected in our snRNA-seq dataset. This might be due to stochastic dropout\nevents, arising from the limited RNA content in single cell nucleus and inefficiencies\nof capture and reverse transcription steps during snRNA-seq experiment (Kharchenko\net al., 2014)2. As a result, transcripts expressed at low abundance or in a transient\nmanner are often not detected, even when they are biologically present in the cell. After\n20E treatment, we observed that Usp gene have a higher expression level than EcR\ngene, as shown in Supplementary Fig. 13a above. This might explain why Usp\nexpression was sufficiently robust to be captured during the experiment, whereas EcR\nmay have been underrepresented due to such stochastic dropout effects.\nComment: 12. In the text, you refer to 12 supplementary tables. These tables are found\nin the excel file. However, in the pdf there are only 3 supplementary tables and the\nformatting of these tables appears to be off.\nResponse: Thank you for your careful work. We sincerely apologize for the confusion\ncaused by the inconsistency between the number of supplementary tables and those\nincluded in the PDF. This has now been corrected: all 12 supplementary tables are\nproperly provided in the updated Excel file, and the PDF version has been revised to\ninclude the appropriate references and formatting. We would also like to note that the\nsupplementary tables were intentionally not converted into PDF format for\nresubmission. Due to the layout and data density of these tables, conversion may cause\ndisplay or formatting issues when generating the integrated manuscript PDF. We have\nexplained this to the editor in advance and clarified that all supplementary tables are\nprovided as separate Excel files to ensure clarity and accessibility of the data. Please\nsee the revised supplementary tables.\nComment: 13. For Figure 3d (and the accompanying supplementary figure), it wasn’t\nclear what was meant by relative information flow.\nResponse: We thank the reviewers for pointing out this issue. We have clarified the\nmeaning of “relative information flow” in the revised version of Figure 3d legend and\nthe associated supplementary figure 14a legend. Specifically, it refers to the normalized\nstrength and direction of ligand–receptor–mediated signaling interactions among cell\ntypes, inferred based on gene expression profiles and prior knowledge of ligand–\nreceptor pairs. Please see page 15 lines 361–363 in the revised manuscript.\nComment: 14. Figure S18: What do BF and Trans refer to?\nResponse: We thank the reviewers for pointing out this issue. We have clarified the\nmeaning of this abbreviation in the revised figure legend. Specifically, “BF” has been\nclarified as Bright field, and “Trans” as Transverse section (Fig. S20 and Page 11, lines\n245-251). All relevant instances in the figure have been updated accordingly.\nReference:\n1. Ali MS, Mishra B, Swapon AH, Yamaguchi M. Successively expressed cuticular protein genes\nat the prepupal stage in wing discs of Bombyx mori. The Journal of Basic and Applied Zoology\n80, (2019).\n2. Peter V Kharchenko, Lev Silberstein, David T Scadden. Bayesian approach to single-cell\ndifferential expression analysis. Nat Methods 11, 740–742 (2014)\nReviewer #3\nComment: I co-reviewed this manuscript with one of the reviewers who provided the\nlisted reports. This is part of the Nature Communications initiative to facilitate training\nin peer review and to provide appropriate recognition for Early Career Researchers who\nco-review manuscripts.\nResponse: We sincerely thank the reviewer for their contribution as a co-reviewer and\nfor taking the time to evaluate our manuscript. We truly appreciate your effort and\nthoughtful participation in the review process.\nYour collaborative review, along with the feedback provided by the other\nreviewers, has been highly constructive and has helped us substantially improve the\nclarity, structure, and scientific rigor of our study.\nWe have carefully addressed all comments and incorporated the suggested\nrevisions throughout the manuscript. Please see the detailed point-by-point responses\nabove for the corresponding changes.\nResponse to reviewers\nWe deeply appreciate the respected reviewer for thier very professional\ncomments and suggestions, which help us to improve the quality of the manuscript\nsignificantly!"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed all my previous concerns, and the revisions have notably improved the clarity and quality of the\nmanuscript. The newly added analysis and experiment strengthen the conclusions. It is a valuable dataset resource for\nfuture studies on insect development and offers meaningful insights into hormone-regulated wing disc morphogenesis. I\nhave no further comments.\nI have also read through the comments provided by the other reviewers. Their suggestions were constructive and\nhighlighted several important aspects of data interpretation and presentation. The authors have taken these points seriously,\nand the revision shows clear improvement in areas such as the explanation of analytical methods, the organization of the\nfigures, and the interpretation of key results. After considering these revisions together with the insights raised by the other\nreviewers, I find the current version more coherent and scientifically rigorous.\nOverall, I have no additional concerns and recommend the manuscript for publication.", "reviewer_2": "Thank you for addressing all of our comments so thoroughly."}, "author_response": "Response to reviewers\nWe truly appreciate your insightful comments throughout the review process, which greatly improved the quality of our work. We appreciate the reviewer’s recognition of the value of this dataset and the improvements made in response to the other reviewers’ constructive comments.\nResponse: We thank the reviewer for the positive assessment of our revision. We truly appreciate the reviewer’s constructive comments throughout the review process, which greatly improved the quality of the manuscript."}]

Decision letter