L. crocea (average weight: 50 g) were obtained from a mariculture farm in Ningde, Fujian, China. Fish were acclimated in 100-L fiberglass-reinforced plastic tanks equipped with a recirculating aquaculture system (24 ± 2°C, UV-sterilized seawater) and fed daily with commercial pellets. After a 14-day acclimatization period, healthy individuals were selected for subsequent experiments. All procedures adhered to the Laboratory Animal Guideline (GB/T 35892–2018) and were approved by the Animal Ethics Committee of Fujian Agriculture and Forestry University. All fish were anesthetized using tricaine methanesulfonate (MS-222; MCE) before surgical dissection.

C. irritans tomonts were harvested from naturally infected L. crocea and maintained in sterile seawater at 8°C with daily water renewal (23). The excystation of tomonts was induced by transferring them to 26°C for 72 hours. Freshly released theronts, less than 2 hours post-hatching, were quantified and used to infect fish. The experimental group was challenged with 3,000 theronts per fish through a 2-hour static immersion, while the control group underwent identical handling procedures without parasite exposure. Following the challenge, fish were transferred to separate flow-through tanks maintained at 24 ± 2°C. Gill tissues were collected at 60 hours post-infection (peak trophont maturation phase) for leukocyte isolation.

Gills were aseptically dissected, minced into 1-mm³ fragments, and digested in cold DMEM containing 1 mg/mL collagenase I, 1 mg/mL collagenase IV, 0.1% (w/v) NaCl, and 1× penicillin/streptomycin at 4°C for 2 hours. Tissue homogenates were filtered through a 70-μm cell strainer (BD Biosciences), washed twice with DMEM supplemented with 2% fetal bovine serum (FBS), and separated using a 34% Percoll gradient by centrifugation at 650 × g for 30 minutes at 4°C. Pelleted cells were treated with red blood cell (RBC) lysis buffer, washed, and assessed for viability (>90% by trypan blue exclusion). Cell counts were quantified via flow cytometry.

Cell suspensions from the gills of five fish per group were pooled and adjusted to a concentration of ~1,000 cells/μL. These preparations were processed using the Chromium™ Single Cell 3’ Kit (10× Genomics). Gel bead-in-emulsions generation, cDNA amplification, and library preparation were performed following the manufacturer’s protocols. Libraries were sequenced on an Illumina NovaSeq 6000 platform (150-bp paired-end) by Gene Denovo Biotech (Guangzhou, China).

Raw sequencing data were aligned to the L. crocea reference genome (GCF_000972845.2) using Cell Ranger (v3.1.0). Doublets were removed using DoubletFinder (v2.0.3). The cells were excluded if they met any of the following criteria: <500 or >4,000 detected genes, >20,000 UMI counts, or >10% mitochondrial reads. Post-filtering, data were normalized, integrated, and clustered using Seurat (v4.0.4) with 30 principal components at a clustering resolution of 0.5. Cell clusters were visualized using Uniform Manifold Approximation and Projection (UMAP). Cell identities were annotated based on conserved immune cell markers.

For comprehensive identification of cluster-specific markers across all cellular subsets, the Wilcoxon rank-sum test was implemented through the FindAllMarkers function, comparing each cluster against the rest cells from all other clusters. Significantly upregulated genes were identified based on the following criteria: (1) genes had to be expressed at least 1.28-fold higher in the target cluster compared to other clusters; (2) genes must be detected in more than 25% of the cells within the target cluster; (3) an unadjusted p-value < 0.05. To investigate condition-specific transcriptional changes within defined immune subclusters, intra-cluster differential expression analysis was performed between control and infected groups using a hurdle model in MAST (Model-based Analysis of Single-cell Transcriptomics). Differentially expressed genes (DEGs) were defined as meeting all of the following thresholds: (1) absolute log₂-fold change (log₂FC) ≥ 1; (2) adjusted p-value ≤ 0.05 (Benjamini-Hochberg correction); (3) detection in ≥ 20% of cells within the target cluster.

DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using a hypergeometric test with false discovery rate (FDR) correction (FDR < 0.05). Enriched terms were mapped to biological processes, molecular functions, or pathways.

