Identification of hub genes and key signaling pathways by weighted gene co-expression network analysis for human aortic stenosis and insufficiency

Background: Human aortic valve stenosis (AS) and insufficiency (AI) are common diseases in aging population. Identifying the molecular regulatory networks of AS and AI is expected to offer novel perspectives for AS and AI treatment.

Methods: Highly correlated modules with the progression of AS and AI were identified by weighted genes co-expression network analysis (WGCNA). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed by the clusterProfiler program package. Differentially expressed genes (DEGs) were identified by the DESeqDataSetFromMatrix function of the DESeq2 program package. The protein-protein interaction (PPI) network analyses were implemented using the STRING online tool and visualized with Cytoscape software. The DEGs in AS and AI groups were overlapped with the top 30 genes with highest connectivity to screen out ten hub genes. The ten hub genes were verified by analyzing the data in high throughput RNA-sequencing dataset and real-time PCR assay using AS and AI aortic valve samples.

Results: By WGCNA algorithm, 302 highly correlated genes with the degree of AS, degree of AI, and heart failure were identified from highly correlated modules. GO analyses showed that highly correlated genes had close relationship with collagen fibril organization, extracellular matrix organization and extracellular structure organization. KEGG analyses also manifested that protein digestion and absorption, and glutathione metabolism were probably involved in AS and AI pathological courses. Moreover, DEGs were picked out for 302 highly correlated genes in AS and AI groups relative to the normal control group. The PPI network analyses indicated the connectivity among these highly correlated genes. Finally, ten hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) in AS and AI were found out and verified.

Conclusion: Our study may provide the underlying molecular targets for the mechanism research, diagnosis, and treatment of AS and AI in the future.

Introduction

Heart failure is a severe terminal stage of all kinds of cardiovascular diseases, which include hypertension (1), coronary heart disease (2), myocardial infarction (3), valvular heart disease (4), and cardiomyopathy (5). As the important cause of heart failure, valvular heart diseases consist of stenosis or insufficiency with specific pathophysiology among the four cardiac valves (aortic valves, mitral valves, tricuspid valves and pulmonary valves) (6). In the globe, valvular heart diseases have increasingly become the important contributor to cardiovascular morbidity and mortality according to the epidemiologic studies, which have resulted in serious social burden and economical cost on valvular heart diseases diagnosis and treatment (7). The prevalence of valvular heart diseases gradually increases with age of clinical patients (8). In terms of aortic valve lesion, aortic valve stenosis (AS) and aortic valve regurgitation (AR) are the highly popular valve lesions among various valvular heart diseases (9), the morbidities of AS and AR were 0.7% and 0.2% with the age 55–64 years, 1% and 1.3% for the age-bracket of 65–74 years, and 2% and 2.8% after 75 years old, respectively (8). Thus, due to the large aging population around the world, the aortic valvular heart diseases are still the important public health problem.

AS and AR are a kind of common aortic valve diseases, characterized by aortic valve opening area reduction or aortic valve insufficiency (AI), respectively. Currently, numerous studies have reported the pathological features and molecular mechanisms about aortic valve damage (10–13). The well-known etiologies for AS include aortic valve degeneration, rheumatic aortic stenosis, congenital valve defects, systemic inflammatory diseases, endocarditis, and many other conditions (10). Whereas, the major causes of AI are made up of various pathological changes of aortic valves, such as leaflet abnormalities, rheumatic fever, myxomatous degeneration, infective endocarditis, etc (12). Although, the current available reports have uncovered the molecular mechanisms for pathological processes of AS and AI, which include but not limit to fibro-calcific remodeling, osteogenic differentiation, lipid accumulation, inflammation, angiogenesis and hemorrhage, disorganization and remodeling of the valvular extracellular matrix (ECM) (10, 13). The vital molecules for the regulation and indication of AS and AI pathological courses still need to be further investigated. Hence, using high throughput sequencing techniques and identifying the key regulatory or indicative molecules for AS and AI may provide a feasible strategy for the diagnosis and treatment of AS and AI from the microscopic molecular viewpoints.

In this study, we analysed the expression profiles of human aortic valve samples of aortic valve stenosis (AS) and aortic insufficiency (AI) by systematic bioinformatics approaches of weighted gene co-expression network analysis (WGCNA). We constructed the gene co-expression modules by WGCNA algorithm and screened the highly correlated modules with the degree of AS, degree of AI, and heart failure, which included orange, steelblue, darkgreen, and grey60 modules. Furthermore, we selected highly correlated genes in indicated modules and performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The differentially expressed genes (DEGs) were identified and intersected with those 30 genes possessing high connectivity among the highly correlated genes from various modules to screen out hub genes. Finally, the mRNA expression values of ten hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) were validated by analyzing results of high throughput RNA sequencing from AS and AI aortic valve samples and by examining the mRNA expression levels of human AS and AI aortic valve tissues. These results may provide an avenue for the diagnosis and treatment of AS and AI in the future.

Materials and methods

High-throughput data acquisition and preprocessing

The high-throughput RNA-sequencing datasets were acquired from the public Gene Expression Omnibus (GEO) database with the accession number GSE153555, which contained the gene expression data from 5 human normal control (NC) aortic valves, 5 human aortic stenosis (AS) aortic valves, and 5 human aortic insufficiency (AI) aortic valves, each individual contained 2 biologically repeated aortic valve high-throughput RNA-sequencing results. The gene expression Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values and count values were analysed by R software (Version 4.1.2). Clinical traits, including age, sex, body mass index (BMI), degree of AS and AI, left ventricular ejection fraction (LVEF), and disease history (diabetes, hypertension, coronary heart disease, and heart failure), for each sample were collected from the Series Matrix Files in the GEO database with GSE153555 number (14). The average gene expression FPKM values of 30 samples were calculated and ranked by size, and the top 6,000 genes with the highest average expression were screened out and used for weighted gene co-expression network analysis (WGCNA) computation. The FPKM values of the 6,000 genes from 30 samples were subjected to log₂(FPKM + 1) conversion followed by samples hierarchical clustering to eliminate 2 outlier samples (GSM4647040 and GSM4647041) using the hclust function in the R software (Version 4.1.2).

