The details of overall analysis flowchart are shown in Figure 1. In the beginning, clinical annotation and gene expression profiling by array data of GSE119717 (14) were retrieved from Gene Expression Omnibus (GEO) database. This dataset belonged to Illumina^® platforms (GPL10558), including PVAT around the AAA or NAA (30 vs. 30 samples). The raw data were processed and normalized using the “limma” (version 3.46.0) (23) R package. Gene annotations for each probe set were derived from the Ensembl database. If multiple probe sets represent the same gene, the probe set with the highest mean intensity across all samples was retained.

Weighted gene coexpression network analysis (WGCNA) was applied to summarize gene expression data into coexpression modules. WGCNA has an obvious advantage to reveal the complex relationships among modules, networks, and phenotypes. In detail, based on gene expression profiles and GPL annotation data, a total of 14,016 genes were identified from each sample in this study. The top 5,000 genes were selected according to the median absolute deviation (MAD) of each gene. Then, the top 5,000 genes were used to construct a coexpression network by “WGCNA” (version 1.70–3) R package (24). To reach the goal of scale-free network distribution, the Pearson's correlation coefficient cor (i, j) was used to acquire the similarity matrix by calculating the link strength between nodes i and j from coexpression network and selecting the optimum soft threshold β. Next, the weighted adjacency matrix was constructed as follows: a_ij = [0.5 × (1 + cor (i, j))]^β. Subsequently, the adjacency matrix was transformed into topological overlap matrix (TOM), representing the overlap of network neighbors, and (1-TOM) described the dissimilarity between genes to identify hierarchical clustering nodes and modules. Eventually, the hierarchical clustering results were partitioned by dynamic tree algorithm (deep-split = 2; merge module height = 0.30; minimum module size = 30; cut tree height = 0.99). To evaluate and investigate the coexpression similarities of all modules, the dissimilarity of module eigengene was calculated. Finally, modules with a correlation of 0.70 were merged by “mergeCloseModules” function with essential parameters.

The MEs were used for the component analysis of each module, and modules with similar expression profiles showed highly correlated eigengenes. The correlation between coexpression modules and clinical traits was estimated based on the phenotypic information of adipose tissues around the abdominal aorta. The gene module with the highest correlation coefficient and p < 0.05 was considered the most relevant module to ischemic events and was defined as the key module.

Functional enrichment was assessed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The process was complete by “clusterProfiler” (version 3.18.1) R package (25). The terms with false discovery rate (FDR) < 0.05 were considered significant.

Protein–protein interaction (PPI) is a protein network analytic algorithm, which was used to better understand molecular function and cellular mechanisms, whereas network analysis is an effective method to reveal all kinds of PPI networks. Here, we submitted all genes within the key module to the STRING (version 11.0; https://string-db.org/) database with a minimum level of confidence >0.4 (26, 27). Then, PPI networks were constructed by Cytoscape software (version 3.8.0; http://www.cytoscape.org) (28, 29). Based on maximal clique centrality algorithm, significant models and hub genes were calculated and selected by molecular complex detection (MCODE) which was a plugin for Cytoscape; parameters setting for MCODE: maximum depth = 100, node score cut ≥2, K-core ≥2, degree cut ≥2.

We used random forest (RF), least absolute shrinkage and selection operator (LASSO) logistic regression, and support vector machine-recursive feature elimination (SVM-RFE) to identify key biomarkers for PVAT. RF is a nonparametric, supervised method of classification (30). It is a combination of tree predictors and with the same distribution for all trees in the forests. The RF algorithm was applied with the “RandomForestSRC” package. LASSO regression has a contractile penalty function on variables that induces sparsity of predictors in the expression profile. To achieve this purpose, the “glmnet” package was applied (31). Furthermore, SVM-RFE is used to find the best variables by deleting SVM-generated eigenvectors, which is a machine learning method based on support vector machine. To further identify the diagnostic value of biomarkers in PVAT, SVM module was developed by the “e1071” package (32). Ultimately, we combined the genes from the three machine learning algorithms for further analysis. Bilateral p < 0.05 were considered as statistically significant. The diagnostic efficacy of biomarkers was assessed by receiver operating characteristic (ROC) curve, area under the ROC (AUC), and concordance index (C-index).

This project was approved by the Ethics Committee Board of The First Affiliated Hospital of Zhengzhou University. According to a previous study (33), the eight-week-old C57BL/6 male mice (n = 10) were fed with normal diet. A total of 0.2% 3-aminopropionitrile fumarate salt drinking water was provided to mice 2 days before surgery until the end of study. A total of 10% chloral hydrate (0.4 mL/100 g) was injected intraperitoneally. The infrarenal segment of the abdominal aorta was exposed after complete anesthesia in mice, and the abdominal aorta and inferior vena cava were carefully separated. Instill 20 μL elastase (Sigma-Aldrich; 10.3 mg protein/mL, 5.9 U/mg protein) was around the exposed aorta adventitia. After incubating for 20 min, the abdominal cavity was flushed with normal saline three times, and the abdominal incision was sutured and disinfected. Following 28 days, mice were sacrificed. The abdominal aorta (vessel wall) was anatomized and divided into AAA and NAA groups and stored at −80°C until use.

A total of 55 peripheral venous blood samples (AAA vs. NAA = 40:15) were collected from patients in the First Affiliated Hospital of Zhengzhou University. Baseline clinical data and patient characteristics are shown in Supplementary Table S1. For all the study participants, venous peripheral blood samples were collected on admission, and the plasma supernatant was stored at – 80°C after centrifugation. This project was approved by the Ethics Committee Board of The First Affiliated Hospital of Zhengzhou University. All patients signed informed consent documents.

Total RNA was isolated from human serum and mice abdominal aorta vessel wall tissues using RNAiso Plus (Takara, Dalian, China) according to the manufacturer's protocol. The RNA quality was evaluated using a NanoDrop One C (Thermo Fisher Scientific, Waltham, USA) ultramicro UV spectrophotometer, and the RNA integrity was assessed using agarose gel electrophoresis. Reverse transcription was performed using the PrimeScript RT reagent Kit (Takara, Dalian, China) with genomic DNA (gDNA) Eraser.

After removing the gDNA at 42°C for 2 min, the tissue RNA was reverse-transcribed into cDNA under the following conditions: 37°C for 15 min and 85°C for 5 s. Serum RNA was reverse-transcribed into cDNA using a RevertAid H Minus First-Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, USA) under the following conditions: 25°C for 5 min, 42°C for 60 min, and 70°C for 5 min. The product was immediately stored at −80°C until use.

