Modification of m6A mediates tissue immune microenvironment in calcific aortic valve disease.

Background: Several human diseases are associated with aberrant expression of regulators involved in N6-methyladenosine (m6A) RNA modification. However, their role in aortic valve calcification (AVC) is largely unknown. The aim of this study was to determine the general expression pattern and potential function of m6A regulators in AVC by bioinformatics methods.
Methods: We obtained AVC datasets from the Gene Expression Omnibus (GEO). The identification of m6A-related differentially expressed genes (DEGs) and the Consensus Clustering method was performed to type AVC individuals based DEGs. Then, we quantified the effect of typing by principal component analysis (PCA). Next, we performed the weighted gene co-expression network analysis (WGCNA) and identified the main modules as well as functional analysis. Additionally, the key genes were screened by protein-protein interaction network (PPIN) analysis and identifying important genes of important modules. We again typed AVC individuals by the same method using key genes. Finally, we evaluated the link between key genes and immune infiltration.
Results: We discovered that METTL14, ZC3H13, FTO, FMR1, HNRNPA2B1, HNRNPC, LRPPRC, YTHDC1, YTHDC2, and YTHDF1 expression levels decreased considerably in AVC tissues. Based on 10 genes, we typed 240 AVC samples as clusters A and B. We assessed the immune cell content in 240 samples using Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) and found that B cell memory, CD8 T cells, T follicular helper cells, monocytes, M0 macrophages, resting dendritic cells (DCs), and interleukin-10 (IL-10) were concentrated in the cluster A group. Additionally, based on the important WGCNA modules, we identified 7 key genes. Next, 240 samples were retyped based on 7 key genes; we found that T cells CD8, T cells CD4 memory activated, T cells follicular helper, and macrophages M1 were significantly increased in gene cluster-1. Finally, we performed functional enrichment of gene cluster-typed samples, showing potential functional differences between different types. Conclusions: Our study provides a review of the m6A regulators' expression pattern and functional importance in human AVC. The data from this study might serve as a significant resource for future mechanistic and therapeutic investigations into the role of critical m6A regulators in AVC. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Aortic valve calcification (AVC) is the most common cause of surgical aortic valve replacement across the globe (1,2). Individuals aged above 65 years are more likely to be affected by AVC (2–7% of the population) in developed countries (3), and it is also a frequent consequence of chronic kidney disease (4). The incidence of AVC has been growing over the last decade, a trend which is expected to continue in the future (5-7). The lack of effective conservative treatment of AVC means that surgery remains the only readily available choice (8-12).

In AVC, the valvular interstitial cells (VICs) undergo transdifferentiation and become osteoblast-like cells, leading to the mineralization of the tissue (8,13,14). A range of pathogenic stimuli, including increased phosphate (15,16), may promote the translocation of Runx2 (a transcription factor) to the nucleus, triggering the cell to change to an osteoblast phenotype. Pro-inflammatory cytokines, namely, tumor necrosis factor (TNF) and interleukin-1 (IL-1), are associated with AVC, and inflammation is a characteristic of AVC (17-19). Rajamannan et al. (2) at the National Heart, Lung, and Blood Institute have shed light on AVC through a variety of innovative methods. They discovered that AVC is a dynamic and intricate osteogenic process that is initiated and promoted by inflammation. As previously documented, abnormal hemodynamic forces experienced by valve leaflets, including shear stress, hypertension, or increased strain, may result in tissue remodeling and inflammation which can advance to calcification, stenosis, and heart valve failure (2,8). Furthermore, according to the findings of an immunohistochemistry investigation conducted on the heart valves of 285 AVC patients, 28.4% of the affected tissues showed considerable leukocyte infiltration and high TNF-expression due to chronic inflammation (18). In a mouse model with IL-1 receptor antagonist (IL-1Ra) deficiency, inflamed aortic valve stenosis was identified, and TNF played a role in the initiation of valvular calcification (20). However, more research is required to determine the potential mechanism of the advancement of AVC.

In eukaryotes, N6-methyladenosine (m6A) is a common and well-studied non-cap reversible marker found on messenger RNAs (mRNAs) and long non-coding RNAs (lncRNAs) (21,22). The methylation of m6A is a process in which many proteins are involved. These proteins may broadly be divided into 3 categories: writer complexes, m6A demethylases, and function executions (readers). Within these categories, WTAP, RBM15, RBM15B, METTL3, METTL14, METTL16, and KIAA1429 are among the writer complexes; FTO and ALKBH5, the 2 evaluated m6A demethylases, are among the erasers; HNRNPs, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, YTHDC3, EIF3A, IGF2BP1, IGF2BP2, and IGF2BP3 are among the function executions (readers) (23-26). Recent studies have shown that m6A is required for a wide range of biological activities, including viral infection (27,28), stress (29), heat shock (30), and DNA damage (31). The m6A modification has an impact on molecular functions (MFs) such as RNA-protein interaction (32), RNA stability (33), and translation efficiency (34).