Co-expression module construction of AS and AI by WGCNA algorithm

The WGCNA co-expression module construction was conducted as previously described (15). The soft threshold β power value for WGCNA module construction was computed by pickSoftThreshold function in the WGCNA program package. The adequate β value 16 was picked out once Scale Free Topology Model Fit, signed R² value was ≥0.8. Then, the adjacency matrix and topological overlap matrix (TOM) were constructed using the power value 16. The co-expression modules were constructed and merged for those modules with similar expression profiles by step-by-step network construction methods. The correlation analysis among the indicated modules were performed by calculating the eigengenes, which were defined as the principal component 1 (PC1) of principal component analysis (PCA) for gene expression values in the indicated modules.

Correlation analysis between co-expression modules and clinical traits

The correlation analyses between modules and clinical traits were performed using the eigengenes of the corresponding modules and clinical traits data to screen out the highly correlated modules with indicated clinical traits. To screen the highly correlated genes in AS and AI pathological courses, these modules relevant to degree of AS, degree of AI, and heart failure were further analysed. These modules with correlation coefficient more than 0.6 and P value less than 0.05 were regarded as highly correlated modules. The correlation coefficients of 6,000 genes expression values and indicated clinical traits were defined as gene significance (GS). For the associations of each module with the genes involved in the WGCNA process, module membership (MM) was defined as the correlation of module eigengenes and gene expression levels. The scatterplot of Gene Significance vs. Module Membership in the indicated module was plotted. These scatterplots with correlation coefficients more than 0.5 and P value less than 0.05 were selected. These 302 genes in indicated scatterplot with MM > 0.8 and GS > 0.8 were regarded as the highly correlated genes with corresponding module or trait, respectively, and were selected for subsequent analysis.

Gene ontology (Go) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses

We performed GO and KEGG enrichment analyses for the highly correlated genes identified by the above procedures using the clusterProfiler program package. Firstly, the ENSEMBL number for each gene was transformed into the ENTREZID number using the bitr function. GO analyses including biological processes (BP), molecular functions (MF), and cellular components (CC) terms were conducted. These terms in GO and KEGG enrichment analyses with P value less than 0.05 were screened out and considered as significant terms in AS and AI pathological processes. The top ten terms in GO and KEGG enrichment analyses were selected for visualization and further analyses.

Identification of differentially expressed genes (DEGs)

The DEGs analysis using the high throughput RNA-sequencing data was performed to evaluate the gene expression situation among the indicated groups using the DESeqDataSetFromMatrix function of the DESeq2 program package. Firstly, the results of DEGs analyses for total genes in high throughput RNA-sequencing data were obtained from the AS and AI group. The indicated DEGs results for the 302 highly correlated genes were screened out. These genes in AS or AI group with |log₂FoldChange| ≥ 1 and adjust P value <0.05 relative to the normal control group, were considered as the DEGs. These genes with log₂FoldChange ≥ 1 were defined as upregulated (UP) genes. These genes with log₂FoldChange ≤ −1 were defined as downregulated (DOWN) genes. These genes with −1 < log₂FoldChange < 1 were defined as unchanged (NOT) genes. The volcano plots and heat maps were plotted to visualize DEGs.

Protein-protein interaction (PPI) network analysis

