Transcriptomics Analysis Revealed Key Genes Associated with Macrophage Autophagolysosome in Male ApoE^−/− Mice Aortic Atherosclerosis

Introduction

Atherosclerosis (AS) corresponds to the formation of fibrofatty lesions in the arterial wall and is the pathological basis of myocardial infarction, stroke, and disabling peripheral arterial diseases.¹ The AS plaque formation process involves endothelial cell injury, oxygen free radical damage, lipid aggregation, foam cell formation, platelet adhesion aggregation, inflammatory cell migration and infiltration, vascular smooth muscle cell calcium overload, pathological angiogenesis, and local thrombus formation; however, its pathogenic mechanism has not been fully elucidated.² Recent studies have shown that the autophagolysosomal system plays an important role in macrophage lipid metabolism in AS plaques.^3,4 Macrophages play a central role in the progression of atherosclerosis. They develop progressive autophagic dysfunction in the developing plaques. Autophagy dysfunction triggering involves a critical factor that appears to affect the functional st atus of macrophage lysosomes.^5,6 Oxidized low-density lipoprotein (LDL) and cholesterol oxidation products (ChOx) can damage the macrophage lysosomal membrane, leading to the inactivation or relocalization of some lysosomal enzymes, which in turn leads to the leakage of lysosomal contents and macrophage death. Lipid-induced lysosomal dysfunction in plaque macrophages is an essential factor in autophagy defects.⁷ Stimulation of autophagy-lysosomal biogenesis can reduce atherosclerosis.

High-throughput gene mapping has become an effective tool to reveal the pathogenesis of cardiovascular diseases, which provides a new direction for the study of cardiovascular diseases in the future. But at present, most of these studies have focused on comparing the differential gene expression between human atherosclerotic plaques and normal tissues. ApoE^−/− mice readily form lesions similar to human atherosclerosis, even without any dietary challenge, and are currently the most commonly used mouse model for studying atherosclerosis, being successfully used to study the occurrence and progression of atherosclerosis and for drug discovery.⁸ In this study, we used bioinformatics methods to identify the core genes involved in the progression of atherosclerosis, hoping that our findings can provide a new therapeutic target for atherosclerosis.

Materials and Methods

All data and materials are publicly available at the NCBI Gene Expression Omnibus (GEO). The datasets can be accessed at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE232472.

Animal Experiments

The procedures involving animals and their care were conducted in accordance with institutional guidelines in compliance with national policies [Regulations on the Administration of Experimental Animals, Decree No. 676 of the State Council of the People’s Republic of China (2017); Guidelines for Ethical Review of Experimental Animal Welfare, GB/T35892-2018] and followed the principles of ethical animal research outlined in the Basel Declaration and the ethical guidelines of the International Council for Laboratory Animal Science. This study was approved by the Ethics Committee of Shandong Provincial Qianfoshan Hospital (Jinan, China) [2018; ethical approval number: S0065].

Feeding and Grouping

Six-week-old (19 ± 2 g) ApoE^−/− mice were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China). To prevent the possible impact of the hormonal changes of female mice on the results, only male mice were employed in the study. After two weeks of adaptive feeding, the 8-week-old male mice were divided into four groups (n = 12), two of which were fed a Western diet containing 21% anhydrous milk fat and 0.15% cholesterol (Keao Xieli Feed Co., Ltd., Beijing, China) for 6 and 22 weeks, respectively, and which were used as atherosclerosis models (AS6W and AS22W), and the other two groups were fed a normal purified diet (AIN93M, Keao Xieli Feed Co., Ltd., Beijing, China) for 6 and 22 weeks, respectively (Sham6W and Sham22W). All mice were housed in the SPF experimental animal center of Shandong Provincial Qianfoshan Hospital at a temperature of 22 ± 1°C, an ambient humidity of 50% ± 5%, a 12-hour light/dark cycle, and received food and water at will.

Ultrasound Examination