Moreover, m6A dysregulation is linked to many diseases, including various forms of cancer, such as leukemia and glioma (35,36). Increasing evidence suggests that dysregulated m6A regulators are intimately linked to various heart-related complications, including cardiac failure, cardiac hypertrophy, vascular calcium deposition, and pulmonary hypertension (37,38). A recent study conducted by Zhou et al. (14) showed that METTL3 enhances osteogenic differentiation in human aortic VICs via an m6A YTHDF2-dependent cascade. The findings above revealed a unique method for the modification of the m6A genes in AVC; nevertheless, the overall expression characteristics of all m6A regulators in AVC, as well as the downstream and upstream molecular processes for important m6A regulators, remain unknown.

According to recent findings, m6A modification is related to immune infiltration (associated with diseases) and inflammation. Meanwhile, research has shown a link between m6A modification-related diseases and aberrant immune regulation (39,40), where upregulated METTL3 stimulated dendritic cells (DCs) activation and maturation, whereas downregulated METTL3 suppressed T cell activation and aggregation. In addition, Han et al. (41) found that suppressing YTHDF1 increased CD8+ T cell tumor invasion and immunotherapy efficacy in a mouse tumor model. The exact mechanism by which AVC-related m6A modification leads to immune cell infiltration is unknown. To better understand immune regulation in AVC, a detailed investigation of the diverse immune infiltration patterns regulated by the m6A regulators is required. These patterns of m6A modification were found via the analysis of genomic data from AVC samples in this study. On the basis of these findings, we further investigated the association between the m6A modification pattern and immune infiltration in AVC, as well as the roles and underlying mechanisms of m6A modification patterns in the immune microenvironment (IME) in AVC. We found that m6A modification patterns were associated with diseases as well as immune cell infiltration. By investigating differentially expressed genes (DEGs) in distinct m6A modification patterns, we showed that the m6A phenotype-related hub genes influence disease features and immune cell infiltration characteristics in AVC patients. In addition, we found that the m6A phenotype-related hub genes also have high diagnostic efficiency for AVC patients. A heat map was constructed to depict the pathways in which different gene clusters might contribute. Additionally, because AVC-induced stenosis is typically only treatable surgically, we sought potential therapeutic targets for early intervention. We present the following article in accordance with the STREGA reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3627/rc).

Methods

Calcific aortic valve disease data set source and preprocessing

The Gene Expression Omnibus (GEO) database was queried for public data and complete clinical annotation. We collected the GSE102249 (42) (240 diseased samples), GSE51472 (43) (5 normal samples and 10 AVC samples), and GSE12644 (44) (10 normal samples and 10 AVC samples) datasets for further investigation (Table S1). For the GSE51472 dataset, 5 sclerotic aortic valve samples were removed because the possible other lesions of the sclerotic aortic valve affected the analysis results. To begin, we used the platform’s annotation file to perform gene name conversion for probe IDs separately. For multiple probes corresponding to the same gene, we used the average expression of the gene. Each data was normalized for subsequent analysis. Since they shared a similar platform and might be used to merge data, GSE51472 and GSE12644 were integrated into a merged dataset. The batch effect was also eliminated using the combat function of the ‘SVA’ package. We observed the results of batch effect removal using principal component analysis (PCA). A total of 21,655 target gene expression profiles were generated. Next, 23 m6A methylase regulators (i.e., writers, erasers, and readers) in 3 types were used with reference to a previous study (45) and use these 23 m6A methylase regulators for follow-up analysis (Table S2). Finally, we found that 19 m6A-related genes are expressed in samples of merged data set.

Differential expression of m6A-related genes and immune cells

Based on the Wilcoxon test, we further analyzed the differences of 19 m6A methylase regulators and 6 immune cells (IL-10, IL-11, IL-6, IL-15, IL-16, TNF) between normal and AVC samples. A P value <0.05 was considered a significant difference.

Immune cell assessment and differential expression

The combined data set was calculated based on the immune cell content of Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) samples, and the difference in immune cell content between normal and AVC samples was assessed by the Wilcoxon test.

AVC samples of GSE102249 dataset classification based on 23 m6A methylation regulators

Based on 23 m6A methylation regulators, we found that 20 m6A-related genes are expressed in samples of GSE102249 data set. We performed Consensus Clustering analysis on the 240 samples of the GSE102249 dataset, and divided them into m6A clusters A and B. The optimal k value was determined by evaluating the point of inflection for the sum of the squared errors. The decline slowed down after k=i, and k=i was selected. We used PCA to quantify the effect of typing.

Generation of the gene co-expression network and evaluation of the vital modules

