Imbalance of stem-like and effector T cell states in children with early type 1 diabetes across conventional and regulatory subsets

Domain: neuroscience
Published: 2025-12-11
Editorial outcome: major_revision

Abstract
Type 1 diabetes (T1D) is an autoimmune disease caused by the loss of self-tolerance toward insulin-producing pancreatic β-cells. Its etiology remains incompletely understood but involves dysregulated T cell responses. Here, we perform single-cell transcriptomic analysis of peripheral blood T cells from children newly diagnosed with T1D, the same children after one year, and healthy donors. We observe that children with diabetes show diminished effector and cytotoxic programs and enhanced stemness-associated gene signature across diverse T cell subsets, especially at diagnosis. In parallel, we detect signs of impaired regulatory capacity in regulatory T cells and regulatory TR3-56 cells. These findings are supported by flow cytometry analysis of the same cohort and reanalysis of publicly available datasets. Overall, our results suggest that T1D is associated with impaired T cell effector differentiation and regulatory T cell dysfunction, both of which may contribute to immune imbalance and loss of self-tolerance.

Introduction
Type 1 diabetes mellitus (T1D) is an autoimmune disease characterized by the destruction of pancreatic β-cells, leading to insulin deficiency, hyperglycemia, and impaired metabolic homeostasis. As there is currently no cure for T1D, people with T1D rely on insulin replacement therapy for the rest of their lives. Unlike most autoimmune diseases, T1D typically manifests in childhood. Both genetic and environmental factors contribute to the development of T1D. The genetic factors include specific HLA haplotypes1, polymorphisms in the insulin gene (INS)2 and several immune-related genes (e.g., CTLA4, PTPN22, IFIH1, and CD226)3. The importance of the environmental factors is highlighted by an increase in T1D incidence in high-income countries4–6. It has been proposed that specific viral infections, such as enteroviruses, can trigger T1D7,8. On the contrary, the hygiene hypothesis proposes that the reduction of the infection burden in children may partly explain the rising incidence of T1D and other autoimmune diseases4,9–11. Despite some controversies in this area of research, it is well-established that a genetically and/or environmentally altered state of the immune system is a significant risk factor for T1D and may contribute directly to its onset. However, the etiology of T1D is still largely unknown.

T cells are extensively studied in the context of T1D because of the genetic evidence of HLA involvement, their clear role in the development of T1D in animal models12,13, and their presumed role in β-cell destruction via cell-mediated cytotoxicity14,15. Moreover, it has been proposed that T1D susceptibility may be induced by insufficient self-tolerance mediated by FOXP3+ regulatory T cells (Treg) and non-classical regulatory T cell populations, such as CD3+ CD56+ TR3-56 cells16. Despite numerous studies in the field17–23, the T1D-associated alterations in the T cell compartment are still not fully understood.

Traditionally, the immune cells in T1D are studied using targeted flow cytometry panels or bulk transcriptomics on whole blood, peripheral blood mononuclear cells (PBMC), and/or sorted leukocyte subsets24–28. Whereas flow cytometry is limited to the pre-selected protein markers, bulk transcriptomic analysis lacks flow cytometry’s single-cell resolution. These limitations can be overcome by multi-parameter approaches on the single-cell level, such as cytometry by time-of-flight or single-cell RNA sequencing (scRNAseq). Indeed, pioneering studies established these techniques as the method of choice to capture the differences between the immune cells in children with and without T1D21,26,29,30. However, these studies also revealed potential caveats, such as a trade-off between the sufficient size of the cohort to address reproducibility and the number of sequenced cells per donor to achieve sufficient resolution26,29,30. In parallel efforts, big consortia such as The Network for Pancreatic Organ Donors with Diabetes (nPOD)31 and Human Pancreas Analysis Program (HPAP)22,32 are generating large atlases of multimodal data from tissues of deceased donors, including scRNAseq, flow cytometry, immune repertoire profiling, and imaging data. Although these efforts provide invaluable datasets, including samples from solid organs that are hardly accessible, their disadvantage is their inability to recruit donors specifically based on features such as their age or time after diagnosis.

In this study, we compare the T cell compartments in children with and without T1D using single-cell transcriptomics, and we monitor longitudinal changes during the first year after diagnosis. T cells from newly diagnosed children display a bias toward naïve-like and stem-associated states, which is partially normalized after one year. In addition, Tregs show impaired functional differentiation with a gene signature of dysfunction. Together, these findings highlight a disrupted balance between effector and regulatory programs in T cells, providing new insights into mechanisms that may underlie loss of self-tolerance in T1D.

Results
Generation of a T cell atlas of children with T1D and healthy donors
We performed a prospective study on a cohort of 43 pediatric donors with T1D or without T1D (henceforth referred to as healthy donors) using scRNAseq to identify potential gene expression patterns and heterogeneity in blood T cells associated with T1D. Specifically, we searched for differences in the T cell compartment between the healthy donors and children with T1D, and for changes occurring during the first year after the diagnosis of T1D. We based our analysis on the scRNAseq, including T cell receptor sequencing (VDJ profiling) and flow cytometry analysis of T cells from the peripheral blood of our Czech cohort, complemented by the analysis of openly accessible data from other cohorts.

Blood samples from 30 children with T1D were taken 4–24 (median 9) days after their diagnosis (time T0) and again at 359–423 (median 377) days after diagnosis (time T1). For comparison, we collected blood from 13 healthy donors (Fig. 1A). The T1D and healthy donor cohorts were comparable in age distribution, as confirmed by descriptive statistics (Supplementary Table 1, Supplementary Fig. 1A) and standardized mean difference (SMD) analysis (Supplementary Fig. 1B). However, the sex distribution differed, with a higher proportion of males in the T1D cohort (67%) compared to the healthy donors (33%) (Supplementary Table 2). We characterized the donors by their age, sex, partial clinical remission status at T1, T1D-associated autoantibody levels, HLA haplotypes (Fig. 1B), as well as blood concentrations of glycated hemoglobin, random and fasting C-peptide, and leukocyte levels (Supplementary Fig. 1C–E). At T0, children with T1D showed a higher proportion of lymphocytes (mean 48.1%, SD 12.9%) and a lower proportion of neutrophils (mean 38.9%, SD 13.1%) compared to healthy donors (lymphocytes: mean 40.4%, SD 11.0%; neutrophils: mean 46.5%, SD 12.6%) and to their own levels at T1 (lymphocytes: mean 37.9%, SD 11.2%; neutrophils: mean 49.2%, SD 11.8%) (Supplementary Fig. 1E).

First, we performed an initial scRNAseq analysis of CD4+ and CD8+ lymphocytes from 16 samples, which revealed that the majority of T cells in our child donors had a naïve CD45RA+ CD45RO- phenotype (Supplementary Fig. 2A–D), determined by CD45 splicing-sensitive analysis tool IDEIS33. To enable a more detailed characterization of antigen-experienced (non-naïve) T cells in the final scRNA-seq analysis, we enriched for these cells by sorting at a fixed 1:5 ratio of naïve to antigen-experienced cells for both CD4+ and CD8α+ lymphocytes (Supplementary Fig. 2E–H). Details on sample counts, donor numbers, and enrichment strategies for both initial and final experiments are summarized in Supplementary Table 3.

In total, we processed 176 samples, with an average of 992 cells per sample. On average, each cell contained 2689 non-redundant transcripts and 1464 detected genes (Supplementary Table 4). We processed the data using our standard pipeline34 including the removal of dead cells and doublets (Supplementary Fig. 3A, B). Because the library preparation was performed in six batches, we used integration by STACAS to remove the batch effect35,36. Integration successfully removed the unwanted variability between experiments, but importantly, it preserved biologically relevant differences, such as the difference between initial and final experiments, which contained non-enriched versus enriched cells (Supplementary Fig. 4A–H). Comparison of data from initial and final experiments showed that, as expected, our approach reduced the frequency of naïve cells from ~70–75% to 20–25% (Supplementary Fig. 4I–J). All subsequent analyses were restricted to the final experiments, i.e., enriched samples, to ensure consistency.

In the next step, we generated CD8+ and CD4+ lymphocyte atlases using the processed scRNAseq data. The CD8+ lymphocyte compartment consisted of three clusters: CD8αlow NK cells, unconventional CD8+ T cells represented by semi-invariant MAIT cells and γδ T cells, as well as conventional CD8+ αβ T cells (Supplementary Fig. 5A–D). Based on the typical markers, we identified that conventional CD8+ T cells contain naïve, memory (Tmem), unconventional KLRC2+ memory, also previously described as KILR-like cells37,38, which we abbreviate as Tmk (T cell memory KLRC2+), terminally differentiated CD45RA+ (Temra), and proliferating subsets (Fig. 1C, D and Supplementary Fig. 5E). Temra and to a lesser extent, proliferating and Tmem cells, but not Tmk cells, were clonally expanded (Fig. 1E). These subsets could be further subdivided into 41 smaller subsets (Supplementary Fig. 5F). The identity of cell subsets was validated using a previously published RNAseq dataset of sorted populations of CD8+ T cells39 (Supplementary Fig. 5G, H).

The CD4+ lymphocyte compartment consisted of a subset of unconventional T cells including iNKT cells expressing their signature TCRα (TRAV10, TRAJ18; CDR3: CVVSDRGSTLGRLYF) (Supplementary Fig. 6A, B) and conventional CD4+ T cells, which could be further divided to naïve, central memory (Tcm), Treg, Th1/Th17, Th2, Temra, CD4+ T cells with a high expression of interferon signaling associated genes (ISAGhi)40, proliferating cells, and CD4+ T cells enriched for genes involved in signal transduction pathways (Tsig) (Fig. 1F, G and Supplementary Fig. 6C–F). Temra cells were the only CD4+ T cell subset showing a substantial level of clonal expansion (Fig. 1H). These subsets could be further split into 38 smaller subsets (Supplementary Fig. 6G). The identity of cell subsets was validated using a previously published RNAseq dataset39 (Supplementary Fig. 6H, I).

