Shared gene characteristics and molecular mechanisms of macrophages M1 polarization in calcified aortic valve disease

Background CAVD is a common cardiovascular disease, but currently there is no drug treatment. Therefore, it is urgent to find new and effective drug therapeutic targets. Recent evidence has shown that the infiltration of M1 macrophages increased in the calcified aortic valve tissues, but the mechanism has not been fully elucidated. The purpose of this study was to explore the shared gene characteristics and molecular mechanisms of macrophages M1 polarization in CAVD, in order to provide a theoretical basis for new drugs of CAVD. Methods The mRNA datasets of CAVD and M1 polarization were downloaded from Gene Expression Omnibus (GEO) database. R language, String, and Cytoscape were used to analyze the functions and pathways of DEGs and feature genes. Immunohistochemical staining and Western Blot were performed to verify the selected hub genes. Results CCR7 and GZMB were two genes appeared together in hub genes of M1-polarized and CAVD datasets that might be involved in the process of CAVD and macrophages M1 polarization. CCR7 and CD86 were significantly increased, while CD163 was significantly decreased in the calcified aortic valve tissues. The infiltration of M1 macrophages was increased, on the contrary, the infiltration of M2 macrophages was decreased in the calcified aortic valve tissues. Conclusion This study reveals the shared gene characteristics and molecular mechanisms of CAVD and macrophages M1 polarization. The hub genes and pathways we found may provide new ideas for the mechanisms underlying the occurrence of M1 polarization during CAVD process.


Introduction
As the most prevalent valvular heart disease, calcific aortic valve disease (CAVD) is a major health problem with risk of severe morbidity and mortality (1). It ranges from valve thickening and mildly calcified to severe valve calcification with impaired leaflet motion and vast blood flow obstruction (2). In developed countries, calcific aortic stenosis is the second-most frequent cardiovascular disease after coronary artery disease and systemic arterial hypertension with a prevalence of 0.4% in the general population (3) and 1.7% in the population over 65 years (4). Currently, there is no effective medical treatment other than surgical or transcatheter aortic valve replacement (5,6).
Currently, immune cells are also believed to play important roles in aortic valve calcification (12)(13)(14). Zhou's study (15) showed that immune cells accounted for about 5.4% in the aortic valve tissue, and macrophages were the most dominant immune cells. Another study showed that (16) in the calcified aortic valves, the infiltration of M1 macrophages was increased, while the infiltration of M2 macrophages was decreased. Therefore, M1 macrophages may play an important role in the CAVD process (17).
At present, there are few studies on the mechanism of macrophages M1 polarization in CAVD. However, exploring the mechanism of macrophages M1 polarization in CAVD may help to prevent the occurrence of CAVD or delay the progression of CAVD. Microarray technology can effectively screen the biomarkers of CAVD. Gene Expression Omnibus (GEO) 1 is a public database, and contains a large number of genes for various diseases. In this study, many bioinformatics analysis methods and microarray technology in GEO database were combined to screen and analyze the shared gene characteristics and molecular mechanism of macrophages M1 polarization in CAVD from the transcriptome level. This will be beneficial to deepen our understanding of CAVD and promote the research of therapeutic drugs for CAVD.

Acquisition of GEO datasets
We use the keywords "calcified aortic valve disease" and "M1 macrophages" in GEO database (see text footnote 1) to search the mRNA datasets of CAVD and M1-polarized macrophages. The datasets we chose included mRNA microarray data of CAVD samples and healthy controls, as well as mRNA microarray data of M1 and M0 macrophages.
The GEO datasets were numbered as:

Microarray data processing
In this study, microarray datasets GSE49240 and GSE61298 were first used to identify differentially expressed genes (DEGs) during M1 polarization. The "limma" package in R language (version 4.1.3) was used to screen out DEGs in the process of M1 polarization. The cut-off values were |log2(foldchange)| > 1and p-value < 0.05. Then two other microarray datasets (GSE51472, GSE83453) were used to identify the feature modules in CAVD by WGCNA analysis. The overlapping DEGs between two M1 macrophage datasets and feature genes between two CAVD datasets were implemented by the "VennDiagram" package.