To identify DEGs among m6A cluster patterns, 8,464 genes and their associated gene expression profiles were identified among the m6A cluster patterns. The DEGs were defined as those with an adjusted P<0.05 and |log2fold change (FC)| ≥1.

The “WGCNA” R-package was utilized to evaluate the top 3,000 DEGs according to the P value and modules with a high degree of connectivity. Each module was constructed by evaluating the soft threshold power β for all components. Following the establishment of the soft threshold power value, the weighted gene co-expression network analysis (WGCNA) algorithm was utilized to build the co-expression module. The information on adjacent genes was converted into topological overlap, and the network connection of the genes was determined. Based on the dissimilarity of the topological overlap matrix (TOM), hierarchical clustering was carried out to categorize genes with identical expression profiles into modules (the minimum size of gene dendrogram =20; merged highly similar modules with height =0.85). The association between clinical characteristics and module eigengenes (MEs) was utilized to identify the modules that were related to the clinical traits. The relationship between gene expression profiles and the MEs was utilized to calculate the quantitative measurement of module membership in this study. It was determined that genes were significant based on the absolute value of the connection between genes and clinical traits. Then, relevant modules with an elevated level of value for clinical characteristics were evaluated (Table S3).

Functional enrichment analysis of genes in the key modules

The R package ‘clusterProfiler’ was employed to evaluate the functional enrichment analysis on key module genes. The findings of the Gene Ontology (GO) including MF, cellular component (CC), and biological function (BF) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed. A P value <0.05 was considered statistically significant.

Identification of the key genes in the key module

To identify the interacting genes, we entered genes from important modules into the Search Tool for the Retrieval of Interacting Genes/Proteins [STRING (3); score >0.4] that predicted the protein-protein interaction network (PPIN) (Table S4). Important genes were determined as those having a significance greater than 0.3 and module membership greater than 0.6 in the modules. Finally, the hub genes of the disease-associated module and those of the PPIN were intersected, generating 7 disease-associated key genes.

AVC samples of GSE102249 dataset classification based on 7 key genes

Based on 7 key genes, we performed Consensus Clustering analysis on the 240 samples of the GSE102249 dataset again and divided them into gene cluster-1 and 2. The optimal k value was determined by evaluating the point of inflection for the sum of the squared errors. The decline slowed down after k=i, and k=i was selected. We used PCA to quantify the effect of typing.

Immune cell assessment and function enrichment

The combined data set was calculated based on the immune cell content of CIBERSORT samples, and the difference in immune cell content between 2 subgroups was assessed by the Wilcoxon test. The R package ‘clusterProfiler’ was employed to evaluate the functional enrichment analysis on key module genes. We performed GO findings on gene cluster-1 and 2, including MF, CC, and BP, and KEGG analysis, respectively, and created a thermal figure display. A P value <0.05 was considered statistically significant.

Statistical analysis

The R software (The R Foundation for Statistical Computing, Vienna, Austria) was employed to evaluate statistical analyses. The P values (bilateral) <0.05 were considered statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

The landscape of m6A regulators and the tissue IME in AVC

Our study used datasets from the GEO (GSE51472 and GSE12644, Table S1). We used PCA analysis to observe the combined dataset batch effects after removing the combined results (Figure S1). We identified 23 m6A regulators based on previously published data (45), including 5 “writers”, 2 “erasers”, and 12 “readers” (Table S2). However, the gene profile of the combined dataset contained only 19 m6A methylation regulators. illustrates the positions of m6A regulators on the chromosomes. To identify the impact of the m6A regulators on the development of AVC, the transcriptional level of m6A regulators in AVC and healthy tissues was determined. The data indicated that the expression of “readers” (i.e., FMR1, HNRNPA2B1, HNRNPC, LRPPRCYTHDC1, YTHDC2, and YTHDF1) was significantly lower in the AVC tissues than in the normal tissues (P<0.01) (). The expression of “writers” (i.e., METTL14) and “erasers” (i.e., FTO) was also significantly lower in the AVC tissues than in the normal tissues (P<0.01) (). The underlined findings indicated that m6A regulators showed a high degree of heterogeneity in their expressional landscapes between normal and AVC patients, implying that differential m6A regulator expression might influence disease progression. Furthermore, we evaluated the characteristics of immune cell infiltration in AVC tissues and identified inflammatory factor IL-10 (P<0.05) () and M0 and M2 macrophages (P<0.05) (). Compared to normal tissues, AVC tissues were enriched in IL-10 and M0 macrophages, which indicated an upregulation of the inflammatory factor. Conversely, M2 macrophages were more abundant in normal tissues, implying that infiltration of M2 macrophages may inhibit the progression of calcification. These findings suggested that changes in the aortic valve tissue’s IME might be the predominant pathogenic factor for AVC, and the infiltration of local M0 macrophages and the inflammatory factor IL-10 might promote AVC progression.