The qRT-PCR was performed on a QuantStudio 5 Real-Time PCR System (Applied Biosystems, Foster City, USA) using a Hieff qPCR SYBR Green Master Mix (Yeasen, Shanghai, China). qRT-PCR assays were performed in triplicate with the following conditions: (1) 95°C for 5 min and (2) 40 cycles of 95°C for 10 s and 60°C for 30 s. The relative expression of mRNA was calculated using the ΔCT (Ct mRNA-Ct GAPDH) method. The relative quantification values for mRNA were calculated by the 2^−ΔΔCt method. The qRT-PCR primer sequences were provided in Supplementary Table S2. GAPDH was used as an endogenous control for normalization.

Before GSEA, we used the “limma” package (23) to calculate the value of log2 (fold change) between PVAT around AAA and NAA. GSEA was used to decipher the underlying biological mechanisms of the adipose tissues around AAA by KEGG and GO terms (Molecular Signatures Database, version: c2.cp.kegg.v7.4.symbols.gmt and c5.go.v7.4.symbols.gmt). To calculate immune cell infiltration abundance, we uploaded the gene expression matrix data to CIBERSORT. After that, we drawn a correlation heatmap by “corrplot” package to display the correlation of immune cells. Furthermore, “ggplot2” package was used to draw boxplots to visualize the differences in PVAT immune cell infiltration around AAA and NAA. A two-dimensional principal component analysis (PCA) clustering map was drawn based on the PCA clustering analysis on immune cell infiltration data (which was calculated by “factoextra” package).

The gene expression matrix of GSE119717 was obtained (14,016 genes) after data preprocessing. WGCNA analysis was carried out on top 5,000 genes which were selected based on the MAD. In detail, the soft-thresholding power in WGCNA was selected based on scale-free R², and we chose six as the optimal soft threshold by the “pickSoftThreshold” function (Figures 2A,B). Subsequently, the coexpression modules in the network were identified by the “cutreeDynamic” function, and eleven modules were identified based on dynamic tree clipping and average hierarchical clustering (Figures 2C,D). The number of genes in every module is shown in Supplementary Table S4. Ultimately, the grey60 module was highly related to adipose tissues around AAA (Figures 3B,C); this module was selected as a key module for further analysis. The heatmap depicted TOM among 400 genes (which was randomly selected from total genes) in WGCNA (Figure 3A), and the scatterplot displayed correlation between PVAT and 75 genes in the grey60 module (Figure 3D).

To obtain a better understanding of the functional features of 75 genes in the grey60 module, GO and KEGG enrichment analyses were performed by “enrichGO” and “enrichKEGG” function. Terms with p < 0.05 (GO and KEGG) were visualized with R language as a bubble plot. As shown in the results, terms such as ERK1 and ERK2 cascade, mononuclear cell migration, cellular response to interleukin-1, and monocyte chemotaxis were enriched in GO biological processes (Figure 4A). Furthermore, toll-like receptor signaling pathway and lipid and atherosclerosis were also significantly enriched (Figure 4B).