We compared the proportions of particular subsets using both frequentist (Supplementary Fig. 7A–D) and Bayesian (Supplementary Fig. 7E, F) statistics. These analyses did not reveal substantial differences between children with T1D and healthy donors, with the exception of CD8+ Temra cells, which were enriched in the healthy donors in comparison to the children with T1D (Supplementary Fig. 7A, E). The Bayesian statistics revealed apparent differences in ISAGhi and Tsig subsets (Supplementary Fig. 7F), which, however, originated from a limited number of outlier participants (Supplementary Fig. 7C). Within the T1D group, the frequency of MAIT cells among CD8+ T cells and the frequency of Th2 cells among CD4+ T cells correlated with fasting C-peptide level, which we used as a proxy for residual β-cell activity (Supplementary Fig. 7G, H). However, because of the prior enrichment for antigen-experienced T cells, the relative abundance of the particular subsets in the scRNAseq data does not correspond to the real frequency in the original sample.

The onset of diabetes is associated with low effector T cell signatures
Next, we compared gene expression between healthy donors and children with T1D at T0 and T1, separately in CD4+ and CD8+ T cells (Fig. 2A, B). We visualized the consistency of gene expression changes by plotting log-fold changes in T1D T0 vs. healthy donors against T1D T1 vs. healthy donors for each gene, showing the key differences at both timepoints (Fig. 2A).

At both timepoints, we found an enrichment of effector and cytotoxic genes (GZMA, GZMB, GZMH, PRF1, GNLY, IFNG, TNF, CCL5, CCL4) in healthy donors compared to T1D. In contrast, genes associated with naïve, quiescent, and/or stem-like T cells (TCF7, LEF1, BACH2, IL7R, CXCR4, ZFP36L2) were upregulated in children with T1D. Hereafter, we refer to the latter set of genes as a stemness-associated gene signature, reflecting transcriptional features shared by naïve, quiescent, and stem memory T cells, rather than being strictly limited to canonical naïve cells. This bias toward a stemness-associated signature was most pronounced at T0 (Fig. 2B).

T cells from T1D T0 samples also showed higher expression of several regulators of intracellular activation signaling pathways (DUSP1, TSC22D3, TNFAIP3), which might be associated with the quiescent state (Fig. 2B). In contrast, samples from healthy donors showed a higher expression of genes associated with calcium signaling and NFAT activation (ORAI1, NFATC2, NFATC3), potentially supporting a more activated, effector-like state.

To further support these observations, we performed gene set enrichment analysis (GSEA) using signatures of genes upregulated or downregulated in naïve versus effector or effector memory T cells from previously published datasets (GSE965041, GSE1105742). This analysis confirmed the enrichment of the naïve (stemness-associated) signature in T cells from children with T1D and an effector signature in T cells from healthy donors (Fig. 2C).

To validate these findings at the sample level, we performed differential expression analysis using DESeq2 on pseudobulk expression profiles aggregated per donor. This analysis confirmed a strong correlation of log-fold changes with the original per-cell analysis (Supplementary Fig. 8A, B). Key genes representing effector-, stemness-, and NFAT-associated signatures were identified by both types of analysis (Supplementary Fig. 8C–F).

In addition, for each donor, we quantified the similarity of their gene expression to reference contrasts of naïve versus memory cells from a previously published dataset (GSE17961339) using decoupleR43, which, unlike GSEA, explicitly accounts for the sign and weight of network interactions. This comparison revealed enrichment of naïve and stem memory signatures in T1D samples at T0, and enrichment of Temra signatures in healthy donors (Fig. 2D). These patterns were observed across different age and sex groups, confirming that these differences are not driven by age or sex (Supplementary Fig. 9A, B).

In addition to the bulk CD4+ and CD8+ T cell compartments, we also analyzed differential expression within individual T cell subsets (Supplementary Figs. 10A, B, 11A, B). This revealed that while most genes followed similar patterns across subsets, certain differentially expressed genes (DEG) (e.g., MYOM2, IFI27, CD38, GZMK, FOS, FOSB) showed distinct, subset-specific regulation. To quantify these overlaps, we assessed shared top DEGs across subsets, showing high similarity among naïve and memory populations (Supplementary Figs. 10B, 11B).

To check whether the downregulation of effector genes and upregulation of naïve genes in children with T1D is associated with the disease severity, we carried out the comparison of differential expression between children with T1D at T0 who suffered from ketoacidosis and those who did not, and between T1D who fulfilled or not the criterion of partial clinical remission (insulin dose adjusted HbA1c < 944,45) (Supplementary Fig. 12A–D). We did not observe a pronounced bias towards stem-like or effector-like gene expression in either of these subgroups (Supplementary Fig. 13A, B).

For the interpretation of the previous results, it should be noted that the ratio of naïve and antigen-experienced T cells was normalized for the scRNAseq analysis. For this reason, we analyzed the samples by flow cytometry to reveal that both antigen-experienced CD8+ and CD4+ T cells were less abundant in children with T1D than in healthy donors at T0, but not at T1 (Fig. 3A). Accordingly, the expression of a key cytotoxic protein Granzyme B (GZMB) is lower in children with T1D in CD8+ and CD4+ T cells than in the healthy donors at both timepoints (Fig. 3A). The frequencies of populations analyzed by scRNAseq before enrichment correlated with the frequencies obtained by flow cytometry, indicating the consistency of the two methods (Supplementary Fig. 14A).

Overall, the data suggested that the T cell compartment has a more stem-like and less effector-like state in children with T1D than in healthy donors. This stem-like state is partially normalized within one year after the diagnosis and the associated introduction of the insulin therapy.

To extend our findings to unrelated T1D cohorts, we addressed the question of the stem-like vs. effector signature in people with T1D in previously published data. Using the HPAP database32, we re-analyzed scRNAseq data obtained from T cells from the spleens of 53 children and adult donors with or without T1D. In this cohort, the effector signature genes were enriched in the control samples, and the stemness-associated (naïve-like) signature genes were enriched in donors with T1D (Fig. 3B). In addition, we re-analyzed flow cytometry data from the HPAP and observed a reduction of effector CD8+ and CD4+ T cells in the blood of donors with T1D (Fig. 3C and Supplementary Fig. 14B), consistent with the findings in our pediatric cohort.

Next, we collected transcriptomic data from previously published studies26,29,30,46,47 (Supplementary Table 5) and re-analyzed it using the same pipeline. Genes consistently upregulated across T1D datasets included those connected to the naïve, stem-like, and quiescent T cell states such as LEF148, FOXP149, BACH250, NELL251, IKZF152 and a negative regulator PI3 kinase signaling INPP4B53 (Fig. 3D). In contrast, genes associated with cytotoxicity and effector functions, such as NKG7, GNLY, CCL4, CCL5, CX3CR1, TNFRSF4 (OX40), TNF, and TNFRSF9 (4-1BB) were upregulated in healthy donors. Altogether, our results showing a T cell bias towards a naïve-like, stem-associated, and quiescent state in T1D are supported by the majority of the published datasets.

BTN3A2 expression depends on HLA haplotypes
It was shown that BTN3A2 and to a lower extent, BTN3A3, encoding butyrophilin family members, are upregulated in T cells of children with T1D compared to healthy donors26. However, we observed the opposite in our dataset, i.e., the higher expression of BTN3A2 and BTN3A3 in healthy donors than in children with T1D (Fig. 2A, B). Since BTN3A2 is a part of the extended MHC region54, we hypothesized that its expression may depend on specific HLA haplotypes rather than T1D status.

In our cohort, most children with T1D carried known susceptibility alleles, specifically HLA-DR3-DQ2 (haplotype HLA-DRB1 03:01–DQA1 05:01–DQB1 02:01) and/or HLA-DR4-DQ8 (HLA-DRB1 04:01–DQA1 03:01–DQB1 03:02)1, while healthy donors reflected the general HLA distribution in the Czech population55 (Supplementary Fig. 15A). We used our data and publicly available data to estimate the role of the T1D status, particular haplotype, and study on BTN3A2 expression. We observed that the study and HLA haplotypes, but not the T1D status, are the factors significantly influencing BTN3A2 expression (Supplementary Figs. 15B, 16). This indicates that the selection of healthy donors (e.g., unmatched as in our study or partially matched as in ref. 26) might have a great impact on the results concerning BTN3A2 expression.

Naïve-like regulatory T cells in children with T1D
It has been proposed that the dysfunction of Treg cells might contribute to the loss of self-tolerance and the development of autoimmune diabetes in humans and animal models56,57. We identified four Treg subclusters labeled Treg1 to Treg4 (Fig. 4A). The gene expression analysis revealed that Treg1 represents Tregs with the most naïve/stem-like phenotype, with high expression of stemness-associated genes such as CCR7, TCF7, SELL and relatively low expression of Treg signature genes such as FOXP3, CTLA4, and IL2RA (Fig. 4B–D). However, the expression of these signature molecules in Treg1 cells was higher than in naïve T cells suggesting that they are bona fide Tregs (Fig. 4D). In contrast, the Treg4 subset represents Tregs with the most activated phenotype, which predispose them for potent regulatory functions58,59.

The naïve-like Treg subset was more abundant among Tregs in the children with T1D than in healthy donors (Fig. 4E, F and Supplementary Fig. 17A) at T0 and T1. Treg1 cells expressed high levels of CD226 and low levels of TIGIT (Fig. 4C), which both bind to common ligands CD112 and CD155, but have the opposite functions. CD226 is associated with Treg dysfunction, whereas TIGIT in Tregs is linked to the self-tolerance in the context of T1D60–63.

To further understand the relationship among Treg subsets, we performed pseudotime analysis64, which revealed a trajectory consistent with progressive Treg activation and differentiation (Fig. 4G). Pseudotime strongly correlated with a Treg activation score (Fig. 4H), and median pseudotime values indicated that Tregs were more activated in healthy donors compared to children with T1D (Fig. 4I). To formally test for differences in trajectory progression, we applied the condiments package65, which confirmed significant differences in pseudotime dynamics between groups (Fig. 4J–L).

