We used published genome-wide association studies (GWAS) summary data for 91 circulating inflammatory cytokines, 731 immune cells, and aortic aneurysm (AA). In our study, multiple two-sample MR analysis was carried out in order to investigate causality between inflammatory cytokines and AA and between immune traits and AA. Meanwhile, the reverse MR method was conducted to explore the potential reverse causal effects between AA and inflammatory cytokines. The MR analysis was based on three key assumptions: 1. the genetic variants selected as IVs were truly associated with the exposure; 2. the genetic variants had to be independent of any confounding factors; 3. the genetic variants affected the outcome only through exposure. Bulk RNA sequencing data and single-cell RNA sequencing data were utilized to validate the MR findings. The study design overview is presented in Figure 1 .

GWAS summary statics for AA were obtained from the GWAS Catalog (accession number GCST90018783). A total of 3,230 European ancestry cases, 475,964 European ancestry controls (Finngen and UK Biobank), 1,155 East Asian ancestry cases, and 173,601 East Asian ancestry controls (Biobank Japan) were included (8).

GWAS summary statistics for the level of 91 circulating cytokines were derived from the largest and latest available genome-wide protein quantitative trait locus (pQTL) study, which covers 14824 participants from 11 independent cohort studies (9).

GWAS summary statistics for each immune trait are publicly available from the GWAS Catalog (accession numbers from GCST0001391 to GCST0002121). A total of 731 immunophenotypes including absolute cell (AC) counts (n = 118), median fluorescence intensities (MFI) reflecting surface antigen levels (n = 389), morphological parameters [MP] (n = 32) and relative cell (RC) counts (n = 192) were included. The original GWAS on immune traits was performed using data from 3,757 European individuals and there were no overlapping cohorts (10).

Firstly, we set a threshold of p<5×10^-8 to screen SNPs strongly related to inflammatory cytokines, immune traits, and AA, however, quite a few SNPs were screened out. Then we applied a genome-wide significance threshold of p<5×10^-6. Secondly, we clumped SNPs with the criteria where kb = 10000, and r² = 0.001 to avoid linkage disequilibrium (LD). Thirdly, the bottom line of F-values for each SNP was set to 10 to avoid weak instrumental bias (11).

MR analysis is conducted utilizing the “TwoSampleMR” R package (Version 0.5.6) of the RStudio (version 4.2.3). Inverse variance weighting (IVW) was employed as the main method, with MR Egger regression, weighted median, and simple median as complementary methods. MR results were ultimately identified as statistically significant only if they met with the following conditions: 1. the 95% confidence intervals (CI) derived from IVW did not cross over 1, 2. p<0.05, 3. IVW results met the consistent trend of effects derived from the other three MR analyses.

We employed three sensitivity analysis methods to evaluate the sensitivity of MR results, including the heterogeneity test, pleiotropy test, and leave-one-out sensitivity test. Cochran’s Q test was used to determine the heterogeneity, and p>0.05 was taken as no heterogeneity. The intercept of the MR Eggar analysis results was used to evaluate the horizontal pleiotropy, with p>0.05 considered as no pleiotropy. Leave-one-out sensitivity was used to investigate whether the causal association was influenced by SNP shared between exposure and outcome.

Transcriptional profiles of AA were obtained from the GEO dataset with the accession number GSE141032 containing 12 AA samples and 12 organ donor samples utilized as controls. The raw counts from the original matrix were processed with R package “DESeq2” to obtain the normalized counts and significantly differentially expressed genes (DEGs) with p<0.05 and log2 fold change (FC) >0.95.

Previously published single cell sequencing (scRNA-seq) data from control (n = 3) and ascending thoracic aortic aneurysm cases (n = 8) were included for reanalysis (3). Cells were removed under the condition of expressing fewer than 200 genes or greater than 20% mitochondrial genes. R package Seurat (version 4.2.0) were used to perform dimensional reduction of scRNA-seq data. The “NormalizeData” and “ScaleData” function from Seurat was used for normalization and the scale factor was set to 100,000, then followed by “FindVariableFeatures” with default parameters to calculate highly variable genes for each sample. “FindIntegrationAnchors” in the Seurat package were applied to remove batch effect and merge samples. Top 2 000 genes with the highest expression and dispersion from each sample were determined as the integration anchors and used for data integration. Then, cells were projected in 2D space using Uniform Manifold Approximation and Projection (UMAP). To identify differentially expressed genes (DEGs) among each cluster, the “FindAllMarkers” function from Seurat was used and non-parametric Wilcoxon rank sum tests were set to evaluate the significance of each individual DEG. DEGs with adjusted P value less than 0.05 were thought to be significant and used in downstream analysis. Hierarchical clustering and heatmap generation were performed for single cells on the basis of normalized expression values of marker genes curated from the literature or identified significant DEGs. “DimPlot” and “VlnPlot” were used to visualize the expression of individual genes, cells were grouped by their cell type as determined by analysis with Seurat.

