Transcriptome sequencing data for GSE55492 and GSE148219 were downloaded from the European Nucleotide Archive.¹ The mRNA microarray dataset GSE12644 was downloaded from the GEO database.² The dataset GSE55492 was carried out on GPL11154 (Illumina HiSeq 2000), including 10 noncalcified aortic valve samples and nine calcified tricuspid aortic valves. Series GSE148219 was performed on GPL16791 (Illumina HiSeq 2500) and contained six noncalcified aortic valves and seven calcified tricuspid aortic valves. Dataset GSE12644 was performed on GPL570, which included 10 normal and 10 aortic stenosis valves. Supplementary Table 1 shows the demographic information of patients in these three GEO datasets.

The quality control of these two transcriptome datasets was controlled by fastqc (version 0.11.9). The reference genome, Homo sapiens GRCh38.104, was downloaded from Ensembl.³ The index for the reference genome was built from HISAT2 (version 2.2.1). Thereafter, sequenced reads were aligned to the reference genomes using HISAT2 (version 2.2.1). Samtools (version 1.13) was used to convert the sequence alignment map (SAM) format derived from aligned RNA-Seq samples into binary alignment map (BAM) format. Finally, the read counts matrix for each gene was generated using featureCounts (version 2.0.1). The fragments per kilobase million (FPKM) and transcripts per kilobase million (TPM) matrices were also calculated according to a previous study (8). The read counts data were uploaded into R (software version 4.0.5), after which the differentially expressed gene (DEG) analysis was performed using the voom function in the Limma package (version 3.48.3). The threshold for DEGs was P-value < 0.05 and logFC > mean [abs (logFC) + 2 * SD (logFC)]. According to the gene types in the reference genome, differentially expressed LncRNAs and mRNAs were filtered for further analysis.

The heat map, volcano plot, and principal component analysis (PCA) picture were plotted using pheatmap (version 1.0.12), ggplot2 (version 3.3.5), and FactoMineR (version 2.4) in R, respectively.

Based on the TPM matrix data of GSE55492 and GSE148219, a step-by-step network construction function was applied to establish the mRNA co-expression network and to filter out key modules using the “WGCNA” package (version 1.70-3) in R. First, a proper soft thresholding power was selected to form an adjacency matrix and turn it into a topological overlap. Thereafter, the hierarchical clustering function was executed to plot the clustering trees. The genes in the TPM matrix were clustered into various modules by using a method called average linkage hierarchical clustering. The least gene number in every module was 50. The threshold to merge similar modules was set at 0.3. Positive or negative correlation modules with P-value < 0.05 were candidates for further analyses.

An intersection between positive correlation modules of WGCNA in GSE55492 and GSE148219 and the differentially expressed mRNAs of GSE55492 and GSE14821 were illustrated using Venn diagrams. The same was done using the negative modules. The union set of these two Venn diagrams was considered the DEGs. The intersection results of differentially expressed LncRNAs in GSE55492 and GSE148219 were selected to construct the LncRNA-miRNA-mRNA regulatory network.

The clusterProfiler package (version 4.0.5) (9) was used to perform the Gene Ontology (GO) analysis, including biological processes (BPs), cellular components (CCs), molecular functions (MFs), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. P-values < 0.05 indicated statistical significance.

Metascape⁴ (10) was used to conduct and visualize pathway and process enrichment. DisGeNET (11) was used to exhibit the relationship between input genes and human diseases. Transcriptional regulatory relationships unraveled by sentence-based text-mining (TRRUST) database (12) integrated into Metascape were also applied to show the transcription factors of input genes.

Differentially expressed genes were uploaded into STRING⁵ to build a PPI network with a confidence score > 0.4, which was visualized by Cytoscape (version 3.8.2). A sub-network of the PPI network was then created using the plug-in, MCODE. The genes in subnetworks were picked up and further validated by plotting ROCs in the GSE12644 dataset. Genes with an area under the ROC curve (AUC) greater than 0.7 were considered hub genes associated with CAVD.

The miRNA targets of differentially expressed LncRNAs were predicted using mircode⁶ and LncBase⁷. The miRNAs binding to the hub genes were predicted by four online databases: miRTarBase,⁸ miRDB,⁹ TargetScan,¹⁰ and mirDIP.¹¹ After that, LncRNA-miRNA and miRNA-mRNAs pairs were uploaded and visualized by Cytoscape.