WGCNA analysis
Weighted gene co-expression network analysis (WGCNA) was an analytical method for analyzing gene expression patterns in multiple samples. It could cluster modules with similar expression and analyze the correlation between specific modules and disease phenotypes. In this paper, we used WGCNA analysis to obtain modules that were highly related to CAVD. We selected the top 5,000 genes according to their variance and analyzed them by the "WGCNA" package in R language. Before analysis, the "Hclust" function was first used for hierarchical cluster analysis to exclude the outlier samples. Then the "pickSoftThreshold" function in "WGCNA" package was used to select an appropriate power value (ranging from 1 to 30, with R 2 ≥ 0.85) according to the criteria of the scalefree network. Then we constructed a hierarchical clustering dendrogram, clustered similar gene expressions into different modules, and calculated the correlation between each module and the disease phenotype. In this study, the soft threshold β of dataset GSE51472 was 19. The soft threshold β of dataset GSE83453 was 13. The other parameters were as followed: networkType = "unsigned, " minModuleSize = 30, mergeCutHeight = 0.25 and deepSplit = 2.

Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed for Gene function and pathway analysis. GO enrichment analysis analyzed protein function at three levels: biological process, cellular component, and molecular function. GO and KEGG enrichment analysis were implemented by the "clusterProfiler" package. ClueGO enrichment analysis was a powerful complement to GO and KEGG enrichment analysis and was implemented by the ClueGO plug-in in Cytoscape (version: 3.7.1).

Protein-protein interaction network and hub genes screening
The construction of the protein-protein interaction (PPI) network was implemented by the online analysis tool String 2 and Cytoscape. The cytoHubba plug-in in Cytoscape was used to evaluate the scores of gene nodes and screened out the top 20 hub genes (Top 20 genes were ranked by MCC method).