One week before the end of the dietary treatment, after an overnight fast, the mice were anesthetized using 2% isoflurane (RWD, Shenzhen, China). The left carotid artery of each ApoE^−/− mouse was scanned from the proximal end to the bifurcation, and then longitudinally scanned along the common carotid artery after transverse scanning using a small animal ultrasound machine, VISUALSONICSVevo3100 (FUJIFILM Japan), a 55-MHz scanning probe, and a 4.5 mm focus, and the B-ultrasound, M-ultrasound, and PW images were preserved. The stored images were imported into the VevoLAB3.2.6 system to measure the maximum carotid systolic flow velocity (PSV), end-diastolic velocity, velocity time integral, and vascular resistance index (RI) (Figure 1A). All the measurements were performed twice at the same site.

Figure 1 Carotid artery ultrasonography and aorta sampling in ApoE−/− mice. (A) B-ultrasound and Doppler images of the left carotid artery of mice were collected and imported into the VevoLAB3.2.6 system to calculate the vascular resistance index (RI). The red arrow points to the left carotid artery of the ApoE−/− mice under ultrasound. (B) Carotid artery ultrasonography showed that RI in AS22W was significantly higher than that in the other three groups (P<0.0001) (ns p≥0.05, ****p<0.0001). (C) Stripping the hearts and aortas. The red arrow points to the aortic atheromatous plaque of the ApoE−/− mice in AS22W group under direct vision.

Blood and Tissue Harvesting

At the end of the dietary treatments, after an overnight fast, the mice were euthanized under general anesthesia using 2% isoflurane (RWD, Shenzhen, China). Blood was collected from the retro-orbital plexus into tubes containing 0.1% (w/v) EDTA and centrifuged in a microcentrifuge for 20 min at 4°C at 2000 rpm. Plasma was stored at −80°C until the analyses were performed. The chest of each mouse was opened quickly and blood was removed via perfusion using 1× PBS. The aorta was rapidly detached from the aortic root to the iliac bifurcation and the periadventitial fat was removed (Figure 1C). Then, each aorta was snap-frozen in liquid nitrogen for the RNA-seq (n = 4) and qRT-PCR (n = 6), or externally fixated with 4% paraformaldehyde (n = 2) for oil red O staining. The hearts with attached aortic roots were obtained immediately, fixed overnight in 4% paraformaldehyde, embedded in optimum cutting temperature (OCT) compound, and then cut into serial 6-µm cross sections for hematoxylin and eosin, oil red O, Masson (n = 6), and immunofluorescence staining (n = 3). Aortic sinus histological analysis was performed according to the American Heart Association recommendations.⁹

Plasma Lipid Measurements

The plasma total cholesterol (TC), triglyceride (TG), and low-density lipoprotein cholesterol (LDL-C) levels were measured using GPO-PAP enzymatic methods (A110-1-1, A111-1-1, A113-1-1, Nanjing Jiancheng Bioengineering Institute) and a multimode detection platform (SpectraMax i3x Molecular Devices, LLC., USA).

RNA-Seq Analyses

The mRNA quality was tested before performing RNA-seq. The sequencing platform used was an Illumina Novaseq 6000 in sequencing mode PE150 (LC-Bio Technology Co., Ltd., Hangzhou, China).

RNA-Seq Data Processing and Visualization

Data Processing of the Differentially Expressed Genes (DEGs)

Network analyst¹⁰ is an online visualization platform for gene expression analysis and meta-analysis, which can be used for comparison, quantitative analysis, gene expression difference analysis and enrichment analysis, protein interaction analysis, and integration analysis of multiple datasets, among others. Here, differential expression analysis was performed using Network analyst to determine the RNA-seq counts. The DEGs were grouped based on cut-off values of adj.P < 0.05 and |logFC| > 1.0. Then, R heatmap and ggplot2 packages were used to draw the heatmap, volcano, and PCA plots to visualize the DEGs.

Protein-Protein Interaction (PPI) Network and Hub Gene Identification

STRING database was used to construct the protein-protein interaction network of the DEGs between the AS22W and Sham22W groups. The maximum neighborhood component (MNC)¹¹ was then calculated using the CytoHubba plug-in, and the genes with the top nine MNC values were selected as the hub genes.