To delineate the immune landscape of L. crocea gills, we conducted comparative scRNA-seq profiling of leukocytes isolated from control and C. irritans-infected specimens ( Figure 1A ). Sequencing achieved an average of 31,660 reads/cell (997 genes/cell) and 37,067 reads/cell (980 genes/cell) in the control and infected samples, respectively. Following stringent quality control ( Supplementary Figure S1 ), 16,832 high-quality transcriptomes (control: 9,028; infected group: 7,804) were retained for further analysis. Unsupervised clustering delineated 25 transcriptionally distinct cell states, validated through UMAP visualization and cluster-specific marker heatmaps ( Supplementary Figure S2 ).

The 25 clusters were annotated into eight immune lineages and non-immune contaminants ( Figures 1B, C ; markers in Supplementary Table S1 ). T lymphocytes occupied clusters 0, 2, 3, 4, 6, 7, and 10, defined by the expression of canonical T-cell markers cd3ξ, zap70, and lck. B-cell compartments comprised clusters 5 (B-5) and 13 (B-13), identified by pax5, cd79a, and cd79b. Myeloid populations, including macrophages (Cluster 12 with the expression of csf1rb ⁺/cfp ⁺/macro ⁺/cd68 ⁺), neutrophils (Cluster 14: mpx ⁺/csf3r ⁺/mmp9 ⁺), DCs (Cluster 19: cd74a ⁺/xcr1 ⁺/zbtb46 ⁺), and cpa5 ⁺ cells (Cluster 21: cpa5⁺ ), were identified. Notably, cpa5 is specifically expressed in the zebrafish mast cell lineage (24). However, multiple mast cell markers (including kita and tps) were not detected (data not shown) in cpa5 ⁺ cells. These suggest that these cells can be classified as cpa5 ⁺ granulocytes and potentially represent a teleost-specific granulocyte subset with partial functional overlap with mast cells. Innate lymphoid diversity spanned NK-like cells (Clusters 17-18: nkl ⁺/gzma ⁺/gzmb ⁺), and ILC2-like cells (Cluster 16: gata3 ⁺/il4/13a ⁺/il4/13b ⁺). A total of 3,762 non-immune cells (Clusters 1, 8, 9, 11, 15, 20, and 22), which lacked cd45 expression, were excluded from immune analyses. Transcriptomic analysis identified the top 20 differentially expressed genes per immune lineage, revealing both evolutionarily conserved and teleost-specific immune signatures ( Supplementary Table S2 ). Strikingly, T cells constituted 71.6% (9,360 cells) of immune cell population ( Figure 1D ; Supplementary Table S3 ), establishing their central role in piscine mucosal immunity. In mammals, T cells comprised functionally distinct subsets, which mediate specialized roles in response to the parasite infection (25). This dominance provides the cellular substrate for exploring the diverse T cell subset responses (proliferation, differentiation, and cytokine production) during C. irritans infection.

Our analysis partitioned T cells from two biological samples into 14 clusters (designated T0–T13) using a clustering resolution of 0.5 ( Supplementary Figure S3A ). Cluster validation through visualization of the top five variably expressed genes revealed a clear demarcation between all subsets in the heatmap analysis ( Supplementary Figure S3B ). Differential gene expression profiles across clusters, integrated with typical T-cell lineage markers, allowed for annotation of definitive subsets ( Figure 2 ; DEGs in Supplementary Table S4 ; markers in Supplementary Table S5 ). The cluster T7 was classified as CD8⁺ cytotoxic T lymphocytes (CTLs) by expressing cd8a, cd8β, nkl, gzmb, and prf1 ( Figure 2A ). CD4⁺ T cell populations segregated into four functionally distinct subgroups, including regulatory T cells (Tregs), Th17 cells, Th1 cells, and naïve T cells ( Figures 2A, B ). The results found that the cluster T8 was annotated as Tregs (expressing foxp3 and tgfb1), T6 as Th1 cells (expressing lineage-defining transcription factor tbx21 and effector cytokines ifng and tnfa), T1 and T4 as Th17 cells (expressing characteristic markers rorc (RORγt), il17a/f1, il1r1, and il23r), T9 as γδ T cells (expressing tcrgc and a transcription factor sox4 (26), lacking CD4 and CD8 receptors), T2 as CD4^-CD8^- double-negative T cells (DNTs; lacking CD4 and CD8 expression), T5, T11, and T12 as proliferating T-cell subset (expressing DNA replication genes pcna, cdk1, and top2a). The ten most statistically significant subset-defining genes across all T cell populations showed the potential markers for T cell subsets in fish ( Supplementary Table S6 ).

