In the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database, we chose GSE122063 which included the autopsy samples of VaD for analysis. GSE122063 was based on GPL16699 which used Agilent-039494 SurePrint G3 Human GE v2 8 × 60 K Microarray to detect the expression of genes. The microarray data includes eight VaD patients, 12 AD patients, and 11 controls postmortem frontal and temporal cortex samples. Each sample was run with at least two technical replicates. Data from AD patients were excluded from analysis and VaD sample 1063 was removed due to poor data quality according to the clustering result. The raw expression matrix was directly downloaded from the website, and the SOFT format file was downloaded and parsed by the GEOquery package (Davis and Meltzer, 2007). Then we used GPL1699 to transit ID into gene names and gene symbols using merge function in R. In addition, we checked if the data need log transformation or normalization. After pre-processing, a normalized expression matrix was constructed. The group matrix was constructed based on clinical information. All bioinformatic analyses and visualization were processed based on R.

The WGCNA package (Langfelder and Horvath, 2008) was used to create a gene co-expression network. By median absolute deviation (MAD), the top 5,000 ranking genes were selected at first. Then a soft-thresholding power β was calculated by using the “pickSoftThreshold” function. A suitable power value was defined as the first number reaching which the degree of independence was at least 0.9. The gene expression matrix was then converted into a topological overlap matrix (TOM), and the genes were divided into several gene modules, each represented by a distinct color. Next, a hierarchical clustering analysis was performed by using the hclust function. Except for the WGCNA package, the gplots package (Warnes et al., 2020) was used for visualization. In addition, the top 100 networks sorted by weight were exported to Cytoscape software for visualization.

In WGCNA, gene significance (GS) was used to describe the relationship between gene and phenotype. Module membership (MM) was calculated to evaluate the importance of a gene in the module by using the cor function. In this study, genes with both GS > 0.3 and MM > 0.9 was defined as hub genes among the candidate gene modules (Jin et al., 2021). The correlation relationship of hub genes was explored by using the gpairs package (Emerson and Green, 2020).

We first used lmFit and eBayes functions in the limma package (Ritchie et al., 2015) to identify the DEGs between VaD and control groups. The statistical method to calculate false discovery rate (FDR) was the Benjamini–Hochberg method. Then a threshold of adjust-p < 0.05 and the absolute value of log₂ fold change (log₂FC) > 1 were set, and the significant DEGs between the VaD and controls were filtered out. A volcano plot was presented by using EnhancedVolcano (Blighe et al., 2018). The distribution shape of TLR2 was shown in the violin plot by using the ggpubr package (Kassambara, 2020).

A GO enrichment analysis was run to annotate the functions of the significant DEGs with GO terms. The GO enrichment analysis could explain the features of changed genes from the following three structural networks of terms: biological processes (BP), cellular components (CC), and molecular functions (MF). The KEGG pathway analysis was performed to investigate the pathway that the significant DEGs might be involved in. The org. Hs.eg.db package (Carlson, 2021) was used for transition from gene symbols to Entrez ID. Then the clusterProfiler package (Yu et al., 2012; Wu et al., 2021) was used for the enrichment analysis. At last, the ggplot2 (Wickham, 2016) package was used for visualization. The aforementioned analysis results enabled us to discover the biological pathways of the altered genes in the VaD group.