To identify the hub genes of these genes related to AAA in the grey60 module, we constructed a PPI network by STRING database and Cytoscape software (Figure 4C).

Subsequently, key modules in the PPI network were identified by MCODE plugin. The key module (score = 8.667, nodes = 13, edges = 52) consisted of 13 target genes, including KLF6, MYC, FOS, BCL3, DUSP1, KLF4, IL6, MUC1, DUSP5, JUN, PTGS2, EGR1, and BTG2 (Figure 4D).

To identify diagnostic markers for AAA, a variety of machine learning methods was utilized. First, we used the LASSO regression algorithm to identify six genes from the grey60 module as diagnostic markers for PVAT around AAA, including EGR1, DUSP1, CYR61, FOSB, KLF4, and CREB5 (Figures 5A,B); eight genes were obtained from the grey60 module by the SVM-RFE algorithm as diagnostic markers, including APOLD1, PLA2G2A, NR4A2, EGR1, KLF4, CCL8, CCL3L3, and FOS (Figure 5C); we also used the RF algorithm to identity eight genes from the grey60 module as diagnostic markers, including FOS, EGR1, DUSP1, CYR61, FOSB, JUN, KLF4, and CREB5 (Figure 5D). Subsequently, we overlapped the three sets of diagnostic genes and finally obtained two diagnostic markers (EGF1 and KLF4) (Figure 5E). To further examine the diagnostic efficacy of EGF1 and KLF4, ROC curve analysis was used and the AUC values were 0.916, 0.926, and 0.948, respectively, for EGF1, KLF4, and two-gene combined, respectively (Figure 5F) and the C-indexes all >0.9 (Supplementary Figure S1A).

As illustrated in Figure 6A, we successfully constructed the AAA mice model. Subsequently, the evaluation of diagnostic efficacy of biomarkers was performed by ROC curve and C-index. The AUCs of EGF1, KLF4, and two-gene combined respectively were 0.870, 0.800, and 0.910 in the mice cohort (Figure 6B) and 0.840, 0.845, and 0.898 in the inhouse cohort (Figure 6C). The C-indexes of these diagnosis makers in the two cohorts were also larger than 0.8 (Supplementary Figures S1B,C), which EGF1 and KLF4 have high diagnostic and stability. By further analysis, we found significant differences in the expression levels of EGR1 and KLF4 among the GSE119717, mice, and inhouse cohorts (Figures 6D–I).

To characterize the relevant differences in gene expression between PVAT around AAA and NAA, we performed GSEA analysis based on the rank log2 (fold change) of all gene list. The results displayed that the significant GO terms mainly focus on extracellular matrix and fiber constitute (Figures 7A,B), such as “collagen-containing extracellular matrix” (normalized enrichment score (NES) = 2.324, FDR = 0.007), “extracellular matrix” (NES = 2.320, FDR = 0.007), contractile fiber (NES = −2.471, FDR =0.017), and structural constituent of ribosome (NES = −2.579, FDR = 0.017). Furthermore, the significant KEGG pathway (Figures 7C,D) mainly focus on “MAPK signaling pathway” (NES = 1.863, FDR = 0.005), PI3K-Akt signaling pathway (NES = 1.790, FDR = 0.012), and amyotrophic lateral sclerosis (NES = −2.002, FDR = 0.009).

The immune cells' score was calculated by CIBERSORT in two groups. Correlation heatmap (Figure 8A) of the 22 types of immune cells displayed that T cells gamma delta and macrophages M1, B cells memory, and regulatory T cells had a significant positive correlation. In contrast, T cells CD memory resting and regulatory T cells, T cells CD native, and macrophages M2 had a significant negative correlation. PCA revealed that there was a significant difference in immune cell infiltration between PVAT around AAA and NAA samples (Figure 8C). The boxplot of the immune cell infiltration difference displayed that, compared with NAA group, T cells follicular helper, NK cells activated, monocytes, dendritic cells activated, and neutrophils infiltrated were more abundant, whereas NK cells resting infiltrated less (Wilcoxon test, Figure 8B).

Correlation analysis revealed that EGR1 was positively correlated with mast cells activated (r = 0.451, p = 0.012) and plasma cells (r = 0.360, p = 0.02) and negatively correlated with T cells CD memory resting (r = −0.252, p = 0.007) (Figure 8D). KLF4 was positively correlated with mast cells activated (r = 0.279, p = 0.028) and negatively correlated with macrophages M0 (r = −0.355, p = 0.003) (Figure 8E).