To investigate the T cell-mediated immune responses in gill tissue during pathogen challenge, we first assessed the redistribution of T cell subsets after C. irritans infection. At 60 hours post-infection, the number of proliferating T cells (T5, 11, and 12) increased from 8.28% to 11.59% ( Figure 3A ). Concurrently, an upregulation of pcna (encoding the proliferating cell nuclear antigen that enhances DNA polymerase processivity during replication) was observed in infected fish ( Figure 3B ). These data demonstrate the induction of T cell proliferation in gill tissue by C. irritans infection. Although overall T cell number decreased, mainly due to a sharp increase in myeloid cells, Treg numbers increased 2.89-fold (from 1.74% to 5.03%) after infection ( Figure 3C ). This Treg amplification was supported by significant induction of lineage-defining markers foxp3 and tgfb1 following pathogen exposure ( Figure 3D ). Furthermore, the proportion of naïve T cells (T0 and T3) decreased significantly, from 39.71% to 26.67% of total T cells ( Figure 3E ). Notably, Th1 and Th17 subsets remained stable after infection ( Figure 3F ), supported by the unchanged expression of their signature transcriptional regulators (tbx21 and rorc) and effector cytokines (ifng and il17a/f1). These results suggest preferential differentiation of naïve T cells into immunosuppressive Tregs rather than effector subsets during C. irritans infection.

DNTs exhibit functional plasticity through cytokine secretion (IL-4, IL-17, IFNγ, and TNFα) mirroring the activities of CD4⁺ T helper cells, termed helper DN T cells or Th-like DN T cells (27, 28). These functional DNTs demonstrate context-dependent protective or pathogenic roles across infection models (27, 29). Our results revealed that L. crocea DNTs (LcDNTs) accounted for 12.48% of total gill T cells under homeostasis, increasing to 14.57% after C. irritans challenge ( Figure 4A ). Comparative transcriptomic profiling across CD4⁺, CD8⁺, and DNT populations through integrative analysis identified subset-defining molecular signatures ( Supplementary Figure S4A ). Gene ontology analysis revealed DNT-specific enrichment of regulatory pathways, particularly metabolic process modulation and systemic biological regulation ( Supplementary Figure S4B ). LcDNTs were prominently engaged in cytokine-cytokine receptor and C-type lectin signaling pathways, as evidenced by differential expression of key mediators including cytokine receptors (il2rb, il1r1, and il17rα), death receptors (dr4), chemotaxis regulators (ccr6 and ccl20a.3), inflammatory cytokines (il22 and tnf), and signaling components (fcrγ, mapk1, pik3r1, and ikba) ( Figure 4B , Supplementary Table S7 ). Notably, LcDNTs displayed Th17-polarized features, characterized by the elevated expression of lineage-defining transcription factor (rorc) and (molecules il17a/f1, il1r1, and il23r) compared to conventional CD4⁺ and CD8⁺ T cells ( Figures 4C , 2A ). Parasitic challenge induced the expression of il17a/f1 production ( Figure 4D ), as well as cytotoxic effector (gzma) and immunomodulatory factors (ccl5l, ccl20, il8, tnf, and tnfl) in LcDNTs ( Figure 4D ; Supplementary Table S8 ), confirming a functional Th17-like programming and dual cytotoxic and regulatory capacities, respectively. Systematic analysis of alternative Th-like DNT subtypes (e.g., Th1 and Th2) demonstrated minimal expression of lineage-specifying transcription factors (tbx21 and gata3), supporting a predominant Th17-like polarization phenotype.