Figure 1

The landscape of m6A regulators and IME in AVC. (A) The GSE51472 and GSE12644 cohorts were used to map the position of m6A regulators on 24 chromosomes. (B,C) The expression of 24 m6A regulators in AVC and healthy tissues. (D,E) The variations in immune cell infiltration between AVC and healthy tissues. (F,G) The expression of inflammatory factors in healthy and AVC tissue panel. *P<0.05; **P<0.01; ns, not significant. AVC, aortic valve calcification; m6A, N6-methyladenosine; IME, immune microenvironment.

The m6A methylation modification patterns are facilitated by 19 m6A regulators in AVC patients

The m6A modification patterns were characterized using the R package ‘ConsensusClusterPlus’ based on the 20 m6A regulators. The fraction of unclear clustering assessment and the consistency of m6A regulator expression were used to determine the optimal clustering stability (k=2–10), and k=2 was obtained (). We categorized the modification patterns of 240 AVC samples in GSE102249 using unsupervised clustering and identified 23 m6A regulators, including 5 “writers”, 2 “erasers”, and 9 “readers” (Table S2). Finally, 2 varied modification patterns were found, namely, 124 cases in cluster A and 116 cases in cluster B (Table S5). The expression of these m6A regulators could be used to differentiate the AVC samples (). To investigate the relationship between “writers”, “erasers”, and “readers”, we analyzed pairwise correlations between the expression of these m6A regulators in AVC and discovered that positive correlations were infrequent (). We found a considerable association among “writers”, “erasers”, “readers”, and m6A regulators with the same role showing considerable linkage in expression. Cases of AVC with high expression of HNRNPA2B1 (a “reader” gene) exhibited downregulation of YTHDF1 (a “reader” gene), YTHDF2 (a “reader” gene), LRPPRC (a “reader” gene), FTO (an “eraser” gene), RBMX, and IGFBP1, and upregulation of WTAP (a “writer” gene) and YTHDC1 (a “reader” gene), and did not affect the expression of other genes. Cases of AVC with elevated expression of YTHDC2 (a “reader” gene) showed downregulation of YTHDF2 (a “reader” gene) and IGFBP1, and upregulation of METTL14 (a “writers” gene), YTHDC1 (a “reader” gene), HNRNPC (a “reader” gene), FTO (an “eraser” gene), and RBMX, but it did not affect the expression of other genes. The expression of the “writer” genes METTL14, RBM15, and RBM15B presented a similar trend as that of FTO (an “eraser” gene). Furthermore, variations in FMR1 and ZC3H13 expression had no impact on the expression of other genes. Hence, crosstalk between distinct m6A regulators could be critical for generating m6A modification patterns in AVC patients. We evaluated the expression of 20 m6A regulators in various m6A subtypes and disease subtypes post unsupervised clustering () and found that the expression of YTHDC2, IGF2BP1, RBMX, and FTO increased in the cluster B compared to their expression in the cluster A, and the expression of HNRNPA2B1 increased in the cluster A compared to its expression in the cluster B (P<0.05).

Figure 2

The methylation modification patterns (mediated by m6A) were associated with clinical features and the tissue IME in AVC. (A-C) CC evaluated 2 subgroups based on the expression of m6A regulators. (D) Evaluation of the pairwise correlations for the expression level of 24 m6A regulators. (E,F) The expression of the 24 m6A regulators in various m6A and subtypes of the disease. (G-J) The variations in tissue inflammatory factor panel (H) and immune cell infiltration panel (I) between m6A cluster 1 and m6A cluster 2. *P<0.05; **P<0.01; ***P<0.001; ns, not significant. CDF, cumulative distribution function; PC, principal component; m6A, N6-methyladenosine; IME, immune microenvironment; AVC, aortic valve calcification; CC, cellular component.

Regulation of tissue IME by m6A modification patterns in AVC

We used CIBERSORT to analyze the variations in immune cell infiltration between the 2 m6A alteration patterns. We found that cluster A was strongly related to immunological activation. Subsequent analyses of immune infiltration revealed that IL-10 expression was considerably elevated (P<0.001) in cluster A (), and cluster A was considerably enriched in infiltrating innate immune cells, including memory B cells (P<0.01), CD8 T cells (P<0.05), T follicular helper cells (P<0.05), monocytes (P<0.05), M0 macrophages (P<0.01), and resting DCs (). Meanwhile, naïve B cells (P<0.01), M2 macrophages (P<0.001), activated DCs (P<0.05), and neutrophils (P<0.05) were significantly abundant in cluster B. Additionally, we observed increased expression of TNF (P<0.01) and IL-11 (P<0.05) in cluster B. Significant differences in immune cell and inflammatory factor infiltration between the 2 m6A clusters suggested that m6A methylation modifications might affect the immunological milieu of the AVC tissue.

