All the experimental procedures involved were performed strictly in accordance with the protocols (code 2021 KT-53) approved by the Ethics Committee of Xiangya School of Basic Medical Sciences, Central South University, Changsha, China. The protocols of housing, breeding, and care of the animals followed the ethical requirements of the government of China.

Six M. fortis and six Kunming mice were used in this study. The Chinese mainland strain of S. japonicum, maintained by serial passage through O. hupensis snails and inbred Chinese Kunming mice in our laboratory, was used to infect the animals. For each species, three animals were infected with 100 cercariae, and three others were kept as a control group. The animals with free access to food and water were housed under specific-pathogen-free conditions at the experimental animal center. The peripheral blood was collected from uninfected and infected animals at 10 dpi. For sequencing cost reduction and processing a larger number of cells, we pooled the biological replicates to perform scRNA sequencing.

The PBMCs were isolated by density gradient centrifugation using Ficoll-Paque Plus medium (GE Healthcare) and washed with Ca/Mg-free PBS. 2 mL GEXSCOPE^® red blood cell lysis bufier (RCLB, Singleron) was added at 25 °C for 10 min to remove the red blood cells. The solution was then centrifuged at 500 g for 5 min and suspended in PBS. The samples were centrifuged at 400g for 5 min at 4 °C, and the supernatant was discarded. After removing red blood cells, the PBMCs were resuspended in PBS and centrifuged at 400 g for 10 min at 4 °C to collect the cells after discarding the supernatant. Finally, the samples were stained with Trypan Blue, and the cell viability was evaluated microscopically.

Single-cell suspensions (2×10⁵ cells/mL) with PBS (HyClone) were loaded onto a microwell chip using the Singleron Matrix^® Single Cell Processing System. Barcoding Beads were subsequently collected from the microwell chip, followed by reverse transcription of the mRNA captured by the barcoding beads to obtain cDNA and PCR amplification. The amplified cDNA is then fragmented and ligated with sequencing adapters. The scRNA-seq libraries were constructed according to the protocol of the GEXSCOPE^® Single Cell RNA Library Kits (Singleton). Individual libraries were diluted to 4 nM, pooled, and sequenced on Illumina Novaseq 6000 with 150 bp paired-end reads.

Raw reads from scRNA-seq were processed to generate a gene expression matrix using CeleScope (https://github.com/singleron-RD/CeleScope) v1.9.0 pipeline. Briefly, raw reads were firstly processed with CeleScope to remove low-quality reads, trim poly-A tail, and adapter sequences with Cutadapt v1.17 (Martin, 2011). STAR v2.6.1a was used to map the reads from the two species to the mouse reference genome GRCm38 (Ensembl version 92 annotation) (Dobin et al., 2013), because the current annotated reference genome of M. fortis (M_Fortis_MF-2015_v1) available in the National Center for Biotechnology Information (NCBI) is incomplete (Li et al., 2020). Each cell’s UMI and gene counts were acquired with feature Counts v2.0.1 software to generate expression matrix files for subsequent analysis (Liao et al., 2014).

For each sample dataset, we filtered the expression matrix by the following criteria: 1) cells with a gene count less than 200 were excluded; 2) cells with a top 2% UMI count were excluded; 3) cells with mitochondrial content > 20% were excluded; 4) genes expressed in less than 5 cells were excluded. The raw count matrix was normalized by total counts per cell and logarithmically transformed into a normalized data matrix. Seurat was used to cluster the merged object into subsets of cells. This workflow includes finding variable genes, running principal component analysis on variable genes, and running Uniform Manifold Approximation and Projection on Principal Components (UMAP) (Mangiola et al., 2021). In the first round of cells clustering (pre-clustering), resolution = 0.3 was used to identify major cell types, and resolution = 0.2 was used for the second round to identify mononuclear phagocytes (MPs), T and B cell subsets. Cell cluster identity was determined by finding DE genes for each cluster using Seurat’s FindAllMarkers () function and by comparing those markers to known cell-type-specific genes in the SynEcoSys database, combined with knowledge from the literature (Zhang et al., 2023) (Zhang et al., 2023). Differential gene expression analysis between conditions was performed using Seurat’s ‘findmarkers’ function based on the Wilcox likelihood-ratio test with default parameters. The genes with an avg_log2FC value greater than 0.25 and adjusted p-value less than 0.05 were considered differentially expressed genes (DEGs).

