Four AD-related datasets were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). Table 1 displays the basic information of these datasets. The RNA-seq dataset GSE153434 containing 10 AD samples and 10 control samples was deployed to pick differentially expressed genes (DEGs) and to construct a weighted gene co-expression network to screen AD-related genes. Microarray datasets GSE98770 and GSE52093 were used to verify the expression and diagnostic efficiency of the hub genes. AD samples (n = 5) and control samples (n = 3) of scRNA-seq dataset GSE213740 were enrolled to verify the expression and function of the hub genes at the single cell level. The workflow of this study is shown in Figure 1.

Hierarchical cluster analysis was performed to detect outliers in the samples derived from GSE153434. The No. 4 and No. 10 samples in the control group were removed by setting the high threshold for 120 to cut tree. Thus, we finally included 8 control and 10 patient samples in our following weighted gene co-expression network analysis (WGCNA) (Figure 2A).

The WGCNA R package (version 1.72-1) was utilized for the creation of a co-expression network and associated analyses (17, 18). First, the scale-free network topology analysis indicated that the soft threshold β = 8 made the R² close to 0.8; meanwhile, the mean connectivity was below 100 (Figure 2B). Next, an adjacency matrix was calculated based on the expression matrix and transformed to the topological overlap matrix (TOM). The TOM was further transformed into the distance matrix to build a tree of gene clusters. Dynamic tree cut algorithm was used to generate different gene modules labeled in different colors. Finally, to increase the capacity of the modules, analogous modules with a cut height of <0.25 were combined.

Along with module merging, the module eigengene (ME) of each module was determined. The Pearson correlation coefficient with the sample grouping variable of each ME was calculated. Gene significance (GS) and module membership (MM) were calculated to investigate the connection between the pink module and patients. The significant genes in the pink module were selected by setting the criteria GS > 0.6 and MM > 0.6.

GEO2R online tool and limma R package (version 3.58.1) were used to analyze DEGs of GSE153434, GSE98770, and GSE52093 datasets (19). Adjust P-value <0.05 and |log-fold change| ≥1.0 was used as the threshold criteria. A volcano plot was plotted by the ggplot2 package (version 3.4.2).

The overlapping genes of DEGs and the significant genes in the pink module were picked by the Venn diagram generated by VennDiagram R package (version 1.7.3).

The overlapping genes were extracted and then entered into the STRING database (https://cn.string-db.org) to construct a protein–protein interaction (PPI) network. The network type was chosen as the full string and the required score was set to 0.4 (medium confidence). Cytoscape software (version 3.7.1) was used for visualization of the network. Five methods, closeness, degree, EPC, MCC, and MNN of cytoHubba plug-in, were used to identify hub genes.

The receiver operating characteristic (ROC) curve was utilized to evaluate the diagnostic efficacy of the hub genes. The pROC R package (version 1.18.4) was used to build the ROC models based on the gene expression data (20). The area under the curve (AUC) was given by the module simultaneously.

Seurat R package (version 5.0.0) was performed to assess the expression of the hub genes in the scRNA-seq dataset GSE213740 (21). The Seurat object was created with the function “CreateSeuratObject” and preliminary filtering parameters “min.cells = 3, min.features = 200” were set. The mitochondrial gene expression ratio was calculated with the function “PercentageFeatureSet.” The cells were selected according to the standard of “nFeature_RNA >500 & nFeature_RNA <5,000 & percent.mt <25.” The raw counts were normalized through the “LogNormalize” method. A principal component analysis (PCA) was conducted with the function “RunPCA.” Samples were integrated with the “IntegrateLayers” function and reciprocal PCA (RPCA) method. The “ElbowPlot” function and the 1–30 dimensions were used to perform “FindNeighbors.” Then cells were clustered with the “FindClusters” function under 0.4 resolution and annotated by the canonical cell markers from the CellMarker (http://xteam.xbio.top/CellMarker/) database and literature. The uniform manifold approximation and projection (UMAP) was generated by using the “RunUMAP” function. The differential genes were calculated using the “FindMarkers” function with the Wilcox method.

Gene Ontology (GO) enrichment analysis with over-representation analysis methods was performed by a ClusterProfiler R package (version 4.8.2) (22). Terms with adjusted P value <0.05 were considered significant. The bar chart was plotted by the ggplot2 package.

Three-week-old male C57BL/6J mice were fed a diet containing 0.4% beta-aminopropionitrile (BAPN) (Sigma-Aldrich, St. Louis, MO, USA) to establish the TAAD model as previously reported (23), and the same aged mice fed with regular diet were served as controls. We chose 10 days to represent the early stages prior to TAAD formation to detect the changes in gene expression. Mice were euthanized by intraperitoneal injection of overload sodium pentobarbital (100 mg/kg). Then the entire aortas were separated for RNA extraction.

Total RNA was extracted by using TRIzol (Invitrogen, Carlsbad, CA, USA). Equal amounts (1 μg) were reverse-transcripted into cDNA by using HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme, Nanjing, Jiangsu, China). ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, Jiangsu, China) was used according to the manufacturer's instructions. All amplification reactions were carried out in a program involving a step at 95℃ for 30 s followed by 40 cycles of 95℃ for 10 s and 60℃ for 30 s. The mRNA expression levels were normalized to those of GAPDH using the 2^(-ΔΔCt) method. The primer sequences for the target genes are listed in Table 2.