Features of key m6A modules

We utilized the “limma” package to identify the DEGs in the 2 m6A clusters. Subsequently, the first 3,000 DEGs with the lowest P value were filtered out. Following that, a WGCNA was built. After evaluating the eigengenes of each module, clustering the modules, and combining adjacent modules into new modules (), 5 modules were acquired (). The brown module showed the strongest association with the m6A cluster classification phenotype (Table S3). The R package clusterProfiler was utilized to assess the genes in the brown module for function and cascade information. As illustrated in , the genes in the brown module were significantly enriched in BP (), including response to estradiol, extracellular matrix (ECM), extracellular structure organization, response to drugs, smooth muscle contraction, and muscle contraction (). The CCs were significantly enriched in myosin filament, collagen-comprising ECM structural constituent, and apical part of the cell (). The MFs were significantly enriched in protein phosphatase and phosphatase regulator activity, growth factor binding, and ECM structural component (). Finally, the KEGG pathways were significantly enriched in vascular smooth muscle contraction and ECM ().

Figure 3

Features of the significant m6A module. (A,B) Hierarchical clustering was conducted to construct a dendrogram of genes. (C,D) The link between the MEs and modification patterns (mediated by m6A). (E) The GO terms BP, CC, and MF, were used to functionally annotate the genes in the vital m6A module. (F) The KEGG pathway terms were used to annotate the genes in the significant m6A module. GEO, Gene Expression Omnibus; BP, biological process; CC, cellular component; MF, molecular function; m6A, N6-methyladenosine; MEs, module eigengenes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Seven key genes were linked with the AVC characteristics

To identify the key genes in important module, a PPIN was constructed, followed by analysis of 25 associated genes in the brown module (). After excluding non-interacting genes, the PPIN identified 13 hub genes. Genes with a gene significance above 0.3 and module membership above 0.6 were evaluated as important genes in the brown module, with 10 genes in total (). Next, the 10 important genes and 13 hub genes were intersected to obtain 7 key genes (). We found considerable differences in the expression of the 7 key genes (associated with diseases) between normal and AVC tissues using the previously merged dataset (GSE51472 and GSE12644) and found that the expression of ITGA8, MYH11, and SMOC2 was considerably elevated in AVC tissues relative to that in the normal tissues, yet the expression of the other 4 genes was indifferent (; P<0.05). Furthermore, we found that 7 key genes exhibited significant differences in the m6A cluster categorization based on the GSE102249 dataset. The 7 key genes were expressed at a greater level in cluster A than in cluster B (; P<0.05). Additionally, compared to cluster A, MYH11, SMOC2, and RCAN2 were considerably higher in cluster B, but ITGA8 expression was significantly higher in cluster A (; P<0.05).

Figure 4

The pathological features of AVC were influenced by seven m6A phenotype-associated hub genes. (A) Evaluation of genes in the vital module by constructing the PPIN. (B) The significant genes in the brown module. (C) The intersection of the hub genes of the interaction network with the hub genes of the m6A-associated module. (D,E) The expression of seven genes between healthy and AVC tissues. (F,G) The expression of the seven genes in different m6A cluster classification. *P<0.05; **P<0.01; ***P<0.001. PPI, protein-protein interaction; AVC, aortic valve calcification; m6A, N6-methyladenosine; PPIN, protein-protein interaction network.

Consistent cluster analysis of GSE102249 based on 7 key genes

The AVC patients were categorized into various genomic subtypes using unsupervised clustering analysis based on the expression of 7 key genes. We found that the gene clusters could be separated (; Table S6). The majority of individuals with the cluster B subtype were assigned to gene cluster-1 (). Considerable variations were observed in gene cluster classification for 7 key genes (P<0.001). The expression of RCAN2 was lower in gene cluster B. The remaining 6 key genes were expressed at significantly higher levels in gene cluster-2 than in gene cluster-1 (; P<0.05). Additionally, we evaluated the different infiltration immune cells in the 2 subgroups () and found considerable variations in CD8 T cells, memory CD4 T cells (activated), T follicular helper cells, monocytes, M1 macrophages, M2 macrophages, and eosinophils (P<0.01, 0.001, 0.001, 0.001, 0.01, 0.001, and 0.05, respectively). The M1 macrophages (P<0.01), CD8 T cells (P<0.01), memory-activated CD4 T cells (P<0.001), and T follicular helper cells (P<0.001) were highly abundant in gene cluster-1. However, monocytes, M2 macrophages, eosinophils, and neutrophils (P<0.001, 0.001, and 0.05, respectively) were substantially more abundant in gene cluster-2 ().

Figure 5