TRRUST¹² is a human-annotated transcriptional regulatory network database that includes 8444 and 6552 transcription factor target regulatory relationships for 800 human and 828 mouse transcription factors, respectively. TRRUST was used to predict the hub gene-related transcription regulators to construct mRNA-TF interaction networks. The StarBase database,¹³ in particular, the high-throughput CLIP-seq and degradation group experimental data, which provide a variety of visual interfaces, was searched for miRNA targets. This database contains abundant data on miRNA-ncRNAs, miRNA-mRNAs, RBP-RNAs, and RNA-RNAs. Here, we predicted the RBP combined with hub genes using the StarBase database. The parameters used were the Clade (mammal), Genome (mouse), Assembly (mm10), CLIP-Data (≥ 3), and pan-Cancer (≥ 0).¹⁴ CTD is a public database that links the toxicological information on chemicals, genes, phenotypes, diseases, and exposures. Through this platform, we have predicted chemical drugs related to the hub genes, while retaining the chemical drugs with ≥ 10 interactions. Finally, an interaction network diagram was created using the R package iGraph.

Enrichment Analysis

Gene Ontology (GO) analysis is a method commonly used for large-scale functional enrichment studies, including biological processes (BPs), molecular functions, and cellular components.^15,16 The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database that contains information on genomes, biological pathways, diseases, and drugs. The “clusterProfiler” R package¹⁷ was used to perform GO and KEGG enrichment analyses of the DEGs between the AS22W and Sham22W groups. A P value < 0.05 was considered significant. Moreover, to study the differences in BPs between different groups, gene set enrichment analysis (GSEA) against the DEGs between the AS22W and Sham22W groups was performed using GSEA Java.¹⁸ The gene set “c2.cp.v7.2. symbols.gmt [Curated]” was downloaded from MSigDB¹⁹ for GSEA, and a P value < 0.05 was considered significant. Furthermore, we used R-package gene set variation analysis (GSVA)²⁰ to calculate the scores of the related pathways based on the gene expression matrices of the AS22W and Sham22W samples and performed differential screening for the enrichment functions using the “limma” package.

Evaluation of Immune Cell Infiltration

We used CIBERSORT to evaluate the abundance of immune cells in the samples from the AS22W and Sham22W groups. Moreover, we analyzed the correlation between the hub genes and immune cells.

Verified by the Published Human Datasets

The dataset GSE100927²¹ was downloaded from GEO and used for transcriptome analysis of the atherosclerotic and non-atherosclerotic tissues of the carotid, femoral, and popliteal arteries. The dataset analysis platform used was GPL17077 from the Agilent-039494 SurePrint G3 Human GE v28×60K Microarray 039381 (probe name version). Here, 29 carotid atherosclerotic tissues and 12 normal carotid tissues were selected as the experimental and control groups, respectively, to test our RNA-seq dataset.

Real-Time Quantitative PCR

SteadyPure universal RNA extraction kit (Accurate Biotechnology Co., Ltd., Hunan, China) was used to extract the total RNA from the aortic samples and Hifair III 1st Strand cDNA Synthesis SuperMix for qPCR (Yeasen Biotechnology Co., Ltd., Shanghai, China) was used to reverse transcribe the total RNA into complementary DNA (cDNA). The polymerase chain reaction (PCR) was performed using Hieff qPCR SYBR Green Master Mix (No Rox)(Yeasen Biotechnology Co., Ltd., Shanghai, China) and specific primer pairs (sequences are listed in Table 1). Amplification was performed in duplicate using LightCycler^® 480 II (Roche). The experiment was repeated three times for each gene and β-actin was used as the endogenous control.

Table 1 Gene-Specific Primers Used for Real-Time RT-PCR

Immunofluorescence

Mice were sacrificed using 2% isoflurane (RWD Shenzhen, China), and the heart of each mouse was obtained immediately, fixed in 4% paraformaldehyde overnight, embedded in OCT compound, and serially cut into 6-µm-thick sections for immunofluorescence. Rabbit anti-ATP6V0D2 antibody (bs-12548R, Bioss, Beijing, China) was used for ATP6V0D2 staining. Anti-CD68 rabbit pAb (GB113109, Servicebio, Wuhan, China) was used for macrophage staining.

Statistical Analysis

Data analysis and plot drawing were performed using GraphPadPrism 8.0. The quantitative data were expressed as the mean ± standard deviation, and the normality test was carried out. The data of multiple groups were compared by single factor analysis of variance (ANOVA) (normal distribution) or KruskalWallis test (non-normal distribution). P value < 0.05 was considered significant.

Results

Verification of the Atherosclerosis ApoE^−/− Male Mouse Model

