The medaka (Oryzias latipes) lines were obtained from the National Bioresource Project (NBRP) Medaka and maintained at Utsunomiya University: standard line OK-cab (strain ID: MT830) and cab-Tg (rag1:egfp), a transgenic medaka expressing egfp under the control of an immature lymphocyte-specific rag1 promoter (strain ID: TG848) (34). Fish were maintained in temperature-controlled tanks (26°C) with a water circulation system under a 14 h light/10 h dark cycle. This study was conducted in accordance with the ethical guidelines of the Utsunomiya University Animal Experimentation Committee and the experimental protocols were approved by the committee (approval no. A23-0011). The developmental stages of medaka were designated as previously described (37).

The Alt-R CRISPR-Cas9 crRNA targeting exon 4 of rag1 (ENSORLG00000011969.3), Alt-R CRISPR-Cas9 tracrRNA, and Cas9 protein were purchased from Integrated DNA Technologies (IDT; Singapore). The sgRNA sequences are listed in Supplementary Table S1 . RNA–protein complexes were prepared according to the manufacturer’s instructions. Glass capillaries with 1.0 mm outer diameter were pulled (temperature 62°C, force: two light and two heavy weights provided by the manufacturer) using a micropipette puller PC-10 (Narishige Instruments, Tokyo, Japan). The capillary was filled with crRNA (0.75 mM), tracrRNA (1.5 mM), and 0.25 mg/mL Cas9 protein. Then, 1–2 nL of the mixture was injected into 1-cell Cab embryos using FemtoJet 4i (Eppendorf, Hamburg, Germany). G0 fish were crossed with wild-type (WT) Cab fish to establish the rag1 mutant strains. rag1^del8 was chosen for the phenotypic analysis. The primers used for genotyping are listed in Supplementary Table S1 . The mutant line was deposited in NBRP medaka as rag1 (del8); MT1583.

The 4.7kb of the medaka lck promoter, upstream of the ATG initiation codon situated in exon 2, was amplified from the genomic DNA of the Hd-rR strain. egfp and SV40 polyA were amplified from pEGFP-1 plasmid (Takara Clontech). Primer sequences are listed in Supplementary Table S1 . Amplicons of the lck promoter and egfp+SV40 polyA were cloned into the linearized pUC19 plasmid using In-Fusion HD (Takara Clontech). Accordingly, the lck:egfp construct, including the medaka lck promoter, egfp, and SV40 polyA, was generated. 25 ng/µL lck:egfp construct was injected into fertilized eggs of the cab strain for transgenesis using FemtoJet 4i (Eppendorf, Hamburg, Germany). Embryos around the day of hatching were observed under a fluorescence microscope and those exhibiting EGFP signals in the pharyngeal region as was previously reported (38) were further used to establish lck:egfp transgenic lines. Two stable transgenic lines were generated, one of which was used for further analysis. The transgenic line was deposited in NBRP Medaka as cab-Tg (lck-egfp); TG1582.

Larvae and adult medaka were anesthetized with 0.1% (v/v) 2-phenoxyethanol and immobilized in 3% methylcellulose. Images were captured using a DFC450c digital camera (Leica, Houston, TX, USA) under an M205FA fluorescence stereo microscope (Leica, Houston, TX, USA) and acquired using Leica Application Suite X (LASX) v1.1.0 (Leica, Houston, TX, USA). The area and fluorescence intensity of the images were measured using the ImageJ software (National Institutes of Health, MD).

Whole bodies of the adults were fixed in 4% Paraformaldehyde/PBS overnight at 4°C. Paraffin sectioning was performed as previously described (39). Transverse sections (5 µm thick) were stained with hematoxylin and eosin. Images were captured using a DP73 digital camera (Olympus, Tokyo, Japan) under a BX60 microscope (Olympus, Tokyo, Japan).

