Dynamic single-cell transcriptomics reveals lsamp-guided neural network formation in male S. japonicum driving female reproduction

Domain: biology
Published: 2026-01-22
Editorial outcome: minor_revision

Abstract
Sexual maturation and egg production in schistosomes requires sustained male-female pairing, yet the cellular dynamics underlying this process remain poorly understood. Here, we generate a high-resolution single-cell transcriptomic atlas of 104,671 cells from male and female Schistosoma japonicum across four developmental stages, identifying 76 distinct cell populations and mapping lineage trajectories in reproductive and neural systems. Notably, we discover a male-specific neuron subtype-N4.3 neuron located in the gynecophoric canal, and identify lsamp as a critical gene co-expressed with nrps in N4.3 neurons. Loss of lsamp impairs axonal integrity and disrupts the male ventral nerve network, reducing the targeted delivery of pheromone BATT precursor to N4.3 neurons, thereby hindering female reproductive maturation. This work maps schistosoma male-female interactions at single-cell resolution, revealing a male ventral nerve network that facilitates intersex communication and promotes female reproductive development, offering strategies to disrupt parasite reproduction and combat schistosomiasis.

Introduction
Schistosomiasis is a prevalent and devastating neglected tropical disease caused by parasitic flatworms called Schistosoma spp., which affects over 250 million people in disadvantaged areas1. Currently, no vaccine is available, and treatment relies on a single drug (praziquantel)2,3. Therefore, new therapeutic and preventive strategies are urgently needed.

Eggs produced by mature female schistosomes drive the pathology and transmission of this disease. In a definitive mammalian host, one mature female parasite could produce hundreds to thousands of eggs daily that are deposited in the host tissues, e.g., the liver, where they induce the local granuloma and further lead to tissue fibrosis4,5. Eggs excreted with feces or urine hatch in the water environment, and then the released miracidia invade the intermediate snail host before spreading. Therefore, developing strategies for obstructing the egg production of these parasites will contribute to the control of this disease.

Schistosomes have a unique reproductive feature. The pairing with the male worm initiates the sexual development of a female worm and maintains her mature state6. Unlike females, the male worm can reach sexual maturity independently7. Continued investigations have been carried out to uncover the underlying mechanism of this male-induced female reproduction since the first report of the phenomenon in Schistosoma japonicum nearly a century ago8. Recently, a dipeptide pheromone (β-alanyl-tryptamine) synthesized and secreted by male worms has been identified in S. mansoni that stimulates female sexual development9. The discovery of this neuron-derived key molecule that can be transmitted from males to females highlights the importance of understanding the complex network during male-female interplay.

The direct physical connection between males and females occurs at their late juvenile stage when the male worm clasps the female through the gynecophoric canal to form a worm pair10. Thereafter, the coupled male and female continue to grow and turn to be sexually mature. Besides the extension of body length, a series of previous studies have observed dramatic gene expression changes in males and females upon pairing10–13. These changes are not only restricted to the sexual organs but also to other somatic tissues, such as neurons and muscles. Additionally, labor division is more apparent between the two sexes at their adult stage, which is represented by the distinct sex-associated gene expression patterns10,14–17. Recently, researchers have revealed that the cell types of S. mansoni expanded along with the growth and development of the parasite in the mammalian hosts18–20. Notably, single-cell RNA sequencing showed distinct cell clusters between the adult male and female and between the virgin female and mature female in S. mansoni21, which highly indicates a dynamic sex-related cellular differentiation during the male-female interaction. However, the whole picture of the molecular changes of males and females during their interaction at the single-cell level is still lacking.

In this study, we construct the dynamic single-cell atlas of male and female S. japonicum across 4 timepoints throughout the reproductive development during their interplay. We sequence 104,671 single-cell transcriptomes and characterize 76 molecularly distinct cell populations. Along with the male-female interaction, the cell clusters of the reproductive and nervous systems undergo the most dramatic changes in both sexes. Trajectory inference analysis reveals the detailed lineage of cells related to reproductive development in both sexes. We analyze and predict gene expression patterns and key regulatory factors involved in the maturation of late male and female germ cells. Furthermore, we notice a complex heterogeneity and sex-biased characteristics of the nervous system during the maturation of the parasite. Notably, we identify a lsamp gene expressed within a male-specific neuron population, which triggers female reproductive development by maintaining axonal integrity within the male gynecophoric canal in S. japonicum. Taken together, our datasets provide a comprehensive single-cell atlas of male and female S. japonicum across their maturation and reveal the key regulator of a male-specific neuron that is required for the stimulus of female sexual development, which will further prompt the innovation of strategies against this human parasite.

Results
Dynamics single-cell transcriptomes of male and female S. japonicum during their interplay

S. japonicum primarily undergoes pairing and sexual maturation between 14 and 26 dpi (days post-infection)10. To characterize the dynamic single-cell atlas of S. japonicum throughout the male-female pairing and subsequent maturation process, we collected worms at 4 timepoints (14, 18, 22, and 26 dpi) from a mouse model (Fig. 1a). Whole worms were separated by sex at each stage, except for those from 14 dpi (since their sex could not be distinguished morphologically at this stage). Subsequently, we dissociated these samples and isolated live single cells using fluorescence-activated cell sorting (FACS) (Supplementary Fig. 1a). Utilizing a 10× Genomics Chromium controller, we generated 7 single-cell RNA-seq (scRNA-seq) libraries (Fig. 1a and Supplementary Data 1). Through the analysis of these libraries using Cell Ranger, we obtained a total of 110,204 cells (Supplementary Data 1). After excluding cells with more than 30% mitochondrial reads or fewer than 500 detected genes, the final count was 104,671 cells, which were further clustered into 76 distinct cell clusters (Supplementary Fig. 1b and Data 1). Furthermore, we annotated the identities of these cell clusters by referring to cell markers defined in the single-cell sequencing of S. mansoni18,19,21,22 (Supplementary Data 2).

In total, the S. japonicum single-cell atlas comprised 33 transcriptionally distinct neuron cell clusters, 14 neoblast clusters, 9 tegument cell clusters, 5 muscle cell clusters, 5 parenchyma clusters, 3 germline stem cell (GSC) clusters, 1 vitellaria cluster, 1 late male germ cell cluster, 1 Mehlis’ gland cell cluster, 1 gut cell cluster, 1 flame cell cluster, 1 esophageal gland cell cluster, and 1 unknown cluster (Supplementary Fig. 1b).

Our single-cell atlas provided detailed lineage information for the parasite’s specific tissues. Interestingly, our data unexpectedly revealed many schistosome neural progenitor cells and neural precursor cells. In comparison, the single-cell atlas of S. mansoni depicts neural cells as scattered clusters lacking apparent lineage connections21. Based on their marked heterogeneity, we classified the neural cells into five categories: neural progenitor cells, neural precursor cells, kk7+ neurons, nonciliated neurons, and ciliated neurons (Fig. 1b). The tegument cells are divided into tegument progenitor, tegument progeny 1, tegument progeny 2, syncytial 1, and syncytial 2 (Fig. 1b). In addition to the well-defined cell markers in previous studies on S. mansoni21, we further identified the marker genes as shown in Supplementary Fig. 2 and Supplementary Data 3 and 4 through whole-mount in situ hybridization (WISH) and double fluorescence in situ hybridization (double-FISH), including new markers such as Sjc-0003010 and Sjc-0008393 for GSCs (Supplementary Fig. 3a), Sjc-0001566 (ccdc) for late male germ cells (Supplementary Fig. 3c), Sjc-0004922 (Tubulin alpha-1 chain) for neoblasts (Supplementary Fig. 3e), Sjc-0008047 (klf11) for neural progenitor cells (Supplementary Fig. 3f), Sjc-0005931 (acss2), and Sjc-0007590 (pcdh7) for neural precursor cells (Supplementary Fig. 3g), Sjc-0002588 (slc6a5) for parenchyma cells (Supplementary Figs. 3i and 4a, b), Sjc-0001890 (sgf1) for esophageal gland cells (Supplementary Fig. 3k), Sjc-0008478 and Sjc-0000277 (slc7a14) for syncytial of tegument cells (Supplementary Figs. 3m and 4a, c), Sjc-0008416 (cpb) and Sjc-0000095 (cacnb1) for muscle cells (Supplementary Figs. 3n and 4a, d).

We first analyzed the distribution of cells in each sample, which allowed us to identify developmental differences at the cellular level between male and female worms at various stages (Supplementary Fig. 5). We also examined the relative proportions of cells from different stages within the cell types of the combined S. japonicum dataset (Supplementary Fig. 6a). As expected, the proportion of the reproduction-related cell types gradually increased throughout the sexual maturation process in both sexes, such as GSCs, female vitellaria, and late male germ cells (Supplementary Fig. 6a and Data 5 and 6). To validate this, we performed WISH using the GSCs marker gene eled18. As shown in Supplementary Fig. 6b, the expression region of eled expanded significantly in both female and male sexual organs from 18 to 26 dpi, which was consistent with the single-cell data. The proportion of neural progenitor cells and precursor cells decreased after 14 dpi in both sexes. In male worms, the populations of differentiated neural cells (e.g., nonciliated neurons and ciliated neurons) continuously expanded. In female worms, these populations did not show sustained expansion but rather exhibited an initial increase followed by a subsequent decline after 22 dpi, which is consistent with the weakening of neural activity in females, as revealed by bulk RNA-seq analysis10 (Supplementary Fig. 6a and Data 5).

Detailed lineage analysis of cells related to reproductive development in both sexes

The schistosome germ cells originate from a subset of stem cells at the asexual stage and eventually separate from somatic cell lineage at the juvenile stage18. Upon pairing, the committed GSCs in both sexes proliferate and differentiate into mature gametes. To demonstrate the cellular dynamics of the parasite’s reproductive system and the relationship between germline and somatic lineages, we subset the cells defined as stem cells (neoblasts) and those present in the reproductive organs, including GSCs, late male germ cells, late female germ cells, vitellaria, and Mehlis’ glands, from both female and male samples (18, 22, and 26 dpi), and performed subcluster analysis. Then, we identified 17 clusters with representative markers for each subtype based on well-defined markers from previous studies (Supplementary Fig. 7a and Data 2), such as the neoblast marker nanos-2 (Sjc-0005962)18,23, the vitellaria S1 cell marker vf1 (Sjc-0000108)24, the GSC markers nanos-1 (Sjc-0006622), and eled (Sjc-0009220)18,25 (Supplementary Fig. 7b, c).

We also validated newly identified tissue-specific markers using WISH, including Sjc-0000843 (ribc2, RIB43A-like with coiled-coils protein 2) for late male germ cells, Sjc-0003040 for late female germ cells, esg-2a (eggshell protein gene 2a26, Sjc-0006798) for vitelline cells, and Sjc-0000795 for Mehlis’ gland (Fig. 2a, Supplementary Fig. 7b, c, and Data 7). These markers can effectively characterize the developmental status of each organ. Interestingly, late female germ cells were also observed in male worms, albeit in relatively low numbers (Supplementary Fig. 7a, c). To confirm the presence of late female germ cells in males, we performed WISH using the late female germ cell marker gene Sjc-0004512 (clec, C-type lectin domain-containing protein) and observed a signal in the testes of male worms (Supplementary Fig. 7d). To further characterize clec expression within the testis, we conducted double-FISH with Sjc-0000843 (ribc2), a marker of late male germ cells. Based on testicular morphology visualized by DAPI staining (outlined by circles), both clec and ribc2 were expressed in the testes. Notably, weak ribc2 signals were also observed in clec⁺ germ cells (Supplementary Fig. 7e). Interestingly, although schistosomes are the only dioecious species within the class Trematoda, various forms of hermaphroditism have been documented, including males developing ovarian-like structures and females exhibiting testicular characteristics27–29. Our identification of cells expressing female germ cell markers in phenotypically normal male worms provides additional molecular insight into the evolution of sexual differentiation in schistosomes.

By performing RNA velocity analysis30, we inferred the possible developmental trajectories from GSCs to late germ cells, separately for female and male worms (Fig. 2b). To further investigate the differences in the germline lineage, we performed trajectory inference using Slingshot31, an expression-based method, on the GSCs, GSC progenies, and the germ cells of male and female worms. It revealed the cell trajectory originating from the GSCs population into late germ cells, in both sexes (Fig. 2c, d).

Subsequently, we analyzed the expression patterns of genes involved in the differentiation of late male or female germ cells. Temporally dynamic genes were grouped into six modules for each sex, each showing distinct expression trends and functional enrichments (Supplementary Fig. 8 and Data 8 and 9).

In females, 1596 genes formed modules associated with translation, RNA processing, DNA conformation change, vesicle organization, and muscle contraction (Supplementary Fig. 8a and Data 10). In males, 1,913 genes were grouped into modules linked to translation, protein processing, ERBB signaling, cilium organization, and microtubule-based movement (Supplementary Fig. 8b and Data 11). Shared dynamic genes between sexes were also identified (Supplementary Data 12).

Analysis of known and newly identified germline regulators showed that meiob, nuob, and horm-2 exhibited similar biphasic expression in both sexes—rising early then declining—indicating conserved roles in early-to-mid gametogenesis; In contrast, Sjc-0003040 (female) and ribc2 (male) were upregulated during late gametogenesis, suggesting sex-specific functions (Fig. 2e, f). Temporal profiling of established germline regulators revealed that nanos-1 was highly expressed in early stages, consistent with its role in germ stem cell maintenance; Conversely, oc-1 and boule promoted spermatogenic differentiation by suppressing self-renewal, while eled inhibited sperm differentiation32 (Fig. 2g). Together, these findings provide a pseudotemporal framework of germ cell development in S. japonicum, revealing that distinct gene modules coordinate sex-specific and stage-specific processes during gametogenesis.

Vitellaria is the largest reproductive organ in adult female schistosomes that occupies two-thirds of the body’s length. We found that both vitellaria S1 cells and GSCs express a certain level of eled (Sjc-0009220) (Supplementary Fig. 7b, c), suggesting a potential shared origin from a population of eled+ stem cells, as previously reported18. Vitellaria maturation is crucial for egg production in schistosomes. Previous studies have shown that vitellaria development can be divided into stages S1–S4, all of which express vf1 (Sjc-0000108)24. In the single-cell transcriptomic analysis of S. mansoni, the vitellaria development process was further categorized into five stages: S1, S1 progeny, early vitellocytes, late vitellocytes, and mature vitellocytes21. However, the detailed differentiation trajectory of vitellaria cells has not been fully elucidated at the single-cell level. Here, we performed subclustering of vitelline cells and identified eight distinct subtypes, including S1 (S1-0, and S1-1), S1 progeny (S1p-0, S1p-1, and S1p-2), early vitellocytes (EVi), and late vitellocytes (LVi-0, and LVi-1) (Supplementary Fig. 9a, b and Data 13). Trajectory analysis with Slingshot revealed two differentiation branches (Lineage 1 and Lineage 2; Supplementary Fig. 9c, d), among which only Lineage 1 showed high differentiation potential according to CytoTRACE33 (Supplementary Fig. 9e). Gene expression dynamics along Lineage 1 indicated that, in addition to the known marker vf1, p48 and the Sjc-0006798 (esg-2a, identified in this study) gradually increased in expression (Supplementary Fig. 9f), supporting Lineage 1 as the primary vitelline developmental pathway. This trajectory, annotated in the RNA velocity plot (Fig. 2b), traces the progression from S1 (vf1+ esg-2a−) to S1 progeny (vf1+ esg-2a+), and subsequently to EVi and LVi (Supplementary Fig. 10). Pseudotime-associated genes by tradeSeq34 (p < 1 × 10−5) were grouped into six co-expression modules (Supplementary Fig. 9g and Data 14) reflecting stage-specific processes: early–mid modules (C1, C3, C6) were enriched for DNA replication and cell cycle functions, whereas late modules (C2, C4) were enriched for eggshell formation–related pathways, including protein localization to the ER and purine metabolism (Supplementary Data 15).

Mehlis’ gland, a specialized organ unique to trematodes and cestodes, has been long speculated to play an essential role in eggshell formation and the release of eggshell granules35,36. However, research on its function in schistosomes remains limited. In our study, we identified candidate marker genes for the Mehlis’ gland, enriching its molecular characterization (Supplementary Data 7). Additionally, we observed that Mehlis’ gland cells emerge earlier than vitelline cells, appearing in large numbers by 18 dpi, while vitelline cells do not proliferate in large numbers until 22 dpi (Supplementary Fig. 11a, b). These findings were further validated through WISH experiments using the marker gene Sjc-0000795 for Mehlis’ gland, the marker gene esg-2a for vitellaria and the marker gene Sjc-0000278 (rreb1) for both vitellaria and Mehlis’ gland (Supplementary Fig. 11c).