PPI network analysis was performed using an online network tool STRING (https://cn.string-db.org/, version 11.5) (16). The 302 highly correlated genes were imported into STRING. The TSV file contained 302 highly correlated genes was downloaded and the PPI network was visualized by Cytoscape software (version 3.9.0) (17). The CytoHubba plug-in was used for identifying genes with high connectivity ranked by Betweenness (18). The top 30 genes with highest connectivity were picked out for further analysis. The DEGs of AS and AI group were intersected concurrently with the top 30 genes from PPI network analysis to identify these genes with significantly changed expression and high connectivity. These genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) that met the above conditions were identified as the hub genes.

Human aortic valves sample collection and grouping

Human aortic valve samples were obtained from patients with pure AS or AI. These aortic valves samples from heart transplantation receptors or aortic dissection patients without definite lesions in aortic valves were used as normal control (NC) and mild aortic valve samples. A total of 35 aortic valve samples including 5 samples from normal control (NC) patient, 15 samples from aortic stenosis (AS) patient and 15 samples from aortic insufficiency (AI) patient were included in this study. Doppler echocardiography was used for evaluation of AR or AS severity (19, 20). The specimens were classified into 4 groups containing normal control, mild, moderate, and severe. The characteristics of aortic valves used in this study are shown in Table 1. All procedures involving human aortic valves samples conformed to the principles outlined in the Declaration of Helsinki. Exemption from informed consent for patients and human sampling procedures were approved by the Human Research Ethics Committees of Tongji Hospital of Huazhong University of Science and Technology (21).

TABLE 1

Table 1. Characteristics of clinical samples used in this study.

Real-time PCR

Real-time PCR analyses of mRNA levels in human aortic valve samples were performed as previously described (21). Briefly, the total RNA was extracted from 5 normal control aortic valve samples, 15 AS aortic valve samples, and 15 AI aortic valve samples using TRIzol reagent (15596018, Thermo Fisher Scientific) and reverse transcribed into cDNA using the Transcriptor HiScript III RT SuperMix for qPCR (+gDNA wiper) (R323-01, Vazyme). The quantitative analyses of ten hub genes mRNA expression levels were determined by real-time PCR assay using SYBR (Q311-02, Vazyme). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used for the internal reference. Primers used in this study were listed in Table 2.

TABLE 2

Table 2. Primers for the real-time PCR assays in this study.

Statistical analysis

Statistical analysis was conducted by SPSS 23.0 software. All data are presented as the mean ± SD. Non-parametric Kruskal-Wallis H test was used for comparisons among multiple groups. P value <0.05 was considered to be statistically significant.

Results

Construction of WGCNA co-expression modules

The data analysis process used in this study was depicted in the flow diagram (Figure 1A). High throughput RNA-sequencing data were preprocessed using the R software (Version 4.1.2). The average FPKM expression values for a total of 48,162 genes were calculated and ranked by size. The top 6,000 genes with the highest average FPKM expression values in the datasets of the 30 samples were chosen for WGCNA computation. Sample hierarchical clustering was performed with hclust function and the height 60 was set as the threshold to screen outlier samples (Figure S1A). Two outlier samples GSM4647040 and GSM4647041 were identified and eliminated from all samples (Figure S1A). Before network construction and module detection, the clinical traits related to the sample dendrograms were visualized as the heatmap in Figure S1B. Finally, 28 samples with 6,000 genes were selected for WGCNA module construction. To screen out the suitable soft threshold power value used for WGCNA algorithm, we set an indicated range for power values and the power value 16 was picked out by pickSoftThreshold function in WGCNA package (Figure 1B). The WGCNA module construction was conducted, and those modules with similar expression profiles were merged. As depicted by the gene dendrograms, total eight modules were finally constructed, which included darkgrey (1,978 genes), lightcyan (1,158 genes), darkorange (242 genes), orange (170 genes), steelblue (35 genes), darkgreen (706 genes), grey60 (649 genes), and grey (1,062 genes) modules (Figure 1C). Those genes uncorrelated with other modules were assigned to grey module. The eigengene dendrogram and eigengene adjacency heatmap were plotted to exhibit the associations among the modules by eigengenes (Figure 1D).

FIGURE 1

Figure 1. Construction of WGCNA co-expression modules. (A) The flow diagram for our WGCNA modules construction procedure in this study. (B) Analysis of scale-independence index (left panel) and mean connectivity for various soft threshold β power values, the most suitable β power value 16 was screened out. (C) Clustering dendrograms of 6,000 genes and module construction and merging according to the similar expression profiles by WGCNA algorithm. Total of eight co-expression modules were obtained including darkgrey, lightcyan, darkorange, orange, steelblue, darkgreen, grey60, and grey modules. The grey module is reserved for unassigned genes. (D) Cluster analysis for eigengenes dendrogram (Top) and correlation degree heatmap (Bottom) of eigengenes in each co-expression modules.

Identification of highly correlated modules connected with AS or AI

PCA for gene expression values in the indicated modules were computed. The PC1 for indicated genes was defined as eigengene. The correlation analyses for modules and clinical traits were performed by eigengenes from various modules and clinical traits data of 28 samples (Figure 2A). To find the highly correlated genes with AS and AI pathophysiologic mechanisms, these modules (correlation coefficient >0.6 and P value <0.05) associated with the degree of AS, degree of AI, and heart failure were subjected to further analyses (Figure 2A).

FIGURE 2

Figure 2. Identification of highly correlated modules connected with AS or AI. (A) The correlation heat map for various modules and clinical traits (age, sex, BMI, degree of AS, degree of AI, LVEF, diabetes, hypertension, coronary heart disease, and heart failure). Each cell contains the corresponding correlation coefficients and P value. (B–G) Scatterplot of Gene Significance (GS) for indicated clinical traits vs. Module Membership (MM) in the indicated modules.

To further study the relationship between modules and clinical traits, the associations of gene expression values with degree of AS, degree of AI, and heart failure were analysed, the correlation coefficients were defined as gene significances (GS). The correlations of gene expression values with the modules (orange, steelblue, darkgreen, and grey60 modules) screened out in Figure 2A were calculated, which were marked as module memberships (MM). The scatter plots for MM and GS in the indicated modules were plotted, these plots (module membership vs. gene significance) with correlation coefficients more than 0.5 and P value less 0.05 were regarded as highly correlated for MM and GS (Figures 2B–G). These plots were used to screen highly correlated genes in AS and AI pathogenesis. These genes with GS > 0.8, MM > 0.8, and P value <0.05 were singled out in the plots. There were 30, 108, 133, 29, 77, and 5 genes in indicated modules screened out by the pre-set parameters (Figures 2B–G). Finally, 302 highly correlated genes were identified after removing duplicates.

GO and KEGG enrichment analysis for highly correlated genes

In order to parse the molecular functions involved by these 302 highly correlated genes, gene function enrichment analyses of GO and KEGG were conducted using clusterProfiler program package. The results of GO enrichment analyses including BP, MF, and CC were obtained. These terms of GO and KEGG enrichment analyses were ranked by ascending P value and descending gene counts. We selected top 10 highly correlated terms with AS and AI as significantly enriched terms in GO and KEGG analyses.

Among these terms in BP analysis, our results showed that genes highly correlated with AS and AI were mainly enriched in collagen fibril organization (GO:0030199), extracellular matrix organization (GO:0030198), extracellular structure organization (GO:0043062), external encapsulating structure organization (GO:0045229), cell-substrate adhesion (GO:0031589), response to oxidative stress (GO:0006979), negative regulation of cellular protein localization (GO:1903828), negative regulation of interleukin-2 production (GO:0032703), pentose metabolic process (GO:0019321), and regulation of cell morphogenesis (GO:0022604) (Figure 3A). For MF analysis, the enrichment analysis results mainly included extracellular matrix structural constituent (GO:0005201), extracellular matrix structural constituent conferring tensile strength (GO:0030020), and kinds of molecular binding (GO:0005518, GO:0048407, GO:0019001, GO:0032561, GO:0005525, GO:0032550, GO:0030246, GO:0001883) (Figure 3B). In CC enrichment analysis, the significantly enriched terms mainly contained extracellular matrix (GO:0062023), endoplasmic reticulum lumen (GO:0005788), collagen (GO:0005583, GO:0098643, GO:0005581, and GO:0098644), focal adhesion (GO:0005925), cell-substrate junction (GO:0030055), cell division site (GO:0032153), and actin filament bundle (GO:0032432) (Figure 3C). These results indicated that the highly correlated genes mainly functioned in these courses.

FIGURE 3

Figure 3. Go and KEGG enrichment analysis for highly correlated genes. (A–C) The GO enrichment analyses were conducted using the 302 highly correlated genes, the top 10 terms of from BP (A), MF (B), CC (C) were depicted. (D) The KEGG enrichment analyses were performed, the 10 significantly enriched pathways were identified. The circle size represents gene counts in each enriched term. The different color means significance for each enriched term.

To investigate the enriched pathways of these highly correlated genes, the KEGG enrichment analyses were conducted. As shown in Figure 3D, the significantly enriched pathways included protein digestion and absorption (hsa04974), glutathione metabolism (hsa00480), Th17 cell differentiation (hsa04659), Th1 and Th2 cell differentiation (hsa04658), ferroptosis (hsa04216), phagocytosis (hsa04658), mitophagy (hsa04137), valine, leucine and isoleucine degradation (hsa00280), hippo signaling pathway (hsa04392), and pentose phosphate pathway (hsa00030). These results revealed that these pathways may be important in the development of AS and AI.

Identification of DEGs for highly correlated genes

To understand the expression of the 302 highly correlated genes in the AS and AI groups, we performed DEGs analyses using the DESeq2 package. These genes with |log₂FoldChange| ≥ 1 and adjust P value <0.05 were considered as the DEGs. As depicted in Figures 4A–D, there were 31 downregulated genes, 229 unchanged genes, and 42 upregulated genes in the AI group and 35 downregulated genes, 190 unchanged genes, and 77 upregulated genes in the AS group in comparison to the control group. These differential expression profiles among AS, AI, and control group indicated the different molecular mechanism in AS and AI pathogenesis.

FIGURE 4

Figure 4. Identification of DEGs for highly correlated genes. (A–D) The 302 highly correlated genes with AS and AI were used for DEGs expression analyses. These genes with |log₂FoldChange| ≥ 1 and adjust P value <0.05 were regarded as the DEGs in AS and AI group in comparison to the normal control group. The volcano plots and heatmaps were used for visualization of the DEGs for indicated group. The blue dots of volcano plots represent the downregulated genes (DOWN) in AS (A) and AI (C) group. The gray dots of volcano plots represent the unchanged genes (NOT) in AS (A) and AI (C) group. The red dots of volcano plots represent the upregulated genes (UP) in AS (A) and AI (C) group. The heatmaps of DEGs from (A,C) of AS (B) and AI (D) group were depicted, the color scale bar of heatmap represents the scale (from −2 to 2) for the expression levels of genes presented in the heatmaps with a breakpoint of zero represented by white.

Hub genes identification by PPI analysis

Based on the above results, we next explored the important hub genes with high connectivity among these highly correlated genes. The interaction network for the 302 highly correlated genes with AS and AI pathogenesis was constructed and visualized utilizing the PPI analysis tool STRING and Cytoscape, respectively (Figure 5A). To screen out these genes with highest connectivity, the recognized plug-in Cytohubba of Cytoscape software was used to select the top 30 genes with the highest connectivity among the 302 genes by Betweenness button (Figure 5B). To identify the functionally crucial hub genes in mediating AS and AI common pathological processes, we overlapped the DEGs from AS group, AI group, and the top 30 highly connected genes (Figures 5C,D). Finally, ten hub genes including CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB were obtained according to above-mentioned analytical methods (Figures 5C,D). These findings indicate that AS and AI shared the common molecular regulatory network.

FIGURE 5

Figure 5. Hub genes identification by PPI analysis. (A) The PPI network analysis for the 302 highly correlated genes was conducted using the STRING online tool and visualized by Cytoscape software. Each node represents the gene, the different quantity of connecting line among these nodes represents connectivity. (B) The PPI network shown the top 30 genes with highest connectivity from (A) screened out with Cytohubba plug-in by Betweenness button, the node size represents different connectivity. (C) Venn diagram exhibited the overlapping of the DEGs from group and the top 30 genes with highest connectivity from (B) to screen out the hub genes. (D) These hub genes were exhibited in the gene list.

Validation of hub genes expression levels

For further verifying the results acquired from above analyses, we first investigated the expression levels of the ten hub genes in AS and AI group by analyzing expression values in high throughput RNA-sequencing data. The expression levels of ten hub genes were re-analysed. The mRNA expression levels of CD74, COL1A1, CCND1, COL5A1, SERPINH1, FOS, and JUNB were obviously upregulated, whereas those of TXNRD1, BCL6, and ITGA10 were evidently downregulated in AS and AI aortic valves relative to the normal controls (Figures 6A,B). Meanwhile, we selected 5 normal aortic valves, 15 AS aortic valves, and 15 AI aortic valves. These diseased aortic valves were divided into mild, moderate, and severe groups according to the degree of aortic valves lesion, respectively. The mRNA expression levels of CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB were detected. As shown in Figures 6C,D, although there was not evident degree of lesion-dependent trend, the mRNA expression levels of the ten hub genes were consistent with the results in Figures 6A,B from high-throughput sequencing. These highly consistent data suggest that the vital function of the ten hub genes in regulating and indicating the disease course for AS and AI.

FIGURE 6

Figure 6. Validation of hub genes expression levels. (A,B) The mRNA expression levels of the ten hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) of AS and AI group were analysed using the expression values from the high throughput RNA-sequencing data. These genes were marked in the indicated volcano plot. (C,D) The mRNA expression levels of the ten hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) of AS and AI group were detected by real-time PCR assay using the human aortic valve samples of AS and AI. The mRNA levels were normalized to those of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (n = 5). *P < 0.05 vs. the normal control group, n.s., no significance.

