This study was a prospective, controlled trial aimed at analysing the possible immune mechanisms in T2DM patients. The study was registered as a Chinese clinical trial on the WHO international clinical trial registry platform (ChiCTR2100049613). All methods were carried out in accordance with the CONSORT statement. All experimental protocols were approved by the Ethics Committee of the First Affiliated Hospital of Anhui University of Traditional Chinese Medicine. The ethics approval number is 2021AH-39. All participants provided written informed consent before participating. Three healthy participants and three patients with T2DM were included.

Patients were diagnosed with T2DM according to the guidelines for the prevention and treatment of T2DM in China (2020 edition). The diagnostic criteria for healthy participants were the absence of a history of systemic disease, such as hypertension, T2DM, or cardiopulmonary insufficiency, and no use of systemic or topical medications. The inclusion criteria for T2DM patients were a glycosylated haemoglobin A1c (HbA1c) level ≤ 7.5%, course of T2DM disease ≤ 3 years, and aged 18–70 years; both sexes were included. The participants were informed of the study procedures and voluntarily signed an informed consent form.

The exclusion criteria for T2DM patients were the presence of T1DM, gestational DM, T2DM requiring insulin therapy, or other special types of DM; acute complications of T2DM; severe cardiovascular and cerebrovascular diseases; severe primary diseases, such as liver, kidney and haematopoietic system diseases; allergy to the known ingredients of the study drug or chronic allergies; pregnancy, lactation, having recently given birth, or planning to become pregnant; long-term alcoholism, drug dependence or mental illness; and participation in another drug clinical trial within one month before the screening period for this study. Finally, all participants were deemed suitable for participation in this clinical study according to the opinion of the investigator.

The collected background information included age, sex, height, weight, etc. Fasting blood glucose (FBG), alanine aminotransferase, aspartate transaminase, blood urea nitrogen, and creatinine levels were determined with a Beckman AU5800 fully automatic biochemical analyser (Beckman, US). The level of HbA1c was determined with a Bio-Rad D-100 HbA1c analyser (Bio-Rad, US). The level of fasting insulin (FINS) was determined with an AutoLumo A2000 Plus fully automatic chemiluminescence analyser (Autobio, Zhengzhou, China). The homeostasis model assessment–IR (HOMA–IR) score was calculated as FINS (μIU/ml) × FBG (mmol/L) ÷ 22.5 (16). Safety indicators included blood pressure, heart rate, routine blood tests, and routine urine tests ( Figure 2 ).

ScRNA-seq was performed by BGI Shenzhen. The 10x Genomics Chromium allows high-throughput single-cell 3’ mRNA quantitative analysis. Afterwards, 5 mL of whole blood containing ethylene diamine tetraacetic acid was added to a 15 mL centrifuge tube with 3 mL of Ficoll lymphocyte separation media. An equal volume of 1X PBS was added to each blood sample. The diluted blood samples were layered carefully in Ficoll lymphocyte separation liquid and then centrifuged at 400 × g at 18–20°C for 30 min continuously. The mononuclear cell layer was transferred to a 15 mL sterile centrifuge tube with a sterile pipette. Three volumes of 1× PBS were added to the lymphocyte layer, which was carefully mixed via pipetting. The samples were subsequently centrifuged at 400 × g, after which the supernatant was discarded. Then, 6 mL of 1× PBS was added to the lymphocyte layer, which was again carefully mixed via pipetting. The mixture was subsequently centrifuged at 400 × g for 10 min, and the supernatant was discarded. The cells were resuspended in the desired volume of 1× PBS and stained with 0.4% trypan blue. Samples with greater than 80% cell viability were used for library construction. The prepared single-cell suspensions were subsequently partitioned into gel beads in emulsion in an automated Chromium Controller, after which the mRNAs were reverse transcribed into cDNAs. The reaction system was configured in sequence for breaking gel beads in emulsion, cDNA amplification, fragmentation, end repair, A-tailing, and adaptor ligation polymerase chain reaction. After reacting at a suitable temperature for a fixed period, the products were separately purified in an appropriately configured reaction system. After library quality control, single-stranded polymerase chain reaction products were produced via denaturation. Single-stranded cyclized products were produced with a circularization reaction system. Single-stranded circular DNA molecules were replicated, and a DNA nanoball that contained multiple copies of DNA was generated. DNA nanoballs of sufficient quality were loaded into patterned nanoarrays and sequenced through combinatorial probe–anchor synthesis.