Cibersort immune infiltration analysis
Cibersort (cell-type identification by estimating relative subsets of RNA transcripts) was an immune infiltration algorithm based on transcriptome. It could calculate the relative content (LM22 gene signature was downloaded from https://www.nature.com/articles/nmeth.3337#MOESM 207) and immune score (perm replacement = 1000, QN quantile normalization = TRUE) of 22 kinds of immune cells by using microarray or high-throughput sequencing data. On this basis, the Spearman correlation between hub genes and immune cells were calculated by the "ggcorrplot" package.

Analysis of single cell sequencing results
The single cell sequencing data in this paper was downloaded from Xu's article (18) 3 and included 9,410 CAVD cells and 3,366 normal aortic valve cells. The raw expression matrix was integrated into Seurat objects and filtered by the "Seurat" package. Then the "NormalizeData" function was used for normalization and the "UMAP" package was used for dimensionality reduction (FindClusters functions with a resolution of 0.5). Finally, cell clusters were annotated by the "singleR" package.

Clinical samples and histological examination
In this study, we collected aortic valve specimens from six patients who underwent aortic valve replacement or heart transplantation at the Department of Cardiovascular Surgery, Changhai Hospital Affiliated to Naval Medical University from April 2022 to July 2022. The study has been approved by the Ethics Review Committee of Changhai Hospital. Among the six patients, three patients who underwent aortic valve replacement were diagnosed with calcified aortic valve, and another three patients who underwent heart transplantation were diagnosed with dilated cardiomyopathy with normal aortic valve, which were identified by echocardiography and Alizarin red S staining. Patients with rheumatic heart disease, myocarditis, and bicuspid aortic valve were excluded. Part of the specimens were embedded in paraffin and the thickness of the paraffin section was 5 µ m.
Alizarin red S staining: Paraffin sections were dewaxed and rehydrate routinely, Alizarin red S staining solution (1% dilution, pH = 4.2) was incubated at room temperature for 5 min, and then the sections were washed with anhydrous alcohol and sealed before observation.
Von Kassa (VK) staining: After dewaxing and hydration, paraffin sections were stained with silver nitrate, irradiated continuously with UV lamp for 1 h and washed with distilled water three times, and then sealed for observation.
Immunohistochemical (IHC) staining: The paraffin sections were dewaxed and hydrated before heat mediated antigen retrieval (Tris-EDTA buffer, pH = 9.0). Endogenous peroxidase activity was blocked by peroxidase blocking solution and non-specific binding sites were blocked by goat serum. Then the slides were incubated with the primary antibody overnight (4 • C). The next day, biotin-labeled secondary antibody was added. Finally, DAB solution and hematoxylin were used for staining and counterstaining, respectively.

Western blot analysis
Aortics valve specimens were lysed in SDS buffer containing a protease inhibitor PMSF (1:100 dilution) on ice for 30 min. Total protein concentrations were evaluated using a protein assay kit (Beyotime Biotech, Shanghai, China). Cell lysate was separated using 10% SDS-PAGE and then transferred to polyvinylidene difluoride membranes. Subsequently, the membranes were blocked for 1 h at room temperature by      the top 50 DEGs ranked by |Log 2 fold change| were shown by heatmap in Figures 1C, D (mRNA expression levels were processed by Z-score). The results showed that IDO1 was the largest increased DEGs and FABP4 was the largest decreased DEGs in both datasets. In addition, the overlapping DEGs in the two M1-polarized datasets were visualized by Venn diagram (Figure 1E), and 277 overlapping DEGs were screened (including 157 up-regulated and 120 down-regulated mRNAs, shown in Supplementary Material 1).

Enrichment analyses of the overlapping DEGs in M1-polarized datasets
The overlapping DEGs were analyzed by GO, KEGG, and ClueGO enrichment analyses. GO enrichment analysis (Figure 2A) showed that DEGs were significantly associated with the biological processes of "Positive regulation of cytokine production, " "Regulation of response to biotic stimulus, " "Immune response-regulating signaling pathway, " "Response to virus, " and "Activation of Immune response." KEGG pathway analysis ( Figure 2B) showed that the overlapping DEGs were mainly enriched in "NOD-like receptor signaling pathway, " "Phagosome, " "Influenza A, " "PPAR Signaling Pathway, " and "Viral protein interaction with cytokine and cytokine receptor." ClueGO enrichment analysis (Figure 2C) showed that the overlapping DEGs were mainly involved in "Interleukin-1 beta production, " "Type I interferon signaling pathway, " and "Response to interferon gamma." The signaling pathway "Interleukin-1 beta production" accounted for 27.27% of the total terms, and 82 genes were involved ( Figure 2D). These results suggested that chemokine-chemokine receptor pathway and inflammatory response may play important roles in macrophage M1 polarization.

PPI analysis of the overlapping DEGs in M1-polarized datasets
We further constructed a PPI network at the level of protein and screened out hub genes. Firstly, the online analysis tool String was used to analyze the interactions among 277 overlapping DEGs (minimum required interaction score = 0.4, 276 nodes and 1,107 edges were included in the PPI network) ( Figure 3A). Then the cytoHubba plug-in in Cytoscape was used to calculate the score of gene nodes and screen out the top 20 hub genes. The results (shown in Supplementary Material 2) showed that CXCL10, TNF, TLR4, CXCL8, TLR7, STAT1, CCL5, IL15, CXCL9, CD40, CD274, ICAM1, GZMB, CASP1, IRF1, CCL20, TNFSF10, CCRL2, IDO1, and CCR7 were identified as hub genes (Figure 3B).

Identification of feature genes in CAVD datasets
A total of five modules were identified by WGCNA analysis in CAVD mRNA dataset GSE51472. Then, a heatmap of module-phenotype relationship (Figure 4A) was drawn according to spearman correlation coefficient to evaluate the association between each module and disease. Two modules were selected as CAVD-related module (blue module: r = 0.68, p = 0.005, containing 1,263 genes; Brown module: r = 0.74, p = 0.002, containing 202 genes). Five modules were also identified in CAVD mRNA dataset GSE83453 (Figure 4B). Only the blue module was selected as CAVD-related module (blue module: r = 0.89, p = 2e −06 , containing 1,195 genes).
CAVD-related modules in both datasets were shown in Supplementary Material 3. The CAVD-related modules in the two CAVD datasets were overlapped by Venn diagram. The results showed that there were 279 overlapping upregulated genes between the two CAVD datasets, which may be associated with the pathogenesis of CAVD (Figure 4C, shown in Supplementary Material 3).

Enrichment analyses of feature genes in CAVD datasets
The feature genes in two CAVD datasets were analyzed by GO, KEGG, and ClueGO enrichment analyses. GO enrichment analysis ( Figure 5A) showed that the feature genes were significantly associated with the biological process of "Leukocyte cell-cell adhesion, " "Regulation of cell-cell adhesion, " "Positive regulation of cell activation, " "Leukocyte migration, " and "Positive regulation of leukocyte activation." KEGG pathway analysis ( Figure 5B) showed that the feature genes were mainly enriched in "Cell adhesion molecules, " "Natural killer cell mediated cytotoxicity, " "Th1 and Th2 cell differentiation, " "Hematopoietic cell lineage, " and "T cell receptor signaling pathway." ClueGO enrichment analysis (Figure 5C) showed that the feature genes were mainly involved in "Positive regulation of leukocyte cell-cell adhesion, " "Myeloid leukocyte migration, " and "T cell differentiation." The signaling pathway "Positive regulation of leukocyte cell-cell adhesion" accounted for 25.4% of the total terms, and 479 genes were involved ( Figure 5D). These results suggested that the cell-cell adhesion pathway may play an important role in CAVD.

PPI analysis of feature genes in CAVD datasets
We further constructed a PPI network at the level of protein and screened out hub genes. Firstly, String was used to analyze the interactions among 279 feature genes (minimum required interaction score = 0.4, 275 nodes and 1,705 edges were included in the PPI network) (Figure 6A). Then the cytoHubba plug-in was used to calculate the score of gene nodes and screen out the top 20 hub genes. The results (shown in Supplementary Material 4) showed that CD86, ITGAM, CD8A, CD2, GZMB, CD28, CD27, IL7R, ITGAX, IL10RA, IL2RB, SLAMF1, TNFRSF4, CCR7, CXCR3, ITGAL, CD48, KLRB1, CD247, and ITGB2 were identified as hub genes (Figure 6B).

Identification and analysis of shared genes between M1-polarized and CAVD datasets 3.3.1. Identification of shared genes between M1-polarized and CAVD datasets
The overlapping DEGs in M1-polarized datasets were overlapped with the feature genes in CAVD datasets. The results showed that CCR7 and GZMB appeared simultaneously (Figure 7A), and both were hub genes in the CAVD datasets and the M1-polarized datasets (Figure 7B). These results suggested that CCR7 and GZMB may play important roles in M1 polarization and CAVD.

Immune infiltration analysis of CAVD datasets
Next, two CAVD datasets were analyzed for immune infiltration. The results showed that compared with the normal aortic valves, the proportion of M0 and M1 macrophages were increased, but the proportion of M2 macrophages was decreased in the calcified aortic valves (Supplementary Figures 1, 2). By calculating the correlation between hub genes and immune cells in CAVD datasets, this study found that CCR7 was positively correlated with M1 macrophages, while GZMB was negatively correlated with M1 macrophages (Supplementary  Figures 3, 4). Therefore, we focused on the relationship between CCR7 and M1 macrophages.

Distribution of CCR7 in aortic valve tissue cells
To investigate cellular heterogeneity in aortic valves microenvironment at single-cell resolution, aortic valve tissue cells were divided into 11 clusters (Figure 8A, FindClusters functions with a resolution of 0.5) by using the "UMAP" package. The annotation of single cell sequencing (Figure 8B) showed that the aortic valve tissue included four types of cells,  Representative Alizarin red S staining and VK staining images. other clusters. Then the distribution of M1 and M2 macrophages were marked by CD86 and CD163, respectively. The results showed that the distribution of CCR7 was consistent with that of CD86, but not consistent with that of CD163 (Figures 8D-G). These results indicated that CCR7 was associated with M1 macrophages.

Experimental verification of clinical specimens 3.4.1. Pathological verification of calcification of clinical specimens
Aortic valves from patients with calcified or normal aortic valves were selected and sectioned. Calcified nodules were identified by Alizarin red S staining and VK staining in calcified aortic valve samples. Representative Alizarin red S staining and VK staining images were shown as follows (Figures 9A, B).

IHC staining and western blot verification of clinical specimens
Finally, the calcified and normal aortic valve tissues were sectioned for CCR7, CD86, and CD163 IHC staining. The IHC staining results showed the expression level of CCR7 and CD86 were significantly increased in the calcified aortic valve tissues compared to the normal aortic valve. On contrary, the expression of CD163 was significantly decreased as shown in Figures 10A-C. We further investigated the expression profile of CCR7 in calcified and normal aortic valves tissues. Western blot assay confirmed that the expression level of CCR7 was increased in calcified aortic valve tissues (Supplementary Figure 5).

Discussion
In this study, we screened out the DEGs in M1 polarization. Enrichment analyses showed that chemokine-chemokine receptor pathway and inflammatory response may play important roles in the process of M1 polarization. Then, the feature genes of CAVD were screened out. Enrichment analyses indicated that cell-cell adhesion pathway may play an important role in the process of CAVD. CCR7 and GZMB appeared in hub genes of M1-polarized and CAVD datasets. These results suggested that CCR7 and GZMB may regulate CAVD by regulating the polarization, chemotaxis, and adhesion of M1 macrophages. Cibersort immune infiltration analysis showed that the infiltration of M1 macrophages was increased and the infiltration of M2 macrophages was decreased during the CAVD process. Single cell sequencing results showed that the distribution of CCR7 was consistent with the distribution of CD86, but not with CD163, suggesting that CCR7 was correlated with M1 macrophages, which may affect the process of CAVD. IHC staining results showed that, compared with the normal aortic valve tissues, the expressions of CCR7 and CD86 were significantly increased in the calcified aortic valve tissues, while the expression of CD163 was significantly decreased. These results indicated that the infiltration of M1 macrophages was increased and the infiltration of M2 macrophages was decreased. Western blot assay further confirmed that CCR7 was increased in calcified aortic valve tissues.
Recent studies have shown that macrophage play an important role in CAVD (19), the infiltration of M1 macrophages was increased and the infiltration of M2 macrophages was decreased in the calcified aortic valves (20). Raddatz (24).
The protein encoded by CCR7 is a member of the G protein-coupled receptor family (25)(26)(27) and is up-regulated in response to inflammation, pathogens, and tissue damage (28). CCR7 chemokine axis is composed of chemokine ligands CCL19 and CCL21, and chemokine receptor CCR7 (29,30). At the cellular level, CCR7-mediated signaling pathways mainly include the activation of integrins leading to cell adhesion and the polarization of actin cytoskeleton (31). Activation of CCR7 can activate T cells, secrete a large amount of IFNγ, and promote Th1 polarization (32). The deficiency of CCR7 leads to Th2 polarization and activation of B cells (33). CCL19 and CCL21 only induce M1 polarization, and CCL19 or CCL21 induced activation of both MEK1-ERK1/2 and PI3K-AKT cascades in M1 but not in M2 macrophages (27). In rheumatoid arthritis, it has been reported that the activation of CCL19/CCL21-CCR7 signaling pathway increases the infiltration of M1 macrophages (34, 35). It has been shown that patients with symptomatic aortic stenosis have higher serum CCL21 level compared to healthy controls (36). Another study reported that O 2 consumption and CO 2 production were higher in CCR7 −/− mice than in wild-type mice. CCR7 −/− mice are protected from diet-induced obesity and subsequent insulin resistance (37), which may protect CCR7 −/− mice from CAVD. This study hypothesized that the CCR7 axis may promote the polarization and adhesion of macrophages, thereby promoting the progression of CAVD.
GZMB encodes a member of the granzyme subfamily of proteins, granzyme B. GZMB-induced apoptotic cells, when phagocytosed by macrophages, can induce the production of TGF-β and thereby influence the Th1/Th2 cytokine and Ig balance (38). Another research reported that GZMB mediates IL-8/macrophage inflammatory protein-2 secretion (39). Further, Guo et al.'s study also reported that GZMB was differentially expressed gene in CAVD datasets (40), which supports our analysis.
To sum up, CCR7 may act as a potential target for the treatment of CAVD. This study still has many limitations. At first, the increased expressions of CCR7 and M1 macrophages in calcified aortic valve tissues were verified only by bioinformatics analyses and clinical specimens in this study. The function of CCR7 needs to be explored by a more comprehensive design, both in vitro and in vivo, to determine its in-depth mechanism leading to CAVD. Secondly, the sample size of the datasets in this paper are too small, which may reduce the reliability of the analysis results. Analysis of datasets with larger sample size and more detailed description of patient characteristics will facilitate the understanding of the results. At last, our grouping of the samples in the dataset GSE51472 is not optimal. Putting mildly calcified and severely calcified samples together may result in false negative results. We will pay more attention to the quality of our datasets and achieve homogeneity in the same group of samples in the future.

Conclusion
In this study, bioinformatics analyses suggested that CCR7 and GZMB might be involved in the process of CAVD. In addition, the infiltration of M1 macrophages was increased and the infiltration of M2 macrophages was decreased during the CAVD process, and the distribution of CCR7 was consistent with that of M1 macrophages. These results suggest that CCR7 may regulate the process of CAVD by regulating macrophage polarization and adhesion.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement
The studies involving human participants were reviewed and approved by the Ethics Review Committee of Changhai Hospital. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions
ZX and GW designed the study. MQ and NL analyzed the data and finished the manuscript. QC, CW, and XX completed