Five neural lineages in S. japonicum

Unlike the scattered neuronal clusters revealed from the previously reported single-cell transcriptomic atlas of adult Schistosoma mansoni21, our datasets exhibit high connectivity among complex neuronal clusters (Fig. 1b and Supplementary Fig. 1b), which is likely to enable cell lineage reconstruction of the whole schistosome nervous system from transcriptomic data. We first reclustered all the neural cells of males and females, including neural progenitor cells, neural precursor cells, nonciliated neurons, and ciliated neurons, identifying 19 subclusters and summarizing them into 5 lineages, which were named N1 lineage, N2 lineage, N3 lineage, N4 lineage, and N5 lineage (Supplementary Figs. 12 and 13a).

Subsequently, we employed Slingshot to confirm and refine the reconstructed cell differentiation lineages and to infer pseudotime (Supplementary Figs. 12 and 13b), which revealed a clear alignment with sample collection times (Supplementary Figs. 12 and 13c). Additionally, we performed RNA velocity analysis on the male and female datasets separately using scVelo algorithms37. The velocity vectors were visualized on the UMAP embedding, revealing potential differentiation trajectories from neural progenitor cells toward mature cell types (Supplementary Figs. 12 and 13d).

To better understand the development trajectories leading up to different neuron populations, we performed CellRank38 analysis on both male and female neuron cells. In male S. japonicum, we identified 7 metastable states in the neural cells, with neural progenitor cells predicted as the initial state and N1.1, N2.2, N3.2, N4.2, N4.3, and N4.4 neuron inferred as six terminal states (Supplementary Fig. 12e). In female S. japonicum, we identified 7 metastable states in the neural cells, with neural progenitor cells predicted as the initial state and N1.1, N2.1, N3.1, N4.1, N4.2, and N5 neuron inferred as six terminal states (Supplementary Fig. 13e). In summary, we defined neural cell differentiation from an early neural progenitor state to several mature states and provided a mechanistic framework for characterizing worm neuronal cell fate decisions.

Male-enriched neuron clusters

To explore the neuronal diversity in S. japonicum, we performed a subcluster analysis of neuronal cells in both male and female parasites (Fig. 3a). From the integrated single-cell atlas of both sexes, we identified three neuron populations enriched in males: N2.2, N3.2, and N4.3 (Fig. 3a). The proportional distribution of neuronal subtypes in both sexes (Supplementary Fig. 14a) and the differential abundance analysis of neuronal cell types across sexes (Supplementary Fig. 14b) further confirmed the existence of these male-enriched neurons.

Using RNA probes targeting specific marker genes (Supplementary Fig. 15 and Data 16), we further validated the spatial distribution of male-enriched neuronal populations—N2.2 (Sjc-0003639+ Sjc-0000814 (nrps)−), N3.2 (Sjc-0001200+), and N4.3 (nrps+)—as well as additional mature neuronal populations, including N2.1 (ano7, Sjc-000798+), N3.1 (Sjc-0006108+), N4.1 (bmp, Sjc-0004182+), N4.4 (fhl2, Sjc-0002320+), and N5 (Sjc-0000744+), in both sexes (Supplementary Figs. 16 and 17 and Data 16). As anticipated, all neurons expressed the pan-neuronal marker 7b2 (Supplementary Figs. 16 and 17). Notably, all three male-enriched neuronal subtypes were located within the gynecophoric canal, a specialized male structure that secures the female during pairing (Fig. 3b and Supplementary Fig. 16). Among these, the N4.3 population was exclusively distributed throughout the gynecophoric canal, whereas N2.2 and N3.2 populations were additionally detected in the head of both male and female parasites (Fig. 3b and Supplementary Fig. 16). To integrate these observations, we generated representations summarizing the spatial distribution patterns of all validated neuronal populations (Supplementary Fig. 18).

Interestingly, the marker gene Sjc-0000814 of N4.3 is nrps (nonribosomal peptide synthetase), which encodes the NRPS enzyme. In S. mansoni, the NRPS enzyme has been reported to synthesize β-alanyl-tryptamine (BATT), a dipeptide pheromone that passes from the males to induce female reproductive development9. Thus, our study clarified that the N4.3 neurons are the key populations that regulate female sexual development by expressing the nrps to produce BATT.

To elucidate the potential functions of these three male-enriched populations, we performed GO (Gene Ontology) enrichment analysis on the highly expressed genes within each neuronal subclass. These enriched genes indicated that N2.2 and N3.2 populations exhibited typical neuronal functions, such as synaptic transmission and neurotransmitter transport. Moreover, the N3.2 population showed more specific enrichments, particularly in ion transport (Fig. 3c, Supplementary Fig. 19, and Data 17 and 18). Unlike them, the N4.3 population emerged as a specialized neuronal subtype, expressing genes closely related to cilium organization, cell projection assembly, and microtubule bundle formation (Fig. 3c and Supplementary Fig. 19). This finding aligns with previous observations in S. mansoni, where SEM studies identified Sm-nrps+ cells with cilium-like protrusions, suggesting a conserved role of ciliated neurons in male-specific functions9. Collectively, these results indicate that these three types of neurons play different roles in male parasites. Considering the importance of the key enzyme coded by nrps in the N4.3 neuron that synthesizes the pheromone required for the female sexual development in schistosomes, we wondered whether there are essential genes maintaining the function of the N4.3 population.

Male lsamp is essential for female sexual development

To identify genes potentially driving the differentiation of N4.3 neurons, we analyzed the driver genes and pathways associated with the N4.3 neuron metastable state by correlating fate probabilities with gene expression using CellRank’s compute_lineage_drivers function38 (Supplementary Data 19 and 20). Among these, genes such as nrps (Sjc-0000814), ddc (Sjc-0006754, Aromatic L-amino acid decarboxylase), and lsamp (Sjc-0000373, limbic system-associated membrane protein) exhibited the highest correlations (Fig. 3d and Supplementary Data 19). Additionally, these three genes are specifically expressed in N4.3 neurons as indicated on the UMAP plot (Fig. 3e). The colocalization of nrps and ddc has been confirmed previously (Supplementary Fig. 16c). We performed double-FISH and found that lsamp is co-expressed with nrps and ddc in male N4.3 neurons (Fig. 3f and Supplementary Fig. 20). Additionally, lsamp and ddc also exhibited colocalization in the large cerebral ganglia (Fig. 3f). Nrps is a male-specific gene, and its expression is induced by pairing9. In S. japonicum, both ddc and lsamp displayed similar expression patterns to nrps, as reported in our previous study10 (Fig. 3g). We therefore reasoned that lsamp or ddc may play an important role in the N4.3 neuron that may further affect the male-female interaction.

By employing RNAi in vitro, we found the male worms treated with lsamp dsRNA lost their ability to initiate female sexual development (Fig. 4a, b). These female worms paired with lsamp RNAi males displayed reduced Fast Blue BB staining in their vitellaria and possessed poorly developed ovaries (Fig. 4b). In addition, they barely laid eggs (Fig. 4c). The qRT-PCR result indicated a strong reduction of the lsamp expression in the male worms following lsamp RNAi treatment (Fig. 4d). In contrast, treatment with ddc dsRNA did not affect female sexual development (Supplementary Fig. 21). The above results suggest that male lsamp is essential for female sexual development.

To determine whether the female phenotype resulting from the male lsamp knockdown was due to a reduction in the number of N4.3 neurons, we performed FISH experiments to quantify these cells. We selected the marker gene ddc to detect N4.3 neurons specifically. However, there was no significant difference in the lsamp silencing group compared to the control group (Fig. 4e).

Lsamp affects acetylated α-tubulin formation and vesicular transport within neurons

Since lsamp colocalized with nrps in N4.3 neurons (Fig. 3f and Supplementary Fig. 20), and nrps has been recognized as the functional gene encoding an enzyme to synthesize the male pheromone9, we therefore wondered whether the lsamp RNAi phenotype is caused by affecting the nrps expression. qRT-PCR results showed that the nrps expression was not affected in the lsamp knockdown males (Fig. 4f and Supplementary Fig. 22). According to Chen et al., the NRPS enzyme uses two simple precursors (β-alanine and tryptamine) to form the dipeptide pheromone BATT9. To explore why knocking down lsamp in male worms blocks female sexual development, we detected the levels of two precursors in worms (D20) using LC-MS/MS but found no changes between the gfp and lsamp RNAi groups (Fig. 4g, h). We then measured the amount of BATT, the product of the NRPS, that is made in males and can be released to the medium in vitro. Surprisingly, BATT contents were significantly reduced in the lsamp RNAi group, both in the male worms and the culture medium (Fig. 4i), even though the nrps level was not changed (Fig. 4f and Supplementary Fig. 22). These findings suggest that lsamp RNAi impairs the ability of male N4.3 neurons to synthesize BATT, consequently failing to stimulate female sexual development, not through affecting nrps expression.

To investigate the underlying mechanism of male lsamp in regulating BATT synthesis, we performed RNA-seq to overview the gene expression changes between paired lsamp RNAi and gfp RNAi males (D20). Silencing lsamp resulted in the downregulation of 50 genes and the upregulation of 19 genes (Supplementary Fig. 23a, b and Data 21 and 22). Our RNA-seq data confirmed that the expression of nrps and ddc remained unaffected by lsamp knockdown (Supplementary Fig. 23c), consistent with our previous results (Fig. 4f, and Supplementary Fig. 22). Since the upregulation of nrps is triggered by physical contact between males and females9, these results suggest that the knockdown of lsamp does not impair male worms’ ability to sense females. Furthermore, GO term analysis of the downregulated genes in lsamp RNAi males revealed significant enrichment in biological processes, including the generation of precursor metabolites and energy, and peptide biosynthesis, as well as molecular functions such as structural constituents of the ribosome and ATPase-coupled ion transmembrane transporter activity (Supplementary Fig. 23d, e and Data 23 and 24). KEGG pathway analysis of the downregulated genes showed that lsamp mainly affects pathways involved in neurodegenerative diseases such as Alzheimer’s disease and Parkinson’s disease (Fig. 5a and Supplementary Data 25), both of which are closely related to the disorder of microtubule acetylation and axon vesicular transport39–43. In addition, GO enrichment on the highly expressed genes in the N4.3 population indicated it as a specialized neuron that expressed genes highly related to cilium organization, cell projection assembly, and microtubule bundle formation (Fig. 3c and Supplementary Fig. 19). Previous studies in mice and nematodes have shown that lsamp homologues play a role in synapse formation, axon guidance, and neuronal stability44–48. We therefore hypothesized that lsamp may affect the processes of microtubule formation and axon vesicular transport in male schistosomes.

Acetylation of α-tubulin, a post-translational modification of microtubules, is associated with microtubule stabilization and regulates cell morphology in nerve fibers39,49. It plays a crucial role in various neuronal processes, including axonal formation, dendritic growth, and neuronal migration39. And acetylated α-tubulin plays a pivotal role in the transport of synaptic vesicles within neurons. It stabilizes microtubules, which serve as “tracks” for motor proteins, facilitating the efficient movement of vesicles and ensuring the proper delivery of neurotransmitters to the synapse. In addition to classical neurotransmitters like dopamine, serotonin, and GABA, synaptic vesicles also transport other small molecules, such as the precursor β-alanine for BATT50–52. To further investigate the specific function of lsamp in schistosomes, we examined the acetylation status of microtubules and synapsin levels in mature male worms. As shown in Fig. 5b, acetylated α-tubulin antibody labeling revealed crisscrossed nerve fibers in the mature male gynecophoric canal. After knocking down the lsamp, the nerve fibers were markedly eliminated without affecting the nrps expression (Figs. 4f and 5b and Supplementary Fig. 22). These data suggest that LSAMP is essential for nerve fiber maintenance, likely by mediating appropriate axon guidance and target recognition. The observed alterations in microtubules may represent a secondary morphological consequence of impaired axon maintenance. Using a synapsin antibody, we detected disrupted synaptic protein distribution along nerve axons (Fig. 5b). Comparative co-staining of acetylated α-tubulin and synapsin revealed severe degeneration of neural network functionality in the lsamp knockdown group (Fig. 5b).

We next examined the spatial relationship between these neural fibers and N4.3 neurons. By labeling N4.3 neurons with three markers—ddc, nrps, and lsamp—we observed that N4.3 neurons were closely associated with neural fibers marked by acetylated α-tubulin (Fig. 5c and Supplementary Fig. 24a, b). Synaptic proteins along neural axons also exhibited substantial colocalization with N4.3 neurons, consistent with possible vesicular trafficking between these structures (Fig. 5d and Supplementary Fig. 24c). Together, these findings indicate that the loss of lsamp in N4.3 neurons disrupts the integrity of the surrounding ventral neural network in male schistosomes.

Lsamp guides the formation of the ventral neural network in males

LSAMP (Limbic system-associated membrane protein) is a glycosylated cell adhesion molecule of the IgLON family47,53. It is a GPI-anchored protein characterized by three immunoglobulin (Ig) domains (Supplementary Fig. 25). Previous studies in mice and nematodes have demonstrated that lsamp plays essential roles in neuronal development and function, including promoting neuronal growth and guiding axon targeting44–48,54. To investigate the potential role of lsamp in schistosomes, we examined the development of the ventral neural fiber network in male worms from 14 to 26 dpi, covering key stages from male-female pairing to full sexual maturity. At 14 dpi, when males begin pairing with females, only a few elongated nerve fibers were detected within the male gynecophoric canal. As the males developed, a characteristic crisscrossed neural fiber network became evident by 22 dpi and was maintained thereafter (Fig. 5e). Notably, the formation of this neural network coincided with a marked upregulation of lsamp expression (Fig. 3g). This temporal correlation is consistent with a potential role for lsamp in axon guidance and neural connectivity.

Overall, our study has profiled the dynamic single-cell atlas of male and female S. japonicum during pairing and sexual development. From a male-specific N4.3 neuron population, we identified the lsamp gene, which is crucial for female reproductive development and operates independently of a previously reported key regulator in males, nrps. Our observations indicate that lsamp expression is associated with the formation of ventral neural fibers in the male gynecophoric canal, which appear to connect to the N4.3 neurons (Fig. 5c–e and Supplementary Fig. 24). Based on these findings, we propose a model in which this nerve network may serve as a conduit, potentially facilitating the axonal transport of the BATT precursor—likely β-alanine rather than tryptamine, given that N4.3 neurons intrinsically express ddc required for tryptamine synthesis—to the N4.3 neurons for pheromone production and subsequent stimulation of female reproductive development (Fig. 5f).

Methods
Ethics
All experiments related to animals were conducted in accordance with the guidelines for the Care and Use of Laboratory Animals of the Ministry of Science and Technology of the People’s Republic of China (2006398) and were approved by the Animal Care and Use Committee of Fudan University (Fudan IACUC 201802158S) to ensure ethical and responsible treatment of the animals.

S. japonicum parasites
The cercariae of S. japonicum were provided by the Chinese Center for Disease Control and Prevention. Single-sex infections of S. japonicum were obtained by infecting mice with male or female cercariae obtained from Oncomelania hupensis snails (Anhui strain) infected with single miracidia.

Animals
For all infection experiments in our study, six-week-old female C57BL/6J mice were purchased from Shanghai Jiesijie Laboratory Animal Co., Ltd (Shanghai, China). Female mice were selected to maintain consistent infection conditions and minimize stress related to aggressive behavior in social housing, which is crucial for reducing experimental variability in the schistosome-mouse model. Standard housing conditions were maintained on a 12-h light/12-h dark cycle, with ambient temperature at 22±1°C and relative humidity at 55±10%. All mice were fed a standard, certified commercial rodent diet and had ad libitum access to sterile drinking water (filtered and autoclaved).

Fluorescence-activated cell sorting
Twenty female C57BL/6J mice were infected percutaneously through the abdomen with 80–200 cercariae. At each time point (14, 18, 22, and 26 dpi), five mice were sacrificed for the recovery of schistosomes10. Approximately 300 (14 dpi, days post infection), 150 (18 dpi), 100 (22 dpi), and 70 (26 dpi) male and female S. japonicum worms were collected respectively, and processed into single-cell suspensions for scRNA-seq.