Discussion

Herein, using the published high throughput RNA-sequencing data (GSE153555) for AS and AI aortic valves, we conducted WGCNA to identify the key co-expression modules and hub genes involved in AS and AI pathogenesis. Total of eight modules including darkgrey, lightcyan, darkorange, orange, steelblue, darkgreen, grey60, and grey modules were constructed. The associations of modules and clinical traits were computed, and these modules associated with the development and outcome of AS and AI were further analysed to screened out the highly correlated genes with AS and AI pathogenesis. The 302 highly correlated genes were obtained from the indicated modules, and GO and KEGG functional enrichment analyses were performed to explore the potential biological processes and signaling pathways involved in these genes. Furthermore, the expression profiles of the 302 highly correlated genes were used for DEGs analyses. There were 31 downregulated genes, 229 unchanged genes, and 42 upregulated genes in AI group and 35 downregulated genes, 190 unchanged genes, and 77 upregulated genes in AS group, respectively. The PPI network analyses were performed and visualized by STRING online tool and Cytoscape software. The top 30 genes with highest connectivity were singled out and intersected with the DEGs in AS and AI group to find out the common functional hub molecules in AS and AI pathogenesis. Ultimately, ten hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) were obtained and the validation of expression levels was performed to elucidate the molecular mechanism of AS and AI pathological processes.