In mammalian systems, γδ T cells represent a functionally heterogeneous lymphocyte population with specialized roles in antimicrobial defense and immunomodulation. These cells exhibit remarkable plasticity, differentiating into cytotoxic effectors and helper-like subsets (Th1/Th2/Th17 analogs) upon activation while producing diverse cytotoxic mediators and cytokines to modulate the immune response (30). Our results identified L. crocea γδ T cells (cluster T9) as multifunctional effectors, demonstrating enriched expression of chemotaxis regulators (e.g., cxcl10, mcp1b, and ccl3), endocytosis and antigen processing components (e.g., cxcr4, chmp4b, arf6, cd74, h2-q9, rt1-b, and hla-drα), cytotoxic effectors (e.g., gzmb), and inflammatory mediators (e.g., tnf, tnfl, and il17a/f3) ( Supplementary Tables S9 , S10 ). Unique adhesion molecules cd2l and cldn4 (claudin-4) were specifically expressed in γδ T cells, potentially mediating intercellular interactions within gill microenvironments ( Supplementary Tables S9 , S10 ). The data revealed type 2 immune polarization in L. crocea γδ T cells, as evidenced by elevated expression of gata3, ap-1, cxcr4, and type 2 cytokines (il4/13a and il4/13b, homologous to mammalian il13 and il4, respectively) ( Figure 5A ). Parasitic challenge triggered a significant increase in the number of γδ T cells, ranging from 1.98% to 3.61% of total T cells ( Figure 5B ). Results also found that γδ T cells in the gills of controls exhibited minimal production of cytokines (such as tnf, il17a/f3, and tnfl). In contrast, the upregulation of cytotoxic effectors (gzma and gzmb), Th1/Th17-polarizing cytokines (tnf, tnfl, faslg, and il17a/f3), and Th2-associated transcriptional regulators (gata3) and cytokines (il4/13a and il4/13b) was found in infected fish ( Figure 5C ).

Furthermore, the CD8⁺ CTL population (cluster T7) was found to express cytolytic effectors (nkl, gzmb, gzma, and prf1), immunoregulatory cytokines (il34 and tnfsf14), and chemoattractant cxcl10 ( Figures 2A , 5D ). Although the proportion of CD8⁺ CTLs decreased from 5.48% to 4.13% after infection ( Figure 5B ), these cells exhibited enhanced expression of cytolytic genes (nkl, gzmb, gzma, and prf1) and chemokines (ccl5l, ccl20, mcp1b, and ccl4l), while decreased the expression of il34 ( Figure 5E ). This dichotomy suggests CD8⁺ CTLs balance direct cytotoxicity (gzmb, gzma, and prf1) with immunoregulatory chemokine signaling during anti-parasitic responses.

A set of genes highly or specifically expressed in LcILC2-like cells (LcILC2s) was identified across samples, establishing a core transcriptomic signature for ILC2 populations in teleost gills ( Supplementary Table S11 ), including transcription factor gata3 (essential for ILC2 lineage development) (31) and key effector cytokines il4/13a and il4/13b ( Figures 1B , 6A ; Supplementary Table S12 ). Notably, these cells expressed the st2 (il1rl1, encoding the IL-33 receptor), a hallmark surface marker for ST2⁺ ILC2s, and batf, a crucial regulator of IL-25-responsive migratory ILC2 cell fate and function ( Figure 6A ) (32, 33). Conservation of intracellular signaling was demonstrated through enrichment of MAPK pathway components (flt3, dusp5, mapk1, myd88, mknk1, and mras), which control proliferation and cytokine production in mammalian ILC2s ( Figure 6B ) (34–36). Intriguingly, typical mammalian ILC2-activating receptors (il-7r, il-9r, tslpr, and il25r) (36) were undetectable in LcILC2s, potentially reflecting teleost-specific adaptations. Instead, this population uniquely expressed endocytic receptors (cd209d, cd163l, and cd163), along with cytokine receptor (il21r) ( Figure 6A ). Following parasitic challenge, LcILC2s underwent an 8.7-fold expansion (from 0.43% to 3.74% of gill leukocytes; Figure 6C ). The infection activated the upregulation of cytokine expression (il4/13a: 2.6-fold; il4/13b: 5.4-fold), and downregulation of il1b and il8 ( Figure 6D ).