In the next step, we addressed whether Tregs from T1D resemble functionally impaired Tregs from conditions in which their dysfunction is well established and characterized. First, we compared T1D Tregs from our cohort to Tregs from patients with a disease called Immune dysregulation, polyendocrinopathy, enteropathy X-linked syndrome (IPEX), which is caused by loss-of-function mutations in FOXP3. Tregs from children with T1D upregulated genes specific for defective Treg-like cells observed in IPEX patients66,67 (Fig. 5A). We obtained the same results when we compared Tregs from children with T1D to Treg-like cells in Foxp3-/- mice68,69 and bona fide human Tregs after their key transcriptional regulators FOXP3, HIVEP2, IKZF2 (Helios), or SATB1 were knocked-out70 (Supplementary Fig. 17B, C). Moreover, we observed an inverse correlation between IL4R expression in Tregs and residual insulin production measured as fasting C-peptide levels (Fig. 5B), which points to a previously proposed inhibitory role of IL-4 signaling in Treg differentiation71,72. In contrast, the expression of IL10RA in Tregs correlated with fasting C-peptide levels, which highlights a possible role of IL-10 signaling for Treg-mediated tolerance in T1D (Fig. 5B). Fasting C-peptide levels also correlated with the frequency of Treg1 cells (negative correlation) and the frequency of Treg4 cells (positive correlation) (Fig. 5C), suggesting that the Treg maturation status might be involved in the disease severity. Overall, this evidence suggests that Tregs in children with T1D might be less mature and/or less functional than their counterparts in healthy donors.

To address our findings in independent cohorts, we analyzed data from three publicly available scRNAseq datasets studying peripheral immune cells of donors with T1D and healthy donors: a recent publication30, HPAP dataset32, and a reference dataset published by Parse Biosciences (henceforth ParseBio). In each dataset, we identified Treg cells based on canonical markers, then extracted and reclustered these cells to define four clusters (A–D), which roughly corresponded to our Treg1–4 clusters based on their gene expression profiles (Fig. 5D). Similar to our dataset, cluster A represented the most naïve and D represented the most mature Treg state (Fig. 5E, F). Cluster A, which showed high expression of naïve genes and IL4R, was overrepresented in Tregs from donors with T1D in comparison to healthy donors, whereas cluster D, which showed high expression of Treg effector genes, was overrepresented in healthy donors in all three datasets (Fig. 5G, H).

In the next step, we examined FOXP3+ Tregs in our cohort using flow cytometry (Supplementary Fig. 18A, B) and analogous data from the HPAP database (Supplementary Fig. 18C, D). In both cases, there were slightly more FOXP3+ cells among CD4+ T cells in children with T1D than in healthy donors. However, we did not see consistent differences in the expression of the naïve or effector signature genes between children with T1D and healthy donors. A possible reason is the low expression of FOXP3 in the naïve Treg subset (Fig. 4C, D), which makes it difficult to identify these cells by the conventional flow cytometry panels, which highly depends on FOXP3 detection.

Low frequencies of effector-phenotype unconventional T cells in children with T1D
It has been proposed that a loss of unconventional CD3+ CD56+ T cells (named TR3-56) with regulatory properties might contribute to the development of T1D16. In our cohort, we observed that children with T1D at T0, but not at T1, exhibited slightly lower expression of the TR3-56 signature than healthy donors, both in CD4+ and CD8+ T cells (Fig. 6A). We observed similar results in two unrelated cohorts with publicly available data (Supplementary Fig. 19A).

In the next step, we aimed to identify TR3-56 cells in our scRNAseq data, starting with the CD8+ T cells atlas, consisting of naïve and antigen-experienced conventional T, NK, γδ T cells, and MAIT cells (Fig. 6B, C and Supplementary Fig. 19B). Among these populations, TR3-56 cells were predominantly annotated within the MAIT and γδ T cell subsets (Fig. 6B, C and Supplementary Fig. 19C). Accordingly, these clusters contained the most cells coexpressing CD3 and CD56 (Supplementary Fig. 19D). Furthermore, TR3-56 cells in these clusters exhibited the highest annotation scores (Fig. 6D), supporting their overlap with MAIT and γδ T cells.

To identify unique markers of TR3-56 cells independent of their parental subsets, we performed differential gene expression analysis comparing TR3-56 cells and non-TR3-56 cells in the CD8+ lymphocyte clusters (Fig. 6E). This analysis revealed that TR3-56 cells upregulate transcription factors associated with cytotoxicity and T cell activation (RUNX3, BHLHE40, NR4A3), as well as genes associated with homing to tissues (CCR6, ITGB1). Finally, we assessed the regulatory potential of TR3-56 cells by comparing their gene expression to a canonical Treg signature, and found that TR3-56 cells exhibited the highest similarity among all CD8+ lymphocyte subsets. (Fig. 6F).

Re-clustering of the unconventional T cells identified two clusters of γδ T cells and four clusters of MAIT cells (Fig. 6G, H). We observed that the γδ T1 cluster was enriched in healthy donors in comparison to children with T1D at T0 and T1 (Fig. 6I). This cluster corresponded to effector-phenotype γδ T cells (Fig. 6H), suggesting that the state of γδ T cells in T1D essentially phenocopies conventional T cells (Fig. 2A–D). In contrast, γδ T2 cells, enriched in some children with T1D (Fig. 6I), expressed lower levels of effector genes (Fig. 6H). Accordingly, the total unconventional CD8+ T cells from children with T1D express lower levels of cytotoxic effector genes such as GZMB, GZMH, GNLY, and CCL5 (Supplementary Fig. 19E).

Using flow cytometry, we observed that CD8low T cells are enriched for γδ TCR+ and CD56+ cells (Fig. 6J and Supplementary Fig. 20A). Thus, gating for CD8low cells can be used as a proxy for CD8+ γδ T cells or CD8+ CD3+ CD56+ T cells. In our flow cytometry analysis, we found that percentages of γδ T cells and CD8low T cells among CD3+ cells were low in children with T1D at T0, which was partially reverted at T1 (Fig. 6K and Supplementary Fig. 20A). Moreover, a higher percentage of CD8low T cells expressed GZMB in healthy donors than in children with T1D, essentially confirming the scRNAseq data (Fig. 6K). A higher percentage of CD8low T cells among CD3+ T cells from healthy donors compared to donors with T1D was observed also in the HPAP flow cytometry data (Supplementary Fig. 20B, C).

In the CD4+ T cell compartment, the Temra cluster contained most TR3-56 phenotype cells (Fig. 7A–C) and had the gene expression profile most similar to TR3-56 cells (Fig. 7D, E). Re-clustering of the Temra cluster revealed four subclusters (Fig. 7F). The Temra2 cluster, corresponding to the cells with the strongest cytotoxic phenotype and the highest expression of GZMB (Fig. 7G), was more abundant in healthy donors than in children with T1D at T0 (Fig. 7H).

The analysis of T cell receptor repertoire in children with T1D
Alongside the single-cell transcriptomics, we analyzed T cell receptor sequences in T cells from the children with T1D and healthy donors. We observed only minimal clonal overlap between individual donors regardless of their disease status (Fig. 8A). However, there was a significant overlap of TCRα and TCRβ sequences in the same children with T1D at T0 and T1 (Fig. 8A), indicating the persistence of these clones in the blood during the one year period. We did not identify any T cell clones previously associated with T1D. The usage of variable, diversity, and joint segments did not significantly differ between the children with T1D and healthy donors (Supplementary Figs. 21, 22).

It has been reported previously that the CDR3 segments of TCRβ chains in people with T1D are generally shorter than in healthy individuals73. However, we did not observe any differences in the CDR3β length in CD4+ or CD8+ T cells in our cohort or in TCR repertoire profiling data from the HPAP database (Fig. 8B–D).

Finally, we analyzed the biochemical properties of the CDR3β TCR sequences in our cohort and in the HPAP cohort. We observed that the CDR3β sequences from people with T1D showed significantly higher Boman index and H Moment and a slightly lower Hydrophobicity index (Fig. 8E). This effect was apparent also on the average values per donor (Fig. 8F, G), but with a relatively low significance. A large cohort is required for eventual confirmation of this observation.

Methods
Collection of clinical metadata
The samples were collected between February 2021 and February 2023 at a setting of a tertiary center for pediatric diabetes (Motol University Hospital, Prague, Czechia). After the consent of parents/caregivers and the study participants was granted, two study samples were obtained from each participant. The first sample was obtained at a median of 9 (4–24) days after the clinical diagnosis of T1D, the second sample at a median of 377 (359–423) days after diagnosis. At both time points, fasting C-peptide, HbA1c, and the standard blood count with the differential count were performed. At T0, the measurements of the T1D-specific autoantibodies (a-GAD, IAA, a-IA2, a-ZnT8) were performed. Data on diabetic ketoacidosis status at first contact (pH below 7.3 and/or bicarbonate < 15) were obtained from individual medical records. To assess partial clinical remission at T1, we employed the insulin dose adjusted HbA1c (IDAA1c) as a marker. IDAA1c was calculated as HbA1c (%) plus [4 times insulin dose (units per kilogram per 24 h)]. IDAA1c < 9 was considered as the presence of partial clinical remission112.

The samples were obtained from 30 children with T1D at T0 and from 29 of these children at T1. One T1D participant withdrew from the study before the one-year follow-up visit. Their sample was included in all analysis except for paired statistical testing. A cohort of 13 healthy donors was also enrolled. Healthy donors were frequency-matched to the T1D cohort by age. Healthy donors were not HLA-matched to participants with T1D. The parents/caregivers and the participants signed a written consent.