By exploiting WGCNA module construction, we screened seven significant gene modules associated with AS and AI. We selected highly correlated modules with AS and AI for further analysis, and the 302 highly correlated genes were picked out. Our study found that these highly correlated genes were mainly implicated in the regulation of collagen fibril and extracellular matrix (ECM) by GO analyses, which was consistent with the previous studies (22, 23). These results indicate that the disorders of collagen fibril and ECM may largely impact the normal function of aortic valves. According to the results of KEGG analyses, we enriched top 10 pathways for AS and AI affected mechanism, which included nutrient metabolism, T cell differentiation, ferroptosis, phagocytosis, mitophagy, and Hippo signaling pathway. Glutathione metabolism is the pivotal pathophysiological course for anti-oxidative stress and anti-aging (24). We identified that glutathione metabolism pathway was the significantly enriched term, this hinted that oxidative stress response could be the vital molecular mechanism in regulating the development of AS and AI. Our analysis results were similar with those of David R. A. Reyes et al. (25) and Michael Mahmoudi et al. (26), which manifested the highly reliability of our study. Besides, our study uncovered Th17 cell differentiation and Th1 and Th2 cell differentiation pathway were associated with AS and AI pathogenesis. As the distinctly important contributors for orchestrating adaptive immune responses, Th1, Th2, and Th17 cells are responsible to various intracellular or extracellular pathogens as well as organ-specific autoimmunity, which were activated by a series of cytokines (27, 28). Further, Immune Cell Abundance Identifier (ImmuCellAI) (29) is introduced for precisely estimating the abundance of immune cell types from the high throughput RNA-sequencing data (GSE153555). As shown in Supplementary Figure S2, the box plots demonstrated that AS and AI patients had a higher level of cytotoxic T cells, gamma delta T cells (γδ T), iTreg, Th2 and Tr1 and a lower level of macrophages, neutrophils, and Th17. Our findings indicated that tackling immune responses may become a possibility for harnessing the pathogenesis of AS and AI. Another, we identified the other significantly enriched pathways as the AS and AI underlying mechanisms, such as ferroptosis, phagocytosis, mitophagy, and Hippo signaling pathway. However, the detailed and direct functions in AS and AI for these identified pathways still needed for investigation deeply in the future.

To find out the important regulatory and indicative molecules for AS and AI, we focused on the 302 genes for further analysis. By DEGs analysis and screening these genes with high connectivity, we found out ten hub genes including CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB. In terms of the expression levels for the ten hub genes, our results manifested that CD74, COL1A1, CCND1, COL5A1, SERPINH1, FOS, and JUNB were significantly up-regulated. Whereas those of TXNRD1, BCL6, and ITGA10 were obviously down-regulated. Although, transcriptional profiles of AS and AI pathological processes have been well parsed by Christina L. Greene et al. (14), the identification of highly correlated genes with the progression of AS and AI by systematic computerized algorithm remains unimplemented. Greene et al. used DEGs analysis only to screen for genes of interest, our analysis strategy may be more comprehensive and diverse. Our study may provide the relatively reliable molecular markers for mechanism researches, diagnosis, and treatment of AS and AI.