Within innate lymphoid populations, two cytotoxic NK cell subsets (NK-1 in cluster 17 and NK-2 in cluster 18) were characterized. Both subsets expressed core effector molecules (nkl, faslg, and prf1), chemotactic factors (e.g., regakine-1 and cxcl10), and immunoregulatory properties (e.g., timd4l, ptpn22, btlal, and crtam) ( Figures 1B , 6E ). NK-1 mainly expressed gzma and Fc receptor-related genes (fcrl5, fcgr1l, and fcgr2l), whereas NK-2 expressed gzmb and cytokines (mcp1b, tnf, and tnfl) ( Figure 6E ). This transcriptional divergence suggests compartmentalized functional roles between the cell subsets. The infection increased the proportion of both NK-1 (from 1.69% to 2.06%) and NK-2 (from 1.02% to 2.42%) ( Figure 6C ) and upregulated the pro-apoptotic effector (faslg) and chemokines (ccl5l and ccl4l) ( Figure 6F ). Subset-specific responses included an elevated expression of gzma and traf5 in NK-1 compared with prf1 and tnf upregulation in NK-2 ( Figure 6F ). These differential profiles suggest distinct cytotoxic mechanisms and compartmentalized immunoregulatory roles for NK-1 and NK-2 during C. irritans infection.

Granulocytes are considered to be important effector cells in the protection of the gill mucosa of teleosts (8). However, the specific roles of individual subsets during parasitic infection are not well understood. To address this gap, we compared neutrophils and cpa5 ⁺ granulocytes to elucidate their functions in the gills of L. crocea during C. irritans challenge. Among 17,781 detected genes, 549 DEGs were identified in neutrophils, compared to 243 DEGs in cpa5 ⁺ granulocytes. GO enrichment analysis revealed that neutrophils were predominantly associated with proton transport (GO:0006818), hydrogen transport (GO:0015992), and macromolecular complex assembly (GO:0034622) ( Supplementary Figure S5A ). The cpa5 ⁺ granulocytes exhibited enrichment in processes related to single-organism signaling (GO:0044700), signaling receptor activity (GO:0023052), and cell communication (GO:0007154) ( Supplementary Figure S5B ). KEGG pathway analysis revealed significant enrichment of phosphatidylinositol signaling (ko04070), insulin signaling (ko04910), and C-type lectin receptor signaling pathways (immune recognition) in cpa5 ⁺ granulocytes, while oxidative phosphorylation (ko00190), proteasome activity (ko03050), phagocytosis (ko04145), and endocytosis (ko04144) pathways in neutrophils ( Supplementary Figures S5C, D ). In addition, cpa5 ⁺ granulocytes are involved in the expression of immune-related genes (such as clec4el, il6rβ, and cd97) and immunoregulatory factors (ccl8l, il16, and tgfb1), while neutrophils are related to the expression of chemotaxis receptors (cxcr1), inflammatory mediators (il1β and il8), chemokines (cklf and lect2l), and respiratory burst components (rac2) ( Figure 7A ; Supplementary Table S13 ). Additionally, protease profiling revealed lineage-specific enzymatic arsenals (mmp17l, mmp13, mmp9, cpe, ctsl, and alox15b), further supporting their different roles in immune defense ( Figure 7A ). These findings highlight significant functional differences between the two granulocyte populations, with neutrophils specializing in phagocytosis and oxidative killing and cpa5 ⁺ granulocytes involved in immunomodulation and pathogen recognition.

Granulocytes were present in limited numbers in the gills of healthy fish (neutrophils: 2.31%; cpa5 ⁺ granulocytes: 0.86%) ( Figure 7B ). After C. irritans challenge, neutrophil counts increased to 5.07% (a 2.2-fold increase), while cpa5 ⁺ granulocyte counts remained stable ( Figure 7B ), suggesting the pathogen-specific mobilization of phagocytic defenses. This stability suggests that cpa5 ⁺ granulocytes may act as tissue-resident sentinels or immunomodulators, rather than undergoing rapid recruitment like neutrophils, potentially providing sustained regulatory signals or acting as a reservoir for later effector responses. Infection-responsive transcriptomic shifts revealed metabolic activation in both lineages, with upregulated oxidative phosphorylation (ko00190) The pathway Ribosome, which extends beyond the scope presented in Figure 7C (top 10), is not displayed in Figure 7C , despite the presence of significant enrichment. Therefore, we have removed this pathway, and Proteasome (ko03050) pathways ( Figures 7C, D ). The infection induced upregulation of chemoattractant (cklf) and respiratory burst effectors (rac2 and mpx), while suppressing inflammatory mediators (il8, tlr5, and clec4el) ( Figure 7E ). This dual response may indicate feedback regulation to mitigate tissue damage during late infection phases.