Processing of blood samples
Three to ten mL of peripheral blood were collected into EDTA-coated tubes, kept on ice, and transferred from Motol University Hospital to the Institute of Molecular Genetics within two hours. PBMCs were separated using Ficoll-Paque (GE Healthcare) and immediately frozen following a cryopreservation protocol (10x Genomics). Briefly, 4 ml of Ficoll-Paque was overlaid with blood and centrifuged at 400×g for 30 min at room temperature (brake set to one). The mononuclear cell layer was washed in PBS and resuspended in RPMI medium containing 40% FBS. After adding an equal volume of freezing medium (RPMI, 40% FBS, and 30% DMSO), two to five aliquots of PBMC were frozen and stored in liquid nitrogen. PBMCs were gently thawed by slow, sequential dilution in RPMI medium containing 10% FBS. To minimize cell stress, wide-bore tips were used during the thawing process.

ScRNAseq experimental design
ScRNAseq was performed in two different experiments, hereafter referred to as the initial experiment and the final experiment.

Initial scRNAseq experiment
The initial experiment was performed in two batches: batch 1 – CD8+ T cells from 6 children with T1D and 2 healthy donors in one well, batch 2 – CD4+ T cells from 6 children with T1D and 3 healthy donors, CD8+ T cells from 1 healthy donor in one well, and CD4+ and CD8+ cells from 6 children with T1D and 3 healthy donors in two wells. Cells were incubated on ice for 5 minutes with Human TrueStain FcX (BioLegend #422301) and for additional 30 min with anti-human CD4 and CD8 antibodies (CD8 APC, LT8, Exbio #1A-817-T100; CD4 AF700, MEM-241, Exbio #A7-539-T100) and one of the hashtag antibodies (TotalSeq™-C0251 anti-human Hashtag 1 Antibody, LNH-94 2M2, BioLegend #394661; TotalSeq™-C0252 anti-human Hashtag 2 Antibody, LNH-94 2M2, BioLegend #394663; TotalSeq™-C0253 anti-human Hashtag 3 Antibody, LNH-94 2M2, BioLegend #394665; TotalSeq™-C0254 anti-human Hashtag 4 Antibody, LNH-94 2M2, BioLegend #394667; TotalSeq™-C0255 anti-human Hashtag 5 Antibody, LNH-94 2M2, BioLegend #394669; TotalSeq™-C0256 anti-human Hashtag 6 Antibody, LNH-94 2M2, BioLegend #394671; TotalSeq™-C0257 anti-human Hashtag 7 Antibody, LNH-94 2M2, BioLegend #394673; TotalSeq™-C0258 anti-human Hashtag 8 Antibody, LNH-94 2M2, BioLegend #394675; TotalSeq™-C0259 anti-human Hashtag 9 Antibody, LNH-94 2M2, BioLegend #394677; TotalSeq™-C0260 anti-human Hashtag 10 Antibody, LNH-94 2M2, BioLegend #394679) and with Hoechst 33258 for viability right before sorting. From each sample, 10,000 CD4+ cells and 10,000 CD8+ cells were sorted. All samples were collected to the same collection tube, washed with PBS/0.05% BSA and counted using the TC20 Automated Cell Counter (#1450102, Bio-Rad). The viability of the cells before loading was higher than 85%.

Final scRNAseq experiment
The final experiment was performed in four batches: batch 1 – CD4+ and CD8+ T cells from 8 children with T1D T0, 8 children with T1D T1, 4 healthy donors in 4 wells, batch 2 - CD4+ and CD8+ T cells from 8 children with T1D T0, 8 children with T1D T1, 4 healthy donors in 4 wells, batch 3 - 8 children with T1D T0, 7 children with T1D T1, 3 healthy donors in 4 wells, batch 4 - 6 children with T1D T0, 6 children with T1D T1, 4 controls in 4 wells. Cells from the same donor from time 0 and time 1 were always processed in the same well to minimize batch effect. Cells were incubated on ice for 5 minutes with Human TrueStain FcX (BioLegend #422301) and for additional 30 minutes with anti-human CD4, CD8, CD45RA and CCR7 antibodies (CD4 APC, MEM-241, Exbio #1A-359-T100; CD8 PE, MEM-31, Exbio #1P-207-T025; CD45RA FITC, MEM-56, Exbio #1F-223-T100; CCR7 PeCy7, G043H7, BioLegend #353226) and one of the hashtag antibodies (the same as in the initial experiment) and with Hoechst 33258 for viability right before sorting. Non-naïve cells were enriched in the final sample as follows: from each sample 1000 naïve (CCR7+ CD45RA+) and 5000 antigen-experienced cells were sorted into two tubes, one for CD4+ T cells and one for CD8+ T cells. All samples were collected to the same two collection tubes, washed with PBS/0.05% BSA and counted using the TC20 Automated Cell Counter (#1450102, Bio-Rad). The viability of the cells before loading was higher than 85%. One control sample was removed from the analysis because of the bad quality of the data at the sort.

Cells from both cohorts were loaded onto a 10x Chromium machine (10x Genomics) aiming for a yield of 1500 cells per sample. cDNA libraries were prepared using the Feature Barcode technology for Cell Surface Protein protocol (#CG000186 RevD) with the Chromium Single Cell 5’ Library & Gel Bead and Chromium Single Cell 5’ Feature Barcode Library kits (10x Genomics, #PN-1000014, #PN-1000020, #PN-1000080, #PN-1000009, #PN-1000084) according to the manufacturer’s instructions. Sequencing was performed on the NovaSeq 6000 platform (Illumina), yielding an average of 45,745 reads per cell in the gene expression libraries and 3712 reads per cell in the V(D)J libraries (Supplementary Table 4).

Quality control, normalization, and integration of scRNAseq data
The human reference used to map sequenced reads was taken from Ensembl version 102113 and prepared using 10x Cell Ranger 5.0.1 Software (mkref tool). The count matrices were generated by 10x Cell Ranger 5.0.1 Software (count tool) in either R2-only or paired-end mode. Afterwards, they were pre-processed using the Seurat 4.3.1 package114 on R 4.2.1. Any cell that was not marked by any expected combination of hashtags was removed. All cells with more than 10% of genes mapping to mitochondrial genes, those expressing less than 200 genes and those that were marked as doubles according to the V(D)J (more than 2 TCRα or more than 1 productive TCRβ sequence present) were excluded during initial object creation (code available in GitHub folder code01_raw_to_init).

Mitochondrial genes, ribosomal genes, genes encoding TCR variable segments (any gene symbol containing the TR[AB][VDJ] substring) and genes present in less than 3 cells were removed. Each data set was then normalized (default method and scale factor=1×104), scaled, subjected to dimensional reduction (PCA with 20 principal components followed by UMAP) and Louvain clustering. PTPRC counts were generated by the IDEIS tool and normalized by the centered-log ratio method33.

Preprocessed datasets were merged, normalized, scaled and used for dimensional reduction (PCA with 12 principal components and UMAP) and Louvain clustering. In a second, more stringent quality control step performed on the merged dataset (code available in GitHub folder code02_analysis_of_init, parts 02 and 03), clusters of dead or low-quality cells and individual cells with fewer than 500 detected genes were excluded. Integration of datasets from different batches was performed using STACAS (v 2.1.3)115. In a subclustering step for CD4+ Temra cells, we revealed a contaminating cluster of NK cells (n=80), which only appeared after reclustering. We kept those cells in the final dataset, but removed them from the analysis in Fig. 7.

Analysis and visualization of scRNAseq data
Transcriptional regulatory interactions were estimated using the decoupleR package (v2.10.0) using the CollecTRI database116. GSEA pathways were processed using gene sets from datasets GSE965041 and GSE1105742 with the R package fgsea (v1.20.0)117. Heatmaps were created with the package pheatmap (v1.0.12). Sankey plots for Supplementary Fig. S3 and S4 were created using RAWGraphs (v2.0)118.

DEGs were identified using the two-sided Wilcoxon rank-sum test (Seurat FindMarkers), with p-values adjusted for multiple testing using the Bonferroni correction. To address concerns about pseudoreplication in cell-level differential expression analysis, we performed a sample-level pseudobulk analysis. For each annotated T cell subset, raw gene expression counts were aggregated per donor, generating one pseudobulk profile per donor per subset. Differential gene expression and log fold changes were then calculated using the DESeq2 package (v1.44.0), enabling rigorous comparison between healthy donors and children with T1D at different timepoints. In parallel, we applied the decoupleR package (v2.10.0) to quantify the similarity of each donor’s expression profile to reference gene expression contrasts between naïve and memory T cells (GSE179613). Pseudotime analysis was performed using the slingshot package (v2.12.0) to infer a lineage trajectory among Treg subsets. The analysis was based on UMAP embeddings and cluster annotations. To test for differences in pseudotime progression between groups, we used the condiments package (v1.12.0), which estimates group-specific changes in pseudotime distributions along the inferred trajectory.

HLA Typing and allele frequency analysis
We used the raw fastq files for estimation of the HLA genotype using arcasHLA119. In cases where samples were profiled twice for the initial and final dataset and/or for two time points, fastq files were merged to one fastq file per donor. These fastq files were processed with the commands genotype and merge. Reference data from the global population of the Czech Republic were obtained from The Allele Frequency Net Database55 using the tool HLAfreq120.

Analysis of published data
Count matrices from the RNA sequencing data were obtained from the Gene Expression Omnibus (GEO) (GSE23721847, GSE123658 (unpublished), GSE1058646) and processed with DESeq2121. Raw fastq files from the Single Cell 3’ sequencing of PBMC of four Finnish children at risk of developing Type 1 diabetes and their gender, age and HLA matched control children were obtained from the European Genome-Phenome Archive (EGAD0000100576826) and mapped with Cell Ranger software (10x Genomics, cellranger-7.1.0) to the GRCh38 human reference genome. Count matrices were processed with Seurat similarly to our data. Processed scRNAseq data of PBMCs from 46 T1DM cases and 31 matched controls was obtained from Synapse under the accession code syn5364184929. Processed scRNAseq data of PBMCs from 5 T1M donors and 3 healthy donors was obtained from GEO GSE22129730. Processed scRNAseq data of PBMCs from 12 T1D donors and 12 healthy donors was obtained from ParseBio (Supplemental Table 1). Genes from all datasets were ranked by the fold changes obtained from comparing children with T1D to healthy donors (DESeq LogFC). For comparison between datasets, the ranks were converted to 20-quantiles.

For the genotype-expression analysis of BTN3A2, we used bulk RNA sequencing data from GSE23721847, GSE123658 (unpublished), EGAD00001005767 (Kallionpaa), and scRNAseq data from the current study and from the HPAP repository. In all cases, raw fastq files were processed using arcasHLA as described above to obtain HLA genotypes for all the donors. HLA alleles were typed to the level of two fields, which distinguishes specific HLA proteins. The information about the expression of BTN3A2 for each of the donors was obtained from the count matrices that were generated for each dataset as described above. The expression of BTN3A2 for the donors in each dataset was scaled to a range from −1 to 1 to allow comparison between datasets, i.e., the donor with the lowest expression of BTN3A2 in each dataset had a value of −1, and the donor with the highest expression in each dataset had a value of 1. The effect of study, diabetes status and particular alleles in all MHC-I and MHC-II loci was calculated using a generalized linear model with gaussian distribution. The analysis was performed for each locus separately.

Annotation of cells with previously published signatures
For annotation of cells with the TR3-56 signature, count matrix from RNA sequencing of the flow sorted CD3+ 56+ cells, NK cells, CD3+ CD56-, and CD8+ cells was obtained from the Gene Expression Omnibus (GSE106082)16. CD4 and CD8 datasets were annotated using SingleR122. For annotation of Treg cells, the differentially expressed genes from dataset GSE1565958 were retrieved from the mSigDB (#M3563). The function AddModuleScore from the Seurat package was used to calculate the expression of the genes in the Treg dataset.

Analysis of scRNAseq data from the HPAP database
From the HPAP repository, we used a collection of scRNAseq samples from splenocytes and lymph nodes of the deceased healthy or T1D donors that were processed using the HPAP CITEseq: Dual-index 3’ HT scRNAseq with Antibody Derived Tags and Hashtag Oligos protocol. The full protocol can be accessed at the HPAP database. For our analysis, raw fastq files were downloaded from the HPAP web and mapped to the GRCh38-2020-A human reference genome using cellranger-7.1.0. Count matrices were merged and subjected to normalization (scale factor=1×104), scaling and dimensional reduction (PCA with 15 principal components on 2000 variable features and UMAP) and Louvain clustering. Clusters representing dead cells or contaminating cell types other than T or NK cells were removed. For analysis of the Treg subpopulations, splenic Tregs were extracted from the whole dataset and subjected to new to normalization (scale factor=1×104), scaling and dimensional reduction (PCA with 10 principal components on 1000 variable features and UMAP) and Louvain clustering. The batch effect of cells from different experiments was removed using STACAS (v 2.1.3).

Analysis of TCR repertoires
TCR repertoires of the CD4+ and CD8+ T cells from initial and final experiments were profiled using 10x Chromium Single Cell V(D)J Enrichment Kit, Human T Cell (#PN-1000005) and mapped by the 10x Cellranger 5.0.1 sotware (vdj tool) to human reference obtained the International ImMunoGeneTics Information System (IMGT)123 for the immune receptor repertoire profiling. Additional V(D)J sequences were extracted from the gene expression library using the MiXCR 3.0.12 software124. For the clonal expansion analyses, cells with the same nucleotide sequences of the CDR3α and CDR3β were considered clones. In cases where we detected two productive rearrangements of CDR3α, the cell was considered the same clone if the nucleotide sequences of CDR3β and both CDR3α were shared, or if the nucleotide sequences of CDR3β and one CDR3α were shared while the second CDR3α was not detected.

For the repertoire overlap analysis, unique amino-acid sequences of the CDR3α or CDR3β were extracted for each donor. Then, the percentage of overlapping sequences in one participant with all of the other participants was calculated and quantified in the following comparisons: (i) healthy donors with other healthy donors, (ii) healthy donor with a T1D donor, or iii) two T1D donors. In addition, the Self T1–Self T0 comparisons were calculated as the overlap between CDR3α or CDR3β sequences in the same donor in the two time points.

For the analysis of the length of the CDR3β, unique amino-acid sequences of the CDR3β were extracted for each donor. The lengths of the sequences were quantified as the frequency of each length in a particular donor and visualized in children with T1D and healthy donors.

For analysis of the biochemical properties of the CDR3 sequences, we used the R package Peptides125. For each of the biochemical properties available in the package, the values for all the unique CDR3β sequences were calculated and compared between healthy donors and children with T1D. Each sequence was taken into the analysis, the number of times it occurred in the dataset, i.e., a CDR3β sequence representing an expanded clone that was detected in three donors, with 50 occurrences in donor 1, 10 occurrences in donor 2 and 1 occurrence in donor 3 would be counted 61 times. To see the global biases of the repertoires in children with T1D and healthy donors not influenced by the clonal expansion, we calculated the biochemical properties for all the unique CDR3β sequences in each donor, and averaged these values per donor.

For all the analyses except for clonal expansion, the analysis was performed on conventional cells only, i.e., MAIT cells, NK cells, and γδ T cells were excluded from the CD8+ dataset, and iNKT cells were excluded from the CD4+ dataset.

The analyses of CDR3β length and biochemical properties were performed also on samples from splenocytes of the non-diabetic or T1D donors from the HPAP dataset. In the HPAP repository, the TCR libraries were prepared as described previously126. Briefly, T cell receptor beta chain gene rearrangements were bulk sequenced from gDNA using custom primers. Typically, 100ng of input DNA was amplified per replicate. Sequencing libraries were prepared using 2×300 bp paired-end kits (Illumina MiSeq Reagent Kit v3, 600-cycle, Illumina Inc., San Diego, CA). The raw sequencing data (fastq files) were downloaded from the HPAP database and mapped to the using Mixcr v4.5.0 to the human reference obtained from the International ImMunoGeneTics Information System (IMGT)123.

Flow cytometry analysis
The list of used antibodies is provided in Supplementary Table 6. Flow cytometry was performed in two batches. In each batch, aliquoted PBMCs stored in the liquid nitrogen were gently thawed by slow, sequential dilution in RPMI medium containing 10% FBS. Cells were counted, and a comparable number of cells (1.2 million in batch one, 2 million in batch two, depending on sample availability) were used for staining. The counted cells were divided into three equal parts. The first part was used for extracellular staining of live cells, the other two parts were used for intracellular staining of fixed cells. To prevent unspecific binding to Fc receptors, Human TrueStain FcX (BioLegend #422301) was added to all staining mixes. For staining of extracellular markers (Panel 1: CD4 OKT4 BV421, CD16 BV510, TCR Vα24-Jα18 BV605, CD45RO BV650, CD197 BV711, CD45RA FITC, CD3 UCHT1 PE-Cy5, HLA-DR PerCP/Cy5.5, CD183 PE, CD25 PE/Fire700, CD196 PE/Dazzle 594, TCR gamma/delta PE-Cy7, CD8 APC, CD56 APC-R700), cells were incubated with diluted antibodies and LIVE/DEAD Fixable Near-IR Dead Cell Stain Kit (Invitrogen #L34975) for 30 min on ice immediately after isolation. For staining of intracellular markers and transcription factors, cells were first stained for 30 min on ice with the mix containing LIVE/DEAD Fixable Near-IR Dead Cell Stain Kit (Invitrogen #L34975) and extracellular markers (CD3 MEM-57 AF700, CD4 SK3 BV421, CD8 APC, CD25 PE/Fire700, CD45RA FITC, CD45RO BV650, CD45R/B220 BV510) and then fixed and permeabilized using the eBioscience Foxp3/Transcription Factor Staining Buffer Set (Invitrogen #00-5523-00) according to the manufacturer’s instructions. The cells were washed, and stained overnight with antibodies for intracellular markers in 4°C (Panel 2: extracellular markers and CD127 BV605, Ki-67 BV750, HLA-DR PerCP/Cy5.5, CD226 PE, FOXP3 PE/Dazzle 594, CD137 Pe/Cy7; Panel 3: extracellular markers and CD127 BV605, Ki-67 BV750, IκBα AF488, GZMK PerCP/Cy5.5, CD184 PE, Eomes PE/CF594, GZMB PE/Cy7).

Flow cytometry was carried out with a Cytek Aurora flow cytometer, configuration 4L 16V-14B-10YG-8R (Cytek). Data were analyzed using FlowJo software v10.10 (BD Life Sciences). Anomalies in flow rate, signal acquisition, and dynamic range were removed using the FlowJo plugin FlowAI127.

HPAP flow cytometry
Flow cytometry data (fcs files) from PBMC of donors with T1D and healthy donors were downloaded from the HPAP. These samples were collected from deceased donors on the date of death and analyzed either fresh or after cryopreservation. For our analysis, we used only PBMC samples from donors with T1D and healthy donors. Two samples were excluded due to suspected duplicity, and six samples were excluded based on the lack of compensation metadata in the fcs file. The list of used and excluded samples is provided in Supplementary Table 7. Samples were analyzed in groups based on three different panels for staining and manually compensated when needed. We analyzed 41 samples in Panel 1 – CD4 phenotyping panel, 19 samples in Panel 2 – CD8 phenotyping panel and 26 samples in Panel 3 – CD8 phenotyping panel, focused primarily on antigen-specific cells. In cases when the same patient had measurements in multiple panels, these values were averaged. For analysis, we excluded donors younger than 5 years as appropriate age-matched healthy controls were not available.

Statistical analysis
The statistical analysis was performed using tests indicated in the Figure legends using R v4.2.1. Quantitative scRNAseq and flow cytometry data were tested using the nonparametric Mann–Whitney test without a correction for multiple comparisons or, for a comparison of the same donors at two different timepoints, using the paired Wilcoxon signed rank test without a correction for multiple comparisons. Two-tailed tests were performed. The number of biological replicates (cells and donors) is indicated in the respective Figure legends.

Bayesian statistics
Differences in population composition were assessed by a Bayesian generalized linear model with a negative binomial response, with an offset term to normalize for the total number of cells in a sample. The input data consisted the counts of cells for each Level 3 population, and the model included a fixed effect for the type of sample (healthy, T1D at time 1, T1D at time 0) and the Level 3 population. Additionally, random effects of the type of sample were added, grouped by both Level 2 and Level 3 populations. The shape parameter was also predicted, with a fixed effect for Level 2 populations and a random intercept for Level 3 populations. The model was fitted with the brms package128. We used the default priors in the brms package, specifically: (i) The intercept for the mean has student-t prior with 3 degrees of freedom and scale 2.5, with location guessed by brms from the overall range of the data (1.81 for CD4 data and 1.97 for CD8 data), (ii) Student-t prior with 3 degrees of freedom and scale 2.5 centered at 0 for the intercept for overdispersion, (iii) Improper flat prior for all fixed effects, (iv) Half Student-t prior with 3 degrees of freedom, scale 2.5 and location 0 for standard deviation of random effects, v) LKJ(1) prior for correlations between random effects.

Ethical consent
The study was approved by the institutional Ethics Committee (EK-819/20) of the Motol University Hospital. Written informed consent was obtained from all the participants and their legal guardians.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "It is important to profile T cells in type 1 diabetes (T1D) since T1D is thought as an autoimmune disease. The study by Veronika et al. investigated the T cell subsets and their features in T1D by comparing with healthy donors, revealing a disrupted balance in T cell phenotypes across several T cell subsets in T1D patients. While the topic is intriguing, the study's persuasiveness is somewhat diminished by the presentation of the results. The authors should conduct more robust experiments or perform a more rigorous analysis to substantiate their findings.\nMajor concerns:\n1. The authors performed batch correction on six batches from initial scRNAseq dataset and final scRNAeq dataset. However, there seems to be a significant batch effect between the initial and final datasets after batch effect correction (Figure S2J), and the batch effect still could significantly bias the results. It is much better to project all the cells on the UMAP plot to examine the integrated results, which could help we understand whether the T cell subsets are common or sample specific.\n2. The authors examined the fraction of T cell subsets in healthy donor and T1D patients and their change mainly using FACS sorting (Fig 2D, 2F, 3E-F, 4F-H, 5H, 6F-G). However, the authors did not examine the fraction of T cell subsets in scRNAseq data. Whether the FACS-results and scRNAseq results are consistent?\n3. The authors identified DEGs and signals of CD8+T cells (or CD4+T cells) between T1D patients and healthy donors without considering T cell subsets (Fig 2A). The author further identified the DEGs of each T cell subset between T1D patients and healthy donors (Figure S6B). However, it is still unclear whether the identified DEGs of these T-cell subsets are consistent with each other or T cell subset specific? The authors could make a more explicit figure to clearly show DEGs and signals are shared by multiple T cell subsets or specific T cell subsets.\n4. The feature, repeatability and functions of these Treg subsets in T1D are not well demonstrated. E.g. the authors claim that Treg1 has phenotypes of naïve Treg. However, Treg1 also highly expressed IL4R, Il7R, LAG3, GZMK and CD226 (Fig 3B) which are memory T or effector T specific genes. How about the fraction of Treg1 in each sample based on scRNAseq data? Whether Treg1 are present in each sample? Is Treg1 more likely to differentiate into Treg4 in healthy donors than in T1D using trajectory analysis?\n5. It seems that TR3-56 cells are also located on conventional CD8+ T cells (Fig4B). What are the differences of the TR3-56 cells on conventional CD8+ T cells and MAIT/γδ T?\nMinor concerns:\n1. The expression of effector genes (NKG7, GNLY, CCL4, CCL5, CX3CR1, TNFRSF4, TNF, and TNFRSF9) in EGAD00001005768 are different the other data. The authors should examine and discuss this inconsistency.\n2. Contaminating NK cluster (Fig. 5F) should be removed at the quality control or cell annotation rather than the final analyses.\n3. Some errors should be checked and corrected, such as\n1) The term \"NK T cells\" should be clarified. Is it referring to NK cells or Natural Killer T (NKT) cells? (Figures S3A-D)\n2) The expression of CD79A, a B cell marker, by iNKT cells is puzzling and requires an explanation.\n3) Legends should be added to figures for clarity (e.g., Figs 6B-6C, S9A, F,).", "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": "The manuscript by Niederlova et al describes the analysis of T-cell population from children with and without T1D at two time points a year apart. The manuscript points to some novel aspects of human T cell biology relating to T1D. Several parts of the manuscript could be improved. These are listed below.\nExcerpts from the text are in italics, my comments are in regular text\nMajor points:\n1. The title not clear. I assume it is the cells, not the phenotypes, that are in imbalance.\n2. To compare the T-cell compartments in healthy donors and T1D children… It is more respectful of people with T1D to refer to people, or children in this case, with and without T1D.\n3. The composition and the matching of the individuals in the cohorts is critical to this study. 30 T1D vs 13 HC, but not HLA matched? What about age and sex?. This is shown diagrammatically in Figure1, but needs to be stated explicitly and clearly in the Results section.\n4. First, we performed an initial scRNAseq analysis of CD4+ and CD8+ lymphocytes from 16 samples…Were these children with, or without T1D? How were these 16 samples chosen? Did the CD45R0/RA+ correlate with the subject’s age?\n5. Can the authors validate that impact of batch specific effects did not impact their analysis and interpretation of the data? Supplementary materials show plots, but can this reduction in batch effects be quantified more objectively?\n6. Have the authors stratified for the age of the participants? Could the differences described be attributable to the age of the participants more than their diagnosis with, or without, T1D?\n7. The first paragraph of the discussion is irrelevant and should be deleted. If the authors feel this information should be included it should be in the Introduction.\nMinor points\n1. In my opinion it is not true to say that there has been a ‘steep increase in T1D incidence…’ (Introduction paragraph 1). This should be reworded to something less hyperbolic.\n2. (Introduction paragraph 2) Despite numerous studies in the field…please cite some of these studies\n3. The Network for Pancreatic Organ Donors with Diabetes. Please include the acronym nPOD.\n4. We used a balanced medium-sized cohort…it is not clear what ‘balanced’ refers to here. Medium-sized is also subjective, state the number of individuals in each cohort.\n5. Blood of 30 donors with T1D was taken 4-9 (median 9) days…this is not possible. If the range is 4-9, the median cannot be 9.\n6. The T1D patients showed higher overall blood lymphocytes and lower neutrophils at T0 in comparison to T1 and healthy donors (Fig. S1C). Give the data and the statistics in the text.\n7. We compared the proportions of particular subsets using both frequentists…please correct typo.\n8. …which could be further divided to naïve, central memory (Tcm), Treg, Th1/Th17, Th2, Temra, ISAGhi. Please define ISAGhi and Tsig in the text,.\n9. How many participants were included in the HPAP study which the authors reanalyzed?\n10. First, we observed that our healthy donors generally followed the common distribution of HLA alleles in the Czech population 45 with the exception of overrepresented HLA-A*02:01, whereas most T1D children carried the respective susceptible alleles. This sentence is very unclear. The non-T1D subjects were more frequently A2+ than expected?...Even though there were only 13? What is meant by ‘respective susceptibility alleles- HLA-DR4-DQ8, DR3-DQ2?", "reviewer_4": "- What are the noteworthy results?\nNiederlova et al have produced an interesting piece of research examining differences in T cell populations between children without T1D, children newly diagnosed, and those one year after diagnosis. They describe in great detail changes in effector T cell populations. A strength of the work is confirmation of the findings by re-examination of previously published scRNAseq and flow cytometry data sets.\n- Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references.\nThis work is useful to the field as we attempt to understand phenotypes of both autoantigen- specific and non specific T cells in T1D. Other papers have performed scRNAseq and flow cytometric analysis of T cells before (and these are well described in the introduction of this manuscript), but the focus on children at point of diagnosis and after one year adds value in this study.\n- Does the work support the conclusions and claims, or is additional evidence needed?\nIn general yes, but see major comments below on additional analysis required\n- Are there any flaws in the data analysis, interpretation and conclusions? - Do these prohibit publication or require revision?\nYes please see the comments below that require revision\n- Is the methodology sound? Does the work meet the expected standards in your field?\nYes\n- Is there enough detail provided in the methods for the work to be reproduced?\nFor the wet lab work- yes.\nFor the data analysis -yes, but see comments on code\nMajor comments\n“Processed data required to generate figures for the manuscript are available on Zenodo: DOI: 10.5281/zenodo.14222418\nRaw data were not deposited to protect the identity of study participants.” Under FAIR data principles, I would expect to have access to the raw data upon point of publication e.g. fastq files as much information is lost by running the cellranger pipeline and conversion to a Seurat object. I don’t understand what how this could compromise the identity of study participants? Files could be renamed if they include study ID or participant names. The authors should discuss this with the editor.\nThroughout there is a problem with the definition of naïve T cells. Frequently in the text they are described as “naïve, quiescent, and/or stem-like T cells (LEF1, BACH2, IL7R, CXCR4, ZFP36L2)” or variations thereof. However, this is then being abbreviated to “naïve” in other places. The authors need to do further analysis to characterise these cells because there can be substantial differences (or overlap) between naïve/stem T cells with quiescence being a property that either may have and this is crucial to interpreting their findings. See https://www.cell.com/trends/immunology/fulltext/S1471-4906(24)00004-8 for example. In particular, I would like to know whether all of these markers are co-expressed within the same cells, or whether there are distinct clusters of “naïve” cells with different expression patterns. Also, do the “naïve” cells show clonal expansions? The expression patterns in the naïve cells are very different at T0 and T1 in fig 2b, suggesting that they are not truly naïve? Oddly, although naïve cells are mentioned in the title they are not in the abstract.\nFrom Fig 1B I do not understand how the participants were age and sex matched, and this is crucial to understanding the differences between controls and children with T1D. Please provide a table listing each donor’s characteristics and summarising mean age, sex %. If there are significant differences in the ages of the controls and children with T1D I would consider this a major flaw. In this case, the authors should produce figures showing how their key findings vary by age in controls and children with T1D.\nLikewise the manuscript states “This indicates that the selection of healthy controls (e.g., unmatched as in our study or partially matched as in 20) might have a great impact on the results concerning BTN3A2 expression.” Does this difference in HLA apply to any of the other DEG you identified?\nMinor comments\n“On the contrary, the so-called hygiene hypothesis explains the recent outburst of T1D and other autoimmune diseases by the reduction of the infection burden in children 4, 9, 10, 11 .” Please could the authors temper this statement?\n“the criterion of partial clinical remission (insulin dose adjusted HbA1c <9)” Please provide a reference for this statement\ndiabetic children Throughout please replace with “children with diabetes”. Likewise “T1D patients” to “people with T1D”\n“In the next step, we addressed whether Tregs from T1D ressemble functionally impaired Tregs from conditions in which their dysfunction is well established and characterized.” Please make it clear in the results section where you obtained the samples/data described in this paragraph and how they were integrated with your data.\n“We obtained the same results when we compared Tregs from the T1D patients to Treg-like cells in Foxp3-/- mice 56, 57 (Fig. 3H)” I am dubious of the value of including a direct comparison between humans and a mouse model and would prefer that this section is removed.\n“The first sample was obtained at a median of 9 (4-24) days after the clinical diagnosis of T1D” This does not match what is said in the results section\nInitial scRNAseq experiment: Please put the antibodies/hashtags into a table for clarity\nPlease state the number of reads/cell for the 10x sequencing\nHPAP flow cytometry :please provide a table of exactly which samples were used/excluded\nFig 6- please use asterisks just for significant p values, as including all p values makes this hard to interpret", "reviewer_5": "This manuscript presents the results of an analysis of single-cell gene expression in T cells from healthy donors and children with T1D. In the T1D group, repeated measures were taken close to the time of diagnosis, and then after 1 year. The authors found reduced expression of genes related to effector and cytotoxic T cells in the children with diabetes, which were somewhat normalized by 1 year. Strengths of this manuscript include the relatively large sample size for an scRNA-seq analysis, with repeated measures within a participant. However, there are a number of problems with the analytical methods. My comments on the paper are as follows:\n1. The DEG analysis using FindMarkers in Seurat treats each cell as an independent observation and therefore is subject to type I error inflation due to pseudoreplication. In addition, there is no correction for multiple testing, which further inflates the type I error (although the legend for Figure 2A describes a Bonferroni adjustment, but that is not described in the methods, and would also be extremely conservative for an scRNA-seq analysis).\n2. The analysis methods do not account for the correlation between T1D and controls induced by performing matching.\n3. No power calculation appears to have been performed. This makes it more difficult to interpret findings that were not statistically significant.\n4. The study design would be stronger if samples were obtained in the people with T1D prior to diagnosis, so gene expression in the period of time leading up to development of diabetes could be studied.\n5. The difference between experiment 1 and experiment 2 is unclear.\n6. Genes present in less than 3 cells were removed. It is more common to see this threshold set to 5-10% of cells, in order to remove effects of rare genes.\n7. In the methods, it is stated that cells with <200 genes are removed, and in another place it is stated that cells with <500 cells are removed.\n8. What criteria were used to select 12 PCs?\n9. What is the justification for using the GSE11057 and GSE9650 databases for GSEA rather than a more standard database?\n10. Why were fastq files merged prior to HLA typing?\n11. No justification is provided for the focus on BTN3A2.Why was this gene scaled differently than the others? The use of a range of -1 to 1 for scaling BTN3A2 seems arbitrary.\n12. In the methods for flow cytometry, the authors state that the same number of cells was used in each batch, but they also state that the same number of cells were used and these numbers look quite different.\n13. The paper is written in a somewhat non-standard way, with text relevant to the methods included in the results, and very few summary statistics in the text of the results which makes it difficult to evaluate the statements made.\n14. Which prior distributions were used in the Bayesian modeling?\n15. In the Figure 1D and 1G heatmaps, it could be useful to present the genes in the same order. Dot plots that show both the average level of expression and the percent of cells in which the genes were expressed could be informative. How were these subsets of genes selected?\n16. Figure 2a and 2b appear to be a volcano plot that incorporates both the T1D T0 vs Healthy comparison and the T1D T1 vs Healthy comparison, but this figure is not explained very well. Figure 2c could also use additional explanation.\n17. The authors claim that this study is the first analysis of repeat measures of single cell transcriptomics in T cells in newly diagnosed patients. At least one other study (PMID: 38967999) seems to meet these criteria?", "reviewer_1_attached": "It is important to profile T cells in type 1 diabetes (T1D) since T1D is thought as an autoimmune disease.\nThe study by Veronika et al. investigated the T cell subsets and their features in T1D by comparing with\nhealthy donors, revealing a disrupted balance in T cell phenotypes across several T cell subsets in T1D\npatients. While the topic is intriguing, the study's persuasiveness is somewhat diminished by the\npresentation of the results. The authors should conduct more robust experiments or perform a more\nrigorous analysis to substantiate their findings.\nWe thank the reviewer for recognizing the importance of our study and for the constructive comments.\nAs detailed below, we addressed these points by clarifying our analyses, revising figures and text, and\nperforming additional supporting analyses to improve the rigor and clarity of the manuscript.\nMajor concerns:\n1. The authors performed batch correction on six batches from initial scRNAseq dataset and final\nscRNAeq dataset. However, there seems to be a significant batch effect between the initial and final\ndatasets after batch effect correction (Figure S2J), and the batch effect still could significantly bias the\nresults. It is much better to project all the cells on the UMAP plot to examine the integrated results, which\ncould help we understand whether the T cell subsets are common or sample specific.\nWe thank the reviewer for this important point. We are sorry that we were not able to explain what the\ndata in Fig. S2J represent in a clearer way, which led to this misunderstanding. The initial scRNAseq\nexperiment (batches 1-2) contained all CD4+ or CD8+ T cells, whereas in the final experiment (batches 3-6)\nnon-naïve T cells were enriched in the 1:5 naïve:non-naïve ratio. For this reason, it is expected that the\nintegration does not normalize the differences between initial and final experiment. The fact that all the\nbatches of the final experiment are well corrected, but still distinct from the initial datasets actually\nshows that the integration method worked well and did not induce over-integration. We did our best to\nexplain this in the revised version of the manuscript: we added a new supplementary Fig. S3 showing\nthat the integration effectively removed differences between batches, but not between enriched and\nnon-enriched samples. Moreover, we added a new Supplementary Table S3 specifying the counts of\nsamples in individual batches and we revised the text of the Results and Methods sections accordingly.\nWe added UMAPs with labeled cells based on the batches showing the enrichment of non-naïve T cells in\nthe final experiment according to the expectation (new Fig. S3I).\n2. The authors examined the fraction of T cell subsets in healthy donor and T1D patients and their change\nmainly using FACS sorting (Fig 2D, 2F, 3E-F, 4F-H, 5H, 6F-G). However, the authors did not examine\nthe fraction of T cell subsets in scRNAseq data. Whether the FACS-results and scRNAseq results are\nconsistent?\nWe are thankful for these comments. As we decided to enrich non-naïve cells in the final experiment, the\nfractions of the T-cell subsets are altered significantly (see newly added Fig. S3J) and therefore the\ninterpretation of their abundance is difficult. Nevertheless, we quantified cluster frequencies and\nevaluated the changes by frequentist and Bayesian statistics (original Fig. S5 which is now Fig S6). To\nshow that in the initial scRNAseq experiment (non-enriched samples), the FACS and scRNAseq data are\nconsistent, we added a new Fig. S10A which shows the correlations of the frequencies of naïve CD4+ T\ncells, CD4+ Treg cells, naïve CD8+ T cells and CD8+ γδTgd cells measured by FACS and scRNAseq (new Fig.\nS10A).\n3. The authors identified DEGs and signals of CD8+T cells (or CD4+T cells) between T1D patients and\nhealthy donors without considering T cell subsets (Fig 2A). The author further identified the DEGs of\neach T cell subset between T1D patients and healthy donors (Figure S6B). However, it is still unclear\nwhether the identified DEGs of these T-cell subsets are consistent with each other or T cell subset\nspecific? The authors could make a more explicit figure to clearly show DEGs and signals are shared by\nmultiple T cell subsets or specific T cell subsets.\nWe thank the reviewer for this important point. We agree that the original Figure S6B was difficult to\ninterpret, largely because we tried to summarize DEGs across all 14 subpopulations.\nTo address this, we revised and expanded Figure S6 by splitting it into two new supplementary figures:\nFig. S8 focusing on the contrasts between T1D patients and healthy donors and Fig. S9 focusing on the\nDEG between T1D patients with/without ketoacidosis and T1D patients with/without partial remission.\nThis allowed us to allocate more space to analysis of DEG in subpopulations of CD8+ T cells (Fig. S8A-B)\nand CD4+ T cells (Fig. S8C-D).\nIn the new Fig. S8: Panels A and C show heatmaps highlighting the top DEGs for each subpopulation of\nCD8+ T cells (Fig. S8A) or CD4+ T cells (Fig. S8C). We focused on the top 5 DEGs per subset to ensure that\nboth shared and subset-specific signals are clearly visible. For example, TNF is upregulated in healthy\ndonors versus T1D T0 in CD8+ Temra cells; IFI27 and CD38 are upregulated in T1D T1 versus T0 in\nproliferating cells (Fig. S8A); and GZMK is upregulated in T1D T0 versus healthy in CD4+ Temra cells (Fig.\nS8C).\nPanels B and D show heatmaps quantifying the overlap of the top 100 up- and down-regulated genes\nacross subpopulations, illustrating which signatures are shared and which are subset-specific. This\nanalysis shows that naïve and memory subsets largely share similar gene expression patterns, while\nsome effector or unconventional subsets show unique patterns.\nWe believe these new figures address the reviewer’s request by explicitly visualizing shared and subset-\nspecific DEGs.\n4. The feature, repeatability and functions of these Treg subsets in T1D are not well demonstrated. E.g.\nthe authors claim that Treg1 has phenotypes of naïve Treg. However, Treg1 also highly expressed IL4R,\nIl7R, LAG3, GZMK and CD226 (Fig 3B) which are memory T or effector T specific genes. How about\nthe fraction of Treg1 in each sample based on scRNAseq data? Whether Treg1 are present in each\nsample? Is Treg1 more likely to differentiate into Treg4 in healthy donors than in T1D using trajectory\nanalysis?\nWe thank the reviewer for this insightful comment. We agree that in conventional T cell subsets, markers\nsuch as IL4R, CD226, GZMK, and LAG3 are not typically associated with a naïve phenotype. However, we\nwould like to highlight that Treg1 cells also express classic naïve markers, including TCF7, LEF1, SELL, and\nCCR7. IL7R is indeed expressed on naïve T cells and upregulated on central memory cells, but\ndownregulated in effectors, which supports our characterization. We have clarified in the revised\nmanuscript that Treg1 cells co-express genes associated with a naïve phenotype as well as genes linked\nto Treg dysfunction (e.g., IL4R, CD226, GZMK), while lacking key Treg functional markers such as FOXP3,\nCTLA4, IL2RA, and TIGIT.\nRegarding their distribution, Treg1 cells are present in each sample, as shown in the frequency plots in\nFig. S12A.\nTo further investigate the relationship between Treg subsets, we performed a pseudotime analysis using\nthe slingshot package (PMID: 29914354) and validated that the inferred trajectory reflects Treg\ndifferentiation (new Fig. 3G-L). We additionally used the condiments package (PMID: 38280860) to test\nfor differences in progression along the trajectory between groups. This analysis revealed that Tregs are\nmore activated in healthy donors compared to individuals with T1D, suggesting differential activation or\ndifferentiation states.\nWe have updated the manuscript text accordingly to reflect these findings.\n5. It seems that TR3-56 cells are also located on conventional CD8+ T cells (Fig4B). What are the\ndifferences of the TR3-56 cells on conventional CD8+ T cells and MAIT/γδ T?\nWe thank the reviewer for this question. To clarify this point, we extracted TR3-56 cells from each\nparental cluster (conventional CD8+ T cells, MAIT cells, γδ T cells, and NK cells) and compared their TR3-\n56 scores (new Fig. 4D). This analysis showed that TR3-56 cells originating from MAIT and γδ T cell\nclusters had the highest TR3-56 scores, indicating that these might be the bona fide TR3-56 cells. We also\nperformed DEG analysis to study the differences between TR3-56 cells coming from the different parental\ncluster, which revealed the signature genes of the parental clusters, such as CD8A and CD8B in\nconventional CD8+ T cells, IL7R, KLRB1 in MAIT cells and TRDC, TRGV2 in Tgd cells (see Response Figure 1\nbelow).\nResponse Figure 1. Heatmap of differentially expressed genes between TR3-56 cells\ncoming from clusters of Conventional CD8 T cells, gd T cells and MAIT cells.\nMinor concerns:\n1. The expression of effector genes (NKG7, GNLY, CCL4, CCL5, CX3CR1, TNFRSF4, TNF, and\nTNFRSF9) in EGAD00001005768 are different the other data. The authors should examine and discuss\nthis inconsistency.\nWe thank the reviewer for this comment. We agree that the expression of certain effector genes in the\nEGAD00001005768 dataset, as well as in the syn53641849 dataset, differs from our data and other\npublic datasets. A direct comparison of results from different studies is crucial for general understanding\nof the problem. Yet, this is not routinely done because of several challenges such as technical differences\nand biases among particular studies. For this reason, we compared our results with all relevant\ntranscriptomic studies with a similar design. This analysis confirmed that our results are consistent with\nthe majority of these studies. An unanimous consensus across studies was not anticipated because of the\nabove-mentioned discrepancies. We believe that a detailed exploration of the reasons for this\ndiscrepancy is outside the scope of our current study.\nWe looked at the genes which showed differences in the studies mentioned by the reviewer in a closer\ndetail by performing a PCA analysis of normalized log fold changes of all genes across studies (see\nResponse Figure 2 below). This analysis shows that these genes cluster together (A), indicating that the\ndifferences are not specific to individual genes but rather reflect broader technical variability in sample\nprocessing across studies (B). We defer to the editor’s decision regarding whether further discussion of\nthis point is necessary.\nResponse Figure 2. PCA of gene expression fold changes across datasets. (A) PCA showing the distribution of genes based on\ntheir normalized log-fold changes across studies. (B) PCA biplot showing the contribution of individual datasets to gene-level\nvariability, with arrows indicating the direction and strength of dataset-specific effects.\n2. Contaminating NK cluster (Fig. 5F) should be removed at the quality control or cell annotation rather\nthan the final analyses.\nWe removed the NK cluster prior to the final analyses.\n3. Some errors should be checked and corrected, such as\n1) The term \"NK T cells\" should be clarified. Is it referring to NK cells or Natural Killer T (NKT) cells?\n(Figures S3A-D)\nWe thank the Reviewer for spotting this typo. It is NK cells. This is corrected in the main text.\n2) The expression of CD79A, a B cell marker, by iNKT cells is puzzling and requires an explanation.\nWe thank the reviewer for raising this interesting point. Although expression of CD79A by iNKT cells is\nindeed puzzling, it has been reported previously in a subset of iNKT cells from human peripheral and cord\nblood (PMID 40305288 and PMID 32528956). As a further validation, we found increased CD79A\nexpression in iNKT cells in RNAseq data from the Immgen database (see figure below). Thus, we believe\nthat expression of CD79A in some iNKT cells is real. The biological relevance of this is not known and is\nquestionable for two reasons: 1) Jayasinghe et al. (PMID 40305288) showed that CD79A protein can not\nbe detected on the surface of iNKT cells, 2) compared to B cells, the expression is ~ 100 times lower than\nin B cells (see Response Figure 3 below). Nonetheless, we believe that CD79A is still useful as a specific\nmarker that distinguishes iNKT cells from conventional T cells, thus we decided to keep it in the heatmap\nin the Fig. S5C.\nResponse Figure 3. Expression of Cd79a by subsets of T cells, B cells and NKT cells. Data from the ImmGen database\nhttps://www.immgen.org/) were visualized using the Gene Skyline application with reference data from ImmGgen ULI RNASeq.\n3) Legends should be added to figures for clarity (e.g., Figs 6B-6C, S9A, F,).\nWe thank the reviewer for noticing this. We added the labels to the panels mentioned.", "reviewer_author_response": "We thank the reviewer for recognizing the importance of our study and for the constructive comments. As detailed below, we addressed these points by clarifying our analyses, revising figures and text, and performing additional supporting analyses to improve the rigor and clarity of the manuscript."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "We are generally satisfied with the authors' responses, except for their response to minor concern #2. The authors identified a contaminating NK cluster within the Temra cluster (Fig. 5F) in the original manuscript. We suggest the authors remove the NK cluster (Fig. 5F) during quality control rather than in the final analyses because the NK contamination may bias all the results in fig 5. Surprisingly, the authors only removed the NK cluster in Fig 5F while Fig5A-5G remained unchanged. We still believe that all panels in fig 5 should be remade after removing the NK cluster.", "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": "The authors have worked hard to improve the quality and clarity of the manuscript. In this process they have addressed most of my concerns.\nOne outlying issue that I feel needs to be mentioned is the HLA-matching. While I agree that 'perfect' HLA matching is not possible, it it possible to match for the HLA alleles of known significance for T1D (HLA-DR4, DQ8, DR3 and DQ2). I also accept that this 'matching' is not perfect it is nonetheless an attempt to avoid being mislead by difference attributable to HLA alleles. Not doing this leaves open the possibility that the differences observed between individuals with and without T1D may simply be attributable to the HLA alleles within each group, not their T1D status.", "reviewer_4": "Thank you for your revisions which have answered all of my queries except for the question below on raw data availability. Although I would still expect the raw data to be made available I am also happy to defer to the editor on what is the minimum acceptable to the journal.\n“Processed data required to generate figures for the manuscript are available on Zenodo: DOI: 10.5281/zenodo.14222418 Raw data were not deposited to protect the identity of study participants.”\nUnder FAIR data principles, I would expect to have access to the raw data upon point of publication e.g. fastq files as much information is lost by running the cellranger pipeline and conversion to a Seurat object. I don’t understand what how this could compromise the identity of study participants? Files could be renamed if they include study ID or participant names. The authors should discuss this with the editor.\nWe thank the reviewer for this comment. The privacy challenges associated with sharing next-generation sequencing data have been discussed recently in multiple studies and reviews (e.g."}, "author_response": "We thank the reviewer for their positive feedback and we apologize for the misunderstanding regarding the contaminating NK cluster. We agree with the reviewer that the original Figures 5A-G may have been affected by the presence of the NK-cell contaminants and thank them for highlighting this point. In the revised manuscript, we have removed the contaminating NK cluster from all analyses of the CD4⁺ Temra subclustering figure (now Fig. 7; corresponding to Fig. 5 in the original submission). The entire figure (panels A–G) was recomputed after excluding these NK cells. We have also updated the Methods to explicitly describe this step. We note that we did not remove NK cells from the global dataset used for the remainder of the manuscript, as doing so would have required reprocessing all analyses and figures without altering any conclusions outside of the Temra subclustering."}]

Decision letter