CD74 (MHC class II invariant chain, Ii), is a kind of type II transmembrane glycoprotein (30). CD74 functions in multiple biological processes and disease types, including lung adenocarcinoma (31), kidney disease (32), spondyloarthritis (33), colitis (34), etc. Collagen type I alpha 1 chain (COL1A1) and collagen type V alpha 1 chain (COL5A1) are the members of collagen family (35), mainly involved in various courses of tumor development, such as hepatocellular carcinogenesis and metastasis (36), immune infiltration in mesothelioma (37), metastasis of lung adenocarcinoma (38), tumor progression in ovarian cancer (39). Thioredoxin reductase 1 (TXNRD1) is a member of the thioredoxin system, regulating hepatocellular carcinoma (40), epilepsy (41), osteosarcoma (42). Cyclin D1 (CCND1) functions as a regulator of CDK kinases and regulates the cell-cycle during G1/S transition (43). Serpin family H member 1 (SERPINH1) is a member of the serpin superfamily of serine proteinase inhibitors and binds specifically to collagen, has been identified acting in gastric cancer metastasis (44) and proliferation and migration of retinal endothelial cells (45). B cell lymphoma 6 (BCL6) is a recognized sequence-specific transcriptional repressor and critical for regulating germinal centers homeostasis (46). Integrin subunit alpha 10 (ITGA10) is a receptor for collagen, has been supposed to be the prognostic biomarker for skin cutaneous melanoma and ovarian cancer (47, 48). Fos proto-oncogene, AP-1 transcription factor subunit (FOS) has been uncovered adjusting cell proliferation, differentiation, and transformation (49). JunB proto-oncogene, AP-1 transcription factor subunit (JUNB) is a member for AP-1 complex, has been identified as the cell proliferation inhibitor and senescence inducer (50), which involves in the regulation of oral squamous cell carcinoma (51) and osteoarthritis (52), etc. However, the molecular functions for these ten hub genes identified by our study in aortic valve were still unclear. Our study provided the feasible molecular bases for the mechanism research or clinical diagnosis and treatment targeting AS and AI, whereas the specific role for these molecules in AS and AI should be further confirmed using in vitro or in vivo experimental models of AS and AI. Collectively, we identified the key signaling pathways and hub genes (CD74, COL1A1, TXNRD1, CCND1, COL5A1, SERPINH1, BCL6, ITGA10, FOS, and JUNB) in AS and AI pathological processes, which may become the potential indicative biomarkers or important regulatory targets in AS and AI pathogenesis. Our findings probably provided the vital theoretical foundation for AS and AI study in the future.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving humans were approved by Human Research Ethics Committees of Tongji Hospital of Huazhong University of Science and Technology. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

LC, YY, JF, and XW conceived the study and designed the experiments; YY, ZP, XF, BX, YC, and QW analyzed the data; LC drafted the manuscript; YY, BX, PZ, and XF prepared the human aortic valve samples. LC and YY performed the experiments; LC, YY, PZ, XF, BX, YC, QW, JF, and XW reviewed and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (NO. 82070488, NO. 81873458, NO. 81800256).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcvm.2023.857578/full#supplementary-material

References

1. Messerli FH, Rimoldi SF, Bangalore S. The transition from hypertension to heart failure: contemporary update. JACC Heart Fail. (2017) 5(8):543–51. doi: 10.1016/j.jchf.2017.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Flaherty JD, Bax JJ, De Luca L, Rossi JS, Davidson CJ, Filippatos G, et al. Acute heart failure syndromes in patients with coronary artery disease. J Am Coll Cardiol. (2009) 53(3):254–63. doi: 10.1016/j.jacc.2008.08.072

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bahit MC, Kochar A, Granger CB. Post-Myocardial infarction heart failure. JACC Heart Fail. (2018) 6(3):179–86. doi: 10.1016/j.jchf.2017.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gill H, Chehab O, Allen C, Patterson T, Redwood S, Rajani R, et al. The advantages, pitfalls and limitations of guideline-directed medical therapy in patients with valvular heart disease. Eur J Heart Fail. (2021) 23(8):1325–33. doi: 10.1002/ejhf.2097

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Seferović PM, Polovina M, Bauersachs J, Arad M, Gal TB, Lund LH, et al. Heart failure in cardiomyopathies: a position paper from the heart failure association of the European society of cardiology. Eur J Heart Fail. (2019) 21(5):553–76. doi: 10.1002/ejhf.1461

CrossRef Full Text | Google Scholar

6. Fleury MA, Clavel MA. Sex and race differences in the pathophysiology, diagnosis, treatment, and outcomes of valvular heart diseases. Can J Cardiol. (2021) 37(7):980–91. doi: 10.1016/j.cjca.2021.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Benjamin EJ, Muntner P, Alonso A, Bittencourt MS, Callaway CW, Carson AP, et al. Heart disease and stroke statistics—2019 update: a report from the American heart association. Circulation. (2019) 139:10. doi: 10.1161/CIR.0000000000000659

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Nkomo VT, Gardin JM, Skelton TN, Gottdiener JS, Scott CG, Enriquez-Sarano M. Burden of valvular heart diseases: a population-based study. Lancet. (2006) 368(9540):1005–11. doi: 10.1016/S0140-6736(06)69208-8

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Iung B, Vahanian A. Epidemiology of valvular heart disease in the adult. Nat Rev Cardiol. (2011) 8(3):162–72. doi: 10.1038/nrcardio.2010.202

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Akahori H, Tsujino T, Masuyama T, Ishihara M. Mechanisms of aortic stenosis. J Cardiol. (2018) 71(3):215–20. doi: 10.1016/j.jjcc.2017.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bing R, Cavalcante JL, Everett RJ, Clavel M, Newby DE, Dweck MR. Imaging and impact of myocardial fibrosis in aortic stenosis. JACC Cardiovasc Imaging. (2019) 12(2):283–96. doi: 10.1016/j.jcmg.2018.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Akinseye OA, Pathak A, Ibebuogu UN. Aortic valve regurgitation: a comprehensive review. Curr Prob Cardiol. (2018) 43(8):315–34. doi: 10.1016/j.cpcardiol.2017.10.004