In the GSE122063 datasets, GSEA was used to explore distinct GO terms and KEGG pathways that may be associated with TLR2. All genes were included in the analysis. Gene sets were directly downloaded from the website (http://www.gsea-msigdb.org/gsea/downloads.jsp). Except for the VaD and control groups, we set the median expression level of TLR2 as the cutoff value to divide patients into TLR2-high and TLR2-low expression groups. The org. Hs.eg.db package (Carlson, 2021) was used for Entrez ID transition, and the clusterProfiler package (Yu et al., 2012; Wu et al., 2021) was used for the enrichment analysis. Furthermore, the gseaplot2 function in the enrichplot package (Yu, 2021) was used for visualization of enrichment results.

We used the STRING online database (https://string-db.org/) to construct a PPI network. Significant DEGs were uploaded to the STRING website. After being filtered by the “no more than 50 interactors” and “k-means clustering” options, the PPI network was exported into a TSV file. At last, the analysis and visualization of the interaction network were achieved by Cytoscape software. The function of network analysis function in the Cytoscape software calculated the degree which was utilized as the continuous mapping of nodes both in size and fill color (from blue to red). The combined score exported directly from the string database was used for the continuous mapping of edges both in width and stroke color (from blue to red). Larger size and bluer nodes indicated the higher degree, while wider and bluer lines indicated the higher combined scores.

Using the expression matrix, WGCNA was used to determine the main module which was most linked with VaD. At first, we chose the top 5,000 genes sorted by MAD in the GSE122063 microarray assay for analysis. According to the calculation result, the soft-thresholding power β was 2 as the plot showed, with the scale-free topology R2-value achieving 0.9 (Figure 1A). To visualize the weighted network, a heat map was plotted. The gene co-expression network was created, and the genes were clarified into five modules represented by distinct colors including grey, turquoise, blue, brown, and yellow. This is called a cluster dendrogram, and it was presented along the axis. a network heat map of all 5,000 genes was shown by using the TOMplot function in Figure 1B. Each row and column in the heat map represented the same gene, and thus the network heat map is a symmetric plot. The genes with strong correlations were clustered into modules, which were represented as dark sections symmetrically distributed along the diagonal in the heat map, corresponding to the cluster dendrogram. The biggest grey module included 2,783 genes, and the smallest yellow module included 400 genes. As shown in the module–trait relationships plot, the yellow module was most positively associated with VaD (correlation coefficient = 0.57, ***p < 0.001; Figure 1C) and was chosen as the key module. The functional annotation of three significantly related modules (blue, turquoise, and yellow) are shown in Supplementary Figure S1. The yellow module was most related to immunity and inflammation.

Among the 400 genes in the yellow module, genes with MM > 0.9 and GS > 0.3 were sorted out as hub genes. The red dotted lines represent the thresholds value of MM > 0.9 and GS > 0.3 set for hub genes and separated an area in the upper right corner. The correlation analysis between yellow module memberships and gene significance showed statistical significance (correlation coefficient = 0.65, ***p < 0.001). In total, 21 hub genes were identified (TLR2, CD163, VSIG4, SLAMF8, C1QB, CD16a, CD32, ALOX5AP, integrinβ2, EBI3, HCLS1, CD14, LAIR-1, CD300a, IFI30, LCP1, C1orf162, γ-parvin, ALOX5, SLA, and CMTM7). According to MM or the chooseTopHubInEachModule function, TLR2 was the top hub gene in the yellow module (Figure 2). Furthermore, we found that TLR2 shows a strong positive correlation with other candidate genes, which indicated that changes in TLR2 expression might cause changes in these genes (Supplementary Figure S2).

The gene expression levels of the samples were distributed at the same baseline after normalization. Compared to the control group, significant DEGs were identified in the VaD group by setting the threshold value as adjust-p < 0.05 and |log2FC| > 1. The expression of the genes was displayed as a volcano plot in which the size of the dot reflects |log2FC| of the gene (Figure 3A). There were 456 significant DEGs between the VaD and control groups among the 23,320 genes detected in microarray, including 198 upregulated ones and 258 downregulated ones. TLR2 was one of the significant DEGs and was marked out in the volcano plot. Specifically, the expression level of TLR2 in the VaD and control groups was shown in the violin plot (***p < 0.001, Figure 3B). TLR2 was significantly differentially expressed between the two groups.

Significantly upregulated and downregulated DEGs were enriched in BP, CC, and MF terms and the KEGG pathway, respectively. The horizontal axis represents −log10 (p-value), while the color indicated the change direction. In detail, BP, Toll-like receptor signaling pathway was enriched, which was consistent with our previous result. Other BPs such as negative regulation of immune system process, antigen processing, and presentation and regulation of B, T, and NK cells were examples of significantly enriched upregulated GO terms (*p < 0.05, Figure 4A), while CCs, including azurophil granule, endocytic vesicle, and secretory granule membrane are shown (*p < 0.05, Figure 4B). Upregulated MFs, such as scavenger receptor activity and RAGE receptor activity, were significantly enriched. Neuropeptide hormone activity, neuropeptide receptor binding, and signaling receptor activation activity were downregulated (*p < 0.05, Figure 4C). Most enriched KEGG pathways did not reach statistical significance in which we observed a trend in Toll-like receptor signaling pathway and neuroactive ligand–receptor interaction was significantly downregulated (Supplementary Figure S3).

GSEA was analyzed in the disease group versus control as well as groups divided by the expression level of TLR2. When comparing the VaD group with the control group, the Toll-like receptor pathway was enriched, which was the same as the results from DEGs. Other immunity and inflammation-related processes were also enriched which indicated the representativeness of the data and complemented evidence for the role of TLR2 in neuroinflammation. The results are shown in Supplementary Figure S4. When comparing the TLR2-high group with the low group, the results showed that BPs, such as cytokine-mediated signaling pathway and defense response to other organism, were significantly enriched in the TLR2-high samples (*p < 0.05, Figure 5A). CCs, such as synapse, vacuole, and cell surface granule, and MFs, such as immune receptor activity and molecular transducer activity, were significantly enriched in the TLR2-high samples, shown in Figures 5B, C, respectively (*p < 0.05, Figures 5B,C). When it comes to the KEGG enrichment analysis, pathways such as antigen processing and presentation, ribosome, and cytokine–cytokine receptor reaction were significantly over-represented in TLR2-high samples (*p < 0.05, Figure 5D). The similar enrichment results in VaD and control groups, as well as in the TLR2-high and low groups further demonstrated the important role of TLR2 in VaD. Moreover, high expression level of TLR2 was related to many genes, including myeloid differentiation factor 88 (MyD88), nuclear factor kappa B (NF-κB), protein kinase B (AKT), glial fibrillary acidic protein (GFAP), ionized calcium-binding adapter molecule 1 (Iba1), and many cytokines according to the expression matrix and the KEGG pathway.

With the combined use of STRING and Cytoscape, the PPI network of the significant DEGs was created. The size and color reflected the degree of nodes in which the more edges connected to this node, the greater its degree. The larger size and bluer node indicated the higher degree. The width and color reflected the combined score of edges in which the combined scores were positively related to the interaction relationships between the two proteins. The wider and bluer line indicated higher combined scores. The overall network of DEG-correlated proteins is shown in Figure 6A. TLR2 got a relatively high degree in this overall network which suggested that TLR2 played a crucial role in the network. Considering the complication of the network, a new network centered on TLR2 was further constructed and amplified. TLR2 was most associated with Complement C5a Receptor 1 (C5AR1), Heat Shock Protein Family A Member 1 A (HSPA1A), cluster of differentiation (CD14), and cytochrome B-245 Beta Chain (CYBB) (Figure 6B).