We searched for related gene expression or high-throughput sequencing datasets using 1) “atherosclerotic” and “type 2 diabetes” as keywords, 2) the test specimens in datasets derived from human tissues, and 3) the largest possible sample size. Finally, two high-throughput sequencing datasets (GSE159984 and GSE164416) and two microarray datasets (GSE100927 and GSE28829) were obtained from the National Center for Biotechnological Information (NCBI) Gene Expression Omnibus https://www.ncbi.nlm.nih.gov/geo/) database. GSE159984 (including 28 patients with T2DM and 58 controls) and GSE100927 (comprising 29 patients with AS and 12 controls) were used to screen for DEGs, while GSE164416 (including 39 patients with T2DM and 18 controls) and GSE28829 (comprising 16 patients with advanced carotid plaque and 13 with early carotid plaque) were used as external validation datasets. Table 1 summarizes the information for the datasets selected.

DEG analysis was filtered using the “edgeR” or “limma” package (8, 9) in R (version 4.3.1), and the results were visualized using the “ggplot2” package. We obtained DEGs with |log2fold change≥1 and p.adj<0.05 in T2DM and AS diseases, respectively. Overlapping DEGs in the same direction between T2DM and AS were identified using the online Venn diagram tool (https://bioinfogp.cnb.csic.es/tools/venny/). KEGG and GO enrichment analyses of the overlapping DEGs were performed using the “clusterProfiler” package (10).

The overlapping genes were imported into the STRING database (http://string-db.org) (11) to construct a protein–protein interaction (PPI) network with complex relationships (interactions combined score >0.4), and this network was visualized in Cytoscape 3.82 (version 3.8.1). Four algorithms (MCC, MNC, Degree, and Closeness) take intersection to identify hub genes using the cytoHubba plug-in (12). Enrichment analysis and co-expression networks of hub genes were performed using “clusterProfiler” package and GeneMANIA (http://www.genemania.org/) (13), respectively. The molecular complex detection technology (MCODE), a plugin in Cytoscape, was used to deconstruct the functional modules. The selection criteria were set as degree cutoff=2, K-core=2, node score cutoff=0.2, and maximum depth=100.

The GSE164416 and GSE28829 datasets were used for external verification. Significance was calculated using the Student’s t-test, and p < 0.05 was considered statistically significant. Then, statistically significant genes were selected for further quantitative real-time PCR (qRT-PCR) analysis in animal models. Male C57BL/6J and ApoE^−/− mice were purchased from Weitonglihua Corporation (Beijing, China). To induce diabetic atherosclerosis models, ApoE^−/− mice were given streptozocin (50 mg/kg/day) by intraperitoneal injection for 5 days consecutively and fed on high-fat diet for 16 weeks (14). C57BL/6J mice were fed on normal-chow diet for 16 weeks after injected with vehicle used in the control group, pool of two groups, five mice per group. The total RNA of thoracic aortas was isolated with RNA extraction reagent (G3013, Servicebio, China) and reverse transcribed with SweScript All-in-one RT SuperMix (G3337, Servicebio, China). Reactions were run using CFX Connect (Bio-Rad) with 2× Universal Blue SYBR Green qPCR Master Mix (G3326, Servicebio, China). GAPDH was deemed as an internal control, and the results were determined with the 2^−ΔΔCt method.

Similarly, the level of significance used was 0.05.

Primers sequences are listed in Supplementary Table S1 . All animal care and experimental procedures were approved by the animal ethics committee of the Ethics Committee of Xiyuan Hospital, China Academy of Chinese Medical Sciences.

Based on the results of qRT-PCR analysis, genes with significant statistical differences (p < 0.05) were regarded as core genes. Thus, core genes were identified successfully through multiple bioinformatics mining, external datasets validation, and qRT-PCR experimental verification. To determine the value of each core genes and multiple genes in diagnosis of T2DM complicated with AS, we executed receiver operating characteristic (ROC) curve analysis, respectively. The diagnostic capacity of core genes was quantified using the area under ROC curve (AUC) in GSE164416 and GSE28829 datasets. The “pROC” R package was used to generate ROC curves (15). The greater the AUC value, the more superior the discriminatory ability of the model. An AUC closer to 1 indicates better prediction, and an AUC>0.7 indicates good diagnostic efficacy.

Transcriptional regulatory relationships unraveled by sentence-based text mining (TRRUST) were used to obtain candidate transcription factors (TFs) that regulate core genes (16). This database contains abundant information about TFs associated with target genes and their regulatory relationships with TFs. Statistical significance was defined as an adjusted p<0.05. We constructed a TF mRNA regulatory network and visualized it using Cytoscape.

We performed immune infiltration analysis to reveal the underlying immune pathogenesis of T2DM complicated with AS. QuanTIseq is a validated deconvolution-based algorithm that estimates the absolute proportions of relevant immune cell types. Thus, we used quanTIseq algorithm to obtain the immune cell infiltration differences between normal group and disease group. Based on “IOBR” packages, relative percentage, different immune cell types were analyzed using GSE164416 and GSE100927 datasets. The results were output as bar graphs and violon plots, respectively. Subsequently, the correlations between core genes and IICs were conducted using Spearman analysis, and the results were visualized using “ggplot2” package.

The study’s flow chart is shown in Figure 1 . In total, 116 and 803 DEGs were obtained from the GSE159984 and GSE100927 datasets, respectively ( Figures 2A, B ). After taking the intersection between the two datasets, 27 overlapping DEGs with the same expression trends (26 upregulated and 1 downregulated) were identified from the two datasets ( Figures 2C, D ). The list of the differential gene expression is included as Supplementary Data . The DEGs lists have been included as Supplementary Tables 2 – 4 .

To explore the potential biological function, GO enrichment and KEGG pathway for the overlapping DEGs were performed using R. According to the GO analysis findings, positive regulation of cytokine production, activation of immune response, and negative regulation of leukocyte activation were significantly enriched in the biological process (BP) entries, specific granule, collagen trimer, and specific granule membrane were significantly enriched in the cellular component (CC) entries, and oxidoreductase activity, growth factor receptor binding, and phosphotyrosine residue binding were significantly enriched in the molecular function (MF) entries ( Figure 3A ; Supplementary Table S5 ). In addition, complement and coagulation cascades, chemokine signaling pathway, and IL-17 signaling pathway were significantly enriched in the KEGG entries ( Figure 3B ; Supplementary Table S6 ). These results indicate that immune-related processes and chemokines may be crucial in the development of T2DM complicated with AS.

The PPI network of overlapping DEGs contained 26 nodes and 77 interaction pairs. The topological parameters of nodes in this network were are in Supplementary Table S7 . Four algorithms (MCC, MNC, Degree, and Closeness) in the CytoHubba plugin were used to identify the hub genes ( Supplementary Table S8 ). The middle part of Figure 3C represents the intersection of the four algorithms. The schematic diagram of PPI network deconstruction and the network constructed with hub genes is shown in Figure 3D . A total of 12 hub genes were identified, and these genes were TYROBP, IL1B, C1QA, C3AR1, FGR, CCR5, LAPTM5, C1QC, CD52, MSR1, HAVCR2, and ALOX5AP. To further investigate the biological characteristics of these genes, we analyzed the related functions of hub genes’ co-expression network using GeneMANIA database. The biological functions are associated with immune and inflammatory-related processes, such as the regulation of humoral immune response, complement activation, humoral immune response, interleukin-2 production, and negative regulation of immune system process ( Supplementary Figure S1 ).

Molecular complex detection (MCODE) plugin was used to screen out an important subnetwork in hub genes network, and this subnetwork included 7 nodes and 18 pairs. The seven nodes were TYROBP, IL1B, C1QA, C3AR1, CCR5, C1QC, and MSR1, and the scores are shown in Supplementary Table S9 . Then, the enrichment analysis revealed that the nodes in MCODE network were significantly enriched in inflammatory response and cytokine transcript regulation, specifically, the results of BP enrichment in activation of immune response, leukocyte-mediated immunity, and positive regulation of cytokine production, and the KEGG pathway enrichment in complement and coagulation cascades, efferocytosis, cytokine–cytokine receptor interaction, and type I diabetes mellitus ( Figure 4 ; Supplementary Tables S10 , S11 ). Together with the preceding results, we suggest that changes in the immune microenvironment affected by cytokines and inflammatory responses may be a common mechanism of the T2DM complicated with AS.

To identify reliable core genes, we carried out both external validation and animal experiments according to genes in MCODE network. The expression levels of seven genes, except for C3AR1, in the T2DM dataset (GSE164416) were significantly higher than that in the control samples ( Figure 5A ). Similarly, the expression levels of the seven genes were upregulated in atherosclerotic samples compared with control (GSE28829, Figure 5B ), and p<0.05 was considered statistically significant. Then, we detected the expression levels of TYROBP, IL1B, C1QA, CCR5, C1QC, and MSR1 in thoracic aortas of mice. The results demonstrated that the expression of IL1B, C1QA, CCR5, and MSR1 was increased (p<0.05) in the model mice of T2DM complicated with AS compared with control groups ( Figure 5C ). Therefore, combining the above results, we identified four core genes, which are IL1B, C1QA, CCR5, and MSR1.

Subsequently, ROC curves were generated to further evaluate the diagnostic value of the validated core genes. The AUCs of IL1B, C1QA, CCR5, and MSR1 were 0.722, 0.778, 0.668, and 0.749, respectively, in the T2DM-related validation dataset ( Figure 5D ). Moreover, the AUC values of all validated core genes were >0.7 in the AS-related validation dataset, with AUC of 0.721, 0.938, 0.962, and 0.856 for IL1B, C1QA, CCR5, and MSR1, respectively ( Figure 5E ). At the multigene expression level, after linear fitting of all validated core gene expression models, the AUC value of the multigene combined diagnosis of T2DM and AS was 0.802 and 1.0, respectively ( Figures 5F, G ). These results reveal that core genes possess good discriminatory ability, and multigene combined diagnosis has a significantly higher predictive power than single gene.

As described above, the result of the enrichment analysis suggested that the immune response might play a crucial role in the course of T2DM complicated with AS. We examine how the immune system works by immune infiltration analyses. This analysis revealed that T2DM diseases were infiltrated by several immune cells, of which regulatory T cells (Tregs), myeloid dendritic cells, and M2 macrophages occupied the top 3 most counts of immune cell subpopulations ( Figure 6A ). Subsequently, we further analyzed differences in immune cell subgroups between the T2DM samples and the control pools ( Figure 6B ). The number of M1 and M2 macrophages and CD4+ T cells in the T2DM group was significantly higher than that in the control group (p<0.05 or p<0.01), while the number of B cells and Tregs was lower (p<0.01).

For the case of AS, the immune cell compositions of the AS and control group are shown in Figure 6C . M1 and M2 macrophages and myeloid dendritic cells occupied the top 3 most counts of immune cell subpopulations. Patients with AS had significantly higher numbers of M2 macrophages and CD4+ T cells than the controls (p<0.001), and there were fewer B cells, monocytes, and Tregs than those in the control group (p<0.01) ( Figure 6D ).

In order to investigate whether or not core genes were linked to IICs, the correlations were conducted using Spearman analysis. In T2DM, there is a positive relationship between C1QA, CCR5, MSR1, and M1 macrophages, and a negative relationship between C1QA, MSR1, and CD8+ T cells, Tregs, and myeloid dendritic cells (p<0.05 or p<0.01). In AS, there have a positive relationship between IL1B, C1QA, CCR5, MSR1, and M1 macrophages, M2 macrophages, CD4+ T cells, CD8+ T cells, and myeloid dendritic cells, and a negative relationship between MSR1 and B cells and Tregs (p<0.05 or p<0.01) ( Figure 6E ). The details of differential IICs selection results and correlation analysis are shown in Supplementary Tables S12 , S13 , respectively.

TRRUST is a TF–target interaction database that shows regulatory regulation between TF and target genes. According to TF binding site information provided in TRRUST, potential key regulators for core genes were selected, a total of 10 associations between five TFs (SPI1, RELA, NFKB1, YY1, and JUN) and four core genes (IL1B, C1QA, CCR5, and MSR1). As shown in Supplementary Table S14 , SPI1 regulated two genes (IL1B and MSR1), YY1 regulated two genes (IL1B and CCR5), JUN regulated two genes (IL1B and MSR1), RELA regulated two genes (IL1B and CCR5), and NFKB1 regulated two genes (CCR5 and IL1B). Based on this result, we constructed a regulatory TF-mRNAs network using Cytoscape software ( Figure 7 ). We use different shapes or colors to distinguish different types of mRNA and TFs. This figure shows the potential pathological regulation process found in this study.