Aortic valve tissues were collected from six severe patients with aortic stenosis and six patients having aortic valve regurgitation without obvious calcification during aortic valve replacement surgery. Patients with coronary heart disease, rheumatic aortic valve disease, congenital bicuspid aortic valve, infective endocarditis, hyperthyroidism, chronic renal dysfunction, smoking, and diabetes mellitus were excluded from this study. Supplementary Table 2 shows the clinical information of these 12 cases. One set of tricuspid aortic valves was then preserved at −80°C for qRT-PCR. Two sets of tricuspid aortic valves without obvious calcification were used to isolate aortic valve interstitial cells. Briefly, the aortic valves were first washed with PBS solution and subjected to 2 mg/ml collagenase II (Sigma, C6885) for 15 min. Epithelium cells in the valves were removed using a sterile cotton swab. After being cut using scissors, the valve tissues were digested again using collagenase II for another 2 h. The cell suspension was filtered, centrifuged, and resuspended with cell medium, including Dulbecco’s modified eagle medium (DMEM) (HyClone, SH30285.FS), 10% of fetal bovine serum (HyClone, SH30396.02), and 1% of penicillin streptomycin combination (HyClone, SV30010). All the experiments were conducted with 3–6 generations of cells. Supplementary Figure 1 shows the purity of valve interstitial cells by immunofluorescence. To induce osteogenic differentiation of aortic valve interstitial cells, the normal cell medium was replaced with an osteogenic medium, including 2 mmol/L of NaH₂PO₄ (Sangon Biotech, A100571-0100), 50 μg/ml of ascorbic acid (Sigma, A4403), and 10^–7 mol/L of insulin (Novo Nordisk) for 7 days. Supplementary Figure 2 shows alizarin red staining for normal valve interstitial cells (VICs) and osteoblast-induced VICs on Day 7.

The total RNA was extracted from aortic valve tissues and aortic valve interstitial cells with Trizol (Takara, Japan), quantified by Epoch Microporous plate spectrophotometer (BioTek, United States), and then reversely transcribed into cDNA by the PrimeScript™ RT Reagent Kit (Takara, Japan). A qRT-PCR was carried out with TB Green™ Premix Ex Taq™ (Takara, Japan) on an LC480 fluorescence quantitative PCR instrument (Roche, Switzerland). The relative expression level of genes was analyzed by the 2^–ΔΔCT method. Table 1 shows all the primers bought from Sangon Biotech (Shanghai, China).

The TPM matrix data of GSE55492 was uploaded to the CIBERSORT.¹² The distribution of the 22 immune cells distribution in each sample and the expression level of immune cells between CAVD and control groups were exhibited by ggplot2 package in R. The Cor function in the Limma package was applied to analyze the correlation coefficients between hub genes and immune cells. The threshold for the correlation coefficient was 0.4.

In GSE55492, we filtered 370 upregulated mRNAs, 257 downregulated mRNAs, 19 upregulated LncRNAs, and 17 downregulated LncRNAs based on P-value < 0.05 and |log₂FC| ≥ 1.20. In GSE148219, 365 upregulated mRNAs, 322 downregulated mRNAs, 7 upregulated LncRNAs, and 19 downregulated LncRNAs were identified according to P-value < 0.05 and |log₂FC| ≥ 1.37. Supplementary Tables 3, 4 show the differentially expressed mRNAs and LncRNAs in GSE55492 and GSE148219, respectively. The volcano plot, heatmap, and PCA plot of differentially expressed mRNAs in GSE55492 and GSE148219 are shown in Figures 2A–F, respectively.

The sample dendrogram and trait heatmap in GSE55492 are presented in Supplementary Figure 3A. The soft-thresholding power of 11 was selected (scale-free R² = 0.92; Supplementary Figures 3B,C). The cluster dendrogram of the different samples is presented in Supplementary Figure 3D, whereas the gene dendrogram and module colors are provided in Supplementary Figure 3E. In GSE148219, the sample dendrogram and trait heatmap are shown in Supplementary Figure 4A. The soft-thresholding power of 12 was selected (scale-free R² = 0.75; Supplementary Figures 4B,C). The cluster dendrogram of the different samples is presented in Supplementary Figure 4D, whereas the gene dendrogram and module colors are provided in Supplementary Figure 4E.