Wilcoxon Mann–Whitney test was used to compare gene expression between two groups. P value <0.05 was considered to be significant. Pearson correlation coefficients were calculated to assess the correlation between genes. Genes with P value <0.05 and |correlation coefficient| ≥0.6 were considered to be significant.

We used WGCNA to identify the genes linked to AD pathogenesis and gene co-expression modules were built according to the process described above. The dynamic tree cut algorithm produced 14 modules. By combining the similar modules, nine modules were harvested (Figure 2C). We subsequently determined the correlation between each module and AD patients, except the gray module in which genes were not co-expressed. The results showed that the green module had the lowest negative correlation coefficient with AD (r = −0.93, P = 3e-08), while the pink module displayed the strongest positive correlation with AD (r = 0.92, P = 7e-08) (Figure 2D). The module–trait relationship was further identified by the correlation seen in the GS-MM plot in the pink module (Figure 2E). A total of 89 significant genes in the pink module were selected for further analysis.

DEGs of the GSE153434 dataset were analyzed as a complement to the co-expression analysis. There were 419 DEGs, including 177 upregulated genes and 242 downregulated genes (Figure 3A). A total of 71 overlapping genes were acquired by intersecting the DEGs with the significant genes in the pink module (Figure 3B).

All the overlapping genes were used to build the PPI network to screen the hub genes. We got a PPI network with 45 nodes and 74 edges after eliminating the genes that had no connections to others (Figure 3C). Figure 3D shows the number of the interaction genes for each gene. We used five different algorithms to calculate the hub genes respectively (Table 3). By integrating the results of the five methods, we identified the top five hub genes, including hypoxia-inducible factor 1 subunit alpha (HIF1A), hepatocyte growth factor (HGF), heme oxygenase 1 (HMOX1), integrin subunit alpha 5 (ITGA5), and integrin subunit beta 3 (ITGB3).

We evaluated the expression and diagnostic effectiveness of the five hub genes for AD in two testing datasets GSE98770 and GSE52093. The expression of HIF1A and HGF increased significantly in the patient group compared to the control (normal) group in both datasets (Figures 4A,B). All these genes in both datasets had AUC values above 0.55, as shown by ROC curves in Figures 4C,D. These results indicate that these genes, especially HIF1A and HGF, have good diagnostic efficacy.

The scRNA-seq dataset GSE213740 was used to investigate the expression pattern of the hub genes. We first gained 23 clusters under 0.4 resolution from 8 samples and 86,131 cells (Figure 5A,B). Then each cluster was annotated with the canonical cell markers. As a result, eight main cell types, including SMCs, endothelial cells, fibroblasts, macrophages, T lymphocytes, B lymphocytes, natural killer cells, and proliferation cells, were identified (Figure 5C). Figure 5D shows the expression distribution of selected markers mapped on the UMAP plot.

Then we analyzed the single-cell location of five hub genes. We found that HIF1A was expressed in almost all the cell types with the highest expression level in macrophages. HGF was expressed with a certain level in SMCs, endothelial cells, and fibroblasts but with low expression in other cell types. ITGA5 had the highest expression level in endothelial cells, followed by macrophages, SMCs, fibroblasts, and proliferation cells. HMOX1 was highly expressed in macrophages. In contrast, the expression of ITGB3 was quite low in all cell types (Figure 6A).