CC for m6A phenotype-associated hub genes linked with AVC patient characteristics. (A) CC was used to cluster the samples. (B) Two subgroups were determined based on the expression of the 7 m6A phenotype-associated central genes. (C) Alluvial diagram indicating the variations in m6A clusters, gene clusters, and disease subtype. (D,E) The expression of 7 hub-genes (associated with m6A phenotype) in gene cluster-1 and gene cluster-2. (F,G) The immune cell infiltration in gene cluster-1 and gene cluster-2. *P<0.05; **P<0.01; ***P<0.001. CDF, cumulative distribution function; PC, principal component; m6A, N6-methyladenosine; AVC, aortic valve calcification; CC, cellular component.

Diagnostic efficacy and pathway analysis of gene cluster subgroups

To investigate the diagnostic abilities of the 7 key genes for AVC, we performed receiver operating characteristic (ROC) curve analysis and found that SMOC2 might be the core gene with the best diagnostic effectiveness (). Also, the pathways in which distinct gene clusters may contribute were identified. The KEGG pathway heat map revealed that ubiquitin-mediated proteolysis, lysine degradation, and aminoacyl-tRNA biosynthesis were significantly enriched in gene cluster-1; olfactory transduction, primary bile acid biosynthesis, and arachidonic acid metabolism were significantly enriched in gene cluster-2. The GO pathway heat map showed that the sin3 type complex, ATPase complex, GOBP covalent chromatin modification, and protein methylation were significantly enriched in the significant gene cluster-1; cysteine-type endopeptidase inhibitor, endopeptidase regulator, and peptidase regulator activity were significantly enriched in the significant gene cluster-2 ().

Figure 6

Diagnostic efficacy of the m6A phenotype-related hub genes, a heat map of the pathways involved in different gene clusters, and prediction of small-molecule drugs. (A) Diagnostic efficacy of m6A phenotype-related hub genes in the merged data set. (B,C) Heat map of the related pathways of different gene clusters. TPR, true positive rate; FPR, false positive rate; AUC, area under the curve; KEGG, Kyoto Encyclopedia of Genes and Genomes; CC, cellular component; BP, biological process; MF, molecular function; GO, Gene Ontology; m6A, N6-methyladenosine.

Discussion

As most studies on disease-related m6A modification have focused on many tumors, the impact of m6A regulators on the IME in AVC tissue may help us to understand disease development and disease-related immune infiltration and to establish more effective therapeutic regimens. Xu et al. (46) found that IL-18 stimulates M1 macrophage-mediated erythrophagocytosis and erythrocyte breakdown by activating HO-1 and FPN through p38 and Erk1/2 in a calcific microenvironment. Additional, recent findings demonstrate an adaptive immune response functions in AVC by activating circulating CD8+ cells, their clonal proliferation and differentiation to the memory-effector subtype, and exchanges of these T-cell clones between the peripheral blood and valve (47). In addition, IFN-γ from CD8+CD28 T cells may also promote calcification (47). Activated effector T cells and their cytokine production may impair osteoclast activity, increasing valvular calcium buildup. There is strong evidence that m6A regulates transcription and protein expression through splicing, translation, degradation and export, affects the activity of cancer cells and/or stromal cells and characterizes the IME (48). At present, the relationship between m6A and the IME is mainly studied in tumors, while the study on AVC is very limited. Therefore, in this study, we can further explore new therapeutic targets by analyzing the relationship between m6A-related genes and immune cells in tissues by bioinformatics.

