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).

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.

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.

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.

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.

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 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).

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.

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.

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).

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).

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.

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.

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.

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.

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.

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.