Ultrasonic Evaluation of Carotid Atherosclerosis

The most commonly used ultrasonic parameter to assess arteriosclerosis is the intima-media thickness (IMT), which shows part of the blood-intima and media-adventitia interfaces in the B-mode image program. However, its disadvantage is the large variability within and among observers, particularly for blood vessels with a small inner diameter. The RI indirectly reflects the degree of atherosclerosis through hemodynamic values, which are equivalent to the IMT. The RI is easier to obtain, with fewer inter- and intra-observer variations and side-to-side differences. The RI of the ApoE^−/− mouse carotid artery was significantly higher in the AS22W group than in the other three groups (P < 0.05), which is consistent with our histological results and indirectly confirmed atherosclerosis progression (Figure 1B).

Histopathological Examination of the Aorta

After 6 weeks of NLD, no atherosclerotic plaques were visible in the Sham6w group, whereas initial atherosclerosis lesions were observed in the Sham22W group; plaques size was comparable to that of plaque measured after 6 weeks of WD. As expected, 22 weeks of WD dramatically accelerated lesion progression and increased atherosclerosis in the AS22w group. Limited plaque progression was observed in ApoE^−/− mice, with no dietary influence and only a time challenge. The dietary challenge (WD) accelerated plaque development, with approximately two-times larger plaques being observed in the AS22W group than in the Sham22W group (Figures 2 and 3A–E).

Figure 2 Aortic gross oil red O staining.

Figure 3 Histological and immunohistochemical characterization of plaques in the aortic sinus. (A) Representative micrographs of oil red O staining, MASSON staining, HE staining and immunofluorescence double staining of CD68 (green particles) and α-SMA (red particles) in frozen sections of aortic sinus of ApoE-/- mice; Objective magnification 4×. (B and C) The plaque area ratio of the aortic sinus stained by oil red O and HE was calculated by Image J. (D and E) Quantification of extracellular matrix and macrophages, using ImageJ to calculate the percentage of positive area to total plaque area. (*p<0.05, ***p<0.001, ****p<0.0001).

Plasma Lipids

The total plasma cholesterol levels were affected by the dietary therapy. LDL-C, triglycerides, and cholesterol accumulated significantly more in the WD-fed AS group than in the ND-fed sham group (Figure 4A–C).

Figure 4 Blood Lipids. (A) the levels of TC; (B) the levels of TG; (C) the levels of LDL-C. Due to the effect of WD, the levels of LDL-C, TC, and TG in the AS22W group were significantly higher than those in the other three groups (** p<0.01,****p<0.0001).

RNA-Seq Analysis

Identification of the DEGs

There were 220 DEGs in AS6w group compared with Sham6w group, including 123 upregulated genes and 97 downregulated genes (Figure 5A). The expression heat map of top 5 high-expressed genes including Myb, Tenm1, Serpinb11, Bex1, Galnt9, and top 5 low-expressed genes including Hoxa10, Hoxc5, Atp4b, Hoxc9, Slc2a5 is shown in Figure 5B. As shown in Figure 5D, the volcanic map of difference analysis between the AS22w group and the Sham22W group showed 45 differentially expressed genes, 32 with high expression and 13 with low expression. The expression heat map of top 5 high-expressed genes including Fabp3, Mmp13, Insig1, Tg, Mmp12, and top 5 low-expressed genes including Ppbp, Nipal4, Tubb1, Fmo3, Saa3 is shown in Figure 5E. PCA figures are shown in Figure 5C and F respectively.

Figure 5 Difference analysis. (A) Difference analysis of volcano maps of AS6w and Sham6W groups. (B) Heat maps of top 5 differentially high/low expression genes in AS6w and Sham6W groups. (C) PCA of AS6w group and Sham6W group; (D) Difference analysis volcano maps of AS22w and Sham22W groups. (E) Heat maps of top 5 differentially high/low expression genes in AS22w and Sham22W groups. (F) PCA of AS22w and Sham22W groups.

PPI Network and Modular Analysis