The expression profile of common m6A regulators in calcified and normal aortic tissues was first elucidated. In the AVC samples, 2 “writers” (METTL14 and RBM15B) and 1 “reader” (HNRNPC) exhibited statistical significance, indicating their important role in AVC. The m6A methyltransferases contribute to the catalysis of the m6A modification (37). Inflammatory development is mediated by the m6A regulators. However, only 1 study found a link between m6A regulators and AVC disease development, and it was reported that METTL3 (a “writer” gene) is highly expressed in human calcified aortic valves relative to normal valves and promotes osteogenic differentiation of human aortic VICs via the m6A-YTHDF2-dependent pathway (14). A homolog of METTL3, METTL14 is a novel family of m6A “readers” (49), which is an RNA-binding scaffold that can recognize substrate RNAs (50) and enhance the proliferative and invasive potential of atherosclerotic vascular endothelial cells (51), which decrease vascular healing function (52). Histopathological analysis revealed that the initial phase of AVC is characterized by an aggressive inflammatory response similar to that seen in atherosclerotic lesions, with processes ranging from lipid deposition to infiltration of inflammatory cells including macrophages and T cells, and ultimately, destruction of the basement membrane (53). Therefore, METTL14 may affect the progression of AVC through an atherosclerosis-like pathway. Additionally, calcification of valve tissue might be the result of osteogenic differentiation of VICs and lead to osteogenic calcification and the formation of bone-like structures (54). Tian et al. showed that the m6A modification in the formation of AVC might be a key regulator involved in the osteogenic or osteoclastic processes of bone. For example, significantly elevated expression of m6A methyltransferases (METTL3 and METTL14) was found in bone mesenchymal stem cells undergoing osteogenic induction (55). Thus, METTL14 might be closely associated with the aortic valve osteocalcification phase. However, Zhou et al. (14) found no significant difference in METTL14 between normal and disease tissues in the progression of AVC. Interestingly, in our analysis, METTL14 showed low expression in AVC tissues; thus, the mechanism of action of METTL14 is not yet known, and it may be related to the regulation of various factors, which needs to be further explored. Additionally, Zhou et al. (14) found that by suppressing TWIST1 via an m6A-YTHDF2-dependent mechanism, METTL3 enhances osteogenic differentiation of human aortic VICs. In this study, we found that YTHDC1, YTHDC2, and YTHDF1 differed significantly between normal and AVC tissues; YTHDF2 and the genes mentioned above are all members of the YTH family. Thus, the YTH family might also contribute to the progression of AVC. In vitro, FTO efficiently eliminated the m6A modification of mRNA (56). Recently, some studies have found strong associations between FTO and hypertension (57), diabetes (58,59), obesity (60,61), non-alcoholic fatty liver disease (62,63), atherosclerosis (64), myocardial infarction (38,65), ischemic reperfusion injury (65), cardiac hypertrophy (65), and heart failure (65). Thus, FTO considerably contributes to the progression of cardiovascular diseases. It might be that FTO is also involved in the process of AVC. Our study revealed that the expression landscape of m6A regulators in normal and AVC patients is highly heterogeneous. The precise molecular mechanisms by which m6A regulators modulate AVC progression, however, remain unknown.

As a reversible RNA methylation process, m6A is inserted by the “writer” and deleted by the “eraser”, suggesting that they may perform antagonistic functions (66). However, our study indicated a positive correlation between the m6A writer and eraser and the frontiers of cell and developmental biology. The FTO group contains the “author” genes RBM15, RBM15B, METL14, and METL16 (a kind of “eraser” gene). Previous studies have indicated that the m6A writers that increase the amount of RNA m6A are oncogenes (67,68). The targets of “erasers” and “writers” described in various studies are distinct, suggesting that disease progression involves several processes and regulatory factors other than the “writer” and “eraser”. Furthermore, the post-transcriptional role of m6A “readers” might partly explain why the “eraser” and “writer” have contradictory roles. Future studies should determine the influence of m6A “writers” and “erasers” on the same RNA to evaluate the relative contribution of other regulatory factors in m6A modification.

In recent decades, m6A regulators have been identified in disease-associated immune cell filtration. By removing the m6A alteration, ALKBH5 reduces the antiviral natural response (69). Additionally, METTL3 might be used as an anti-inflammatory target to methylate the STAT1 mRNA and promote M1 macrophage polarization (70). The m6A reader YTHDF2 modulates the inflammatory reaction induced by lipopolysaccharides in macrophages (71). Knockout of FTO induces the reduction of PPAR-γ and STAT1 transcripts by triggering mRNA degradation (YTHDF-mediated) and attenuating M1/M2 activation (72). The immuno-infiltration analysis identified that M0 macrophages were expressed at a higher level in AVC tissues than in healthy tissues, that demonstrating that m6A regulators might affect inflammatory processes through macrophage recruitment. Furthermore, Hu et al. showed that the absence of FTO resulted in the conversion of macrophages to M1-type pro-inflammatory macrophages (73). These results suggested that FTO agonists might be potential anti-inflammatory targets. We found that FTO was expressed at low levels in AVC tissues and highly expressed in M0 macrophages, which is consistent with the results of Hu et al. Additionally, the data from AVC samples were typed by m6A expression matrix and then subjected to immuno-infiltration analysis. We found that cluster A, with elevated expression of m6A regulators, was considerably more abundant in innate immune cell infiltration, including memory B cells, CD8 T cells, T follicular helper cells, monocytes, M0 macrophages, and resting DCs, whereas, cluster B, with decreased expression of m6A regulators, was enriched in the infiltrative pattern of innate immune cells, including naïve B cells, M2 macrophages, activated DCs, and neutrophils. Thus, m6A methylation modification can considerably elevate the expression of innate immune cells linked to inflammatory processes. Cluster A had higher IL-10 expression, whereas cluster B had higher TNF and IL-1 expression. This indicated that modifying m6A methylation might affect tissue inflammatory factors in AVC. Our results on m6A regulators, however, are restricted to bioinformatics, and more research is needed to determine the exact molecular cascade through which the m6A regulator triggers immune cell infiltration in AVC. From the perspective of a biological process, GO enrichment analysis revealed enrichment in response to estradiol, ECM organization, extracellular structure organization, drug response, smooth muscle contraction, and muscle contraction. It was primarily concentrated in myosin filaments, collagen-containing ECM, and the apical region of the cell for CCs. It was primarily concentrated in phosphatase regulator activity, protein phosphatase regulator activity, growth factor binding, and ECM structural components for MFs. Further, AVC was found to be closely related to vascular smooth muscle contraction and the ECM-receptor interaction signaling pathway in the KEGG pathway analysis, which was consistent with the results of another study (74).