The association between gene modules and CAVD status were also analyzed. There were 25 modules in GSE55492 (Figure 3A). The dark-gray module (r = 0.75, P = 3e−04), light-green module (r = 0.64, P = 0.004), dark-green module (r = −0.49, P = 0.04), and cyan module (r = −0.63, P = 0.005) were considered the four associated modules with P < 0.05. The GSE148219 had 13 modules (Figure 3F). The antique white 4 module (r = 0.67, P = 0.01), orange red 4 module (r = 0.58, P = 0.04), thistle 2 module (r = −0.73, P = 0.005), sky blue 1 module (r = −0.73, P = 0.005), bisque 4 (r = −0.8, P = 0.001), and maroon module (r = −0.86, P = 0.0002) were identified as the six associated modules. The correlation between modules and gene significance were also shown as dark-gray (correlation coefficient = 0.55, P = 1.8e-69; Figure 3B), light-green (correlation coefficient = 0.34, P = 1.4e-28; Figure 3C), dark-green (correlation coefficient = 0.29, P = 1.6e-06; Figure 3D), cyan (correlation coefficient = 0.45, P = 7e-82; Figure 3E), antique white 4 (correlation coefficient = 0.61, P < 1e-200; Figure 3G), orange red 4 (correlation coefficient = 0.59, P < 1.2e-16; Figure 3H), thistle 2 (correlation coefficient = 0.78, P = 3.2e-27; Figure 3I), sky blue 1 (correlation coefficient = 0.81, P = 2.4e-18; Figure 3J), bisque 4 (correlation coefficient = 0.77, P < 1e-200; Figure 3K), and maroon (correlation coefficient = 0.84, P < 5.9e-124; Figure 3L).

After considering the intersection among differentially expressed mRNAs in GSE55492, GSE148219, positive, or negative modules resulting from WGCNA, we screened out 126 differentially expressed mRNAs (Figures 4A,B), which were considered DEGs in CAVD. Four differentially expressed LncRNAs, namely TRHDE-AS1, LINC00092, LINC01094, and LINC00702, overlapped in these two datasets (Figure 4C).

Based on the clusterProfiler package, GO and KEGG were used to elucidate the biological function of the genes that were picked up. As shown in Figure 4D, the BP for differentially expressed mRNAs in GSE55492 were most enriched during extracellular matrix organization, extracellular structure organization, T-cell activation, second-messenger-mediated signaling, and positive regulation of cell adhesion. The CC terms were mostly focused on the regulation of dendritic cell dendrite assembly, collagen-containing extracellular matrix, external side of the plasma membrane, endoplasmic reticulum lumen, and membrane raft. The five enriched MF annotations included receptor ligand activity, signaling receptor activator activity, extracellular matrix structural constituent, glycosaminoglycan binding, and endopeptidase activity. The KEGG pathway analysis revealed that differentially expressed mRNAs in GSE55492 were significantly enriched in the PI3K-AKT signaling pathway, cytokine-cytokine receptor interaction, human papillomavirus infection, neuroactive ligand-receptor interaction, and focal adhesion (Figure 4E).

In GSE148219, the five most enriched BP terms were neutrophil degranulation, neutrophil activation involved in immune response, T-cell activation, leukocyte cell-cell adhesion, and positive regulation of cytokine production. Collagen-containing extracellular matrix, external side of the plasma membrane, secretory granule membrane, secretory granule lumen, and cytoplasmic vesicle lumen were the CC terms with the most enriched genes. The enriched MF were mainly involved in extracellular matrix structural constituent, carbohydrate-binding, G protein-coupled receptor binding, phosphoric ester hydrolase activity, and glycosaminoglycan binding (Figure 4F). The KEGG pathway analysis shows that the PI3K-Akt signaling pathway and cytokine-cytokine receptor interaction were the most enriched pathways, followed by phagosome, focal adhesion, and human papillomavirus infection (Figure 4G).

Functional enrichment analyses were also performed in DEGs for CAVD. As shown in Figure 5A, the top five BP terms were enriched in extracellular matrix organization, extracellular structure organization, lymphocyte differentiation, regulation of leukocyte differentiation, and T-cell differentiation. The top five CC terms were enriched in the collagen-containing extracellular matrix, external side of the plasma membrane, dendritic spine, neuron spine, and coated vesicle. The top five enriched MF terms included extracellular matrix structural constituent, immune receptor activity, extracellular matrix structural constituent conferring tensile strength, IgG binding, and immunoglobulin binding. The KEGG pathways were mostly enriched in PI3K-Akt signaling pathway, ECM-receptor interaction, hematopoietic cell lineage, cell adhesion molecules, and focal adhesion (Figure 5B). More details of enrichment results for DEGs are presented in Table 2.

Moreover, Figure 5C shows the top 20 pathway enrichment using Metascape, such as ECM-receptor interaction, extracellular matrix organization, and focal adhesion: PI3K-Akt-mTOR-signaling pathway and cytokine signaling in the immune system. Enrichment analysis by DisGeNET showed that DEGs were significantly associated with sarcoidosis, nephritis, acute infectious disease, lymphoproliferative disorders, and others (Figure 5D). The TRRUST database provides various transcription factors regulating DEGs, including LIF3, NFKB1, EGR1, SP1, RELA, POU2F1, REST, ETS1, and JUN (Figure 5E).