CrossRef Full Text | Google Scholar

13. Ravalli F, Kossar AP, Takayama H, Grau JB, Ferrari G. Aortic valve regurgitation: pathophysiology and implications for surgical intervention in the era of TAVR. Struct Heart. (2020) 4(2):87–98. doi: 10.1080/24748706.2020.1719446

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Greene CL, Jaatinen KJ, Wang H, Koyano TK, Bilbao MS, Woo YJ. Transcriptional profiling of normal, stenotic, and regurgitant human aortic valves. Genes (Basel). (2020) 11(7):789. doi: 10.3390/genes11070789

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

16. Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. (2021) 49(D1):D605–12. doi: 10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. Cytohubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. (2014) 8(Suppl 4):S11. doi: 10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zoghbi WA, Adams D, Bonow RO, Enriquez-Sarano M, Foster E, Grayburn PA, et al. Recommendations for noninvasive evaluation of native valvular regurgitation: a report from the American society of echocardiography developed in collaboration with the society for cardiovascular magnetic resonance. J Am Soc Echocardiog. (2017) 30:303–71. doi: 10.1016/j.echo.2017.01.007

CrossRef Full Text | Google Scholar

20. Zakkar M, Bryan AJ, Angelini GD. Aortic stenosis: diagnosis and management. BMJ-Brit Med J. (2016) 355:i5425. doi: 10.1136/bmj.i5425

CrossRef Full Text | Google Scholar

21. Fang J, Li T, Zhu X, Deng K, Ji Y, Fang C, et al. Control of pathological cardiac hypertrophy by transcriptional corepressor IRF2BP2 (interferon regulatory factor-2 binding protein 2) novelty and significance. Hypertension. (2017) 70(3):515–23. doi: 10.1161/HYPERTENSIONAHA.116.08728

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Di Vito A, Donato A, Presta I, Mancuso T, Brunetti FS, Mastroroberto P, et al. Extracellular matrix in calcific aortic valve disease: architecture, dynamic and perspectives. Int J Mol Sci. (2021) 22(2):913. doi: 10.3390/ijms22020913

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Urban P, Rabajdova M, Spakova I, Sabol F, Mickova H, Lakatosova K, et al. Molecular recognition of aortic valve stenosis and regurgitation. Eur Rev Med Pharmacol Sci. (2019) 23(24):10996–1003. doi: 10.26355/eurrev_201912_19805

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Homma T, Fujii J. Application of glutathione as anti-oxidative and anti-aging drugs. Curr Drug Metab. (2015) 16(7):560–71. doi: 10.2174/1389200216666151015114515

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Reyes DRA, Gomes MJ, Rosa CM, Pagan LU, Zanati SG, Damatto RL, et al. Exercise during transition from compensated left ventricular hypertrophy to heart failure in aortic stenosis rats. J Cell Mol Med. (2019) 23(2):1235–45. doi: 10.1111/jcmm.14025

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Mahmoudi M, Gormaz JG, Erazo M, Howard M, Baeza C, Feelisch M, et al. Early oxidative stress response in patients with severe aortic stenosis undergoing transcatheter and surgical aortic valve replacement: a transatlantic study. Oxid Med Cell Longev. (2019) 2019:1–8. doi: 10.1155/2019/6217837

CrossRef Full Text | Google Scholar

27. Luckheeram RV, Zhou R, Verma AD, Xia B. CD4+T cells: differentiation and functions. Clin Dev Immunol. (2012) 2012:1–12. doi: 10.1155/2012/925135

CrossRef Full Text | Google Scholar

28. Park H, Park H, Lee J, Bothwell A, Choi J. Sex-based selectivity of PPARγ regulation in th1, th2, and th17 differentiation. Int J Mol Sci. (2016) 17(8):1347. doi: 10.3390/ijms17081347

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, et al. ImmuCellAI: a unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy. Adv Sci (Weinh). (2020) 7(7):1902880. doi: 10.1002/advs.201902880

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Su H, Na N, Zhang X, Zhao Y. The biological function and significance of CD74 in immune diseases. Inflamm Res. (2017) 66(3):209–16. doi: 10.1007/s00011-016-0995-1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Fernandez-Cuesta L, Plenker D, Osada H, Sun R, Menon R, Leenders F, et al. CD74–NRG1 fusions in lung adenocarcinoma. Cancer Discov. (2014) 4(4):415–22. doi: 10.1158/2159-8290.CD-13-0633

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Valiño-Rivas L, Baeza-Bermejillo C, Gonzalez-Lafuente L, Sanz AB, Ortiz A, Sanchez-Niño MD. CD74 in kidney disease. Front Immunol. (2015) 6:483. doi: 10.3389/fimmu.2015.00483

CrossRef Full Text | Google Scholar