The raw gene expression matrix generated from each sample was aggregated via Cell Ranger (v5.0.1) (17), which is provided on the 10x Genomics website. Downstream analysis was performed with the R package Seurat (v 3.2.0) (18). Specifically, cells with fewer than 200 genes or with > 90% of the proportion of the maximum genes were filtered. For the mitochondrial metric, the cells were sorted in descending order of the mitochondrial read ratio, and the top 15% of the cells were filtered. Potential doublets were identified and removed via doublet detection (19). Cell cycle analysis was performed with the cell cycle scoring function of the Seurat program. The gene expression dataset was normalized. Uniform manifold approximation and projection (UMAP) was subsequently used for two-dimensional visualization of the resulting clusters. For each cluster, marker genes were identified with the FindAllMarkers function as implemented in the Seurat package (V3.2.0, logFC > 0.25, minPct > 0.1 and Padj ≤ 0.05). The clusters were then marked as known cell types via the scRNA-seq atlas method (20). Differentially expressed genes (DEGs) across different samples were identified with the FindMarkers function (logFC > 0.25, minPct > 0.1 and Padj ≤ 0.05). Volcano plots were created with the ggplot2 package. The threshold for the log fold change was set at 0.2, and that for p values was set at 0.05. Gene Ontology (GO) analysis was performed via the phyper function of the R package (R-3.1.1). Kyoto Encyclopedia of Genes and Genomes (KEGG, V93.0) enrichment analysis results. GO and KEGG pathways with p or Q values ≤0.05 were considered significantly enriched (21).