The trajectory analysis based on gene expression changes among the cell clusters was performed using R package Slingshot (version 2.1.0.0). We infers the key elements of the global lineage structure using a minimum spanning tree method based on clustered cells. The start point of trajectories was defined based on DEGs.

CellChat analysis was performed by R package CellChat (version 1.6.1). Function ‘createCellChat’ was selected to create a CellChat object for the overall gene expression matrix, and then the cell annotation information identified by scMGCA was added through the function ‘addMeta’. CellChatDB.human provided by CellChat was used for the human ligand receptor reference. Finally, the cell-cell communication probability is inferred using the function ‘computeCommunProb’, and the communication probability at the level of each cell signaling pathway is inferred by the function ‘computeCommunProbPathway’. The figure is generated with the function ‘netVisual_aggregate’.

Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to. The R package “org.Hs.eg.db” was used to obtain the Entrez ID for each DEG, and “clusterProfiler” was employed for biological function analyses.

Cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) was utilized to deconvolute gene expression into immune cell fractions (12). Spearman correlation analysis was conducted to determine the correlation between the expression of CX3CR1 and infiltrated immune cells.

Aortic wall tissue was taken from patients undergoing aortic replacement surgery due to AA and the recipient undergoing heart transplantation due to ischemic cardiomyopathy at Wuhan Union Hospital, Wuhan, China. The specimens were fixed in 4% paraformaldehyde immediately after separation from the body. The Review Board of Union Hospital, Affiliated to Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China, approved the collection of human aorta tissues and their use in our study. The Declaration of Helsinki’s guiding principles were followed when conducting the research. Written informed permission was obtained from each recruited patient for the use of their aorta tissue in the study. Immunofluorescence staining was performed as previously described (13).

Under the threshold of p<5×10^-6, fibroblast growth factor 5 (FGF5), stem cell factor (SCF/c-kit), Fractalkine/CX3CL1 and Interleukin-20 receptor subunit alpha (IL20RA) were first screened out of the 91 genetically determined circulating inflammatory cytokines based on the IVW method ( Figure 2A ). Among the four candidates, FGF5, CX3CL1, and IL20RA were associated with an increased risk of AA, and c-kit was associated with a reduced risk of AA ( Figure 2B ).

None candidates showed significant heterogeneity or horizontal pleiotropy (P>0.05) in the sensitivity test except for FGF5 ( Supplementary Table S1 , Supplementary Figure S1 ). CX3CL1 demonstrated a positive correlation solely by the IVW method (OR=1.228, 95%CI=1.029-1.465, p=0.023) ( Figure 3C ). IL20RA exhibited positive correlation with an OR=1.251 (95%CI=1.015-1.541, p=0.035) ( Figure 3E ). SCF was the only negatively correlated cytokine, which was validated by IVW and weight median method (OR=0.870, 95%CI=0.778-0.973, p=0.015; OR=0.813, 95%CI=0.697-0.948, p=0.008) ( Figure 3G ). The scatter plots of SNPs associated with each of the four candidates were also demonstrated in Figure 3 . The funnel plot and forest plot of SNPs were listed in Supplementary Figure S1 .

To further clarify the causal relationships confirmed by the forward analysis, we performed a reverse Mendelian randomization analysis utilizing AA as the exposure. Initial screening using IVW as the primary method suggested that three inflammatory factors, which did not overlap with the four entrants identified in the forward MR analysis, were affected by AA ( Figure 4A ). Of these, monocyte chemoattractant protein-4 (CCL4) was positively regulated by AA (OR=1.031, 95%CI=1.004-1.059, p=0.026), while glial cell line-derived neurotrophic factor (GDNF) and tumor necrosis factor (TNF) were negatively regulated by AA (OR=0.971, 95%CI=0.945-0.997, p=0.032; OR=0.970, 95%CI=0.941-0.999, p=0.046) ( Figure 4B , Supplementary Figure S2 , Supplementary Table S2 ). The results of heterogeneity and pleiotropy analysis were demonstrated in Supplementary Table S2 .

Circulating cytokines generally operate their immune-regulatory function through the interaction with their corresponding receptors. As reported before, CX3CL1 exclusively binds to CX3CR1, which is a 7-transmembrane receptor coupled to heterotrimeric G proteins (GPCR). SCF also cooperates with its exclusive receptor c-kit, which is a tyrosine kinase receptor. The normalized counts of CX3CR1 were significantly up-regulated in AA, while on the contrary, KIT, which encodes c-kit, was significantly down-regulated in AA ( Figures 5A–C ). Their corresponding ligand, CX3CL1, and KITLG showed no significant change between these two groups ( Figure 5A ).

With the criteria of p<0.05, log2FC>0.95, we obtained a total of 1120 DEGs, where 513 DEGs were up-regulated and 607 DEGs were down-regulated in AA ( Figure 5C , Supplementary Table S3 ). Rather surprisingly, we found that CX3CR1 was among the most dramatically up-regulated DEGs in AA ( Figure 5D ). The functional enrichment analysis of the up-regulated DEGs with GO demonstrated that “immune response”, “inflammatory response”, “positive regulation of immune response” and “antigen processing and presentation of exogenous peptide antigen via MHC class II” were among the top 10 enriched terms ( Figure 5E ).