Using the STRING database, we constructed a PPI network of the DEGs between the AS22W and Sham22W groups (Figure 6A). Screening of the hub genes using the CytoHubba plug-in in Cytoscape revealed that the top nine hub genes were Atp6v1d, Atp6v1a, Atp6v1h, Atp6v0c, Atp6v0d2, Atp6v0b, Atp6ap1, Atp6ap2, and Slc38a9 (Figure 6B). The expression of the hub genes matched between the AS6W and Sham6W groups (Figure 7A); Atp6ap2 was significantly overexpressed in the AS6w group (P < 0.05). The expression of hub genes matched between the AS22W and Sham22W groups (Figure 7B); Atp6Ap2, Atp6v0b, Atp6v0d2, Atp6v1a, and Atp6v1d were significantly overexpressed in the AS22w group (P < 0.05). In the GSE100927 dataset (patients with carotid atherosclerosis and normal carotid arteries), the expression of the hub genes was significantly higher in the carotid atherosclerosis group (P<0.001) (Figure 7C).

Figure 6 PPI network construction. (A) PPI networks of DEGs in AS22w and Sham22W groups. (B) Top 9 hub genes based on PPI network.

Figure 7 The expression of Hub genes. (A) Expression of Hub genes in AS6w group and Sham6W group. (B) Expression of Hub genes in AS22w group and Sham22W group. (C) Expression of Hub genes in the GSE100927 dataset of carotid atherosclerotic lesion group and normal carotid artery group (ns p≥0.05;*p< 0.05;***p<0.001).

The interaction network based on the TRRUST database for predicting transcription factors binding to DEGs of the AS22w group and Sham22W group is shown in Figure 8A. The interaction network based on the TRRUST database for predicting transcription factors binding to AS22w group and Sham22W group DEGs is shown in Figure 8B; the interaction network based on the StarBase database for predicting RPB combined with AS22w group and Sham22W group DEGs is shown in Figure 8B; and the network based on the CTD database for predicting chemical drugs related to AS22w group and Sham22W group DEGs is shown in Figure 8C.

Figure 8 Construction of related regulatory network. (A) mRNA-TF regulatory network. (B) mRNA-RBP regulatory network. (C) mRNA-Chemical network.

Enrichment Analysis

We conducted GO annotation and KEGG signaling pathway enrichment analyses based on the DEGs between the AS22W and Sham22W groups (Figure 9, Tables 2–3). The GO annotation results showed that the DEGs were mainly related to BPs, such as regulation of the body fluid levels, regulation of small molecule metabolic processes, regulation of cytokine-mediated signaling pathways, regulation of ion transmembrane transport, and regulation of the response to cytokine stimulus (Figure 9A). At the same time, it is mainly enriched in Cellular Components (CC) such as extracellular matrix, transporter complex, transmembrane transporter complex, ion channel complex, cation channel complex (Figure 9B). In terms of Molecular Function (MF), it is mainly enriched in receptor ligand activity, endopeptidase activity, organic acid binding, proton transmembrane transporter activity, scaffold protein binding, collagen binding (Figure 9C). The KEGG pathway enrichment analysis showed that the DEGs were mainly enriched in the prion disease, cardiac muscle contraction, Alzheimer’s disease, neurodegeneration-multiple diseases, oxidative phosphorylation, amyotrophic lateral sclerosis, and Parkinson’s disease pathways (Figure 9D). Moreover, GSEA enrichment analysis showed that the atherosclerotic samples were significantly enriched in the KEGG fatty acid metabolism, KEGG peroxisome, WP amino acid metabolism, KEGG PPAR signaling pathway, WP ferroptosis, KEGG tryptophan metabolism, KEGG proteasome, and PID integrin two pathways (Figure 10, Table 4).

Table 2 GO Enrichment Analysis

Table 3 GSEA Enrichment Analysis

Table 4 KEGG Signaling Pathway

Figure 9 GO and KEGG enrichment analysis of DEGs. (A–C) GO analysis based on differential genes of AS22W group and Sham22W group. ((A) Biological Process. (B) Cellular Component. (C) Molecular Function). (D) KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis based on differential genes of AS22W group and Sham22W group.

Figure 10 GSEA analysis. (A) KEGG FATTY ACID METABOLISM. (B) KEGG PEROXISOME. (C) WP AMINO ACID METABOLISM. (D) KEGG PPAR SIGNALING PATHWAY. (E) WP FERROPTOSIS. (F) KEGG TRYPTOPHAN METABOLISM. (G) KEGG PROTEASOME. (H) PID INTEGRIN2 PATHWAY.