RNA was extracted using a NucleoSpin RNA kit (Macherey-Nagel, Düren, Germany) with on-column DNase treatment. RNA was reverse-transcribed using a QuantiTect Reverse Transcription Kit (Qiagen). To detect immunoglobulin mu (igm) and t-cell receptor beta chain (tcrb) cDNA after VDJ rearrangement, nested PCR was performed using primer sets targeting the variable (V) and constant (C) regions. Two V regions that were highly expressed in wild-type kidney were chosen for RT-PCR analysis (data not shown). Amplification was performed in a thermal cycler using the following program: 1 min at 95°C; cycles of 10 s at 95°C, 30 s at 59°C, 30 s at 72°C; and 5 min at 72°C. The cycle numbers were as follows: actb (actin beta), 25 cycles; igm (igV_H-Cm) and tcrb (tcrVb-Cb2), 30 cycles plus nested 30 cycles. The primers used for the RT-PCR are listed in Supplementary Table S1 .

Quantitative PCR (qPCR) was performed using the THUNDERBIRD SYBR qPCR Mix (Toyobo, Osaka, Japan) and LightCycler 96 (Roche, Basel, Switzerland). Primers used for qPCR are listed in Supplementary Table S1 .

Whole kidney marrow (WKM) cells and thymocytes were prepared as previously described (36). Cells were obtained by pipetting the organs into 1 mL ice-cold 1% fetal bovine serum (FBS) in 0.9× PBS, washed with 1% FBS in 0.9× PBS by centrifugation, and then filtered using a 40µm-stainless mesh. Analytical flow cytometry (FCM) was performed on a FACSLyric flow cytometer (BD Biosciences). Populations of lymphocytes, myelomonocytes, and precursors in WKM cells in light scatter profiles were defined using previous information in zebrafish and medaka (40, 41). Data analysis was performed using the Flowjo10.6.2 (BD Biosciences).

FACS Aria III (BD Biosciences) was used for cell sorting. The cell populations in the light scatter profiles of WKM cells and thymocytes were defined according to previous studies (34, 40, 41).