KEGG pathway analysis revealed an upregulation of genes related to lysosomal and phagosomal pathways in macrophages, DCs, and B-5 cells (Cluster 5) ( Figure 8A ). Key genes within these pathways included cathepsins (ctsb, ctsk, and ctsz), lysosomal membrane proteins (limp2, npc2, and laptm4a), glycosidases (hexb), vacuolar ATPase subunits (atp6s1, atp6l, and atp6n), NADPH oxidase components (ncf1, nox2, and cyba), and vesicular trafficking regulators (rab5c, rab7, and rab11b) ( Supplementary Table S14 ). These data suggest that B-5 cells, similar to macrophages and DCs, have phagocytic capacity. Macrophages expressed lineage-specific receptors (e.g., marco, mrc1, and csf1rb), complement system components (e.g., c1qbl, c2, and cfb), antigen presentation molecules (e.g., hla-drα, rt1-b, and cd74), and antimicrobial protein gene (bpi), indicating their roles in complement activation, antigen processing/presentation, and pathogen clearance in the gills ( Figures 8B , 1B ). DCs were characterized by transcription factors essential for their development (flt3, batf3, and znf366) (37), pattern recognition receptor (tlr7- a viral ssRNA sensor) (38), chemokines (ccl17 and ccl19l), and antigen presentation machinery (e.g., hla-drα, rt1-b, and cd74) ( Figure 8B ). B-5 cells exhibited antigen presentation activity (e.g., hla-drα, rt1-b, and ciita) and ribosome biogenesis signatures (e.g., rpl23, rpl39, and rpl37a) ( Figure 8C ), similar to B1-like cells in zebrafish (39). In contrast, B-13 cells showed features of plasma cells, with an enrichment in ER stress pathway activation (e.g., pdia4, erlec1, and txndc5), and upregulation of protein export machinery (e.g., srp14, hspa5, and srp19), immunoglobulin synthesis (ighm and ight), and differentiation-driving transcription factors (irf4l and prdm1a) (40) ( Figure 8C ).

After C. irritans infection, the abundance of DCs, B-5, and B-13 cells increased from 1.20% to 1.66%, 6.81% to 11.60%, and 2.08% to 5.65%, respectively, while macrophages decreased (from 3.89% to 3.39%) ( Figure 8D ). Results found upregulation of pathogen-recognition receptors (clec4e, cd209c, cd209e, and tlr2) in macrophages, DCs, or B-5 cells in infected fish compared to controls ( Figures 8E, F ). Macrophages in the gills of infected fish exhibited increased expression of tnfsf14 (a cytokine that promotes T cell proliferation (41), 8.64-fold), ccl19l (a chemokine critical for naïve T cell recruitment (42), 3.37-fold), and txndc17 (a redox homeostasis regulator (43), 4.22-fold), compared to the controls. DCs showed elevated expression of tnfsf14 (15.06-fold), cd81 (MHC-II compartment regulator, 4.58-fold) and mif (a macrophage migration inhibitor, 3.7-fold) in infected fish compared to controls. Both B-5 and B-13 cells exhibited enhanced antibody production (iglcl and ighm) and upregulated expression of the T cell chemoattractant (ccl5l) in the infected fish ( Figure 8F ). Notably, B-5 cells specifically upregulated a lymphoid-organizing chemokine responsible for B cell follicle formation (44) (cxcl13) and a Th1-associated chemokine mediating T/NK cell recruitment (cxcl10) in the infected fish ( Figure 8F ). These results suggest coordinated roles for B cell subsets in adaptive immune regulation and lymphoid tissue remodeling in the gills of L. crocea during C. irritans infection.