Next, we identified 7 disease-related hub genes utilizing key modules and discovered that the expression of ITGA8, MYH11, and SMOC2 was considerably higher in AVC tissues than in the healthy tissues. The transmembrane cell-surface receptor, ITGA8, is a member of the alpha integrin family (75). A strong association between ITGA8 and tumor development has been reported previously. For example, low ITGA8 expression is related to a poor overall survival prognosis in patients with clear cell renal cell carcinoma (76). High expression of ITGA8 causes an epithelial-mesenchymal transition in patients with early relapsed multiple myeloma that results in increased motility and invasive potential of multiple myeloma cells (77). The smooth muscle myosin, MYH11, is represented by the MYH11 gene. As a member of the myosin heavy chain family (78), MYH11 is a contractile protein that slides through actin filaments to contract muscles via adenosine triphosphate hydrolysis (79,80). Deficiency of MYH11 leads to the loss of vascular smooth muscle cells, disorganization, and hyperplasia in aortic tissues, which is one of the mechanisms causing thoracic aortic aneurysms and dissections (81,82). Thus, the association of MYH11 with AVC remains strong. The SMOC2 gene is a SPARC family member that is expressed in many tissues (83). Similarly, SMOC2 plays a role in several biological processes, including adhesion and migration of cells, angiogenesis (84), fibrosis (85), and cell proliferation (86). The involvement of SMOC2 in fibrosis during AVC needs further investigation. Our findings suggest that 7 important genes may be potential therapeutic targets for AVC. These targets may involve the development of new AVC treatments. However, further studies are needed to confirm their role in AVC progression.

The patients were categorized into various genomic subgroups based on the expression of the 7 hub genes. The CD8 T cells, memory CD4 T cells (activated), T follicular helper cells, and M1 macrophages were considerably more abundant in gene cluster-1, whereas M2 macrophages, monocytes, neutrophils, and eosinophils were considerably more abundant in gene cluster-2. The 7 key genes might play a regulatory role for immune cells in different subtypes of AVC. However, further research is needed to follow up. Furthermore, we identified the pathways in which different gene clusters may participate, which suggested that the pathways of enrichment for different subgroups of AVC were inconsistent and might be strongly related to the degree or stage of AVC lesions.

AVC has been studied for years. Ongoing clinical research examine therapies targeting recently revealed pathways. Denosumab, a human monoclonal antibody against RANKL (a member of the TNF cytokine family) that inhibits its binding to RANK (a transmembrane protein expressed on marrow stromal cells and preosteoclasts), mimics the action of osteoprotegerin, presumably limiting osteoblastic development in aortic valve tissue (87). Bisphosphonates inhibit vascular and AVC (88). They may also block AVC-related inflammatory cytokines (89). Repurposed drugs may also target AVC. Choi et al. (90) studied the effect of DPP-4 inhibitors, such as sitagliptin, on VIC calcification. NO depletion in valvular endothelial cells mimicked endothelial dysfunction, activating NF-κB pathway in VICs. AVC is a complex disease involving metabolic factors, immune cells, endothelial cells and mesenchymal cells. Inadequate in vivo model systems and in vitro method heterogeneity hinder confidence in past information. We need new methods with higher reproducibility that will increase our understanding of this disease and lead to new treatment options.

Although our findings identified potential subgroups of patients with AVC based on molecular analyses of m6A regulators and revealed the characteristics of each subtype, gene expression changes are not necessarily genetically driven and may be influenced by both genetic and environmental factors. Several additional limitations should be acknowledged. First, some clinical follow-up information was unavailable for the samples; as a result, we did not consider several important factors, such as the presence of patient comorbidities, when differentiating between subgroups. The consistency of the AVC subgroups should then be confirmed through additional analyses. The results were obtained solely through bioinformatics analysis, and they should be confirmed through additional experiments.

Conclusions

We identified a link between m6A modulators, immune infiltration, and key genes in AVC. The identification of genetic subgroups of patients with AVC has improved our understanding of the pathogenesis of AVC and has facilitated the development of potential methods for disease diagnosis, classification, and treatment.

The article’s supplementary files as