Worms were collected from mice by perfusing the hepatic portal vein with 4°C PBS (Cat# E607008, Sangon Biotech, Shanghai, China) and then rinsed twice in DMEM (5% FBS)10,55,65,66. Freshly perfused male and female worms were separated by incubation in a 0.25% solution of tricaine for approximately 5 min, except for the 14 dpi worms, which were processed as mixed-sex populations because their sex could not be distinguished morphologically under a light microscope. FACS sorting was performed as follows. Worms were washed twice in 8 mL PBS, then triturated for approximately 15 min in 4 mL 0.125% Trypsin-EDTA (in PBS) until the solution became turbid and no large pieces of worms were left21. The trypsin was inactivated by adding 8 mL 10% FBS DMEM. The dissociated worms were then centrifuged at 500×g for 10 min at 4°C. Next, the worm cells were resuspended in 1 mL of PBS with 10 μL of RQ1 DNAse (Cat# M6101, Promega, Madison, WI, USA) and incubated for 10 min at RT. The worm cells were centrifuged again at 500×g for 10 min at 4°C. Then the worm cells were resuspended in 1 mL of PBS containing Fixable Viability Dye eFluorTM 506 (1 μL/mL) (Cat# 65-0866-14, eBioscience, USA) and incubated for 30 min at 4°C in the dark. The worm cell suspension was centrifuged at 500×g for 10 min at 4°C, then washed twice in PBS, and then filtered through a 40 μm cell strainer into a 12×75 mm FACS tube. Filtered cells were then sorted on a MoFlo XDP High-Speed Cell Sorter (BECKMAN COULTER) with 405/460/488/532/561/640 nm lasers. Sorts were performed with a 100 μm nozzle, and cells were sorted into sorting media (0.04% BSA in PBS, pH=7.40). The sequential gating strategy was implemented using Summit Software (Beckman Coulter, v5.2) through four steps: (1) R1: initial events were gated on FSC-H versus SSC-H to identify the cell population and exclude debris; (2) R2: singlets were isolated using FSC-width versus SSC-log_height to exclude cell aggregates; (3) R3: a second doublet discrimination step was performed using SSC-width versus SSC-log_height to ensure a high-purity single-cell population; and (4) R4: viable cells were finally isolated by selecting the eFluor 506-negative population (identified via FL8-log_height versus SSC-log_height). Post-sort cell viability was further validated by trypan blue staining and manual counting under a light microscope to ensure high sample quality for downstream applications.

Single-cell RNA sequencing
FACS-sorted cells were centrifuged at 500×g for 10 min at 4°C and then resuspended in 0.04% BSA in PBS. Libraries were created using a Chromium Controller (10× Genomics) according to manufacturer guidelines and sequenced using a NovaSeq 6000 (Illumina) to generate 150-bp paired-end reads. Raw data were processed and mapped to the S. japonicum genome (v3) using the Cell Ranger 3.0.0 (10× Genomics) pipeline.

Single-cell RNA-seq data processing
The Seurat package (v3.1.5) was used to analyze the raw values of the matrix as follows67. Cells with greater than 30% mitochondrial reads or fewer than 500 genes, were excluded from the analysis. Mitochondrial genes were identified with the prefix “^ND|^C|^ATP”. Each of the 7 individual datasets was scaled and transformed, and variable genes were identified using the SCTransform function with the parameter “variable.features.n=3000”. Integration anchors across samples were identified using the FindIntegrationAnchors function, and the datasets were subsequently integrated using IntegrateData from the Seurat R package. This sample-level integration approach minimizes batch effects while preserving biological variation across samples. Then we performed principal component analysis using variable genes and used the first 100 principal components (PCs) to perform UMAP to embed the dataset into two dimensions. For clustering, we identified 76 clusters using the FindNeighbors and FindClusters functions with the parameter “dims=1:100, resolution=2”.

To determine cellular identity, we used the function FindAllMarkers with the parameter “test.use=“wilcox”, only.pos=T, min.pct=0.25, logfc.threshold=0.25” to find DE genes for each cluster and compared to the known marker genes from the previous datasets18,19,21,68. GO (Gene Ontology) enrichment analysis was performed with the R package clusterProfiler69. The p-value was corrected for multiple hypothesis testing with the Benjamini–Hochberg false-discovery rate procedure (adjusted p-value). GO terms with corrected p-value < 0.05 were considered significant. R packages ggplot2 (https://ggplot2.tidyverse.org) and tidyverse (https://www.tidyverse.org) were used to generate figures and intermediate data preprocessing.

The density plots of cells from each sample over UMAP embedding were generated using a modified treecor_celldensityplot function from TreeCorTreat package (https://github.com/byzhang23/TreeCorTreat).

Differential abundance analysis
We tested for differences in cell-type abundances between female and male worms using the miloR package (v1.4.0)70, which is a method for differential abundance analysis on KNN graphs from single-cell datasets. Briefly, we computed the k-nearest neighbor graph of cells in the integrated datasets based on latent embedding. Subsequently, cell neighborhoods were computed using miloR’s makeNhoods function with parameter “k=25, prop=0.2”. To leverage the variation in the number of cells between 18, 22, and 26 dpi samples, we counted the cells belonging to each sample in each neighborhood using the countCells function. We assigned to each neighborhood a cell-type label based on the majority voting of the cells belonging to that neighborhood. We then used the testNhoods function to test for enrichment of cells from the different sex datasets. Neighborhoods were enriched if the spatial FDR < 0.2 and |log FC| > 1.

Subclustering analysis
To investigate cellular heterogeneity within reproductive system-associated and neuronal populations, we performed subclustering analysis using the Seurat package (v3.1.5)67. Reproductive system-associated cells (including neoblasts, germ cells, vitellaria, and Mehlis’ gland cells) and neuronal cells (including neural progenitor cells, neural precursor cells, nonciliated neurons, and ciliated neurons) were subset using the SubsetData function. Principal component analysis, dimensionality reduction, and clustering were then performed separately on each subset. We used 20 PCs and a resolution of 3 for reproductive system-associated cells, and 30 PCs with a resolution of 1 for neuronal cells, based on the variance explained in the Principal Component Analysis (PCA) elbow plots and the robustness of clustering results. Subclusters were annotated based on lineage-specific transcriptional signatures and established marker genes.

Trajectory analyses
To infer potential developmental trajectories within the reproductive-related or neuronal lineages, trajectory inference analysis was performed using the R package Slingshot (v2.0.1)30. The vitelline, germline, and neuronal cell lineages identified in the integrated dataset were subsetted, and recomputed a UMAP embedding. Slingshot was then applied with default parameters, using neural progenitor cells, S1-0 cells, or GSCs as the starting clusters to infer cell trajectories and calculate pseudotime. The resulting trajectories were mapped onto the UMAP projection for visualization. To evaluate how pseudotime inferred from Slingshot corresponded to the actual sample collection timepoints, pseudotime distributions across samples were visualized using density plots.

In addition, to identify genes dynamically changing along vitelline and germline pseudotime trajectories, tradeSeq34 analysis was carried out on the Slingshot-derived trajectories. For each dataset, the “associationTest” function was used to identify genes whose expression patterns were associated with pseudotime, with a knots value of 6 for vitelline and 7 for germline lineages. A p-value cutoff of 1e−5 was applied to select significantly changing genes, and a more stringent cutoff (1e−10) was used to define a refined set of dynamic genes in both female and male datasets along the germline lineages. Heatmaps of these genes were generated using the R package ClusterGVis (v0.1.4) (https://github.com/junjunlab/ClusterGVis). To determine the optimal number of gene expression clusters, the “getClusters” function in ClusterGVis was applied, which implements the Elbow method to evaluate within-cluster sum of squares across a range of cluster numbers.

RNA velocity and cellrank analysis
To further characterize differentiation trajectories leading to distinct cell populations, we performed RNA Velocity analysis. In brief, we generated loom files of the spliced and unspliced RNA matrices using velocyto (v 0.17.17)30 with our CellRanger output. Then we applied scVelo (v 0.2.2) to estimate transient cell states and RNA velocities using a generalized dynamical model following the recommended workflow. Additionally, diffusion pseudotime71 was computed using scVelo to provide an alternative measure of cellular progression along developmental trajectories based on transcriptomic similarity and RNA velocity.

To further investigate cell fate decision, we used CellRank package (v 1.0.0)38 as described in (https://cellrank.readthedocs.io/en/stable/index.html) to identify initial and terminal states and characterize cell fate probabilities. Specifically, to infer putative driver genes for any of the trajectories, we correlate expression values with fate probabilities using the “cellrank.tl.lineage_drivers” function with the parameter “method=Fischer”.

CytoTRACE analysis
To evaluate the differentiation potential of cells within the vitelline lineages, we applied CytoTRACE33 with default parameters. The UMAP embedding was used to visualize CytoTRACE scores across cell clusters. A higher CytoTRACE score indicates greater stemness and a higher potential for differentiation into mature cell states.

RNA interference (RNAi)
Oligo sequences used to generate dsRNA templates and qPCR primers are listed in Supplementary Data 26. PCR primer synthesis and DNA sequencing were performed by Beijing Tsingke Biotech Co., Ltd. (Tianjin, China). Double-stranded RNA (dsRNA) was synthesized using the T7 High Yield RNA Synthesis Kit (Cat# 10623ES60, Yeasen, Shanghai, China) following the manufacturer’s protocols. Briefly, DNA templates for in vitro transcription were PCR-amplified from S. japonicum cDNA using gene-specific primers flanked with T7 promoter sequences (5’-TAATACGACTCACTATAGGGAGA-3’) at the 5’ ends. The PCR products were purified using a TIANgel Purification Kit (Cat# DP219, TIANGEN, Beijing, China) and utilized as templates for the transcription reaction. The mixture was incubated at 37°C for at least 4 h. The synthesized dsRNA was subsequently purified, resuspended in nuclease-free water, and quantified using a NanoDrop ND-2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The integrity and quality of the dsRNA were verified by 1% agarose gel electrophoresis. The dsRNA of the gfp (green fluorescent protein) fragment was used as a negative control for all RNAi experiments55. Day 0 represents the first day of the experiment. Males were treated with 30 μg/mL dsRNA (on Day 0/2/4/6 with fresh medium) for one week in BM169 and then paired with virgin female parasites from Day 8 to Day 27 in ABC169 (fresh medium was replaced every other day). The eggs produced by the worms are retained when changing the medium. On Day 28, the egg number and the number of female parasites were counted to calculate egg/female/day.

Bulk RNA-sequencing and data analysis
We collected male worms from the gfp RNAi and lsamp RNAi groups on Day 20 for bulk RNA sequencing, with four biological replicates per group. Total RNA extraction, library construction, and RNA sequencing were performed by the Novogene Co., Ltd (Beijing, China). Briefly, RNA was extracted from worms using the AG RNAex Pro Reagent (Cat# AG21101, Accurate Biotechnology (Hunan) Co., Ltd., Changsha, China). The integrity and quality of total RNA were assessed using a Nanodrop ND-2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Pooled RNA was used to construct a library using Illumina TruSeq™ RNA sample prep Kit v2 (Illumina, USA) and the Ribo-Zero Plus rRNA Depletion Kit (Illumina, USA) for rRNA depletion, following the TruSeq RNA Sample Preparation Guide. The library was subsequently sequenced on the Illumina NovaSeq 6000 platforms to obtain 150 bp paired-end (PE) reads.

QC (Quality control) of the raw sequencing data was performed using the FASTQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low-quality reads and adapter sequences were trimmed using the Fastp tool (v0.20.1) (parameters: -q 15 -u 40 -n 5 -l 15)72. The clean reads were mapped to the chromosome-level S. japonicum reference genome (SjV3) using HISAT2 (v 2.1.0) with default parameters73. The gene-level assignment was then performed to estimate the transcript abundance using featureCounts (1.6.4)74. PCA was performed using the prcomp function in the R stats package (v3.6.0). Hierarchical clustering analysis was performed with pheatmap (https://cran.r-project.org/web/packages/pheatmap/index.html). Differential expression analysis was performed using DESeq275, and the differentially expressed genes (DEGs) were defined with adjusted p-value < 0.05 and |log FC| > 2. To identify the possible functions and pathways of DEGs, GO terms enrichment analysis and Gene Set Enrichment Analysis were performed using the clusterProfiler package69.

Parasite labeling and imaging
Colorimetric and fluorescence in situ hybridization were performed as follows. Riboprobes were synthesized using the Fluorescein RNA Labeling Mix (Cat# 11685619910, Sigma-Aldrich, St. Louis, MO, USA) or DIG RNA Labeling Mix (Cat# 11277073910, Sigma-Aldrich, St. Louis, MO, USA), with templates generated using primers listed in Supplementary Data 4. Typically, all probes were applied at a concentration of 100 ng/mL in the hybridization buffer, with the concentrations of certain probes adjusted to optimize localization performance. For WISH, anti-digoxigenin-AP (1:2000, Cat# 11093274910, Lot# 84080920, Sigma-Aldrich, St. Louis, MO, USA) was used, followed by color development with NBT/BCIP (Cat# 11681451001, Sigma-Aldrich, St. Louis, MO, USA). For FISH, anti-digoxigenin-POD (1:1000, Sigma-Aldrich, Cat# 11207733910, Lot# 77798000, St. Louis, MO, USA) or anti-fluorescein-POD (1:1000, Sigma-Aldrich, Cat# 11426346910, Lot# 57545020, St. Louis, MO, USA) were used. Signal amplification was performed using FITC-tyramide (Cat# B40953, Thermo Fisher Scientific, Waltham, MA, USA) or TAMRA-tyramide (Cat# B40957, Thermo Fisher Scientific, Waltham, MA, USA)23,66.

The immunofluorescence assays for detecting Synapsin levels in the parasite were performed using mouse monoclonal anti-Synapsin (1:100, Clone 3C11, Cat# 3C11, Lot# 10/1/21, Developmental Studies Hybridoma Bank, DSHB, Iowa City, IA, USA), and for detecting acetylated α-tubulin levels using mouse monoclonal anti-acetylated α-tubulin (1:500, Clone 6-11B-1, Cat# sc-23950, Lot# G2522, Santa Cruz, Dallas, Texas, USA)76,77. Subsequently, samples were incubated with Goat Anti-Mouse IgG H&L (Alexa Fluor® 488) secondary antibody (1:1000, Cat# ab150113, Lot# 1004120-39, Abcam, Cambridge, UK) or Goat Anti-Mouse IgG H&L (Alexa Fluor® 594) secondary antibody (1:1000, Cat# ab150116, Lot# GR3402482-1, Abcam, Cambridge, UK) for 1 h at room temperature.

For females collected on Day 28, vitellaria were labeled with 1% Fast Blue BB (Cat# F3378, Sigma-Aldrich, St. Louis, MO, USA) in PBSTx for 1 min at room temperature24,55.

All labeled parasites were counterstained with DAPI (1 μg/mL), cleared in 80% glycerol, and mounted on slides with 80% glycerol. Confocal imaging of fluorescently labeled samples was performed on Nikon A1 Laser Scanning Confocal Microscope or Olympus FV3000 Laser Scanning Confocal Microscope. Brightfield images were acquired on Olympus BX51 Microscope and Zeiss AxioZoom V16 Microscope.

qRT-PCR (quantitative real-time PCR)
All qPCR reactions were performed on a LightCycler® 96 Instrument (Roche, Switzerland) using 2× SYBR green qPCR master mix (Cat# 11201ES50, Yeasen, Shanghai, China) according to the manufacturer’s instructions. Each 20 μL qRT-PCR reaction mixture comprised 2 μL of cDNA (1:4), 10 μL 2× SYBR green master, 0.8 μL (5 μM) of each primer, and 6.4 μL ddH2O. The qRT-PCR cycle parameters were as follows: 95°C for 3 min, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s; melt curve analysis ranged from 60°C to 95°C to ensure that the specific product was amplified in each reaction. The 2−ΔΔCt method was used to calculate the relative fold change of the differentially expressed transcripts.

BATT synthesis and detection
The synthesis of BATT was conducted by the Chinese Center for Disease Control and Prevention. Concentrated stock solutions (100 mM) of BATT were prepared in sterile DMSO and stored at −20°C. LC-MS/MS was employed for the detection of BATT in parasites (thirty male worms) and medium samples collected on Day 20. BATT detection was performed on a triple quadrupole mass spectrometer operating in multiple reaction monitoring employing negative ESI interface using QTRAP® 6500+ LC-MS/MS System (Sciex, Framingham, MA, USA), with the transitions set at m/z 232.4→m/z 143.7, m/z 232.4→m/z 185.1, m/z 232.4→m/z 89.1, respectively9,55.

Quantification and statistical analysis
Quantification and statistical analysis were performed using GraphPad Prism software (v9.0, San Diego, CA, USA). The normality of datasets with an adequate sample size (n ≥ 5) was assessed using the D’Agostino-Pearson test (Supplementary Data 27). For data conforming to a normal distribution or with a sample size (n < 5), unpaired two-tailed Student’s t-tests were applied, and the results are presented as the mean ± SEM (Standard Error of the Mean). For data that did not pass the normality test, the non-parametric Mann–Whitney U test was applied (Supplementary Data 28). Statistical significance was defined as *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001; ns, not significant.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Wei et al describe the first single cell sequencing atlas for S. japonicum—the human-infective Asian blood fluke. This is a really impressive study; with approx. 100k cells it’s the largest single cell atlas produced in schistosome research—creating a rich new resource for understanding schistosome development.\nThe high cell coverage and inclusion of intermediate stages of development has enabled the authors to produce high resolution data on the sexual-tissue development, uncover distinct neural lineages and reveal novel biological insights into sexual signalling/development in schistosomes.\nA key highlight of the work is the discovery of a male-specific neural circuit that directly mediates female development. The discovery builds upon previous research in S. mansoni demonstrating that nrps in males produces the pheromone BATT, which is essential for male-stimulated female development. The authors report a neural cell type (“delta-3”) that expresses nrps, and the S. japonicum homologue of lsamp, which expresses a neural surface protein in other species. The authors functionally characterise lsamp, confirming its expression in delta-3 cells and showing that it is required for BATT production but affects neither delta-3 development nor nrps expression. They further show that lsamp expression is closely colocated with the axons and synaptic junctions of a ventral neural network. Knocking down lsamp disrupts this network. The authors suggest that lsamp is required for males to produce or maintain neural fibrils.\nOverall, this is exciting work that deepens our understanding of sexual signalling and development in these obligate dioecious parasites. The work is highly significant within the parasitology field and will be of broader interest for developmental biologists.\nIt is a large body of work and, in places, the manuscript is currently hard to navigate. There are some sections that may benefit from simplifying/streamlining, including some analyses that didn’t seem strictly necessary. We also felt that there were aspects of the functional characterisation, where the authors have over-reached in their interpretations.\n*Specific comments - Major points*\nTo validate clusters (S. Fig 3), some of the hybridization patterns are unclear. Presenting evidence for the more diffuse tissues is difficult. Nevertheless, ctsb, tpm2, syn2 do look very similar.\nDescribing the proportions of cell types is quite prominent in the manuscript (incl. Figs 1 c and d), but how informative is this analysis, beyond an extremely crude look? Sorting, library prep etc. creates bias for certain cell-types (e.g. more round cells like stem cells are favoured by sorting) and can distort recovered cell numbers. In previously published works, actual cell abundance vs cell captured by the scRNAseq process differ wildly, affecting differing cell types unequally. Some of the numbers make sense, like the expansion of vitellaria in females, but others make less sense like the much higher abundance of flame cells than gut cells. There’s a vague/superficial statement, which could be removed to save space, about there being ~20% neural and neoblast cell types when the numbers actually differ quite markedly between timepoints. Abundance across the whole dataset is shown in S Fig 5; yes it does broadly confirm some expectations but it doesn’t really shed any new light. And some of the findings seem counter-intuitive: like the proportion of muscle cells from males.\nThe authors performed sub cluster analysis on sexual tissues and neoblast populations, from which they define 17 sub-types. The rationale for their naming scheme is not clear and abbreviations need to be explained. A discontinuous series from the Greek alphabet, plus some numbering (i.e. alpha, beta, gamma-1, gamma-2, delta, epsilon, kappa, nu) is used. Previously, a similar scheme has been used in published descriptions of S. mansoni stem cells (see doi.org/10.1016/j.pt.2021.03.005). It’s not clear if published delta, epsilon, and kappa sub-types are equivalent to those reported here. Distinguishing between 17 populations in the figures (S Fig 6a and Fig 2) is also challenging. Perhaps a numerical key could be added to help readability?\nRNA velocity analysis was performed on the sex related subsets, revealing trajectories (Fig 2b) that made sense like the GSC -> GSC progeny in both sexes and the S1 vitellocytes ->->late vitellocyte. There are also a lot of arrows converging on the centre (SSC beta?) that aren’t mentioned. Further subsetting (Fig 2c) reveals trajectories from GSC to gametes in males but the equivalent analysis is not at all convincing in females (contrary to the main text). In Fig 2c the four panels show a lot of redundant information and could be reduced to two panels.\nThere’s quite a long section on clustering male and female specific genes. The authors reveal a cluster of genes with opposite temporal patterns between the sexes. It’s also interesting to note that for the majority of sex-specific genes, the temporal pattern is different between the sexes, with male genes reducing and female genes increasing. However, some of the methods and figures in this long section are difficult to navigate and it’s disappointing that the section is left quite open-ended with no further analysis of the 47 C1 genes. The overall conclusion is that the process is complex and involves modules of genes. Without the composition of the modules being described, it doesn’t really provide anything new.\nThe main text (L201+) and legend for Fig 7 really need signposts to the relevant section of Methods because it’s not immediately obvious that trajectories are the input for the sex-specific analyses. Even after the gam R description is found, it’s still unclear exactly what the inputs are and how the lists of 693, 750 and 156 relate — are they mutually exclusive? If not, overlaps should be indicated. In Fig 7b, pseudotime=0 is presented in the centre, so that one panel reads right-to-left and the other reads left-to-right; that seems counter-intuitive. How were the datasets clustered? It looks like the cluster number was set to 5 in each case, but it’s not clear why, when several clusters could have been further divided and some possibly merged. What is the dotted line on each cladogram?\nAfter the module analysis, the narrative returned to trajectory analysis using RNA velocity. The authors were unable to connect stem cell development to vitellaria development but were able to show trajectory from S1 vitellocyte to S1 progeny, with the latter marked by esg-2a expression. The paragraph finishes with a rather bland statement that the work aligned with previous findings and provides new insights. Vf1 expression was known, but the esg-2a expression later in differentiation is new. What else switched on later? Having established a convincing trajectory, was it not possible to mine further information on genes expressed early vs late? Likewise, the analysis of Mehlis gland is interesting but only two strong markers are shown in Supplementary Figure 6. Could no other markers be found?\nThe authors subset the neuronal data and recluster into a highly connected UMAP showing at least 5 lineages. This is an important part of the reference resource created by the project. However, the authors should maybe avoid using nomenclature based on Greek letters, given that it is already used for stem cells, including within the same paper. How about just numbering them? Over time, more functional names will be applied but, for now, it might suffice to name them as neural lineage 1, etc. Or maybe use N1, etc? The more granular names like delta-3 would then become N4.3.\nThe authors use Slingshot to define lineages and derive pseudotime from their neural UMAP. How well does the pseudotime derived from these Slingshot trajectories agree with the sample-time metadata? (does a UMAP coloured by sample show differential enrichment for earlier or later parts of the trajectories?). The subsequent RNA velocity analysis agrees well with their pseudotime predictions, and Cell Rank is used to define terminal states (fully differentiated cells?). Similar to the earlier section on gamete development, the signposting to the methods is not very clear in the text - it’s not immediately clear the purpose of each analysis until reading further into the paper. Cell Rank is not mentioned in the legend to S Fig 10 and “dpt pseudotime” is not explained. The PAGA analysis, resulting in a network of pie charts (Fig. S10e and 11e), isn’t clear. What exactly is this showing? Is it really needed?\nRegarding the section on male-specific neurons: in schistosomes, a large number of genes are present in sex-specific regions (~900 in S. mansoni) and, given their incomplete dosage compensation between sexes, those sex-specific genes can skew the relative abundance of all other transcripts. Did the authors check whether the clustering was influenced by genes in those regions? If it is, there could be distinct clusters that are expressing the same genes but just differ in the expression levels across their sex-specific genes.\nBeta-2, gamma-2 and delta-3 neurons are described as ‘male-specific’ (Line 291), but the markers for two of these populations (beta-2 and gamma-2) were then detected in the heads of females (Line 304). Was this a sensitivity issue of sc-RNASeq vs. WISH? Is it really appropriate to describe beta-2 and gamma-2 neurons as ‘male-specific’? They have also included beta-2 and gamma-2 neurons on the ‘Female anterior’ schematic in Fig. S16, suggesting they are not male-specific in the authors’ model. This needs some clarification in the text.\nNavigating through some of the results on markers for the male specific clusters is confusing due to the different types of terminology. The markers are initially referred to by their cell type and systematic gene id (in the text and Supplementary Figure 13) but in the WISH validation gene names like nrps and ddc appear but without their systematic names. I think it would be simpler to read, if the more human-readable gene names were used in the main text where possible (with systematic IDs in parentheses). In the legend for Supplementary Figure 13, I had to work out some names, like “ddc” = Sjc-0006754 by subtraction. Similarly, I needed to cross-reference Supplementary Fig 14 to learn that Sjc-0000814 = “nrps”.\nIt’s interesting that lsamp was not identified as a marker for delta-3 cells. It had the third-highest correlation across the delta-3 trajectory, but what was its expression across the data set. Was its expression just too low to be picked as a marker by Seurat?\nThe GO analysis of genes affected by lsamp knockdown isn’t particularly informative. The KEGG analysis (Fig 5) and functional information gleaned from other systems enable the authors to hypothesise about the role of lsamp. The authors found fewer nerve fibres in RNAi knockdown but suggested that LSAMP is therefore “essential for nerve fiber growth”. Given that the experiment was conducted on adult worms, would the nerves be fully grown? In which case, the knockdown presumably destabilised the nerve network rather than stopped growth?\nThe colocalisation experiments seem to show lsamp having a similar relationship as nrps and ddc to acetylated tubulin and synapsin but it’s difficult to interpret when such a small part of the worm is shown. Can lower magnification images be included for orientation? The colocalisation isn’t perfect, so it’s difficult to see how this suggests “active vesicular transport between these axons and delta-3 cells”. Based on expression timing, the authors speculate that “lsamp facilitates axon growth and attachment to delta-3” but evidence isn’t very strong for this. They also say that “lsamp is essential for guiding the formation of [nerve fibers]” but this is also not supported by the evidence. We suggest that the speculation be moved to the discussion, to make it much clearer where the evidence gaps remain. In the discussion L530–531, the authors state that “LSAMP functions as a membrane protein”. Although quite likely, it does remain to be shown, so appropriate caveats need to be inserted.\n*Minor points*\nL147–148: The statement about differentiated neuron populations expanding after D14 only seems to apply to males. In females it goes up then down. But is this a real effect or is this due to the differing ability to capture certain cells as the worms in each sample are morphologically different?\nL149 (word choice): The neural progenitor cells presumably don’t “degenerate”, but differentiate into neurons.\nL189 esg-2a is identified as a new marker. This gene encodes eggshell protein (PMID:2172818).\nL210–218 How were the example genes in Fig 7 b-d picked? It looks arbitrary.\nL300. The authors state that all three neurons are primarily located in the gynecophoric canal. According to Fig 3b, this statement is only actually true for delta-3. For the other two neuron types, they are localised in the canal but are expressed more in other areas.\nL310–315 After the male specific neurons part, the narrative shifts to validating the other neural types (Supplementary Figure 15 and 16), before returning to male specific neurons. This seems to interrupt the flow. Maybe this could be moved to earlier in the text, before the male specific part.\nL318–320 “These enriched genes indicated that beta-2 and gamma-2 populations exhibited typical…In contrast, gamma-2 showed…” This isn’t a contrast, you could say “Moreover, gamma-2 showed…” instead.\nL347: This sentence looks like it needs re-phrasing. Do the authors mean “...worms treated with lsamp dsRNA but not ddc dsRNA lost their ability…”\nL372: “ …suggest that lsamp impairs the ability of male delta-3 neurons to synthesise BATT”. We think the sentence, as written, means the opposite of what the authors intended. It should say “lsamp RNAi”\nL409: “nerve fibers were markedly eliminated without affecting nrps expression (Fig 5b)”. The figure doesn’t show nrps expression.\nL410: “The data suggest LSAMP is essential for nerve fiber growth by stabilizing microtubules”. Yes fiber growth is inhibited but you can really only speculate on the mechanism at this point.\nL444: The nerve network functions as an interstate. Maybe use the word “conduit” instead of “interstate”.\nL530–531: “As a glycosylated cell adhesion molecule of the IgLON family, LSAMP functions as a membrane protein”. This is quite an assertive statement about the function of the LSAMP protein. However, I don’t see evidence that the authors have tested these functional predictions in the text, so it would be good to make that clear in this statement.\nFig. 2b The annotation of “late Vi” on UMAP seems to occupy the position annotated as “Mehlis gland” on the same UMAP when used in Supplementary Figure 6. Which one is it? The trajectory looks like it runs S1 vitellocyte -> S1 progeny -> Early Vi -> Late Vi -> Mehlis\nFig 3a and Supp Fig 12a. These figures are very similar and could be combined e.g keep the colours but highlight the three male-specific populations with boxes.\nFig 3b and d use consistent gene naming\nFig 3d is very small given its importance to the lsamp story\nFig 5 (legend) “tight adherence” between lsamp, delta-3 markers, synapsin and acetylated tubulin in axons is incorrect. It shows colocalisation.\nIn Supplementary Figure 12b, the third panel (% composition of clusters by time) shows the three male specific neurons nicely. The other two panels are less useful and panel 2 is mislabelled (should be M not F).\nFig. S12C – on enriched neural cell-types graph (on right) the red arrow for ‘Enriched in female’ would be more intuitive if facing the other way in the legend\nSupp Fig. 22. Is it really needed — the information has already been presented elsewhere in the manuscript.\nSupp Fig. 24b – what do red and blue signify?\n*Methods*\nFor the scRNAseq, QC details are quite sparse. The >30% mitochondrial reads filter seems high. Were there a lot of cells with high %mito?\nNo methods included for the sex and neural subsets of 17 and 19 cell-type populations\nL611–612 ‘All samples were combined using the functions “FingIntegrationAnchor” and “IntegrateData”’. Which features were the datasets integrated on? Timepoint? Sex? Both?\n- Professor Matt Berriman & Dr Frances Blow", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "In this manuscript, You et al. identified a key gene required for male schistosome-triggered female reproductive maturation in Schistosoma japonicum. Specifically, they applied single-cell transcriptomics, RNA velocity and pseudotime analysis, and RNA in situ hybridization to reveal the cellular and molecular dynamics underlying this process. By focusing on four key time points, they uncovered putative developmental trajectories of male and female cells, focusing on germ cells and various neuronal cell types. Notably, the authors found that one of the three male-specific neuron subtypes, delta-3 neurons, relies on lsamp for normal axonal growth and ventral nerve network formation in the gynecophoric canal of male parasites. Through RNA interference experiments, they demonstrated that defects in this male nerve network due to lsamp knockdown reduce pheromone BATT production and subsequently impair female sexual development. Overall, this study provided a comprehensive single-cell transcriptomic dataset of both male and female S. japonicum during the important male-female pairing process and discovered a male nerve network critical for female reproductive development. The paper is well-written with high-quality figures that present solid experimental evidence. Beyond the important discovery of the role of lsamp in BATT signaling and its contribution to understanding this new signaling molecule, this work provides a useful resource for understanding schistosome development and reproductive biology. Thus, it is well suited for publication here. However, there are a few places where results are over-interpreted or alternative interpretations are not considered. Some suggestions for issues the authors should consider addressing in a revised manuscript are listed below.\n1. Some of the language used to describe cell states identified via single-cell sequencing overinterprets what has actually been learned. For example, male gametes are sperm and female gametes are oocytes. Since transcripts defining the “gamete” clusters are not presented in high-magnification FISH images, it is not clear whether these clusters reflect mature gametes or earlier stages along the differentiation pathway. Similarly, the cluster referred to as neural stem cells has not been validated in a way that allows them to be called stem cells. No evidence is presented that they self-renew, one of the defining features of stem cells. This could be a large pool of neural progenitors. In the absence of evidence that these cells are actual gametes or neural stem cells, the authors should soften their language.\n2. The authors conclude that “LSAMP is essential for nerve fiber growth by stabilizing microtubules” (line 410). Since lsamp does not appear to be expressed in the affected neurons, it is unclear how it would act to stabilize their microtubules. Based on what is known about LSAMP function, it seems much more likely that the affected neurons are not able to extend their axons to the appropriate targets (or maintain them there), so any impact on microtubules would be secondary (i.e., there are no axonal microtubules because there may be no axons, not because the microtubules are less stable).\n3. Understanding the phenotype of lsamp RNAi parasites requires the authors to provide much more detail about the RNAi protocol and the subsequent phenotypic analysis. As written, it is unclear when the male parasites treated with dsRNA were initially harvested from the host. If they were harvested at 4 weeks or later, then lsamp is required to maintain the structure of the neural network (because it would have developed already). If they were harvested earlier, it is important to know when, so the reader can place the phenotype in the context of the developmental time course shown in Fig. 5e. The authors should also clarify the time points used for the subsequent analysis of the lsamp RNAi phenotypes. Were all assays conducted at 21 days post-pairing, or were some phenotypes examined at other time points? The graph in Supp. Fig. 22 overinterprets what is known. It is not clear that LSAMP guides the ventral neural network. Their RNAi experiment shows a single undefined endpoint from an undefined starting point.\n4. The conclusion that the affected axons transport BATT precursors is not supported by the data; other interpretations are possible. The authors begin their argument with the claim that “synaptic vesicles also transport other small molecules, such as the precursors for BATT ( -alanine and tryptamine)38, 39.” As far as I could tell, neither of these references mentions  -alanine or tryptamine. It remains an open question whether these compounds are transported via synaptic vesicles, released via solute carrier proteins, or if some other mechanism is at work. Furthermore, although the authors note that ddc (aadc) is required for ovary formation in planarians, they do not mention that it is also required for BATT production and is co-expressed with nrps in planarians, where it is likely producing tryptamine from tryptophan in the BATT-producing cells (PMID: 38889152). The lack of a phenotype in the ddc RNAi parasites could easily be due to insufficient knockdown (Supp. Fig. 18). Because ddc (aadc) is co-expressed with nrps in schistosomes, it seems most likely that the tryptamine is produced locally, rather than transported. If not, why is ddc (aadc) expressed in these cells?\nMinor comments\n1. To help readability and the interpretation of the figures and text, wherever possible, informative gene names should be used in place of the Sjc numbers. To provide just one example, the in situ in Fig. 3b used Sjc numbers that are not shown elsewhere in the figure and are marked in colors that do not correspond to any others in the figure, so it was unnecessarily difficult to keep track of which neuronal subclass they were marking.\n2. The authors may wish to change the nomenclature of their neuronal subtypes to avoid confusion. Stem cell subtypes in schistosomes have already been defined as delta, using the Greek letter (δ) to define them (PMID: 29988015), so this may add unnecessary confusion to the literature. Furthermore, because delta genes are conserved across animals and play critical roles in the development of many organ systems, defining subclasses as delta-# (i.e., without the Greek lettering) may lead readers to expect that they express a specific delta gene.\n3. Line 148-9: “This finding indicates that neural progenitor cells are active during the juvenile stage but gradually degenerate in adulthood.” There is no evidence that these cells degenerate. It seems more likely that they simply generate their progeny.\n4. Since most readers of this journal will be unfamiliar with details of schistosome biology, instead of noting that they find female germ cell markers in males with no additional context (Supp. Fig. 6d), the authors should point out that hermaphroditism is occasionally observed in schistosomes and provide appropriate citations. It is not clear from the images shown in Supp. Fig. 6d that female germ cells are found in the testes because the testes are not visible at the magnification shown. The authors would need to show higher magnification views that show these cells in the testes to support this claim.\n5. Figure 1C, S5, and Data S5 showed changes in cell type compositions across different time points. These measurements are relative, meaning a significant expansion of one cell type may cause a noticeable reduction of other cell types in proportion. Please provide a plot or table with absolute cell counts.\n6. Line 170: eled is referred to as a germline cell marker, but it is also expressed in vitellaria and neoblasts in the posterior growth zone in juveniles. Can the authors comment on the eled+ cells in males that appear to be in vitellocytes at 22dpi (Fig. 1d)?\n7. In Figure 3E, lsamp and ddc expression seems to be partial in the delta-3 neuron cluster. Is this a gene coverage artifact from the scRNAseq or a reflection of the heterogeneity of delta-3 neurons? Quantifying their double FISH studies should be able to clarify this question.\n8. Unpaired two-tailed parametric t-tests were widely used for the statistics. Have all the data passed a normality test?\n9. Besides the top three genes, nrps, lsamp, and ddc, what other genes and associated pathways are specific to delta-3 neurons? Can the authors elaborate on this further?\n10. For IF images, please add quantification and experimental replications in the figure legends.\n11. Line 347: “worms treated with lsamp dsRNA instead of ddc dsRNA lost their ability…” This statement would be more effective omitting mention of ddc dsRNA here. First, state what happens after lsamp dsRNA treatment, then say that ddc dsRNA treatment did not affect female sexual development.\n12. Line 398: “Acetylated alpha-tubulin is a protein that stabilizes microtubules” is an awkward statement. Acetylation of alpha-tubulin (a subunit of the microtubule polymer) is a post-translational modification of microtubules (on alpha-tubulin) associated with stable microtubules.\n13. Line 469: “We… characterize gene regulatory programs…” This work did not characterize gene regulatory programs; it identified changes in gene expression associated with dimorphic sexual development.\n14. Line 507: “…no active neurogenesis at the adult stage.” Single-cell sequencing does not suffice to make this claim. The population of neural progenitors may just be smaller than can be detected, given the number of cells sequenced here.\n15. “Gender” is a social construct, so the authors should omit that word and refer simply to the parasite’s sex.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We would like to thank the reviewers for their thoughtful and helpful critiques of our work. The work is greatly improved upon the incorporation of their suggestions.\nReviewer #1:\nWei et al describe the first single cell sequencing atlas for S. japonicum—the human-infective Asian blood fluke. This is a really impressive study; with approx. 100k cells it’s the largest single cell atlas produced in schistosome research—creating a rich new resource for understanding schistosome development. The high cell coverage and inclusion of intermediate stages of development has enabled the authors to produce high resolution data on the sexual-tissue development, uncover distinct neural lineages and reveal novel biological insights into sexual signalling/development in schistosomes.\nA key highlight of the work is the discovery of a male-specific neural circuit that directly mediates female development. The discovery builds upon previous research in S. mansoni demonstrating that nrps in males produces the pheromone BATT, which is essential for male-stimulated female development. The authors report a neural cell type (“delta-3”) that expresses nrps, and the S. japonicum homologue of lsamp, which expresses a neural surface protein in other species. The authors functionally characterise lsamp, confirming its expression in delta-3 cells and showing that it is required for BATT production but affects neither delta-3 development nor nrps expression. They further show that lsamp expression is closely colocated with the axons and synaptic junctions of a ventral neural network. Knocking down lsamp disrupts this network. The authors suggest that lsamp is required for males to produce or maintain neural fibrils.\nOverall, this is exciting work that deepens our understanding of sexual signalling and development in these obligate dioecious parasites. The work is highly significant within the parasitology field and will be of broader interest for developmental biologists.\nIt is a large body of work and, in places, the manuscript is currently hard to navigate. There are some sections that may benefit from simplifying/streamlining, including some analyses that didn’t seem strictly necessary. We also felt that there were aspects of the functional characterisation, where the authors have over-reached in their interpretations.\nWe sincerely appreciate your recognition of the novelty of our work. In line with your constructive suggestions, we have carefully revised the manuscript to enhance readability by simplifying analyses that were not central to the core narrative and by cautiously refining the interpretations of certain functional characterizations. Our detailed responses are provided below.\nSpecific comments - Major points\n1. To validate clusters (S. Fig 3), some of the hybridization patterns are unclear. Presenting evidence for the more diffuse tissues is difficult. Nevertheless, ctsb, tpm2, syn2 do look very similar.\nThank you for your feedback. We fully agree that the extensive distribution of tissues like muscle and tegument makes it challenging to resolve specific cell identities based solely on WISH patterns. To definitively discriminate between such broadly expressed but distinct cell types, we performed double-FISH in male worms, using published marker genes (such as tal for tegument and tpm2 for muscle) and the novel marker genes (slc7a14 for syncytial, slc6a5 for parenchymal, and cpb for muscle) identified in our study.\nAs shown in the revised Supplementary Fig. 4, the results confirmed that the previously reported tegument marker tal and the muscle marker tpm2, though both widely expressed, labeled distinct cell types (Supplementary Fig. 4a). Similarly, the newly identified parenchymal marker slc6a5 (Supplementary Fig. 4b) and syncytial marker slc7a14 (Supplementary Fig. 4c) did not colocalize with tpm2—a pattern particularly evident in the oral sucker of male worms (Supplementary Fig. 4c, lower panels). Furthermore, the novel muscle-associated marker gene cpb showed complete colocalization with tpm2, confirming its expression in muscle cells (Supplementary Fig. 4d). Together, these findings resolve cellular identities that were ambiguous from WISH data alone and further validate the clustering accuracy of our single-cell RNA sequencing data. These results have now been incorporated into the revised manuscript (908–915).\n2. Describing the proportions of cell types is quite prominent in the manuscript (incl. Figs 1 c and d), but how informative is this analysis, beyond an extremely crude look? Sorting, library prep etc. creates bias for certain cell-types (e.g. more round cells like stem cells are favoured by sorting) and can distort recovered cell numbers. In previously published works, actual cell abundance vs cell captured by the scRNAseq process differ wildly, affecting differing cell types unequally. Some of the numbers make sense, like the expansion of vitellaria in females, but others make less sense like the much higher abundance of flame cells than gut cells. There’s a vague/superficial statement, which could be removed to save space, about there being ~20% neural and neoblast cell types when the numbers actually differ quite markedly between timepoints. Abundance across the whole dataset is shown in S Fig 5; yes it does broadly confirm some expectations but it doesn’t really shed any new light. And some of the findings seem counter-intuitive: like the proportion of muscle cells from males.\nThank you for raising this important point. We agree that inferring absolute cell type proportions from scRNA-seq data can be misleading due to technical biases introduced during cell dissociation, sorting, and library preparation. Although all samples were processed under the same protocol, making relative comparisons feasible, we fully acknowledge that these biases may result in the over- or under-representation of certain cell types. In line with your suggestion, we have now swapped the positions of the previous Supplementary Fig. 5 (showing abundance across the whole dataset) and Fig. 1c-d, removed statements regarding approximate proportions of neural and neoblast cells, and revised the manuscript to tone down and better contextualize descriptions of cell type abundances throughout in lines 149-159.\n3. The authors performed subcluster analysis on sexual tissues and neoblast populations, from which they define 17 sub-types. The rationale for their naming scheme is not clear and abbreviations need to be explained. A discontinuous series from the Greek alphabet, plus some numbering (i.e. alpha, beta, gamma-1, gamma-2, delta, epsilon, kappa, nu) is used. Previously, a similar scheme has been used in published descriptions of S. mansoni stem cells (see doi.org/10.1016/j.pt.2021.03.005). It’s not clear if published delta, epsilon, and kappa sub-types are equivalent to those reported here. Distinguishing between 17 populations in the figures (S Fig 6a and Fig 2) is also challenging. Perhaps a numerical key could be added to help readability?\nApologies for the confusion. To improve clarity in the naming scheme of neoblast subpopulations, we have replaced the Greek-letter designations with a numerical system in the revised Fig. 2b, Supplementary Fig. 7 and 9. The updated names are: SSCs-1 (formerly alpha), SSCs-2 (formerly beta), SSCs-5 (formerly gamma-1), SSCs-6 (formerly gamma-2), SSCs-7 (formerly delta), SSCs-3 (formerly epsilon), SSCs-4 (formerly kappa).\nAs you noted, a previous study on stem cells in S. mansoni asexual (sporocyst) and sexual (juvenile) stages designated certain populations as delta, epsilon, and kappa subtypes (PMID: 29988015). To aid comparison, we examined the expression of orthologous stem cell marker genes in our dataset, which provides a cross-reference between the stem cell populations of S. mansoni and S. japonicum. As shown in the revised dot plots (Supplementary Fig. 7b–c), SSCs-1 to SSCs-4 exhibit strong similarity to δ-cells, characterized by high expression of fgfra, fgfrb, hes2, and tp53. SSCs-5 and SSCs-6 likely correspond to ε-cells, marked by elevated expression of astf and bhlh (PMID: 33893056). Furthermore, to enhance figure readability, we have regenerated the UMAP plots with a clear numerical key (Fig. 2b, and Supplementary Fig. 7a).\n4. RNA velocity analysis was performed on the sex related subsets, revealing trajectories (Fig 2b) that made sense like the GSC -> GSC progeny in both sexes and the S1 vitellocytes ->->late vitellocyte. There are also a lot of arrows converging on the centre (SSC beta?) that aren’t mentioned. Further subsetting (Fig 2c) reveals trajectories from GSC to gametes in males but the equivalent analysis is not at all convincing in females (contrary to the main text). In Fig 2c the four panels show a lot of redundant information and could be reduced to two panels.\nThank you for these valuable suggestions. We agree with the reviewer that the RNA velocity analysis shows arrows converging on the center, likely around SSCs-2 (formerly beta). In the original version, we focused on the vitellaria lineage in females, but these central arrows may indeed represent putative progenitor differentiation events that occur in both sexes. We have added a discussion of this point to the Discussion section in lines 484-490: “Although our RNA velocity analysis of sex-related subsets revealed clear trajectories from GSCs to their progeny, it consistently showed convergence toward a central cell state identified as SSCs-2 in both sexes (Fig. 2b). We hypothesize that this convergence point represents a shared, multipotent progenitor pool that functions as a developmental hub, capable of giving rise to multiple downstream cell types in both males and females. This finding offers new insight into the hierarchical organization of gonadal cells in this organism, suggesting the presence of a common—rather than sex-specific—early developmental regulatory mechanism.”\nFor the observed differences in germline cell trajectories between males and females, we carefully re-examined the RNA velocity analyses. Although the female GSC dataset is of reasonably high quality (as shown in the Response Fig. 1a), phase plots of high-confidence genes revealed that female cells lacked sufficient splicing signal (e.g., Sjc-0008397, Sjc-0003894, Sjc-0006007), whereas male cells displayed the expected almond-shaped distributions (as shown in the Response Fig. 1b). This deficiency in splicing information likely accounts for the noisy and inconsistent velocity directions observed in females. Additional factors may include sparse sampling of intermediate states or potential sex-specific splicing dynamics. Given these limitations, we concluded that RNA velocity is less suitable for inferring female GSC trajectories. Accordingly, the RNA velocity plots have been removed from the figures and manuscript text accordingly. Trajectory inference was subsequently performed employing the more robust expression-based method (lines 678-687, Slingshot (PMID: 30089906), as detailed in the revised Fig. 2c-d and lines 195-198.\nIn the revised Fig. 2c-d, we have retained two panels for each to improve clarity for readers: the left panel displays the cellular subclusters, while the right panel shows the pseudotime trajectories.\nResponse Figure1. Evaluation of sequencing data quality and RNA velocity fit on male and female GSCs. (a) Quality control summary of male and female GSCs. Left: UMI counts per cell in male and female datasets, with a nUMI cutoff of 1,500 indicated by the line. Middle: Number of genes detected per cell, with a nGene cutoff of 1,000 indicated by the line. Right: Transcript complexity in male and female datasets, with a log GenesPerUMI cutoff of 0.8 indicated by the line. (b) Phase portraits of the top 5 genes with the highest likelihood in female (top) and male (bottom) GSCs, as calculated by scVelo.\n5. There’s quite a long section on clustering male and female specific genes. The authors reveal a cluster of genes with opposite temporal patterns between the sexes. It’s also interesting to note that for the majority of sex-specific genes, the temporal pattern is different between the sexes, with male genes reducing and female genes increasing. However, some of the methods and figures in this long section are difficult to navigate and it’s disappointing that the section is left quite open-ended with no further analysis of the 47 C1 genes. The overall conclusion is that the process is complex and involves modules of genes. Without the composition of the modules being described, it doesn’t really provide anything new.\nThe main text (L201+) and legend for Fig 7 really need signposts to the relevant section of Methods because it’s not immediately obvious that trajectories are the input for the sex-specific analyses. Even after the gam R description is found, it’s still unclear exactly what the inputs are and how the lists of 693, 750 and 156 relate — are they mutually exclusive? If not, overlaps should be indicated. In Fig 7b, pseudotime=0 is presented in the centre, so that one panel reads right-to-left and the other reads left-to-right; that seems counter-intuitive. How were the datasets clustered? It looks like the cluster number was set to 5 in each case, but it’s not clear why, when several clusters could have been further divided and some possibly merged. What is the dotted line on each cladogram?\nThank you for these detailed and constructive comments. To address your concerns, we have completely re-analyzed this section and revised both the Methods and Figures to improve clarity and informativeness.\nIn the revised manuscript, we applied the tradeSeq R package (PMID: 32139671) (Methods, lines 688–699) to identify temporally dynamic genes along male and female germline pseudotime trajectories separately. Using a stringent p-value cutoff (1e-10), we robustly defined dynamic gene sets for each sex (revised Supplementary Fig. 8, and Data 8-9). This new analysis replaces the previous version that included the 47 C1 genes and the three gene lists (693, 750, and 156), resolving the previous ambiguities and improving interpretability.\nTo enhance biological interpretation, we have supplemented the description of module composition with detailed functional enrichment analyses of the genes in each module (lines 193-222). We firstly performed GO enrichment analysis on the dynamic gene modules to reveal the biological processes associated with each cluster (revised Supplementary Fig. 8, and Data 10-11). A full list of module genes is provided in Supplementary Data 8-9, 12 to ensure transparency. The number of clusters was determined using the getClusters function in ClusterGVis, which applies the Elbow method to evaluate within-cluster sum of squares (Methods, lines 694-699). We further incorporated a focused analysis of known germline regulators (e.g., oc-1, nanos-1, eled, boule, nuob, meiob, and horm-2), along with newly identified candidates—ribc2 and Sjc-0003040—which are highly expressed in the late germline cells of males and females, respectively (revised Fig. 2e-g). Their pseudotime expression profiles are now shown explicitly, providing finer resolution into the regulatory programs underlying schistosome GSC differentiation (revised Fig. 2c-d). Taken together, we have revised this section to be more concise and included clearer conclusions (lines 193-218).\nTo improve clarity, explicit cross-references have been added in the main text (line 239) and in the legend of Supplementary Fig. 8, directing readers to the relevant Methods section describing the trajectory-based identification of dynamic genes (lines 688-699). The revised Supplementary Fig. 8 focuses on pseudotime trajectories and gene dynamics patterns for each sex separately, with reorganized orientation and removal of the dashed lines automatically generated by the ComplexHeatmap R package. In addition, we identified shared dynamic genes between male and female datasets by taking the intersection of their respective dynamic gene sets, and we now provide a comprehensive list of genes dynamically expressed in both sexes across each cluster as updated Supplementary Data 12.\n6. After the module analysis, the narrative returned to trajectory analysis using RNA velocity. The authors were unable to connect stem cell development to vitellaria development but were able to show trajectory from S1 vitellocyte to S1 progeny, with the latter marked by esg-2a expression. The paragraph finishes with a rather bland statement that the work aligned with previous findings and provides new insights. Vf1 expression was known, but the esg-2a expression later in differentiation is new. What else switched on later? Having established a convincing trajectory, was it not possible to mine further information on genes expressed early vs late? Likewise, the analysis of Mehlis’ gland is interesting but only two strong markers are shown in Supplementary Figure 6. Could no other markers be found?\nThank you for these suggestions. We performed sub-clustering and trajectory analysis of vitelline cells to provide a finer resolution of their differentiation from S1 to late vitelline cells (revised Supplementary Fig. 9). Sub-clustering identified eight distinct subtypes, including S1, S1 progeny, early vitellocytes (EVi), and late vitellocytes (LVi) (Supplementary Fig. 9a, and Data 13). Slingshot trajectory inference revealed two potential differentiation branches (Lineage 1 and Lineage 2; Supplementary Fig. 9c–d). Evaluation of differentiation potential using CytoTRACE (PMID: 31974247) indicated that only Lineage 1 exhibits high differentiation potential (Supplementary Fig. 9e). Using the tradeSeq R package (PMID: 32139671) to analyze the gene expression dynamics along these lineages, we found, in addition to the known marker vf1, genes such as p48 and the newly identified esg-2a display a gradual increase in expression along Lineage 1 (Supplementary Fig. 9f). These findings support the designation of Lineage 1 as the primary vitelline developmental pathway. Furthermore, to address your comment, we have provided a set of significantly dynamic genes identified along the vitelline pseudotime (Supplementary Data 14), encompassing genes exhibiting high expression at distinct stages of vitelline cell differentiation, including both early- and late-stage–enriched genes.\nThe functional composition of each gene module involved in Lineage 1 differentiation was summarized by GO term enrichment (Supplementary Data 15), and a complete marker gene list is provided (Supplementary Data 13) to enable readers to explore candidate regulators in detail.\nFor the Mehlis’ gland, we extended the analysis beyond the two strong markers originally shown in Supplementary Fig. 6c (now named Supplementary Fig. 11c). Additional candidate marker genes have now been identified and are included in the revised Supplementary Data 7, providing a more comprehensive molecular characterization of this organ.\nThe above content has been updated in the manuscript lines 223-253.\n7. The authors subset the neuronal data and recluster into a highly connected UMAP showing at least 5 lineages. This is an important part of the reference resource created by the project. However, the authors should maybe avoid using nomenclature based on Greek letters, given that it is already used for stem cells, including within the same paper. How about just numbering them? Over time, more functional names will be applied but, for now, it might suffice to name them as neural lineage 1, etc. Or maybe use N1, etc? The more granular names like delta-3 would then become N4.3.\nThank you for the suggestion. We have simplified the nomenclature of neuronal populations by adopting a consistent numerical system (e.g., N1–N5). Specifically, Pre. alpha has been renamed N1. p, with its derivatives designated as N1.1 Neuron (alpha-1 Neuron) and N1.2 Neuron (alpha-2 Neuron). Similarly, Pre. beta is now N2. p, giving rise to N2.1 Neuron (beta-1 Neuron) and N2.2 Neuron (beta-2 Neuron). Pre. gamma has been renamed N3. p, with its sub-lineages renamed N3.1 Neuron (gamma-1 Neuron) and N3.2 Neuron (gamma-2 Neuron). For the delta lineage, Pre. delta and related progenitors are now designated as N4. p1–N4. p3, while the corresponding neuronal subtypes are renamed N4.1–N4.4 Neurons (delta-1 to delta-4 Neurons). In addition, epsilon Neuron has been renamed N5 Neuron. All figures (Fig. 3, and Supplementary Fig. 12-19) have been updated to reflect this revised nomenclature.\n8. The authors use Slingshot to define lineages and derive pseudotime from their neural UMAP. How well does the pseudotime derived from these Slingshot trajectories agree with the sample-time metadata? (does a UMAP coloured by sample show differential enrichment for earlier or later parts of the trajectories?). The subsequent RNA velocity analysis agrees well with their pseudotime predictions, and Cell Rank is used to define terminal states (fully differentiated cells?). Similar to the earlier section on gamete development, the signposting to the methods is not very clear in the text - it’s not immediately clear the purpose of each analysis until reading further into the paper. Cell Rank is not mentioned in the legend to S Fig 10 and “dpt pseudotime” is not explained. The PAGA analysis, resulting in a network of pie charts (Fig. S10e and 11e), isn’t clear. What exactly is this showing? Is it really needed?\nThank you for your insightful comments on the neural lineage analysis. To examine how pseudotime inferred from Slingshot trajectories corresponds to actual sample collection times, we visualized the pseudotime distributions across samples using density plots, which were performed separately for male and female samples. The results show a clear alignment, with earlier samples enriched at the beginning of the trajectories and later samples at the terminal ends. The new plots illustrating this concordance have been added (Supplementary Fig. 12-13c).\nAs noted by the reviewer, the terminal states defined by CellRank can indeed be regarded as fully differentiated cells. This interpretation is consistent with the original CellRank publications (PMID: 35027767), where terminal states were shown to correspond to mature cell types in systems such as pancreatic lineage development and human hematopoiesis.\nTo improve clarity for readers, we added further details in the Methods section (lines 709–713) and explicitly mentioned the use of CellRank in Supplementary Fig. 12-13. We also clarified that “dpt pseudotime” refers to Diffusion Pseudotime (DPT) as calculated using scVelo (PMID: 32747759, lines 266-268 and 706-708).\nFinally, to simplify the presentation and avoid potential confusion, we have removed the PAGA analysis and associated pie-chart network figures, as suggested.\n9. Regarding the section on male-specific neurons: in schistosomes, a large number of genes are present in sex-specific regions (~900 in S. mansoni) and, given their incomplete dosage compensation between sexes, those sex-specific genes can skew the relative abundance of all other transcripts. Did the authors check whether the clustering was influenced by genes in those regions? If it is, there could be distinct clusters that are expressing the same genes but just differ in the expression levels across their sex-specific genes.\nThis is an excellent point. To address it, we examined the genomic locations of marker genes for the N4.3 Neuron cluster (formerly delta-3) and found no significant enrichment of Z-specific genes (Response Fig. 2a, Chi-squared test, p = 0.23). We also assessed the expression of Z-specific genes across cell populations using Seurat’s “AddModuleScore” function, which revealed no significant differences in module scores between sex-specific clusters and other populations (Response Fig. 2b). These findings suggest that the male-specific clusters are unlikely to result from dosage compensation between sexes.\nResponse Figure 2. Examination of the enrichment of Z-specific genes in N4.3 neuron populations. (a) Chi-square test results comparing Z-specific vs. non-specific genes among N4.3 neuron marker genes and non-marker genes. (b) Module scores of Z-linked genes across different neuron subtype populations.\n10. Beta-2, gamma-2 and delta-3 neurons are described as ‘male-specific’ (Line 291), but the markers for two of these populations (beta-2 and gamma-2) were then detected in the heads of females (Line 304). Was this a sensitivity issue of sc-RNASeq vs. WISH? Is it really appropriate to describe beta-2 and gamma-2 neurons as ‘male-specific’? They have also included beta-2 and gamma-2 neurons on the ‘Female anterior’ schematic in Fig. S16, suggesting they are not male-specific in the authors’ model. This needs some clarification in the text.\nThank you for raising this important point regarding the apparent discrepancy in the description of beta-2 (now named N2.2) and gamma-2 (now named N3.2) neurons. We have added the clarification in the lines 297-300: “Among these, the N4.3 population was exclusively distributed throughout the gynecophoric canal, whereas N2.2 and N3.2 populations were additionally detected in the head of both male and female parasites (Fig. 3b, and Supplementary Fig. 16).”\nBased on WISH and FISH results—which offer higher spatial resolution and sensitivity for low-expression markers and thus better reflect the in vivo situation—these neuronal populations are predominantly present in males, while the N2.2 and N3.2 neurons also distribute in female anterior. Therefore, we prioritize WISH/FISH localization data for accurate cell type characterization.\nTo accurately reflect these findings and avoid overinterpretation, we have revised our terminology throughout the manuscript. N2.2 and N3.2 neurons are now referred to as \"male-enriched\" rather than \"male-specific\".\n11. Navigating through some of the results on markers for the male specific clusters is confusing due to the different types of terminology. The markers are initially referred to by their cell type and systematic gene id (in the text and Supplementary Figure 13) but in the WISH validation gene names like nrps and ddc appear but without their systematic names. I think it would be simpler to read, if the more human-readable gene names were used in the main text where possible (with systematic IDs in parentheses). In the legend for Supplementary Figure 13, I had to work out some names, like “ddc” = Sjc-0006754 by subtraction. Similarly, I needed to cross-reference Supplementary Fig 14 to learn that Sjc-0000814 = “nrps”.\nApologies for the confusion. We have provided gene names alongside gene IDs throughout the main text to improve readability and clarity across the manuscript.\n12. It’s interesting that lsamp was not identified as a marker for delta-3 cells. It had the third-highest correlation across the delta-3 trajectory, but what was its expression across the data set. Was its expression just too low to be picked as a marker by Seurat?\nThank you. Marker genes in our analysis were defined based on a ranked list using adjusted p-values. When applying the default filtering criteria in Seurat, including the default pct.2 threshold (the percentage of cells expressing the gene in non–delta-3 (N4.3) populations), lsamp did not appear among the top-ranked markers and was instead ranked further down the list. However, when we applied a more stringent filter on pct.2, lsamp emerged as one of the top-ranked genes, reflecting its strong specificity for delta-3 (N4.3) cells despite relatively low expression across the dataset (Response Fig. 3 and Supplementary Data 16). This explains why lsamp was not initially highlighted as a canonical marker in our results.\nResponse Figure 3. Expression levels of lsamp (left) and nrps (right) across the integrated neural cell datasets.\n13. The GO analysis of genes affected by lsamp knockdown isn’t particularly informative. The KEGG analysis (Fig 5) and functional information gleaned from other systems enable the authors to hypothesise about the role of lsamp. The authors found fewer nerve fibres in RNAi knockdown but suggested that LSAMP is therefore “essential for nerve fiber growth”. Given that the experiment was conducted on adult worms, would the nerves be fully grown? In which case, the knockdown presumably destabilised the nerve network rather than stopped growth?\nThank you for this critical insight. The limited informativeness of the GO analysis for genes affected by lsamp knockdown may be attributed to the current quality and coverage of GO term annotations in this system. Alternatively, it is also possible that lsamp knockdown only induces relatively slightly changes in gene expression, resulting in the downregulation of merely 50 genes (Supplementary Data 21), which could limit the statistical power of the enrichment analysis.\nWe also agree that the adult schistosomes used in our RNAi experiments already possess a mature nervous system in which nerve fibers are formed. Thus, the observed reduction in nerve fibers following lsamp knockdown likely reflects a role in maintaining structural stability rather than in developmental growth. Accordingly, we have revised our conclusion to state that: “These data suggest that LSAMP is essential for nerve fiber maintenance, likely by mediating appropriate axon guidance and target recognition. The observed alterations in microtubules may represent a secondary morphological consequence of impaired axon maintenance.” in lines 407-410. And we have clarified this point in the Discussion in lines 543-547: “Based on these observations and homology inference, we propose—though functional validation remains to be carried out—that LSAMP may serve as a selective adhesion molecule on the surface of N4.3 neurons, potentially maintaining axonal integrity and facilitating neuronal connectivity (Supplementary Fig. 27). Further experimental studies are needed to confirm its precise membrane localization and mechanistic role in schistosomes.”\n14. The colocalisation experiments seem to show lsamp having a similar relationship as nrps and ddc to acetylated tubulin and synapsin but it’s difficult to interpret when such a small part of the worm is shown. Can lower magnification images be included for orientation? The colocalisation isn’t perfect, so it’s difficult to see how this suggests “active vesicular transport between these axons and delta-3 cells”. Based on expression timing, the authors speculate that “lsamp facilitates axon growth and attachment to delta-3” but evidence isn’t very strong for this. They also say that “lsamp is essential for guiding the formation of [nerve fibers]” but this is also not supported by the evidence. We suggest that the speculation be moved to the discussion, to make it much clearer where the evidence gaps remain. In the discussion L530–531, the authors state that “LSAMP functions as a membrane protein”. Although quite likely, it does remain to be shown, so appropriate caveats need to be inserted.\nThank you for these insightful suggestions. We have added the lower-magnification images in the revised manuscript to better illustrate the overall distribution and relationship of N4.3 markers with neural axons (revised Supplementary Fig. 24, lines 414-421).\nWe acknowledge that our interpretation regarding “active vesicular transport” may have been overstated. We have therefore revised the corresponding text to focus on the observed spatial proximity in lines 417-419: “Synaptic proteins along neural axons also exhibited substantial colocalization with N4.3 neurons, consistent with possible vesicular trafficking between these structures (Fig. 5d, and Supplementary Fig. 24c).”\nFurthermore, we agree that the precise mechanistic role of lsamp in axon guidance and growth requires further evidence, and we have clarified and expanded the statement that “…. confirming its direct function as a membrane protein in axon guidance and adhesion requires further investigation” in lines 564–567.\nMinor points：\n15. L147–148: The statement about differentiated neuron populations expanding after D14 only seems to apply to males. In females it goes up then down. But is this a real effect or is this due to the differing ability to capture certain cells as the worms in each sample are morphologically different?\nWe thank the reviewer for this observation. We propose the observed dynamics in female neuronal populations—initial increase followed by a decrease after D14—as a likely biological effect, rather than a technical artifact. This interpretation is supported by our previous bulk RNA-seq data, which showed significant downregulation of neuron-related genes in mature female worms (PMID: 28287085). Nevertheless, we acknowledge that additional biological replicates would strengthen the conclusion and further validate this effect. We have accordingly shortened the corresponding statement in lines 149–153 of the manuscript to avoid potential overinterpretation: “For somatic tissues, the proportions of different cell populations also showed marked dynamics. The proportion of neural progenitor cells and precursor cells decreased after 14 dpi in both sexes (Supplementary Fig. 6a, and Data 5). In male worms, the populations of differentiated neural cells (e.g., nonciliated neurons and ciliated neurons) continuously expanded (Supplementary Fig. 6a, and Data 5).”\n16. L149 (word choice): The neural progenitor cells presumably don’t “degenerate”, but differentiate into neurons.\nThanks, corrected.\n17. L189 esg-2a is identified as a new marker. This gene encodes eggshell protein (PMID:2172818).\nThank you. In the revised manuscript, we have included the full name of esg-2a (eggshell protein gene 2a) and cited the relevant reference (PMID: 2172818) in line 176.\n18. L210–218 How were the example genes in Fig 7 b-d picked? It looks arbitrary.\nThank you for raising this point. We originally selected example genes in the previous Supplementary Fig. 7b–d based on the range of predicted values from the GAM model (range_pr, defined as the difference between maximum and minimum fitted values). This criterion alone may not highlight functionally representative genes; thus, we have presented well-established germline regulators (e.g., oc-1, nanos-1, eled, boule, nuob, meiob, and horm-2) as illustrative examples in the revised Fig. 2e-g. We also have provided the full gene lists in the revised Supplementary Data 8-9, and 12. (lines 199-218).\n19. L300. The authors state that all three neurons are primarily located in the gynecophoric canal. According to Fig 3b, this statement is only actually true for delta-3. For the other two neuron types, they are localised in the canal but are expressed more in other areas.\nThank you for this helpful observation. We have revised the annotations regarding the neuronal subtypes in the head region and adjust our wording to more accurately reflect their spatial distribution in lines 297-300: “Among these, the N4.3 population was exclusively distributed throughout the gynecophoric canal, whereas N2.2 and N3.2 populations were additionally detected in the head of both male and female parasites (Fig. 3b, and Supplementary Fig. 16)”.\n20. L310–315 After the male specific neurons part, the narrative shifts to validating the other neural types (Supplementary Figure 15 and 16), before returning to male specific neurons. This seems to interrupt the flow. Maybe this could be moved to earlier in the text, before the male specific part.\nThank you for the suggestion. We have moved this part to an earlier section in the manuscript to ensure a more logical and coherent structure (lines 289-302).\n21. L318–320 “These enriched genes indicated that beta-2 and gamma-2 populations exhibited typical…In contrast, gamma-2 showed…” This isn’t a contrast, you could say “Moreover, gamma-2 showed…” instead.\nThank you. We have revised this as you suggested.\n22. L347: This sentence looks like it needs re-phrasing. Do the authors mean “...worms treated with lsamp dsRNA but not ddc dsRNA lost their ability…”\nRight, we have revised the sentence accordingly. (lines 341-348)\n23. L372: “ …suggest that lsamp impairs the ability of male delta-3 neurons to synthesise BATT”. We think the sentence, as written, means the opposite of what the authors intended. It should say “lsamp RNAi”.\nThanks, corrected.\n24. L409: “nerve fibers were markedly eliminated without affecting nrps expression (Fig 5b)”. The figure doesn’t show nrps expression.\nThank you for pointing this out. We have now included the figures demonstrating unchanged nrps expression after lsamp knockdown in male worms (Fig. 4f, and Supplementary Fig. 22) in the revised manuscript (line 360).\n25. L410: “The data suggest LSAMP is essential for nerve fiber growth by stabilizing microtubules”. Yes fiber growth is inhibited but you can really only speculate on the mechanism at this point.\nThank you. To be more precise, we have revised the sentence as follows: \" These data suggest that LSAMP is essential for nerve fiber maintenance, likely by mediating appropriate axon guidance and target recognition. The observed alterations in microtubules may represent a secondary morphological consequence of impaired axon maintenance.\" (lines 407-410).\n26. L444: The nerve network functions as an interstate. Maybe use the word “conduit” instead of “interstate”.\nChanged.\n27. L530–531: “As a glycosylated cell adhesion molecule of the IgLON family, LSAMP functions as a membrane protein”. This is quite an assertive statement about the function of the LSAMP protein. However, I don’t see evidence that the authors have tested these functional predictions in the text, so it would be good to make that clear in this statement.\nThank you for raising this important point. The reviewer is correct that we did not directly test the membrane localization of LSAMP in this study. Our statement is supported by well-established characterizations of LSAMP and other members of the IgLON protein family in various organisms (PMID: 20537990, PMID: 30555000, PMID: 29605488, PMID: 31043568, and PMID: 28743494). According to these studies, LSAMP is a glycosylphosphatidylinositol (GPI)-anchored protein, representing a classic class of membrane-associated proteins. This GPI anchoring, together with its three immunoglobulin-like domains, is a defining feature of the IgLON family and is consistently annotated across multiple protein databases (e.g., UniProt ID Q13449).\nTo avoid overstating our findings, we have revised the manuscript to clarify that this is a previously established property and have included the appropriate citation. The revised text now reads: \"As a member of the IgLON family, LSAMP is a glycosylated cell adhesion molecule that functions as a GPI-anchored membrane protein.\" (lines 537-539).\n28. Fig. 2b The annotation of “late Vi” on UMAP seems to occupy the position annotated as “Mehlis’ gland” on the same UMAP when used in Supplementary Figure 6. Which one is it? The trajectory looks like it runs S1 vitellocyte -> S1 progeny -> Early Vi -> Late Vi -> Mehlis\nWe apologize for the mislabeling. In the original UMAP, the annotation “late Vi” was mistakenly placed in the region that should have been labeled “Mehlis’ gland.” This has been corrected in the revised Fig. 2b.\nRegarding the vitelline lineage, the trajectory from S1→ S1 progeny → Early Vi → Late Vi is well supported by our data. However, the connection between Late Vi and the Mehlis’ gland is less clear based on the trajectory arrows (Fig. 2b). Instead, the UMAP suggests that the Mehlis’ gland may originate from the light purple population (SSCs-beta, now renamed SSCs-2).\n29. Fig 3a and Supp Fig 12a. These figures are very similar and could be combined e.g keep the colours but highlight the three male-specific populations with boxes.\nThanks for your suggestion. We have revised the Fig. 3a and removed the original Supplementary Fig. 12a. The UMAP plot in Fig. 3a is now colored, and the three male-enriched populations are highlighted with circles.\n30. Fig 3b and d use consistent gene naming.\nThanks, changed.\n31. Fig 3d is very small given its importance to the lsamp story\nThank you for your suggestion. We have now enlarged this panel in the revised Fig. 3d.\n32. Fig 5 (legend) “tight adherence” between lsamp, delta-3 markers, synapsin and acetylated tubulin in axons is incorrect. It shows colocalisation.\nThank you for pointing this out. We have revised the description accordingly.\n33. In Supplementary Figure 12b, the third panel (% composition of clusters by time) shows the three male specific neurons nicely. The other two panels are less useful and panel 2 is mislabelled (should be M not F).\nThank you. We have removed the first and second panels in the revised Supplementary Fig. 14.\n34. Fig. S12C – on enriched neural cell-types graph (on right) the red arrow for ‘Enriched in female’ would be more intuitive if facing the other way in the legend\nThanks, it has been modified.\n35. Supp Fig. 22. Is it really needed — the information has already been presented elsewhere in the manuscript.\nThank you. Deleted.\n36. Supp Fig. 24b – what do red and blue signify?\nWe apologize for the confusion. We performed the sequence alignment using the Multiple Sequence Alignment Viewer from NCBI, and the colors displayed in the results adhere to the conventional ClustalX color scheme, which groups amino acids based on their physicochemical properties to visually represent conservation. Hydrophobic residues (A, V, L, I, M, F, W, P) are colored blue; positively charged residues (K, R) are red; and negatively charged residues (D, E) are magenta. Polar uncharged residues (S, T, N, Q, G, C) are green, while aromatic and special residues (H, Y) are highlighted in pink or orange. The color intensity at a given position indicates the degree of sequence conservation, with darker shades representing higher conservation. Explanations for the two colors have now been added to the legend of the revised Supplementary Fig. 27b (Lines 1084-1089).\nMethods\n37. For the scRNAseq, QC details are quite sparse. The >30% mitochondrial reads filter seems high. Were there a lot of cells with high %mito?\nThank you for raising this concern. We adopted the 30% cutoff to remain consistent with the criterion used in the published single-cell RNA atlas of S. mansoni (PMID: 32973030). To further address your point, we analyzed the data and found that the vast majority of cells had low mitochondrial content: approximately 97% of cells had <10% mitochondrial reads (as shown in the Table 1 below). Therefore, only a small proportion of cells exhibit a high mitochondrial percentage.\nTable 1. Statistics of mitochondrial percentage across 7 samples.\n38. No methods included for the sex and neural subsets of 17 and 19 cell-type populations\nWe have added the related methods in lines 665-676.\n39. L611–612 ‘All samples were combined using the functions “FingIntegrationAnchor” and “IntegrateData”’. Which features were the datasets integrated on? Timepoint? Sex? Both?\nThank you. This approach is designed to mitigate technical batch effects across datasets (in our case, across samples), while not explicitly aligning data on metadata attributes such as timepoint or sex. To make it clearer, we have added more descriptions in this section in lines 633-637: “Integration anchors across samples were identified using the FindIntegrationAnchors function, and the datasets were subsequently integrated using IntegrateData from the seurat R package. This sample-level integration approach minimizes batch effects while preserving biological variation across samples.”\nReviewer #2:\n1. I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely appreciate your thorough review.\nReviewer #3:\nIn this manuscript, You et al. identified a key gene required for male schistosome-triggered female reproductive maturation in Schistosoma japonicum. Specifically, they applied single-cell transcriptomics, RNA velocity and pseudotime analysis, and RNA in situ hybridization to reveal the cellular and molecular dynamics underlying this process. By focusing on four key time points, they uncovered putative developmental trajectories of male and female cells, focusing on germ cells and various neuronal cell types. Notably, the authors found that one of the three male-specific neuron subtypes, delta-3 neurons, relies on lsamp for normal axonal growth and ventral nerve network formation in the gynecophoric canal of male parasites. Through RNA interference experiments, they demonstrated that defects in this male nerve network due to lsamp knockdown reduce pheromone BATT production and subsequently impair female sexual development. Overall, this study provided a comprehensive single-cell transcriptomic dataset of both male and female S. japonicum during the important male-female pairing process and discovered a male nerve network critical for female reproductive development. The paper is well-written with high-quality figures that present solid experimental evidence. Beyond the important discovery of the role of lsamp in BATT signaling and its contribution to understanding this new signaling molecule, this work provides a useful resource for understanding schistosome development and reproductive biology. Thus, it is well suited for publication here. However, there are a few places where results are over-interpreted or alternative interpretations are not considered. Some suggestions for issues the authors should consider addressing in a revised manuscript are listed below.\nWe sincerely thank the reviewer for the thorough review and for recognizing the value of our single-cell transcriptomic dataset and the biological insights into male–female reproductive interactions in Schistosoma japonicum. In the revised manuscript, we have carefully addressed your concerns. Our detailed responses are provided below.\n1. Some of the language used to describe cell states identified via single-cell sequencing overinterprets what has actually been learned. For example, male gametes are sperm and female gametes are oocytes. Since transcripts defining the “gamete” clusters are not presented in high-magnification FISH images, it is not clear whether these clusters reflect mature gametes or earlier stages along the differentiation pathway. Similarly, the cluster referred to as neural stem cells has not been validated in a way that allows them to be called stem cells. No evidence is presented that they self-renew, one of the defining features of stem cells. This could be a large pool of neural progenitors. In the absence of evidence that these cells are actual gametes or neural stem cells, the authors should soften their language.\nThank you. We agree that the terminology should be more precise. Accordingly, we have softened the language throughout the manuscript by replacing “male/female gametes” with more accurate descriptors such as “late male/female germ cells,” consistent with the previously published single-cell atlas of S. mansoni (PMID: 32973030). In addition, we have relabeled “neural stem cells” as “neural progenitor cells”.\n2. The authors conclude that “LSAMP is essential for nerve fiber growth by stabilizing microtubules” (line 410). Since lsamp does not appear to be expressed in the affected neurons, it is unclear how it would act to stabilize their microtubules. Based on what is known about LSAMP function, it seems much more likely that the affected neurons are not able to extend their axons to the appropriate targets (or maintain them there), so any impact on microtubules would be secondary (i.e., there are no axonal microtubules because there may be no axons, not because the microtubules are less stable).\nWe sincerely thank the reviewer for this insightful and critical assessment of our proposed mechanism. We fully agree that LSAMP's established function as a cell-surface adhesion molecule makes an extracellular mechanism—namely, axon guidance and target recognition—far more plausible than a direct, intracellular role in microtubule stabilization. Accordingly, we now interpret the observed reduction in axonal microtubules as a secondary morphological consequence reflecting the failure of nerve fiber maintenance and retraction, not a primary molecular impairment. To address this, we have removed the original statement and revised the text in lines 407–410: “These data suggest that LSAMP is essential for nerve fiber maintenance, likely by mediating appropriate axon guidance and target recognition. The observed alterations in microtubules may represent a secondary morphological consequence of impaired axon maintenance.”\n3. Understanding the phenotype of lsamp RNAi parasites requires the authors to provide much more detail about the RNAi protocol and the subsequent phenotypic analysis. As written, it is unclear when the male parasites treated with dsRNA were initially harvested from the host. If they were harvested at 4 weeks or later, then lsamp is required to maintain the structure of the neural network (because it would have developed already). If they were harvested earlier, it is important to know when, so the reader can place the phenotype in the context of the developmental time course shown in Fig. 5e. The authors should also clarify the time points used for the subsequent analysis of the lsamp RNAi phenotypes. Were all assays conducted at 21 days post-pairing, or were some phenotypes examined at other time points? The graph in Supp. Fig. 22 overinterprets what is known. It is not clear that LSAMP guides the ventral neural network. Their RNAi experiment shows a single undefined endpoint from an undefined starting point.\nWe sincerely thank the reviewer for raising these important points regarding the experimental timeline and interpretation. We have added the detailed information about the RNAi protocol in the Methods (lines 721-730).\nIndeed, the adult male parasites (26 dpi) used in the RNAi assay have already established their ventral neural network. Therefore, our results more likely support the conclusion that LSAMP is required for maintaining the neural network rather than for its initial development, as you suggested. To avoid potential misunderstanding, we have removed the previous Supplementary Fig. 23 and revised the related sections of the manuscript accordingly.\n4. The conclusion that the affected axons transport BATT precursors is not supported by the data; other interpretations are possible. The authors begin their argument with the claim that “synaptic vesicles also transport other small molecules, such as the precursors for BATT (β-alanine and tryptamine)38, 39.” As far as I could tell, neither of these references mentions β-alanine or tryptamine. It remains an open question whether these compounds are transported via synaptic vesicles, released via solute carrier proteins, or if some other mechanism is at work. Furthermore, although the authors note that ddc (aadc) is required for ovary formation in planarians, they do not mention that it is also required for BATT production and is co-expressed with nrps in planarians, where it is likely producing tryptamine from tryptophan in the BATT-producing cells (PMID: 38889152). The lack of a phenotype in the ddc RNAi parasites could easily be due to insufficient knockdown (Supp. Fig. 18). Because ddc (aadc) is co-expressed with nrps in schistosomes, it seems most likely that the tryptamine is produced locally, rather than transported. If not, why is ddc (aadc) expressed in these cells?\nWe thank the reviewer for their insightful comments, which have significantly strengthened our manuscript. We have removed unsupported statements about β-alanine/tryptamine transport via synaptic vesicles. β-alanine is a known substrate for vesicular uptake systems, specifically via the vesicular GABA transporter (VGAT), for subsequent exocytotic release， and relevant citations (PMIDs: 7746489, 8897468, 38103629) have been added in lines 398-402.\nWe agree with your interpretation regarding local tryptamine synthesis, given the conserved co-expression of ddc and nrps in schistosomes, as previously reported in planarians (PMID: 38889152). Accordingly, the weak ddc RNAi phenotype observed in our study is likely attributable to insufficient knockdown, since DDC provides the essential substrate for BATT synthesis (updated Supplementary Fig. 21). We have therefore cited this study (PMID: 38889152) and incorporated a discussion of this point into the revised manuscript in lines 516-521: “The lack of an observable phenotype, despite significant knockdown, could be due to residual enzymatic activity being sufficient to maintain function, or possibly compensation through alternative pathways. We also acknowledge that the incomplete knockdown may have contributed to the absence of a pronounced phenotype. Future studies employing more potent interference strategies or conditional knockout approaches will be necessary to conclusively determine ddc's functional role”.\nConsequently, we have reframed our model: While axonal transport is essential, tryptamine is more likely synthesized locally within nrps⁺ cells, whereas β-alanine—the other precursor of BATT—may be delivered via axonal cargo. Corresponding figures (Fig. 5f, and Supplementary Fig. 28) have been updated.\nMinor comments\n5. To help readability and the interpretation of the figures and text, wherever possible, informative gene names should be used in place of the Sjc numbers. To provide just one example, the in situ in Fig. 3b used Sjc numbers that are not shown elsewhere in the figure and are marked in colors that do not correspond to any others in the figure, so it was unnecessarily difficult to keep track of which neuronal subclass they were marking.\nFollowing your suggestion, we have provided gene names alongside gene IDs throughout the figures and text to enhance readability.\n6. The authors may wish to change the nomenclature of their neuronal subtypes to avoid confusion. Stem cell subtypes in schistosomes have already been defined as delta, using the Greek letter (δ) to define them (PMID: 29988015), so this may add unnecessary confusion to the literature. Furthermore, because delta genes are conserved across animals and play critical roles in the development of many organ systems, defining subclasses as delta-# (i.e., without the Greek lettering) may lead readers to expect that they express a specific delta gene.\nWe thank the reviewer for this helpful suggestion regarding neuronal lineage nomenclature. To avoid confusion with the Greek-letter system already used for stem cell subtypes, we have revised the naming scheme for neuronal populations. In the revised version, neuronal lineages are now consistently numbered (e.g., N1, N2, N3, N4, and N5) with corresponding progenitor and sub-lineage designations. Specifically, Pre. alpha has been renamed N1. p, with its derivatives designated as N1.1 Neuron (alpha-1 Neuron) and N1.2 Neuron (alpha-2 Neuron). Similarly, Pre. beta is now named N2. p, giving rise to N2.1 Neuron (beta-1 Neuron) and N2.2 Neuron (beta-2 Neuron). Pre. gamma has been renamed N3. p, with its sub-lineages renamed N3.1 Neuron (gamma-1 Neuron) and N3.2 Neuron (gamma-2 Neuron). For the delta lineage, Pre. delta and related progenitors are now designated as N4. p1–N4. p3, while the corresponding neuronal subtypes are renamed N4.1–N4.4 Neurons (delta-1 to delta-4 Neurons). In addition, epsilon Neuron has been renamed N5 Neuron. All figures have been updated accordingly (updated Fig. 3, 5f, and Supplementary Fig. 12-19, 28).\nAs for the stem cell subtypes, we have also replaced the Greek-letter designations with a numerical system in the revised Fig. 2b, Supplementary Fig. 7a, and 9. The updated names are: SSCs-1 (formerly alpha), SSCs-2 (formerly beta), SSCs-5 (formerly gamma-1), SSCs-6 (formerly gamma-2), SSCs-7 (formerly delta), SSCs-3 (formerly epsilon), and SSCs-4 (formerly kappa). To aid comparison with the previous study (PMID: 29988015), we examined the expression of orthologous stem cell marker genes in our dataset, which provides a cross-reference between the stem cell populations of S. mansoni and S. japonicum. As shown in the revised dot plots (updated Supplementary Fig. 7b–c), SSCs-1 to SSCs-4 exhibit strong similarity to δ-cells, characterized by high expression of fgfra, fgfrb, hes2, and tp53. SSCs-5 and SSCs-6 likely correspond to ε-cells, which are marked by elevated expression of astf and bhlh (PMID: 33893056). However, due to certain discrepancies between our data and the S. mansoni data in terms of the developmental stages, species, and clustering resolution, a direct one-to-one nomenclature of the stem cell subclusters based on the published classifications is not feasible.\n7. Line 148-9: “This finding indicates that neural progenitor cells are active during the juvenile stage but gradually degenerate in adulthood.” There is no evidence that these cells degenerate. It seems more likely that they simply generate their progeny.\nThank you. We agree with the reviewer that \"degenerate\" is an overinterpretation without strong evidence. We have deleted this sentence.\n8. Since most readers of this journal will be unfamiliar with details of schistosome biology, instead of noting that they find female germ cell markers in males with no additional context (Supp. Fig. 6d), the authors should point out that hermaphroditism is occasionally observed in schistosomes and provide appropriate citations. It is not clear from the images shown in Supp. Fig. 6d that female germ cells are found in the testes because the testes are not visible at the magnification shown. The authors would need to show higher magnification views that show these cells in the testes to support this claim.\nThank you for this suggestion. We have added a discussion of the occasional hermaphroditism observed in schistosomes, along with supporting citations (PMID: 20267997; PMID:16883991; https://www.journals.uchicago.edu/doi/abs/10.1086/417848), in the relevant section of the manuscript in lines 187-192: “Interestingly, although schistosomes are the only dioecious species within the class Trematoda, various forms of hermaphroditism have been documented, including males developing ovarian-like structures and females exhibiting testicular characteristics. Our identification of cells expressing female germ cell markers in phenotypically normal male worms provides additional molecular insight into the evolution of sexual differentiation in schistosomes.”\nWe acknowledge that the magnification used in the previous Supplementary Fig. 6d was too low to unambiguously identify testicular tissue. Accordingly, we have included higher-magnification images obtained by double-FISH in the revised manuscript to substantiate the presence of clec⁺ germ cells—typically found in the female ovary—within the testes (lines 183-187, updated Supplementary Fig. 7e).\n9. Figure 1C, S5, and Data S5 showed changes in cell type compositions across different time points. These measurements are relative, meaning a significant expansion of one cell type may cause a noticeable reduction of other cell types in proportion. Please provide a plot or table with absolute cell counts.\nFollowing your suggestion, we have added a table (updated Supplementary Data 6) reporting the absolute cell counts for each cell type across all samples corresponding to the updated Fig. 1c, Supplementary Fig. 6a, and Data 5.\n10. Line 170: eled is referred to as a germline cell marker, but it is also expressed in vitellaria and neoblasts in the posterior growth zone in juveniles. Can the authors comment on the eled+ cells in males that appear to be in vitellocytes at 22dpi (Fig. 1d)?\nWe appreciate the reviewer’s observation regarding eled expression. The eled+ cells in males that appear to be associated with vitellaria cells at 22 dpi (revised Supplementary Fig. 6b) may reflect either mis-assignment of germline cells by the UMAP algorithm (PMID: 35253834, the Fig1 in this reference illustrates this possibility) or the presence of a previously uncharacterized vitellaria-like population in males. We favor the latter explanation, as we observed some S1 Vitelline-like cells in males at 22 dpi (updated Supplementary Fig. 7c), suggesting that these cells may represent a male-specific vitellaria-like population rather than misclassified germline cells. This interpretation aligns with the work of George Wendt et al., who identified p48- and ataxin2-labeled vitellaria cells in males (PMID: 32973030), thereby confirming the long-standing observation that male parasites harbor low frequencies of vitellaria-like cells (PMID: 7200108).\n11. In Figure 3E, lsamp and ddc expression seems to be partial in the delta-3 neuron cluster. Is this a gene coverage artifact from the scRNAseq or a reflection of the heterogeneity of delta-3 neurons? Quantifying their double FISH studies should be able to clarify this question.\nThank you for this insightful suggestion. Our scRNA-seq data indicate that nrps is expressed in approximately 89.5% of N4.3 neurons, while ddc and lsamp are detected in about 65.5% and 54.2% of these cells, respectively, suggesting potential heterogeneity (updated Supplementary Data 16). However, careful quantification with double-FISH reveals that both lsamp and ddc signals are present in nearly 100% of nrps+ neurons (Supplementary Fig. 20). We therefore propose that the seemingly heterogeneous pattern in the UMAP projection is more likely attributable to the limited sequencing depth of scRNA-seq, which can fail to detect low-abundance transcripts, rather than to true biological heterogeneity within the N4.3 neuronal population.\n12. Unpaired two-tailed parametric t-tests were widely used for the statistics. Have all the data passed a normality test?\nWe sincerely thank the reviewer for raising this critical point concerning data normality and the use of parametric testing. We have reassessed the normality of all relevant datasets using the D'Agostino & Pearson test (Supplementary Data 27). The majority of the data, including those in Fig. 4d–h and the updated Supplementary Fig. 22, were found to be normally distributed. However, the D8 and D20 datasets in the updated Supplementary Fig. 22a and b, as well as the worm BATT data in Fig. 4i, deviated from normality. For these non-normally distributed datasets, we reapplied the non-parametric Mann-Whitney U test (two-tailed) for group comparisons (Supplementary Data 28). The results are fully consistent with the original conclusions. And we have updated the Methods section to include a detailed description of this statistical approach (lines 794-802).\n13. Besides the top three genes, nrps, lsamp, and ddc, what other genes and associated pathways are specific to delta-3 neurons? Can the authors elaborate on this further?\nWe thank the reviewer for this helpful suggestion. In addition to nrps, lsamp, and ddc, we now provide a list of driver genes (Supplementary Data 19) and associated pathways (Supplementary Data 20) are specific to the N4.3 neuron (formerly delta-3).\n14. For IF images, please add quantification and experimental replications in the figure legends.\nThanks, added.\n15. Line 347: “worms treated with lsamp dsRNA instead of ddc dsRNA lost their ability…” This statement would be more effective omitting mention of ddc dsRNA here. First, state what happens after lsamp dsRNA treatment, then say that ddc dsRNA treatment did not affect female sexual development.\nThank you. We have revised the text as suggested to make the statement more concise and effective in lines 341-348: “By employing RNAi in vitro, we found the male worms treated with lsamp dsRNA lost their ability to initiate female sexual development (Fig. 4a-b). These female worms paired with lsamp RNAi males displayed reduced Fast Blue BB staining in their vitellaria and possessed poorly developed ovaries (Fig. 4b). In addition, they barely laid eggs (Fig. 4c). The qRT-PCR result indicated a strong reduction of the lsamp expression in the male worms after the lsamp RNAi (Fig. 4d). In contrast, treatment with ddc dsRNA did not affect female sexual development (Supplementary Fig. 21). The above results suggest that male lsamp is essential for female sexual development.”\n16. Line 398: “Acetylated alpha-tubulin is a protein that stabilizes microtubules” is an awkward statement. Acetylation of alpha-tubulin (a subunit of the microtubule polymer) is a post-translational modification of microtubules (on alpha-tubulin) associated with stable microtubules.\nWe thank the reviewer for this insightful comment. We have revised the sentence to more accurately describe the relationship between α-tubulin acetylation and microtubule stability: “Acetylation of α-tubulin, a post-translational modification of microtubules, is associated with microtubule stabilization and regulates cell morphology in nerve fibers.” (lines 394-395).\n17. Line 469: “We… characterize gene regulatory programs…” This work did not characterize gene regulatory programs; it identified changes in gene expression associated with dimorphic sexual development.\nThanks. We have revised the wording to more accurately reflect our findings: \"In the present study, we used single-cell RNA sequencing to systematically profile gene expression dynamics underlying dimorphic gametogenesis in S. japonicum.\" (lines 469-471).\n18. Line 507: “…no active neurogenesis at the adult stage.” Single-cell sequencing does not suffice to make this claim. The population of neural progenitors may just be smaller than can be detected, given the number of cells sequenced here.\nThanks. We agree and have removed the corresponding statement.\n19. “Gender” is a social construct, so the authors should omit that word and refer simply to the parasite’s sex.\nThank you. The term “gender” has been replaced throughout the manuscript with “sex”.\nReviewer #4:\n1. I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely appreciate your valuable contribution in reviewing our manuscript.\nWe sincerely appreciate the reviewers’ positive evaluation and recognition of our work. We would like to thank the reviewers for their thoughtful and helpful critiques of our work. The work is greatly improved upon the incorporation of their suggestions.\nReviewer #1 (Remarks to the Author)\nThe authors have made extensive changes to the manuscript and have addressed all our previous points comprehensively. The new version of the manuscript has been greatly improved, including a streamlined analysis of vitelline cell developmental trajectories, neural lineages, and germline gene expression dynamics. We remain enthusiastic about its publication and just have some minor suggestions to consider:\nWe appreciate your affirmation that our revisions have comprehensively addressed the previous points and greatly improved the manuscript. We have implemented the minor suggestions as detailed below.\n1. Proportions of cell types\nAnalyses describing the proportions of cell types are still really prominent in the manuscript but we still don’t really feel that they add any new insights. The context for this comment is that the paper is huge with many figures, so we feel that making some reductions would further improve readability. The authors have swapped some figures between supplementary material and main text but we’re unclear how this improves the manuscript. figure 1c does confirm a couple of known things; for instance, it clearly shows that when vitellocyte cells are combined across the datasets, most come from mature females, as expected.\nThe figure (1c) legend text “…showing the relative proportions of key cell types in S. japonicum at different developmental stages” is a bit misleading. Each bar represents all cells of a given type from across multiple stages. It’s not the proportion of that type within the parasite. A better description might therefore be “…showing the relative proportions of cells from different stages within the cell types of the combined S. japonicum dataset”.\nThe composition analysis also doesn’t consistently support the narrative. For instance, it is true that “male … differentiated neural cells (e.g., non-ciliated neurons and ciliated neurons) continuously expanded” but the authors don’t mention that in females this is not the case.\nThank you for this constructive feedback regarding the prominence of cell-type proportion analyses and the overall readability of our manuscript. We have removed the Supplementary Fig 6a, and moved Fig 1c to the revised Supplementary Fig 6a. Correspondingly, we have shortened the related descriptions in the main text to reduce redundancy (lines 138-153).\nWe have also revised the figure legend of the revised Supplementary Fig 6a as suggested: “…showing the relative proportions of cells from different stages within the cell types of the combined S. japonicum dataset.” (in the updated Supplementary figure informations file)\nAdditionally, we have supplemented the description regarding the changes in neural cells in female worms: “In female worms, these populations did not show sustained expansion but rather exhibited an initial increase followed by a subsequent decline after 22 dpi, which is consistent with the weakening of neural activity in females after 22 dpi, as revealed by bulk RNA-seq analysis10 (Supplementary Fig. 6a-b, and Data 5).” (lines 150-153)\n2. Convergent RNA trajectories\nIn discussing the convergent RNA velocity trajectories, the authors now hypothesise that the convergence point is “a shared multipotent progenitor pool that functions as a developmental hub capable of giving rise to multiple downstream cell types in both males and females.” They conclude that “This finding offers new insight into the hierarchical organization of gonadal cells in this organism, suggesting the presence of a common —rather than sex-specific —early developmental regulatory mechanism.” This seems like a lot to infer from RNA velocity alone, without additional validation. Some reference to published precedents from other systems is needed\nThank you for this insightful comment. We agree that inferring shared progenitors from RNA velocity alone requires caution. Our hypothesis was proposed to explain the convergent pattern observed in both sexes, and it is strongly supported by precedent: in S. mansoni, Wang et al. identified δ-neoblasts—somatic stem cells defined by high nanos-2, fgfrA, and fgfrB expression—as multipotent progenitors for both sexes (PMID: 29988015). The convergent SSCs-2 cluster in our data shows a highly similar molecular signature. Therefore, our RNA velocity results likely reflect a conserved, shared multipotent progenitor pool analogous to δ-neoblasts, extending this model to S. japonicum.\nWe have cited the relevant paper (PMID: 29988015) and revised the Discussion section in lines 476-483: “Although our RNA velocity analysis of sex-related subsets revealed clear trajectories from GSCs to their progeny, it consistently showed convergence toward a central cell state identified as SSCs-2 in both sexes (Fig. 2b). SSCs-2 population shares high expression of key markers (e.g., nanos-2, fgfrA, fgfrB) with the multipotent δ-neoblasts, which constitute a shared somatic stem cell pool in S. mansoni18. We therefore hypothesize that this convergence point may represent a similar, shared multipotent progenitor pool, capable of giving rise to multiple downstream cell types in both males and females, which requires further experimental validation.”\n3. The last sentence of the results (“These results provide new insights into male-female interactions during sexual development in schistosomes and may inform future strategies for combating schistosomiasis ” ) is more of a discussion/conclusion statement. It’s also very similar to that last sentence of the Discussion, so I suggest removing it from Results.\nThank you. We have deleted this sentence.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely appreciate your valuable contribution in reviewing our manuscript.\nReviewer #3 (Remarks to the Author):\nThe authors have satisfactorily addressed the concerns of the referees and the revised manuscript is ready for publication.\nWe sincerely appreciate your recognition of our work.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely appreciate your thorough review."}]

Decision letter