We further collected PBMCs from uninfected and infected M. fortis at 10 dpi to perform qPCR for the validation of some potential candidate genes. Total RNA was extracted from M. fortis PBMCs using a Trizol reagent (Rio et al., 2010).

SYBR Green using Master Mix on a CFX96 Touch System (BioRad) was used to determine the expression level of IL2, IL4, Cxcl9, and Fcgrt (Table 1) in uninfected and infected M. fortis. Beta-actin (ACTB) was used as an internal control, and the relative expression was calculated using the 2^-△△Ct method. The data obtained for analysis were reported as mean values and standard deviation (SD), with 95% confidence intervals (CI); P >0.05 was considered statistically significant. All statistical analyses were conducted using GraphPad Prism 10.4.1.

The total number of valid reads obtained in the different pools of samples varied from 349,110,435 to 412,975,949, and the mapping ratio varied from 70.68% to 80.92% (Table 2). After quality control, we recovered 20667 cells in M. fortis and 23147 cells in Kunming mice, respectively (Table 2). Based on the UMAP plot and the conserved marker genes, we grouped the PBMCs into T cells, natural killer (NK) cells, B cells, mononuclear phagocytes (MPs), and erythrocytes in both M. fortis and Kunming mice (Figures 1A–D). T cells were identified by the expression of Il7r and Tcf7, and NK cells were distinguished from T cells by the relatively high expression of Xcl1 and Nkg7 (Figures 1C, D). B cells were identified by the expression of Cd19 and Cd22, MPs by the expression of Irf7, Cd14, and Cebpb, and erythrocytes by the expression of Alas2 and Bpgm (Figures 1C, D). Cell type percentage comparison revealed that the percentage of T cells and MPs decreased in the blood of both M. fortis and Kunming mice at 10 dpi. In contrast, the percentage of B cells increased, suggesting that infection can induce dynamic changes in peripheral immune cell populations (Figures 1E, F).

Next, we considered the MPs to dissect their heterogeneity. Based on the conserved markers, we divided the MPs into monocytes, characterized by the expression of Cd14, and DCs, characterized by the expression of Flt3 (Figures 2A–D). To identify additional markers, we generated a heatmap displaying the top differentially expressed genes between monocytes and DC in the two species. We found that Fn1 and C1qb were highly expressed in M. fortis monocytes, while Kunming mice monocytes highly expressed Ace and Fcgr4. Similarly, M. fortis DC highly expressed Fxyd5, Lsp1, and Napsa, while the DC from Kunming mice highly expressed Ccr2, S100a6, F13a1, and Sell (Figures 2E, F). By comparing the percentage change of the MPs before and after infection, we observed that the percentage of DCs increased after infection, but that of monocytes decreased in both M. fortis and Kunming mice (Figures 2A, B). However, there were more DCs in the PBMCs of M. fortis when compared to Kunming mice, while the number of monocytes was more elevated in Kunming mice (Figures 2A, B).

We further performed DEG analysis in monocytes and DC to identify the genes responding to the invading worm. Interestingly, we found that the Ccr2 was upregulated in both M. fortis and Kunming mice monocytes (Supplementary Files 1, 2), but Cxcl9 was only detected as DEG in M. fortis (Supplementary File 1). Regarding DC, we observed that infection induced the upregulation of Irf7 and Fcgrt in M. fortis DC (Supplementary File 3), but not in Kunming mice, suggesting that the two species exhibit distinct responses to the invading worms (Supplementary File 4).