Furthermore, the GSVA enrichment analysis showed that the atherosclerotic samples were significantly enriched in hallmark oxidative phosphorylation, PI3K-AKT-mTOR signaling, glycolysis, peroxisome, adipogenesis, mTORC1 signaling, fatty acid metabolism, reactive oxygen species pathway, KRAS signaling, and xenobiotic metabolism (Figure 11, Table 5).

Table 5 GSVA Enrichment Analysis

Figure 11 GSVA Enrichment analysis. (A and B) Enrichment pathway in AS22w group and Sham22W difference analysis box map and heat map (ns p≥0.05, *p< 0.05, **p<0.01, *** p<0.001).

Immunoassay

Finally, we used CIBERSORT to evaluate the abundance of immune cells in the AS22W and Sham22W samples (Figure 12A). Significant negative correlations between Atp6v1d and follicular helper T cells, Atp6v1a and T cells, CD8 and regulatory T cells (Tregs), Atp6v0b and macrophages M1, and Atp6ap2 and resting mast cells were observed (Figure 12B)

Figure 12 Immunoassay. (A) Abundance of immune cells in AS22w and Sham22W samples. (B) Heat map of correlation between Hub genes and levels of immune cell infiltration (*p< 0.05).

Corroboration of the Hub Genes Using QT-PCR

PPI network analysis revealed that Atp6ap2, Atp6v0b, Atp6v0d2, Atp6v1a, and Atp6v1d were significantly expressed in the AS22w group. QT-PCR performed on the ApoE^−/− mouse aortas showed that the expression of two core genes (Atp6v0d2, Atp6v1a), especially Atp6v0d2 (P < 0.0001), increased during atherosclerotic plaque progression (Figure 13A and B)

Figure 13 (A and B) Gene expression of Atp6v0d2, Atp6v1a in aorta of mice in each group (**p<0.01, ***p<0.001, ****p<0.0001). (C) The proportion of positive area of co-localization of CD68 and Atp6v0d2 by immunofluorescence double staining of aortic sinus (****p<0.0001). (D) Representative immunofluorescence analysis for detecting the colocalization (yellow particles) of CD68 (specific for monocyte macrophages, green particles) and Atp6v0d2 (red particles) in frozen sections from the aortic sinus of ApoE−/− mice; red arrow: the co-localization of CD68 and Atp6v0d2; Objective magnification 4×, 18×.

Co-Localization of Atp6v0d2 and Macrophages in Atherosclerotic Lesions

Macrophages, endothelial cells, and vascular smooth muscle cells are the main cell types present in atherosclerotic plaques.²² The colocalization of Atp6v0d2 and macrophages in the atherosclerotic lesions in mice was verified via double immunofluorescence of Atp6v0d2 and CD68 (a macrophage marker) (Figure 13C and D).

Discussion

Atherosclerosis is the root cause of cardiovascular and cerebrovascular diseases and is associated with high morbidity and mortality.^23,24 The mechanisms underlying the formation and development of atherosclerotic plaques are complex. In recent decades, the rapid development of high-throughput detection technology and the progress made in bioinformatics have resulted in the development of new methods for studying the pathophysiological mechanisms of atherosclerosis. However, owing to the complexity of its pathogenesis, effective diagnostic and treatment strategies have not been employed in routine clinical practice.