The PPI network, which consists of 386 edges and 126 nodes, was visualized via the Cytoscape software (Figure 6A). A total of four sub-networks of the PPI network were filtered out via MCODE plug-in, and are shown in Figures 6B–E. A total of 29 genes in the four sub-networks were further validated by applying ROC analyses in GSE12644 (Figures 6F–K), including CD2 (AUC = 0.61), CD3E (AUC = 0.65), ITGAX (AUC = 0.61), HLA-DRB1 (AUC = 0.44), CR1 (AUC = 0.73), TREM1 (AUC = 0.75), TLR8 (AUC = 0.71), MATK (AUC = 0.6), SDC1 (AUC = 0.71), CCL19 (AUC = 0.74), IL2RB (AUC = 0.65), NCAM1 (AUC = 0.73), TLR7 (AUC = 0.57), FCGR3B (AUC = 0.7), IL7R (AUC = 0.62), COL6A6 (AUC = 0.71), SPP1 (AUC = 0.88), COL4A4 (AUC = 0.7), COL28A1 (AUC = 0.48), FCGR1A (AUC = 0.68), HLA-DRA (AUC = 0.62), TREM2 (AUC = 0.68), TRAF3IP3 (AUC = 0.64), CD52 (AUC = 0.67), CD3D (AUC = 0.49), FCGR3A (AUC = 0.7), MAPT (AUC = 0.7), CNTN1 (AUC = 0.73), and GPM6A (AUC = 0.75). Based on the results of ROC, 10 genes (CR1, TREM1, TLR8, SDC1, CCL19, NCAM1, COL6A6, SPP1, CNTN1, and GPM6A) with AUC > 0.7 were considered the hub genes associated with CAVD.

Based on the prediction of online databases, a total of 38 LncRNA-miRNA pairs and 181 miRNA-mRNAs pairs were confirmed. The LncRNA-mediated regulatory network was constructed via Cytoscape software and are exhibited in Figure 7. Interestingly, miR-181b-5p was the downstream target of LINC00702, whereas SPP1 was the target gene of miR-181b-5p. According to the competing endogenous (ceRNA) theory, we speculated that LINC00702-miR-181b-5p-SPP1 axis might participate in the process of CAVD development, which deserves further experimental validation.

The relative expression levels of four LncRNAs and 10 hub genes were detected in aortic valve samples and valve interstitial cells by qRT-PCR. Figure 8A reveals that calcified aortic valve tissues had higher expression levels of LINC00702, LINC01094, LINC00092, SPP1, TREM1, TLR8, and SDC1 than noncalcified aortic valves and that calcified aortic valve had lower expression levels of GPM6A and CNTN1 than the noncalcified aortic valve. Meanwhile, no statistically significant differences between the two groups were observed in the expression levels of TRHDE-AS1, CCL19, CR1, NCAM1, and COL6A6.

Osteogenic-induced aortic valve interstitial cells had higher expression levels of LINC00702, LINC00092, SPP1, TREM1, TLR8, and SDC1 than normal aortic valve interstitial cells. Meanwhile, osteogenic-induced aortic valve interstitial cells had lower expression levels of GPM6A and CNTN1 than normal aortic valve interstitial cells. Moreover, no statistically significant differences were observed in the expression levels of LINC01094, TRHDE-AS1, CCL19, CR1, NCAM1, and COL6A6 (Figure 8B).

In conclusion, statistically significant differences in LINC00702, LINC00092, SPP1, TREM1, TLR8, SDC1, GPM6A, and CNTN1 were observed in both clinical samples and valve interstitial cells.

According to the results of CIBERSORT, the distribution of 22 immune cells between 9 calcified and 10 normal aortic valve samples were distinguished. Supplementary Table 5 shows the results of the distribution of immune cells in different samples. Figure 9A exhibits the proportion of 13 immune cells (cells with zero expression in half of the samples were deleted). The percentage of various immune cells in each sample was illustrated by the bar plot (Figure 9B). Notably, M2 macrophages accounted for the largest proportion in each valve tissue. The boxplot demonstrated that calcified aortic valves had higher levels of Tregs, naïve B cells, and M0 macrophages and lower levels of M2 macrophages compared with the normal valve tissues (Figure 9C). The correlation between validated hub genes and immune cells was also revealed, including SPP1 (Figure 10A), TREM1 (Figure 10B), GPM6A (Figure 10C), CNTN1 (Figure 10D), TLR8 (Figure 10E), and SDC1 (Figure 10F).

The current study has several limitations worth noting. First, this study utilized bioinformatic analysis of transcriptome data, and mRNA expression is not equal to protein expression. Further proteomics of samples obtained from humans or animals is therefore recommended. Second, bioinformatics analysis cannot replace experimental validation. Hence, a well-designed in vitro and in vivo exploration should be conducted to ascertain the hub genes and immune cell regulation in CAVD.