We further performed unsupervised clustering on T cell data to identify the T cell subsets present in PBMCs of the two species. Based on the expression of Cd3d, Cd4, Cd8a, Cd44, Maf and Prf1, we grouped the T cells from M. fortis into four T cell subsets, including CD4 naive T cells, CD4 effector T (CD4 Teff) cells characterized by the specific expression of Maf, cytotoxic T cells characterized by the expression of Cd8a and Prf1, and double negative T (DNT) characterized by the lack of Cd4 and Cd8a markers but expressing Cd3d (Figures 3A, C). Similarly, we grouped Kunming mice T cells into 5 subsets, including those identified in M. fortis and CD8 naive T cells expressing Cd8a but not expressing the effector molecule Prf1 (Figures 3B, D). By comparing the T cell subpopulations across conditions, we observed that the percentage of CD4+ T effector cells increased in M. fortis at 10 dpi, but not in Kunming mice. In contrast, the percentage of cytotoxic T cells decreased in both M. fortis and Kunming mice (Figures 3E, F).

We further considered the CD4 Teff cells to decipher their heterogeneity. Based on the UMAP plot and marker genes, we identified 3 distinct clusters in M. fortis: a cluster with high expression of Tnf that we annotated as Th1 cells, a cluster highly expressing Gata3, IL4, and IL5 that we annotated as Th2, and a cluster highly expressing Foxp3 and Ctla4 that we annotated as Treg (Figures 4A, C, E). By performing similar analysis in Kunming mice, we identified 2 stable clusters: a cluster with a relatively high expression of Tnf, IL4 and IL5, but low expression of Gata3 that we annotated as Th cells because the unsupervised clustering and marker genes failed to distinguish the Th cell subsets, and a cluster with high expression of Foxp3 and Ctla4 that we annotated as Treg (Figures 4B, D, F). By comparing the cell type percentages before and after infection, we observed that the percentage of Th2 cells increased, while that of Th1 cells decreased after infection.

We further merged the data from clusters expressing CD4 markers and those expressing CD8a markers to perform a DEGs analysis across conditions, providing a better understanding of the transcriptome profiles of CD4 and CD8 T cells. Interestingly, we found that IL2 and IL4 were upregulated in M. fortis CD4 T cells (Supplementary File 5), while they were not detected as DEGs in Kunming mice (Supplementary File 6). In addition, we observed that Ccr2 and Tcf7 were upregulated in M. fortis CD8 T cells (Supplementary File 5), but not detected as DEG in Kunming mice CD8 T cells (Supplementary File 6).

To determine whether the humoral immune response was significantly activated at 10 dpi, we characterized B cells in the two species. Based on the marker genes, we grouped the B cells into three subsets in both M. fortis and Kunming mice. These include naive B cells characterized by the high expression of Cd19, Cd79a, Pax5, Ms4a1 and Cr2, marginal zone (MZ) B cells characterized by the high expression of Cd44, Cd80 and Cd86 when compared to naive B cells, and antibody secreting cells (ASC) characterized by the high expression of Cd27 and Cd38 and decreased expression of Cd19, Cd79a, Pax5 and Ms4a1 (Figures 5A–D). By comparing the cell type percentage changes across conditions, we did not observe considerable alteration among the B cell subsets in M. fortis; however, the total number and percentage of naive B cells increased in Kunming mice at 10 dpi (Figures 5A, B).

As ASCs are the terminally differentiated B cells that play the key role in humoral immunity, we further screened the DEGs across conditions in this B cell subset. We found that the immunoglobulin genes Igha, Ighg1, and Ighg3 were upregulated in M. fortis at 10 dpi, while IghM was downregulated (Figure 6A). However, those detected in Kunming mice, including Igha, Ighg2c, were downregulated (Figure 6B).

Due to the lack of M. fortis-specific antibodies to perform flow cytometry, we selected 4 potential candidate genes, including IL2, IL4, Cxcl9, and Fcgrt, and examined their expression patterns by qPCR in M. fortis to validate our scRNA-seq analysis. The qPCR results demonstrated that the expression level of the selected genes increased significantly after infection, suggesting the reliability of the scRNA-seq data (Figure 7).