We used high-throughput sequencing to analyze four groups of aortic transcripts of atherosclerosis ApoE^−/− mouse models under different experimental conditions. A total of 265 DEGs were found. GO annotation showed that the DEGs were mainly involved in BPs such as the regulation of body fluid levels, regulation of small molecule metabolic processes, regulation of cytokine-mediated signaling pathways, regulation of ion transmembrane transport, and regulation of the response to cytokine stimuli. GSEA and GSVA showed that the atherosclerosis samples were significantly enriched in fatty acid metabolism, peroxisome, amino acid metabolism, oxidative phosphorylation, PI3K-AKT-mTOR signaling, glycolysis, among other pathways. These signaling pathways are related to the pathogenesis of atherosclerosis.²⁵ Glycolysis is upregulated in atherosclerosis. Atherosclerotic plaques are characterized by local hypoxia. Under hypoxic conditions, HIF-1α activates the glycolytic pathway, thus reducing cell dependence on oxidative phosphorylation (OXPHOS). Oxidized lipids generated by LPO play an important role in inflammatory responses in atherosclerosis and CVD.^26–28 It is reported that the signal pathways of typical cytokines such as interleukin, interferon, tumor necrosis factor, chemokine, and TGF- β are related to atherosclerosis.²⁹ Amino acid metabolism has been identified as a key regulator of vascular homeostasis and immune cell function. In the process of L-arginine (Arg) and L-hyperarginine (HArg), there are circulatory effects among endothelial cells, innate immune cells, and acquired immune cells, and among L-tryptophan (Trp) metabolism, which may have different effects on the development of atherosclerosis.³⁰ It is well known that inhibition of the PI3K-Akt signaling pathway may reduce the vulnerability of atherosclerotic plaques.³¹ PPARs regulate multiple genes related to lipid metabolism, inflammation, and oxidative stress in the cardiovascular system.^32,33 The involvement of fatty acids and their lipid metabolites in cross-linking with inflammation in atherosclerosis is of interest. Fatty acids are involved in atherogenesis in different ways, demonstrating both pro- and anti-atherogenic functions. They may contribute to inflammation and endothelial dysfunction by participating in many signaling pathways.³⁴ Ferroptosis is a non-apoptotic form of cell death, different from cell necrosis and autophagy. It is characterized by high levels of iron-catalyzed free radicals and the accumulation of lipid reactive oxygen species induced by iron, resulting in oxidative stress and subsequent DNA, protein, and lipid damage. It has been proved that it is involved in the occurrence and progression of atherosclerosis through various signal pathways.³⁵ Busnelli et al used high-throughput sequencing technology to analyze the mouse aortic transcriptome. They showed that in addition to the well-known activation of inflammation and immune response, the damage of the phagosome-lysosome system and osteoclast differentiation was also related to the development of atherosclerosis.

To identify the core genes involved in the progression of aortic atherosclerotic plaque in male ApoE^−/− mice, We constructed DEG PPI networks and screened the hub genes involved in atheromatous plaque progression using the STRING database and the CytoHubba plug-in in Cytoscape. Atp6ap2, Atp6v0b, Atp6v0d2, Atp6v1a, and Atp6v1d were significantly overexpressed in the AS22w group. The hub genes were also significantly overexpressed in human carotid atherosclerotic lesions based on GSE100927 data. Finally, qRT-PCR showed that the mRNA levels of two core genes (Atp6v0d2 and Atp6v1a) were upregulated during plaque progression in atherosclerosis ApoE^−/− mouse models; in particular, the RNA levels of Atp6v0d2 were significantly increased.

Atp6v0d2 (V-ATPaseD2 subunit), as a specific macrophage subunit of vacuolar ATPases (V-ATPases), has been shown to promote the formation of autophagy lysosomes in vitro.³⁶ The V-ATPases are large protein complexes with a peripheral V1 domain that hydrolyzes ATP and an integral V0 domain that transports hydrogen ions, acidifying intracellular compartments, including lysosomes, which are necessary to complete autophagy.³⁷ After severe burns, the activity of myocardial lysosomal V-ATPase decreases in rats, which inhibits myocardial autophagy flux and causes myocardial injury.³⁸ ATPases also work through a pathway that does not depend on lysosomal acidification, such as the fusion of lysosomes with other organelles to form larger compartments. This suggests that V-ATPases may play a role in the fusion of autophagosomes and lysosomes.^39,40 V-ATPases subunits in V1 and V0 domains, including D subunits, play a coordinating role in controlling the coupling of proton transport and ATP hydrolysis.⁴¹ Previous studies have shown that in vitro, the D1 subunit collaboratively promotes the acidification of macrophage lysosomes, and the D2 subunit promotes macrophage autophagosome-lysosome fusion by promoting STX17 binding to VAMP8.³⁹ The expression of Atp6v0d2 in liver macrophages was upregulated after liver ischemia-reperfusion, and by gradually promoting the formation of autophagolysosomes to increase autophagy flux to limit the activation of liver inflammation.⁴²