The single-cell RNA-seq library was constructed using the Chromium Controller and Chromium Next GEM Single Cell 3′Reagent Kits v3.1 (10x Genomics) following the standard manufacturer’s protocols. The FACS-sorted single-cell suspension was immediately loaded onto 10x Chromium iX to recover 10000 cells, followed by library construction. The library was sequenced using the Novaseq6000 system (Illumina), according to the manufacturer’s instructions. Sequencing was performed using a 28/90bp paired-end configuration. The obtained dataset was processed with Cell Ranger (v7.2.0). Oryzias_latipes.ASM223467v1.dna.toplevel.fa.gz (Ensembl) and olatipes-HdRr-ensembl-peaks2utr.v1.gtf.gz (https://figshare.com/articles/dataset/Medaka_Gene_Model_version_1/24080463?file=42403659) were used as the reference genome and annotation files, respectively.

Loupe Browser v7.0.1 and v8.1.2 were used to characterize cell types. Barcodes with total Unique Molecular Identifier (UMI) counts of 2048–32768 (2¹¹–2¹⁵) were used for this analysis. Cell characterization was based on gene expression profiles in t-distributed Stochastic Neighbor Embedding (t-SNE) plots. The proportion of clusters per sample was normalized, and principal component analysis (PCA) was performed using Python library scikit-learn (42). Differentially expressed genes were obtained by “compared to entire dataset” mode based on each cluster. The p-values calculated using the differential expression analysis feature were adjusted for multiple comparisons using the Benjamini–Hochberg correction. Python package Scanpy v1.9.3 (43) (https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1382-0) was used to generate heatmaps. Genes used for cluster characterization were chosen based on expression profiles in zebrafish and mice (22, 25–27, 44–50) ( Supplementary Table S1 ).

ighm, sid1, trac, and trdc were not annotated in medaka Ensembl ASM223467v1. Reads mapped to the target region and their cell barcodes were obtained from the bam files generated by Cell Ranger using SAMtools (ver. 1.10). Target areas were defined as ighm (8:14941689-14941969) (34) sid1 (4:24323750-24324064), trac (17:12460290-12460457) (34) and trdc (17:12510953-12511133) (51). sid1 was retrieved from the genome database of Ensembl ASM223467v1 through TBLASTN, using the amino acid sequence of zebrafish sid1 (ENSDARG00000095540) as a query. Zebrafish and medaka sid1 are commonly located next to si:ch211-1a19.3. The frequency of each cell barcode in the generated bam file was counted using the awk command, and a csv file containing the cell barcodes and their frequencies was imported into Loupe Browser’s cluster mode for visualization on a t-SNE plot.

GO enrichment analysis was performed for 109 upregulated genes and 79 downregulated genes in the DP T (1) cluster against the DP T (2) cluster using Metascape v3.5.20240901 (https://metascape.org/gp/index.html#/main/step1) (52). Analysis was performed using annotation datasets of complete GO biological processes and a reference list of human genes. Twelve upregulated and 20 downregulated genes that were uncharacterized or lacked human orthologs were excluded from the analysis.

The datasets generated in the current study are available from the DNA Data Bank of Japan (DDBJ) database (accession numbers DRA013138 and DRA017334).

To establish an immunodeficient genetic model that lacks mature lymphocytes and use it for the analysis of cellular mechanisms of lymphocyte development, rag1 mutant medaka (rag1^–/– ) were generated by genome editing using the CRISPR/Cas9 system, targeting exon 4 of the rag1 gene ( Figure 1A ). A mutant line with an 8 bp deletion in rag1 was established, which was predicted to cause a frameshift and encode 570 appropriate amino acids, followed by three unrelated amino acids and a stop codon ( Figure 1A ). The deduced protein lacked most of the core region including the central domain, which recognizes the recombination signal sequence heptamer and RAG2 (53, 54) ( Figure 1B ), and was considered a null allele. The homozygous form of this frameshift mutant did not exhibit gross developmental defects.

First, we investigated the transcripts of rearranged lymphocyte antigen receptors in the kidney, which is equivalent to mammalian bone marrow as a site of hematopoiesis. Reverse transcription PCR (RT-PCR) analysis did not detect VDJ-recombined igm transcripts in B cells or tcrb transcripts in T cells in the rag1 mutant ( Figure 1C ; Supplementary Figure S1 ), as expected. Crossing the mutant with rag1:egfp transgenic medaka, which expresses EGFP under the control of the immature lymphocyte-specific rag1 promoter (34), revealed a reduction in EGFP-expressing cells in the larval thymus ( Supplementary Figure S2A ). Crossing the mutant with lck:egfp transgenic medaka ( Supplementary Information , Supplementary Figure S3 ), which expresses EGFP under the control of the immature T and NK cell-specific lck promoter, revealed a comparative area of EGFP-expressing cells in the larval thymus ( Figures 1D, E ). However, the fluorescence intensity of the thymus area in the rag1 mutant was significantly reduced compared to that in the wild-type ( Figure 1E ), suggesting a lower lck:EGFP⁺ cell density in the rag1 mutant. Adult rag1 mutant of rag1:egfp and lck:egfp backgrounds exhibited a reduction in EGFP signals in the thymus ( Supplementary Figure S2B ; Figure 1F ). Furthermore, reduction of thymus size in the rag1 mutant was confirmed by tissue sections ( Figures 1G, H ). These results suggest that the rag1 mutant has defects in T-cell development.

Next, B-cell development in the kidney of the rag1 mutant was assessed. Kidney sections from wild-type and rag1 mutants showed a comparable appearance of hematoxylin-rich blood cells ( Figure 2A ). Fluorescence-activated cell sorting (FACS) analysis of WKM cells was carried out for the rag1 mutant of the rag1:egfp background. Light scatter profiles revealed that the proportion of lymphocytes, which are mostly developing B cells (explained later in detail), against myelomonocytes was severely reduced in the rag1 mutant ( Figures 2B, C ). EGFP expression levels in the lymphocyte population were different between the wild type and rag1 mutant; whereas the rag1:EGFP^high population corresponding to relatively immature B cells was increased in the rag1 mutant, rag1:EGFP^low population corresponding to relatively mature B cells was decreased in the rag1 mutant ( Figures 2B, C , explained later in detail). In qPCR analysis of adult WKM cells, the expression of pan B-cell marker cd79a was significantly reduced. Contrary to a mild decrease in the expression of the early B-cell marker ebf1, expression of the relatively mature B-cell marker cd22l was severely reduced in the rag1 mutant, suggesting arrested B cell-development in the kidney. In addition, the expression of the T-cell marker cd8b was significantly decreased in rag1 mutants ( Figure 2D ). In contrast, the expression of the natural killer cell markers nitr17 and nkl.1 was mildly increased in rag1 mutants ( Figure 2D ). Among the other blood populations, the expression of the neutrophil marker mpx was comparable ( Figure 2D ). These results suggest that rag1 mutations affect T- and B- cell development.

To analyze the cellular and molecular mechanisms of T-cell development, thymocytes from two wild-type and two rag1 mutants at five months post-fertilization (mpf) were harvested separately. Cell sorting using light scatters was performed to obtain lymphocytes for single-cell transcriptome profiling ( Supplementary Figure S4 ). It was confirmed that the lymphocyte gate used for cell sorting covered the lck:EGFP⁺ cells ( Supplementary Figure S4 ). 6192 cells out of the 7390 wild-type cells, and 13239 cells out of 14528 rag1 mutant cells had appropriate UMI numbers and were accordingly used for clustering.

t-SNE visualization separated wild-type lymphocytes into 12 clusters, including 11 lck-expressing T/NK-lineage clusters and one cd79a-expressing B-cell cluster ( Figure 3A ). gata1-expressing erythrocytes, which seem to be contaminated by forward scatter (FSC)^low erythrocytes during cell sorting, were not many enough to constitute an independent cluster ( Figure 3A ). Excluding this B cluster and erythrocytes from wild-type and the corresponding rag1 mutant cells, PCA analysis of 5730 cells from two wild-type and totally 12661 cells from rag1 mutant cells based on the proportion of cells in each cluster was performed, confirming genotype-dependent separation ( Supplementary Figure S5A ).

Next, 5730 wild-type cells in T/NK-lineage clusters were separated into 10 clusters and annotated according to subsequent analyses of gene expression patterns and rag1-dependency as described below ( Figure 3B ). A heat map of gene expression levels was constructed using each cluster-specific representative genes (p<0.05) and genes known to be involved in T-cell development in mammals and zebrafish ( Figure 3C ; Supplementary Tables S2 , S3 ). lck and zap70 were constitutively expressed in all T/NK-lineage clusters ( Figure 3C ). Indication of wild-type T/NK-lineage clusters on t-SNE visualization of T/NK-lineage cells of both wild-type (5730 cells) and rag1 mutant (12661 cells) identified rag1-independent clusters (shared by wild-type and rag1 mutant) of early T-cell development and rag1-dependent clusters (absent in rag1 mutant) of late T-cell development requiring V(D)J recombination of tcr genes ( Figure 3D ). In contrast, re-clustering of T/NK cells from a mixture of wild-type and rag1 mutants did not generate any rag1 mutant-specific clusters (data not shown).

Cells of all five rag1-independent early clusters [Early T-cell progenitors (ETP), DN T, and cycling T (1–3)] were highly accumulated in the rag1 mutant ( Figure 3D ), suggesting developmental arrest at these stages. These clusters expressed rag1 and dntt ( Figures 3C, E ). In addition, two rag1-dependent clusters [DP T (1–2)] expressed rag1, cd4, and cd8, but not dntt ( Figures 3C–E ). Genes labeling proliferating cells (pcna, dut, rrm2) were commonly expressed in four clusters [cycling T (1–3), DP T (1)] ( Figures 3C, E ), suggesting cell proliferation before tcra or tcrb gene rearrangement. These results imply that cycling T (1–3) clusters are in the course of proliferation before tcrb gene rearrangement, and DP T (1–2) clusters expressing both cd4 and cd8 are in the more mature stage of T-cell development with ongoing or completed tcra gene rearrangement. These cd4⁺cd8⁺ double-positive (DP) T clusters commonly expressed rorc ( Figure 3C ).

The remaining rag1-independent cd4^–cd8^– double-negative (DN) clusters without proliferative gene expression (ETP and DN T) were assumed to be the most immature T cells. Although ETPs and genes specifically expressed in ETPs have been identified in zebrafish (27), these genes including cebpa, csf1rb, and fli1 were not expressed in these clusters in medaka (data not shown), suggesting that the timing of ETP immigration or gene expression patterns of ETPs is not conserved between zebrafish and medaka. Instead, the ETP cluster, not the DN T cluster, highly expressed some genes including gria1a and olfml2a ( Figures 3C, E ), suggesting that this ETP cluster contains recent thymic immigrants already expressing rag1.

In the t-SNE profile of wild-type T/NK-lineage cells, two clusters expressing cd8a and cd4-1 [DP T (1–2)] were located separately ( Figures 3B, E ). cd8b, whose gene product is assumed to form a heterodimer with the gene product of cd8a, and cd4-2, a teleost-specific paralogous gene of cd4-1, are widely expressed in the DP T (2) cluster. In contrast, cd8b- and cd4-2-expressing cells in the DP T (1) cluster were limited ( Figure 3C ; Supplementary Figure S6A ). Notably, neither cd4-1⁺cd8a^– nor cd4-1^–cd8a⁺ SP clusters were detected ( Figures 3C, E ; Supplementary Figure S6A ). Looking at the expression levels at the individual cell level, these genes were expressed in cells in the DP T clusters in an overlapping manner ( Figure 3F ).

To characterize the two DP T clusters, gene expression profiles were compared between them, and 109 and 79 genes were identified to be highly expressed in DP T (1) and DP T (2), respectively (p<0.05, log₂fc>2) ( Supplementary Table S4 ). GO analysis showed that genes related to pathways for cell division including “DNA replication-dependent chromatin assembly” (rbbp4, chaf1a, and asf1b) and “spindle elongation” (incenp, kif11, prc1, aurkb, kif23, and cdca8) were significantly enriched in DP T (1) cluster ( Supplementary Figures S6B, D, F ), suggesting that cells in DP T (1) cluster are more actively proliferating. On the other hand, genes related to pathways for lymphopoiesis including “T cell differentiation” and “T cell activation” (cd4-1, cd4-2, cd8a, cd8b, egr1, itk, ptger4, tespa1, and ccr9) were significantly enriched in DP T (2) cluster ( Supplementary Figures S6C, E, F ), suggesting that cells in DP T (2) cluster are more mature T cells. Notably, tespa1, which is required for positive selection in mammals (55), was enriched in the DP T (2) cluster. Both DP clusters express trac ( Figure 4A ), confirming that they are at the most mature stage of αβ T-cell development.

In contrast to trac expression in the two DP T clusters, trdc was expressed in a part of the three clusters [non-αβ T (1–3)] of the wild-type T/NK-lineage cells ( Figure 4A ). Similarly, sox13 and blk, which marked thymic γδ T cells in zebrafish (27), were also expressed in some cells in these three clusters ( Figures 4B , 3C ), suggesting that these clusters include γδ T-lineage cells. The expression of lymphocyte activation gene-3 (lag-3) encoding a CD4-related inhibitory receptor (56) was overlapped with the area of γδ T cells ( Figures 4B , 3C ). In-depth observation of these clusters showed that they also contained cells highly expressing zbtb16a, which encodes the natural killer T (NKT)-lineage-specific transcription factor Plzf (57, 58) ( Figures 4B , 3C ). zbtb16a expression did not merge with the area of γδ T cells. In addition, we found a rag1-independent cell mass expressing the natural killer (NK)-cell markers il2rb and nitr2 (25, 26, 59) ( Figures 4B , 3C ). These results suggest that, although there is no clear cluster distinction, the three non-αβ T clusters contain γδ T, NKT, and NK cells.

To assess the cellular and molecular basis of B-cell development, WKM cells from the three wild-type and two rag1 mutants at 2 mpf were harvested separately. Cell sorting using light scatters was carried out to obtain lymphocytes for single-cell transcriptome profiling using Cell Ranger (10x Genomics) ( Supplementary Figure S4 ). 22929 cells out of 25103 wild-type cells and 15472 cells out of 17604 rag1 mutant cells had appropriate UMI numbers, and accordingly were used for clustering.

t-SNE visualization separated wild-type lymphocytes into 24 clusters including 19 cd79a-expressing B-lineage clusters (19804 cells, 86.4%), 2 lck-expressing T/NK cell clusters (1017 cells, 4.4%), and 3 gata1-expressing erythrocyte clusters, which seemed to be contaminated by FSC^low erythrocytes during cell sorting ( Figure 5A ). Excluding these T/NK and erythroid clusters from wild-type and the corresponding rag1 mutant cells, PCA analysis of 19804 B-lineage cluster cells from three wild-type and 12989 cells from rag1 mutant cells based on the proportion of cells in each cluster was performed, confirming genotype-dependent separation ( Supplementary Figure S5B ).

Next, 19804 wild-type cells in B-lineage clusters were separated into 19 clusters and annotated according to subsequent analyses of gene expression patterns and rag1-dependency as described below ( Figure 5B ). A heat map of gene expression levels was constructed using each cluster-specific representative genes (p<0.05) and genes known to be involved in B-cell development in mammals and zebrafish ( Figure 5C ; Supplementary Tables S2 , S5 ). cd79a and pax5 were constitutively expressed in all the B-lineage clusters. Indication of wild-type B-lineage clusters on t-SNE visualization of B-lineage cells of both wild-type (19804 cells) and rag1 mutant (12989 cells) identified rag1-independent clusters (shared by wild-type and rag1 mutant) of early B cell development and rag1-dependent clusters (absent in rag1 mutant) of late B cell development requiring V(D)J recombination of immunoglobulin genes ( Figure 5D ). In contrast, re-clustering of B-lineage cells from the mixture of wild-type and rag1 mutants did not generate any rag1 mutant-specific clusters (data not shown).

Three clusters with high rag1 expression [cycling rag1 ^high B, rag1 ^high B (1–2)] and two clusters with low rag1 expression [rag1 ^low B (1–2)] were identified ( Figures 5C, E ). Whereas the rag2 expression pattern basically coincided with rag1, expression of DNA nucleotidylexotransferase (dntt) and sid1, an ortholog of zebrafish sid1 that is considered to be orthologous to mammalian VPREB1 (60), was observed only in the clusters with high rag1 expression ( Figures 5C, E ). In the rag1 mutant, cells merging with these rag1-expressing wild-type clusters were highly accumulated ( Figure 5D ), suggesting developmental arrest at this stage. Accumulation of the rag1:EGFP^high population in the rag1 mutant in the FACS analysis ( Figures 2B, C ) supports this notion.

Six Clusters [common lymphoid progenitors (CLP) (1–2), cycling rag1 ^– B (1–4)] were commonly present in the wild-type and rag1 mutants and did not express rag1 ( Figures 5C–E ). Among these clusters, genes labeling proliferating cells (mcm5, pcna, dut, rrm2, and mki67) were commonly expressed in 4 cycling rag1 ^– B clusters ( Figures 5C, E ), suggesting that these clusters represent the stage of cell proliferation before heavy chain gene rearrangement. csf1rb marking hematopoietic stem cells (27) was expressed in CLP (1-2) clusters, and cebpa marking hematopoietic stem cells as well as myeloid lineages (27) was expressed in part of the CLP (1) cluster ( Figures 5C, E ), suggesting that these cell populations are the earliest progenitors of the sorted lymphocytes.

Eight clusters [mpeg1 ⁺ B (1–8)] were rag1-dependent without rag1 expression ( Figures 5C–E ), suggesting that they were the most mature B-cell populations after completion of immunoglobulin gene rearrangement. In contrast to the expression of ebf1, myb, and si-dkey261h17.1 in early B cell-populations, these late clusters commonly express cd22l, dusp2, and mpeg1 ( Figures 5C, E ; Supplementary Figure S7 ). ebf1-expressing immature B cells, but not cd22l-high, relatively mature B cells, were present in the rag1 mutant, reflecting the qPCR results ( Figure 2D ).

Unlike zebrafish, in which clusters expressing igD and teleost-specific igT were separated (27, 61), igT has not been identified in medaka (62), and immunoglobulin m heavy chain (ighm)-expressing cells were scattered ( Supplementary Figure S7 ). In zebrafish, several genes including ccr2 and several genes including irf8 were co-expressed with igcz (igT) and igcd (igD), respectively (27). However, cells expressing one of the ccr2 paralogs and irf8 did not form independent clusters in medaka ( Supplementary Figure S7 ), indicating no sign of separation of the relatively mature B-cell populations.

To characterize T and NK cells in the kidney, lck-expressing T/NK clusters in wild-type kidney lymphocytes ( Figure 5A ) were focused and re-clustered. To this end, 1037 cells were separated into five clusters, all of which expressed lck and il7ra ( Figure 6A ; Supplementary Table S6 ). Indication of these wild-type T/NK-lineage clusters on t-SNE visualization of T/NK cells of both wild-type (1037 cells) and rag1 mutant (276 cells) identified one rag1-independent (NK cells) and four rag1-dependent clusters (T cells (1–3), cd79a⁺ T cells) ( Figure 6B ). The heatmap and t-SNE profiles of these clusters showed that rag1-dependent clusters commonly contained cells expressing cd4 paralogs (cd4-1 and cd4-2) or cd8 (cd8a and cd8b) ( Figures 6C, D ; Supplementary Figure S8 ), suggesting that they are T cells. Based on the expression levels at the individual cell level, cd4-1 and cd8a were separately expressed in the cells of T-cell clusters ( Figure 6E ), in contrast to overlapping expression in thymic DP clusters ( Figure 3F ), suggesting that they are mature CD4⁺ or CD8⁺ T cells.

One of these T-cell clusters (cd79a⁺ T cells) also expressed B-cell markers cd79a, pax5, and blnk ( Figures 6C, D ). Among the three cd79a^– T-cell clusters [T cells (1–3)], the T cells (3) cluster showed higher expression of genes including ccr9 and pcna ( Figure 6C ; Supplementary Table S6 ).

The rag1-independent cluster (NK cells) expressed il2rb, nitr2, and nitr17 ( Figures 6C, D ; Supplementary Figure S8 ), whose zebrafish orthologs are known markers of NK cells (25, 26, 59), suggesting that this population is composed of NK cells.

Another T lineage was identified in cd79a-expressing B-lineage cell clusters in the wild-type ( Figure 5A ). Cells expressing the T/NK-lineage markers lck and il7ra, but not cd4/cd8, were present in a portion of the CLP cluster and merged with cebpa-expressing cells ( Figure 7A ; Supplementary Figure S9A ; Figure 5E ). These cells also expressed notch1a, one of the orthologs of mammalian notch1 whose signaling is involved in T-cell fate decision (63, 64) ( Supplementary Figure S9A ).

Some genes including interferon-induced transmembrane protein 5 (ifitm5) and t-cell-specific immunoglobulin domain family 11 (tidf11), were specifically expressed in the identical population in kidney and T cells, except for DP clusters in the thymus ( Figure 7B ; Supplementary Figure S9B ). In the kidney, they were also expressed in mature T cells but not in NK cells ( Figure 7B ; Supplementary Figure S9B ). ifitm5 is one of the orthologs of mammalian ifitm family genes that encode transmembrane proteins involved in helper T-cell differentiation and immune response against cancer and viral infection (65, 66). tidf11 is one of 18 family member genes on chromosome 17 ( Supplementary Table S7 ). At least five other family member genes (tidf3, tidf7, tidf10, tidf12, and tidf16), for which the Ensembl gene ID was assigned, showed similar expression patterns to tidf11 (data not shown). As tidf family genes, except three members (tidf8, tidf9, tidf13) whose start codons have not been identified in this study, are predicted to encode transmembrane proteins with two immunoglobulin domains ( Figure 7C ), this family has been named the t-cell-specific immunoglobulin domain family. These results suggest that this CLP cluster contains cells committed to the T lineage and identified ifitm5 and tidf family genes already expressed at the T-lineage commitment stage and also in mature T cells.