Enlightened by this, we further explored the immune microenvironment utilizing immune infiltration analysis. It turned out that M1 macrophages and monocytes were significantly infiltrated more in AA samples ( Figure 5F ). Another novel finding was that the expression of CX3CR1 was positively correlated with M1 macrophage and CD4⁺ T cell ( Figure 5G ), suggesting its potential role in pro-inflammatory response and adaptive immune response.

Inspired by the putative role of CX3CR1 in regulating the immune response in AA, especially characterized by the significantly up-regulated infiltration of M1 macrophage, we further explored the causal relationship between immunophenotype and AA at the genetic level. The IVW method was used as the main criteria to explore the causal effect of 731 immunophenotypes on AA. With the threshold of p<0.05, we obtained 13 immunophenotypes negatively and 14 immunophenotypes positively correlated with AA ( Figure 6A , Supplementary Table S4 ).

Among the selectees, we noticed that genetically determined CCR2 on CD14+CD16+ monocytes was positively correlated with the onset of AA (OR=1.036, 95%CI=1.002-1.071, P=0.035) ( Supplementary Table S4 ). Furthermore, the MR eggar method demonstrated that CX3CR1 on CD14⁺CD16^- monocytes was positively correlated with AA (OR=1.219, 95%CI=1.091-1.362, P=0.003), although no significance was acquired with IVW method ( Supplementary Table S4 ). Considering the generally recognized classification of monocytes, CD14 and CD16 were used to distinguish between classical (CD14⁺CD16^-) and nonclassical (CD14^-CD16⁺) subtypes, while growing evidence suggested the potential role of CCR2 and CX3CR1 in classifying monocytes in a more detailed manner.

In order to further identify the cell-type localization of the ligand-receptor relationship (CX3CL1-CX3CR1) as indicated above, we re-analyzed the previously published scRNA-seq data from 3 controls and 8 AA samples (3). Similar to the original research, we identified 12 primary cell types, including T cell, smooth muscle cell (SMC), natural killer cell (NK), monocyte, M2 macrophage, fibroblast, mesenchymal stem cell (MSC), M1 macrophage, endothelial cell (EC), B cell, Plasma cell and Mast cell ( Figures 7A–C ). The cell counts of each sample indicated the aggregation of immune cells in AA samples, especially monocytes, M2 macrophages, and M1 macrophages ( Figure 7D ).

Based upon the preliminary clustering, we identified that CX3CL1 was predominantly expressed by ECs and CX3CR1 was mainly expressed by monocytes as demonstrated by the feature plots and violin plots ( Figures 7E, F ). CellChat analysis exhibited an intense signal from ECs to monocytes ( Figures 7G, H ), which indicated that ECs were the major source of CX3CL1 that binds with CX3CR1 on monocytes to leverage downstream effects.

To further identify the heterogenicity of the myeloid cell lineage, we defined 9 distinct clusters from the original monocyte/macrophage population, consisting of M1-like macrophage, CD16⁺ monocyte, monocytic precursor, M2-like macrophage, bipolar macrophage, dendritic cell (DC), neutrophil, remodeling macrophage and proliferating macrophage ( Figures 8A, B ). Remarkably, the signature genes of CD16⁺ monocyte were enriched in GO terms analogous to DEGs between normal and AA tissue, such as “antigen processing and presentation of exogenous peptide antigen via MHC class II” ( Supplementary Figure S3 ).

We noted that the population of cells with high CX3CR1 expression also exhibited high expression of CD14 and CD16 ( Figure 8C ), which was corroborated by colocalized immunofluorescent staining of dilated ascending aorta ( Figure 8D ). Also, we observed co-localization of EC-specific marker CD31 and signaling molecular CX3CL1 ( Figure 8D ). In addition to the normalized count obtained from transcriptomic data ( Figure 5B ), we observed that the expression of CX3CR1 was significantly higher in AA samples compared with normal samples ( Figure 8E ). Trajectory analysis revealed the pluripotency of CD16⁺ monocytes, which could differentiate either to the M1 or M2 polarization status, as well as self-conversion ( Figure 8F , Supplementary Figure S3 ), whereas the driving force for direction might well be different cytokines.

Thus, we explored the CytoSig database which is composed of transcriptomic profiles of targeted genes modulated by cytokines (14). We obtained the processed microarray data of the primary human CD16⁺ monocytes treated with or without recombinant human CX3CL1 ( Supplementary Table S5 ). Through perusal of the sorting protocol, we noticed that these sequenced cells had been positively sorted in advance by CD14. The top 100 up-regulated genes upon CX3CL1 treatment were enriched in the “G-protein coupled receptor signaling pathway”, “adaptive immune response”, “positive regulation of interleukin-8 production” and “positive regulation of cell adhesion” ( Figure 8G ), which all indicated a shift towards the pro-inflammatory state, thus leading to the inflammatory cascade in the pathogenesis of AA.