We further analyzed the expression alterations of the five hub genes in the patient group compared to the control group. We found that HIF1A was significantly upregulated in macrophages, SMCs, and T lymphocytes. HGF was significantly upregulated in SMCs. ITGA5 was significantly upregulated in SMCs, macrophages, and proliferation cells. HMOX1 was significantly upregulated in SMCs. ITGB3 was significantly upregulated in SMCs and macrophages.

Given the high expression of HIF1A in macrophages and its upregulation in AD samples, its function was further analyzed. We first screened the genes related to HIF1A expression in GSE153434 by calculating Pearson correlation coefficients. Two hundred and forty-nine genes were positively correlated to HIF1A while 246 genes were negatively correlated to HIF1A (Figure 7A). GO enrichment analysis showed that the positively correlated genes were mainly focused on cytokine production, macrophage activation, and leukocyte cell–cell adhesion, while the negatively correlated genes were mainly focused on ion transmembrane transport, cell–substrate adhesion, and cyclic nucleotide biosynthesis (Figure 7B). Next, we analyzed the DEGs of macrophages in the patient group compared to the control group (Figure 7C). Considering that the proteins encoded by HIF1A are transcription factors, we determined the downstream functions that HIF1A may regulate directly by predicting the target genes of HIF1A using the TRUST database (https://www.grnpedia.org/trrust/). By intersecting the target genes of HIF1A and the DEGs of macrophages, we got 21 overlapping genes (Figure 7E). Among them, the upregulated genes were related to the inflammation (e.g., TLR2, ALOX5AP, and MIF), glycolysis (e.g., ENO1, LDHA, and GAPDH), tissue remodeling (e.g., PlAU), and angiogenesis (e.g., SERPINE1 and VEGFA).

We also analyzed the function of HGF which was highly expressed in SMCs and significantly upregulated in AD samples in dataset GSE153434. Pearson correlation coefficients showed that there were 186 and 167 genes positively and negatively correlated to HGF, respectively (Figure 8A). GO enrichment analysis showed that the positively correlated genes were mainly focused on cytokine production, the interleukin-1-mediated signaling pathway, and the vascular endothelial growth factor receptor signaling pathway, while the negatively correlated genes were mainly focused on the cAMP biosynthetic process, the cyclic nucleotide biosynthetic process, and actin filament-based movement (Figure 8B). Given that the proteins encoded by HGF are cytokines, we analyzed the receptors to which HGF may bind to. We searched for proteins that interacted with HGF in the STRING database, and 10 proteins were harvested (Figure 8C). Among them, five proteins were receptors, namely, MET, NTRK1, EGFR, IGF1R, and KDR. We assessed the gene expression level of these receptors in different cell types of the ascending aortas (Figure 8D). MET and KDR were specifically highly expressed in endothelial cells. IGF1R had a certain expression level in fibroblasts, endothelial cells, and SMCs. The expression of EGFR was highest in fibroblasts, followed by SMCs (Figure 8D). Based on the above information, we speculated that HGF secreted by SMCs may bind to their own receptors EGFR and IGF1R in an autocrine manner. In addition. HGF may also bind to EGFR and IGF1R in fibroblasts and MET, IGF1R, and KDR receptors in endothelial cells to regulate their biological functions in a paracrine manner. Thus, HGF may play an important role in the cross talk between fibroblasts, endothelial cells, and SMCs (Figure 8E).

To verify the activation of HIF1A and HGF signaling in AD, we analyzed the mRNA expression levels of Hif1a, Hgf, and their potential target genes in a BAPN-induced mouse TAAD model using qRT-PCR. Our results showed a significant upregulation of Hif1a and its several target genes, including Alox5ap, Serpine1, Tlr2, and Plau, in the aortic tissues of the BAPN-treated mice compared to the control group (Figure 9). Similarly, Hgf and its downstream receptor genes, namely, Egfr and Igf1r, were also significantly upregulated in BAPN-treated mice compared to controls (Figure 9). These results indicate the activation of HIF1A and HGF signaling may participate in the occurrence of AD.