The activation of the phagocytic-lysosomal pathway in macrophages is considered a critical event in atherosclerotic plaque progression.⁴ Lysosomes are tiny acidic organelles that have the activity of up to 60 different hydrolases, including proteases, lipases, and nucleases.⁴³ The dynamic process of autophagy usually begins with the formation of two-membrane autophagosomes, engulfs cytoplasmic substances, and ends with the fusion and degradation of autophagosomes and lysosomes.^44,45 Therefore, it is essential to enhance the formation/maturation of autophagosomes and promote the fusion of autophagosome/lysosome for the function of autophagy.^46,47 Rajat Singh first reported in 2009 that autophagy regulates intracellular lipid storage.⁴⁸ Subsequent studies have shown that lysosomal acid lipases mediate the hydrolysis of cholesterol esters, which are transferred from lipid droplets to lysosomes via autophagy.^49,50 Within macrophages, lysosomes degrade extracellular substances, including lipids, through heterophagy and intracellular substances, such as lipid droplets, through autophagy.⁵¹ Some studies have shown that lysosome / autophagy pathway is the key regulator of cell metabolism, and the genetic and environmental factors that affect lysosome homeostasis can affect systemic metabolism, especially lipid metabolism.⁵² Under oxidative stress, low-density lipoprotein, and cholesterol esters are easily oxidized to complex oxidation products through free radical-induced lipid peroxidation (LPO).⁵³ OxLDL is usually transported to cells through receptor-mediated endocytosis, which not only partially inactivates lysosomal enzymes but also destroys the stability of acid vacuolar chambers and destroys lysosomal membranes, resulting in lysosomal enzymes and cathepsin leakage into the cytoplasm, thus inducing cell apoptosis.⁵⁴ Cholesterol oxidation products (ChOx) have been reported to be the significant cytotoxic components of oxidized LDL/LDL-, 7b-hydroxycholesterol (7bOH) and 7-ketocholesterol (7keto) are the main toxic components of oxLDL, which have toxic effects on arterial wall cells, including monocytes, macrophages, and smooth muscle cells. The mechanism may be related to early lysosomal lipid accumulation, activation of the lysosome, cellular oxidative stress, and mitochondrial pathway.^6,55 In addition, the expression of lysosomal cathepsin L was significantly increased in atherosclerotic plaques, which was mainly related to cd68 positive macrophages, apoptosis, and stress protein ferritin. Macrophage apoptosis was significantly correlated with the expression of cathepsin L in the nucleus and cell membrane.⁵⁶ Hyperlipidemia and atherosclerosis affect the autophagy lysosomal system. Oxidative stress and lysosomal iron induced by increased permeability and rupture of lysosomal membrane lead to ferritin induction, accompanied by mitochondrial pathway activation and cell apoptosis. Impaired autophagy and apoptosis of macrophages lead to a vicious cycle that aggravates atherosclerosis.⁵⁷ The negative association between autophagy activity and human atherosclerotic progression suggests that autophagy activity in early human atherosclerotic lesions may be a transient self-defense that then declines with prolonged lipid oxidation and oxidative stress and that the flux of autophagy through lysosomes decreases with plaque progression in advanced human atherosclerotic lesions. Autophagy activity may play a key role in maintaining the stability of atherosclerotic plaques by regulating lipid accumulation and necrotic core formation.⁵⁸ Therefore, improving autophagy damage and promoting autophagy-lysosome fusion is one of the treatment ideas to improve atherosclerosis. Our study shows that Atp6v0d2 is significantly up-regulated in aortic atherosclerosis in male ApoE^−/− mice, which is considered as the compensatory protective response of pressure overload, and its effect on autophagy-lysosome can be used as a target for the study of atherosclerosis mechanism.

Conclusion

In summary, Atp6v0d2 was significantly upregulated during the progression of atherosclerotic plaques in ApoE^−/− mice and the double immunofluorescence confirmed its colocalization with the macrophages. Atp6v0d2 may affect atherosclerosis progression by regulating atherosclerotic macrophage autophagy, and the upregulation of Atp6v0d2 promotes autophagy as a compensatory protective response to pressure overload. This study is only a preliminary screening of the hub genes involved in the atherosclerosis process, and the specific pathophysiological mechanisms in which these hub genes are involved require further research. In a future study, we will clarify the mechanism by which Atp6v0d2 affects autophagy in atherosclerotic macrophages in vitro and in vivo and provide new potential targets for atherosclerosis treatment.