33. Sogkas G, Klose K, Baerlecken N, Schweikhard E, Matthias T, Kniesch K, et al. CD74 is a T cell antigen in spondyloarthritis. Clin Exp Rheumatol. (2020) 38(2):195–202. doi: 10.55563/clinexprheumatol/lezc69

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Hoffman JM, Rankin CR, Pothoulakis C. Here to heal: mucosal CD74 signaling in colitis. Cell Mol Gastroenterol Hepatol. (2020) 10(1):197–8. doi: 10.1016/j.jcmgh.2020.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Ricard-Blum S. The collagen family. Cold Spring Harb. Perspect Biol. (2011) 3(1):a4978. doi: 10.1101/cshperspect.a004978

CrossRef Full Text | Google Scholar

36. Ma H, Chang H, Bamodu OA, Yadav VK, Huang T, Wu ATH, et al. Collagen 1A1 (COL1A1) is a reliable biomarker and putative therapeutic target for hepatocellular carcinogenesis and metastasis. Cancers (Basel). (2019) 11(6):786. doi: 10.3390/cancers11060786

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhang C, Liu S, Wang X, Liu H, Zhou X, Liu H. COL1A1 is a potential prognostic biomarker and correlated with immune infiltration in mesothelioma. Biomed Res Int. (2021) 2021:1–13. doi: 10.1155/2021/5320941

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Liu W, Wei H, Gao Z, Chen G, Liu Y, Gao X, et al. COL5A1 may contribute the metastasis of lung adenocarcinoma. Gene. (2018) 665:57–66. doi: 10.1016/j.gene.2018.04.066

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhang J, Zhang J, Wang F, Xu X, Li X, Guan W, et al. Overexpressed COL5A1 is correlated with tumor progression, paclitaxel resistance, and tumor-infiltrating immune cells in ovarian cancer. J Cell Physiol. (2021) 236(10):6907–19. doi: 10.1002/jcp.30350

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Fu B, Meng W, Zeng X, Zhao H, Liu W, Zhang T. TXNRD1 is an unfavorable prognostic factor for patients with hepatocellular carcinoma. Biomed Res Int. (2017) 2017:1–8. doi: 10.1155/2017/4698167

CrossRef Full Text | Google Scholar

41. Kudin AP, Baron G, Zsurka G, Hampel KG, Elger CE, Grote A, et al. Homozygous mutation in TXNRD1 is associated with genetic generalized epilepsy. Free Radical Bio Med. (2017) 106:270–7. doi: 10.1016/j.freeradbiomed.2017.02.040

CrossRef Full Text | Google Scholar

42. Huang S, Zhu X, Ke Y, Xiao D, Liang C, Chen J, et al. LncRNA FTX inhibition restrains osteosarcoma proliferation and migration via modulating miR-320a/TXNRD1. Cancer Biol Ther. (2020) 21(4):379–87. doi: 10.1080/15384047.2019.1702405

PubMed Abstract | CrossRef Full Text | Google Scholar

43. González-Ruiz L, González-Moles MÁ, González-Ruiz I, Ruiz-Ávila I, Ramos-García P. Prognostic and clinicopathological significance of CCND1/cyclin d1 upregulation in melanomas: a systematic review and comprehensive meta-analysis. Cancers (Basel). (2021) 13(6):1314. doi: 10.3390/cancers13061314

CrossRef Full Text | Google Scholar

44. Tian S, Peng P, Li J, Deng H, Zhan N, Zeng Z, et al. SERPINH1 regulates EMT and gastric cancer metastasis via the Wnt/beta-catenin signaling pathway. Aging (Albany NY). (2020) 12(4):3574–93. doi: 10.18632/aging.102831

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hu L, Liu Y, Wei C, Jin H, Mei L, Wu C. SERPINH1, targeted by miR-29b, modulated proliferation and migration of human retinal endothelial cells under high glucose conditions. Diabetes Metab Syndr Obes. (2021) 14:3471–83. doi: 10.2147/DMSO.S307771

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Leeman-Neill RJ, Bhagat G. BCL6 as a therapeutic target for lymphoma. Expert Opin Ther Targets. (2018) 22(2):143–52. doi: 10.1080/14728222.2018.1420782

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Nurzat Y, Su W, Min P, Li K, Xu H, Zhang Y. Identification of therapeutic targets and prognostic biomarkers among integrin subunits in the skin cutaneous melanoma microenvironment. Front Oncol. (2021) 11:751875. doi: 10.3389/fonc.2021.751875

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhu T, Chen R, Wang J, Yue H, Lu X, Li J. The prognostic value of ITGA and ITGB superfamily members in patients with high grade serous ovarian cancer. Cancer Cell Int. (2020) 20:1. doi: 10.1186/s12935-020-01344-2

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Acquaviva C, Bossis G, Ferrara P, Brockly F, Jariel-Encontre I, Piechaczyk M. Multiple degradation pathways for fos family proteins. Ann N Y Acad Sci. (2002) 973:426–34. doi: 10.1111/j.1749-6632.2002.tb04677.x

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Piechaczyk M, Farràs R. Regulation and function of JunB in cell proliferation. Biochem Soc T. (2008) 36(5):864–7. doi: 10.1042/BST0360864

CrossRef Full Text | Google Scholar

51. Tsunoda M, Fukasawa M, Nishihara A, Takada L, Asano M. Junb can enhance the transcription of IL-8 in oral squamous cell carcinoma. J Cell Physiol. (2021) 236(1):309–17. doi: 10.1002/jcp.29843

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lin Z, Miao J, Zhang T, He M, Wang Z, Feng X, et al. JUNB-FBXO21-ERK axis promotes cartilage degeneration in osteoarthritis by inhibiting autophagy. Aging Cell. (2021) 20(2):e13306. doi: 10.1111/acel.13306

PubMed Abstract | CrossRef Full Text | Google Scholar