The Type 2 Diabetes Knowledge Portal (https://t2d.hugeamp.org/) contains summary data on genetic correlations, genome annotations, bioinformatics results, expertise in T2DM and related traits, blood glucose, etc. The human genetic evidence calculator integrates several kinds of human genetic results to quantify genetic support for the involvement of a gene in a disease or phenotype of interest (22).

The statistical analyses were conducted in SPSS 23.0. Continuous variables are expressed as the means ± standard deviations. Categorical variables are presented as numbers or percentages. Two-group comparisons were conducted via independent-samples t tests or chi-square tests. A p value less than 0.05 was considered to indicate statistical significance.

There were no significant differences in age, sex, or disease course between the two groups, but there were significant differences in body mass index (BMI) or FBG and HbA1c levels between the two groups, as shown in Table 1 .

After quality control and filtering, 13,591 cells remained for analysis. Unsupervised clustering of single-cell data after normalization and aggregation was performed via Seurat 3.2.0. Four cell types were identified ( Figure 3A ). We annotated the cell types according to the expression of classic marker genes, and the classic genes that were differentially expressed in those cells were consistent with the annotations ( Figure 3B ). The percentages of the 4 cell types in patients with T2DM and healthy participants are shown in Figure 3C .

A total of 20,092 genes were identified in the 6 samples. Among the 4 cell types, T cells expressed the most genes ( Figure 4A ). A total of 3188 marker genes in the two groups were annotated to the KEGG metabolic pathway ( Figure 4B ). According to the KEGG pathway term level 2, 7 endocrine and metabolic diseases and 9 signal transduction pathways were screened ( Figures 4C, D ). In contrast, the presence of 11 marker genes, including PKM, MAPK1, MAPK3, PIK3R1, HK1, HK3, INSR, PIK3CD, SOCS1, IRS2 and TNF, was associated with T2DM status. Pathways with significant differences included the NF-κB, HIF-1 and TNF signalling pathways. Notably, the gene set enrichment analysis (GSEA) results also revealed the pathways with the greatest differences in the TNF signalling pathway ( Figure 4E ).

T cells are the main specific immune cells found in patients with T2DM. Unsupervised clustering of T cells revealed two CD4+ T-cell clusters (including 3127 cells) and three CD8+ T-cell clusters (including 3678 cells) ( Figure 5A ). T cells were annotated separately by canonical genes, and the expression of the canonical genes of these cell types was consistent with the annotation ( Figure 5B ). There were 387 marker genes found in CD4+ T cells and 684 marker genes in CD8+ T cells. The genes in CD4+ and CD8+ T cells were analysed, and the transcription signal score was calculated. The results suggested that CD4+ T cells in T2DM patients tended to be in memory and naïve states, whereas CD8+ T cells tended to be in effector and memory states ( Figure 5C ).

The differences in the expression levels of genes between the T2DM group and the healthy group were compared to construct a volcano plot ( Figure 6A ). Among these genes, 58 were upregulated, and 61 were downregulated in T cells in the T2DM group. We then conducted a correlation analysis between the DEGs and clinical characteristics. The expression levels of RPL27, TXN1P and RPL37 were negatively correlated with HbA1c. The MNDA of genes was negatively correlated with FBG levels, and the expression levels of DDX5 were positively correlated with FBG levels. The expression levels of GIMAP7 were positively correlated with HOMA–IR levels ( Figure 6B ). GO analysis revealed that the biological process (BP) terms enriched among the DEGs were related mainly to leukocyte chemotaxis, cytoplasmic translation, positive regulation of the apoptotic signalling pathway, myeloid cell activation involved in the immune response, the T-cell receptor signalling pathway, and the immune response-regulating signalling pathway. The enriched cellular component (CC) terms were associated mainly with the cytosolic large ribosomal subunit, ribosome, cytosolic ribosome, tertiary granule membrane, and cell–substrate junction. The molecular function (MF) terms were specifically related to structural constituents of ribosome, Toll−like receptor binding, RAGE receptor binding, phospholipase inhibitor activity, and cell–cell adhesion mediator activity ( Figures 6C, D ). The first gene set, HALLMARK_TNFA_SIGNALING_VIA_NFKB, was used to generate a GSEA graph via HALLMARK pathway enrichment analysis ( Figure 6E ). The WGCNA results revealed that the genes were divided into 7 modules. According to the correlation analysis between clinical characteristics and the 7 modules, FBG levels were significantly positively correlated with the brown module. FINS and HOMA–IR levels were significantly and positively correlated with the turquoise and blue modules. With respect to the base genes with differences and the brown module, we found that HLA-DRB5, AHNAK, TYROBP and AIF1 were shared genes ( Figure 6F ). The KEGG results revealed that genes in the brown and turquoise modules were involved in the TNF signalling pathway, T-cell receptor signalling pathway and NF-κB signalling pathway ( Figure 6G ). Moreover, TNFRSF1A was the core gene in terms of network interactions in the brown module ( Figure 6H ).

We further analysed the cell cycle stages of T cells from healthy participants and T2DM patients. Compared with those in healthy participants, the T cells in T2DM patients were more likely to be in the G1, S and G2M states, indicating active proliferation ( Figures 7A, B ). The expression of the 50 top genes varied with developmental time. These genes are associated with cytoplasmic translation, cell–cell adhesion mediated by integrins, and regulation of the inflammatory response ( Figure 7C ). The results revealed the dynamic expression of 6 marker genes in CD4+ T cells, and different expression levels were observed at the seven stages of disease progression. The expression of the RPL32, RPS10, RPS12, RPS14 and RPS23 genes tended to increase, whereas S100A4 expression tended to decrease ( Figure 7D ).

Monocytes were the most abundant nonspecific immune cells in our T2DM patient cohort. Unsupervised clustering revealed two dendritic cell clusters (including 233 cells), CD14+ monocyte clusters and FCGR3A+ monocyte clusters (including 2568 cells) ( Figure 8A ). Monocytes were annotated separately by canonical genes, and the canonical expression of genes in these cell types was consistent with the annotation ( Figure 8B ). There were 666 marker genes in monocytes and 776 in dendritic cells.

The differences in monocyte gene expression between the T2DM group and the healthy group were compared to construct a volcano plot ( Figure 9A ). Among these genes, 51 presented upregulated expression, and 124 presented downregulated expression. The expression levels of the CLEC7A, SIGLEC14 and AC018755.4 genes were negatively correlated with HbA1c levels. The expression level of VSTM1 was negatively correlated with FINS and HOMA–IR levels ( Figure 9B ). The GO results revealed that the enriched CC terms included mainly cytosolic ribosome, cell−substrate junctions, focal adhesion and tertiary granules. The enriched MF terms were specifically related to the structural constituents of ribosome, enzyme inhibitor activity, S100 protein binding, cytokine activity and cytokine binding ( Figures 9C, D ). The first gene set, HALLMARK_INTERFERON_GAMMA_RESPONSE, was utilized to generate a GSEA graph via HALLMARK pathway enrichment analysis ( Figure 9E ). The WGCNA results revealed that the genes were divided into 7 modules. FBG levels were significantly and negatively correlated with the red module. FINS and HOMA–IR levels were significantly negatively correlated with the brown module and the turquoise module. With respect to the base DEGs and the brown module, CLEC2B, B2M and MALAT1 were identified as shared genes ( Figure 9F ). According to the KEGG results, the genes in the red module were involved in the chemokine signalling pathway ( Figure 9G ).

We further analysed the cell cycle stages of monocytes from healthy participants and T2DM patients. In healthy participants, monocytes were common in the G1, S and G2M states, indicating more active proliferation in T2DM patients than in healthy participants ( Figure 10A, B ). The top-50 genes whose expression varied with developmental time were associated with the negative regulation of hydrolase activity, neutrophil migration, and polyamine biosynthetic processes and the positive regulation of cytokine production ( Figure 10C ). The results revealed the dynamic expression of 6 marker genes in CD14+ monocytes, and different expression levels were observed at the three stages of disease progression. The expression levels of the C1QA, HES4 and RHOC genes tended to increase, whereas the S100A12, S100A8, and S100A9 expression levels tended to decrease ( Figure 10D ).