Abstract

Background

   Considered as one of the major reasons of sudden cardiac death,
   hypertrophic cardiomyopathy (HCM) is a common inherited cardiovascular
   disease. However, effective treatment for HCM is still lacking.
   Identification of hub gene may be a powerful tool for discovering
   potential therapeutic targets and candidate biomarkers.

Methods

   We analysed three gene expression datasets for HCM from the Gene
   Expression Omnibus. Two of them were merged by “sva” package. The
   merged dataset was used for analysis while the other dataset was used
   for validation. Following this, a weighted gene coexpression network
   analysis (WGCNA) was performed, and the key module most related to HCM
   was identified. Based on the intramodular connectivity, we identified
   the potential hub genes. Then, a receiver operating characteristic
   curve analysis was performed to verify the diagnostic values of hub
   genes. Finally, we validated changes of hub genes, for genetic
   transcription and protein expression levels, in datasets of HCM
   patients and myocardium of transverse aortic constriction (TAC) mice.

Results

   In the merged dataset, a total of 455 differentially expressed genes
   (DEGs) were identified from normal and hypertrophic myocardium. In
   WGCNA, the blue module was identified as the key module and the genes
   in this module showed a high positive correlation with HCM. Functional
   enrichment analysis of DEGs and key module revealed that the
   extracellular matrix, fibrosis, and neurohormone pathways played
   important roles in HCM. FRZB, COL14A1, CRISPLD1, LUM, and sFRP4 were
   identified as hub genes in the key module. These genes showed a good
   predictive value for HCM and were significantly up-regulated in HCM
   patients and TAC mice. We also found protein expression of LUM and
   sFRP4 increased in myocardium of TAC mice.

Conclusion

   This study revealed that five hub genes are involved in the occurrence
   and development of HCM, and they are potentially to be used as
   therapeutic targets and biomarkers for HCM.

   Keywords: hypertrophic cardiomyopathy, HCM, weighted gene coexpression
   network analysis, WGCNA, hub gene, biomarkers, bioinformatics analysis

Introduction

   As a common inherited cardiovascular disease, hypertrophic
   cardiomyopathy (HCM) is still an unsolved clinical problem. Previous
   studies have reported that HCM is caused by more than 1,440 mutations
   in 11 or more genes encoding cardiac sarcomeric proteins.[48]^1–3
   Sarcomere mutations in the two most common genes, β-myosin heavy chain
   (MYH7) and myosin-binding protein C3, have been found in approximately
   70% of the genotyped patients.[49]^4 However, genetic defects that
   cause HCM have only be identified in less than 50% of the clinically
   diagnosed probands.[50]^3 Thus, the complex pathogenesis of HCM has not
   been elucidated yet. HCM is a heterogeneous heart disease, and its
   prevalence is at least 1 in 500 people (0.2%) in the general
   population.[51]^3^,[52]^5 Its symptoms include exercise intolerance,
   angina, dyspnea, syncope, heart failure, atrial fibrillation, embolic
   stroke, and sudden death.[53]^6 The current treatment strategies for
   HCM mainly include improving the symptoms and preventing sudden death,
   and no specific treatment exists.[54]^7^,[55]^8 Discovery of hub genes
   in the occurrence and development of HCM may facilitate the discovery
   of specific therapeutic targets and candidate biomarkers.

   Presently, systems biology has been applied to the study of multiple
   diseases.[56]^9–11 Among them, the weighted gene coexpression network
   analysis (WGCNA) identifies candidate biomarkers and therapeutic
   targets by finding gene clusters having a high correlation with the
   phenotype.[57]^12–15 In this study, we used multiple independent gene
   expression datasets from the Gene Expression Omnibus (GEO) database to
   discover potential therapeutic targets and candidate biomarkers for
   HCM. First, we merged two independent datasets and performed
   differential expression analysis. Subsequently, we utilized WGCNA to
   search the gene coexpression modules highly related to HCM, and we
   performed gene annotation to identify their functions. Finally, we
   identified five hub genes, namely FRZB, COL14A1, CRISPLD1, LUM, and
   sFRP4, in the key module, and we validated the expression changes and
   diagnostic values of hub genes by using another dataset and animal
   model.

Materials and Methods

Data Collection

   The RNA expression profiles of [58]GSE133054, [59]GSE141910, and
   [60]GSE36961 were downloaded from the GEO database
   ([61]http://www.ncbi.nlm.nih.gov/geo/). The expression profiling of
   [62]GSE133054, which includes 8 normal cardiac tissues and 8 HCM
   cardiac tissues, was performed using high throughput sequencing. The
   expression profiling of [63]GSE141910, which contains cardiac tissues
   from health donors and patients with cardiomyopathies, was performed
   using high throughput sequencing. And we downloaded the expression data
   of 28 normal and 28 HCM cardiac tissues. The expression profiling of
   [64]GSE36961, which includes 39 normal cardiac tissues and 106 HCM
   cardiac tissues, was performed using array.

Preprocessing of Data and Screening of Differentially Expressed Genes (DEGs)

   In this study, we used R version 3.6.1 for analysis. The [65]GSE133054
   dataset is the expression matrix of raw gene counts, and the voom
   function in “limma” package was used to normalize the dataset.[66]^16
   The downloaded [67]GSE141910 dataset was also normalized using the voom
   function in limma. The [68]GSE36961 dataset was normalized by Fastlo
   normalization and log2 transformation.[69]^17

   We merged the [70]GSE133054 and [71]GSE141910 datasets and used the
   “sva” package to correct the batch effect.[72]^18 Two-dimensional
   principal component analysis (PCA) cluster plots were used to show the
   sample distribution before and after correction. The merged dataset was
   used for further analysis. DEGs were screened using the “limma”
   package,[73]^16 the cutoff criteria were set as adjusted P < 0.05 and
   |log2-fold change (log2 FC)| ≥ 1. The volcano map and heatmap were used
   to demonstrate the differential expression of DEGs.

Enrichment Analyses of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes
and Genomes (KEGG) pathways

   The “clusterProfiler” package was used to perform functional enrichment
   analysis.[74]^19 Adjusted P < 0.05 was considered statistically
   significant, and the enrichment terms were sorted according to gene
   counts.

Coexpression Network Construction

   The top 25% (3670 genes) of the variance in the merged dataset were
   selected for co-expression network analysis. The “WGCNA” package
   version 1.68 in R software was used to construct the co-expression
   network.[75]^20^,[76]^21 Adjacency coefficient (aij) was calculated as
   follows:

   a[ij] = S^β[ij], S[ij] = |cor(x[i],x[j])|

   Where xi and xj are vectors of expression value for genes i and j; cor
   represents Pearson’s correlation coefficient of the two vectors; and
   aij is adjacency coefficient, which is acquired via exponential
   transform of Sij.

   The WGCNA method takes topological properties into consideration in
   order to identify modules from a gene co-expression network. Therefore,
   this method not only considers the relationship between two connected
   nodes, but also takes associated genes into account. Weighting
   coefficient (Wij) is calculated from aij as follows:

   W[ij] = (l[ij] + a[ij])/ (min{k[i], k[j]} + 1 - a[ij]), l[ij] =
   ∑[u]a[iu] a[uj], k[j] = ∑[u]a[iu]

   Where u represents common genes linked gene i and gene j together; aiu,
   the connection coefficient of gene i and gene u; and auj, the
   connection coefficient of gene u and gene j. Wij considers overlapping
   between neighbor genes of genes i and j. Modules were identified via
   hierarchical clustering of weighting coefficient matrix W.

   Sample clustering to detect outliers in the samples. An appropriate
   soft threshold was selected to ensure scale-free topology (R^2>0.9).
   The topological overlap matrix (TOM) was constructed to measure the
   network connectivity of the genes. Genes with similar patterns were
   clustered into the same modules (minimum size =30) using average
   linkage hierarchical clustering. The relationships between phenotypes
   and modules were calculated to identify highly related modules.
   Finally, the highly correlated modules were analysed to explore their
   potential roles. In addition, the gene expression profiles of the
   highly correlated modules were visualised using the R software.

Identification of Hub Genes

   The intramodular connectivity (IC) represents the connection strength
   between genes in the module, and the genes with higher IC are more
   important in the module. Therefore, we identified the genes with high
   IC as hub genes. Cytoscape version 3.7.2 was used for
   visualization.[77]^22

Validation of the Hub Genes in Datasets and Analysis of the Receiver
Operating Characteristic (ROC) Curve

   The expression data of hub genes extracted from the merged dataset and
   verification dataset [78]GSE36961 were utilized to validate the
   differential expression of hub genes in HCM. The “glmnet” package was
   used to construct a logistic regression model for hub genes.[79]^23
   Additionally, the merged dataset was used as a training set to
   construct the model while [80]GSE36961 was used for external
   validation. Subsequently, we used the “pROC” package to perform ROC
   curve analysis.[81]^24

Transverse Aortic Constriction (TAC)

   A pressure-overloaded cardiac hypertrophy model was constructed using
   TAC. Briefly, 8–10-week-old C57BL/6 mice were anesthetised with 0.3%
   intraperitoneal pentobarbital. A minimally invasive incision was
   performed at the junction of sternum and the upper edge of left second
   rib with the help of a stereo microscope. After separating the adjacent
   tissue around the aortic arch, a 6/0 silk suture was threaded under the
   aortic arch, and it was tied over a 27-gauge needle. When the needle
   was removed, the model was successfully made. The sham group underwent
   the same procedure without the banding of aortic arch. The animal
   procedures were performed under a project license (NO.: SRRSH20201215)
   granted by institutional ethics committee of Sir Run Run Shaw Hospital,
   in compliance with institutional guidelines for the care and use of
   animals.

Histology

   After 4 weeks of TAC, the hearts were harvested and fixed in 4%
   paraformaldehyde. Subsequently, the hearts were cut into 4 mm slices
   after dehydration and paraffin embedding. The Masson’s trichrome
   staining was conducted for observing the collagen deposition in the
   heart.

Total RNA Isolation and Quantitative Real-Time Polymerase Chain Reaction
(qRT-PCR)

   Total RNA from the sham and TAC-operated hearts was extracted using
   ultrapure RNA kit (CWBio, Beijing, China). The obtained mRNA was
   reverse transcribed to cDNA using the PrimeScriptTM RT Master Mix
   (Dalian, China). Following this, each 384-well plate was mixed with 10
   ng cDNA, 5 µL qPCR SYBR Mix (Yeason, China), 3 µL ddH2O, and 1 µL
   primer. The primer sequences were provided in [82]Supplementary Table
   5.

Western Blot

   Proteins were lysed from myocardium of sham and TAC mice. The lysate
   was resolved on 8–12% sodium dodecyl sulfate/polyacrylamide gel
   electrophoresis (SDS/PAGE), transferred to nitrocellulose membrane, and
   immunoblotted with the indicated antibodies as following, anti-FRZB
   (Bioss, bs-16185R), anti-Lumican (Ablonal, A11593), anti-Collagen XIV
   alpha 1 (Affinity, AF0573), Anti-SFRP4 (Proteintech, 15328-1-ap),
   Anti-CRISP10 (Bioss, bs-14060R), anti-RCAN (Sigma D6694).

Statistical Analysis

   Most bioinformatic analyses were performed using R with default
   statistical settings and cutoff values specified in the individual
   method sections. Measurement data were expressed as mean ± standard
   deviation (Mean ± SD) and analyzed by unpaired Student’s t-test. P <
   0.05 was considered statistically significant.

Results

Preprocessing of Data and Identification of DEGs

   First, the [83]GSE133054 and [84]GSE141910 expression matrices were
   merged. The merged dataset included 36 normal cardiac samples and 36
   HCM cardiac samples. Following this, the inter-batch difference was
   removed using the “sva” package. The sample distribution before and
   after eliminating the batch effect was shown in the two-dimensional PCA
   cluster diagram ([85]Figure 1A, B). The results indicate an elimination
   of batch difference. [86]GSE36961 was used as the validation dataset.

Figure 1.

   [87]Figure 1
   [88]Open in a new tab

   Merging of datasets and differential expression analysis. (A, B)
   Two-dimensional principal component analysis cluster plot before and
   after merging [89]GSE133054 and [90]GSE141910. (C) Volcano plot of 455
   DEGs. (D) The number of DEGs filtered using limma according to the
   cutoff criteria [adjusted P value < 0.05, |log2 FC| ≥ 1]. (E) Heatmap
   of all DEGs.

   Abbreviations: Normal, normal myocardial tissue; HCM, hypertrophic
   cardiomyopathy myocardial tissue; DEGs, differentially expressed genes;
   FC, fold change.

   In the differential expression analysis, we identified 455 DEGs,
   including 302 upregulated genes and 153 downregulated genes [adjusted P
   value < 0.05, |log2 FC| ≥ 1, [91]Figure 1C, D]. Hierarchical clustering
   analysis showed that the expression levels of DEGs were significantly
   different ([92]Figure 1E).

Functional Enrichment Analysis of DEGs

   The functional categories of GO terms include biological process (BP),
   cellular component (CC), and molecular function (MF). The most enriched
   BP terms were associated with extracellular structural organization,
   extracellular matrix organization, regulation of neurotransmitter
   levels, regulation of blood circulation, and heart contraction. The
   most enriched CC terms were mainly associated with extracellular
   matrix, collagen-containing extracellular matrix, and endoplasmic
   reticulum lumen. As for MF, DEGs were mainly associated with receptor
   regulator activity, cation transmembrane transporter activity, and
   extracellular matrix structural constituent ([93]Figure 2A). Detailed
   information is listed in [94]Supplementary Table 1.

Figure 2.

   [95]Figure 2
   [96]Open in a new tab

   Enrichment analyses of DEGs. (A) GO enrichment analysis of DEGs. (B)
   KEGG pathway enrichment analysis of DEGs.

   Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology;
   KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process;
   CC, cellular component; MF, molecular function.

   In KEGG pathway enrichment analysis, the DEGs were mainly enriched in
   neuroactive ligand-receptor interaction, cytokine-cytokine receptor
   interaction, phagosome, protein digestion and absorption, and
   renin-angiotensin system ([97]Figure 2B). Detailed information is
   listed in [98]Supplementary Table 2.

Construction of Weighted Coexpression Network and Identification of Key
Modules

   The 72 samples were clustered using the “hclust” function in the
   “WGCNA” package to detect outliers (cutHeight = 110, [99]Figure 3A).
   Therefore, the expression profiles of N14, N18, N29, H18, H22, H33,
   H35, and H36 were excluded from the following analyses. Soft-threshold
   β = 7 was used to construct a gene coexpression network. ([100]Figure
   3B and C).

Figure 3.

   [101]Figure 3
   [102]Open in a new tab

   Construction of a coexpression network. (A) Sample clustering to detect
   outliers; N14, N18, N29, H18, H22, H33, H35, and H36 were excluded. (B)
   Analysis of the scale-free fit index for various soft-threshold powers;
   the red line was set at 0.90. (C) Analysis of mean connectivity for
   various soft-threshold powers. (D) Cluster dendrogram of genes in the
   coexpression network.

   Abbreviations: N, normal myocardial tissue; H, hypertrophic
   cardiomyopathy myocardial tissue.

   Next, we used a one-step method to construct the coexpression network.
   Genes with similar expression patterns were clustered into the same
   module. In this study, we identified a total of 11 gene modules
   ([103]Figure 3D). Subsequently, the correlation between Module
   eigengenes (ME) and HCM was calculated, and the blue (r = 0.78; P =
   3×10^−14) and magenta (r = −0.83; P = 4×10^−17) modules were identified
   as significant modules having the highest correlations with HCM
   ([104]Figure 4A). To further evaluate the correlation between HCM and
   the modules, the module membership and gene significance were
   calculated. The blue module (cor = 0.74; P = 1.8×10^−82) showed a high
   positive correlation with HCM ([105]Figure 4B), whereas the magenta
   module (cor = −0.74; P = 1.3×10^−19) showed a high negative correlation
   with HCM ([106]Figure 4C). Meanwhile, the genes in the blue module were
   significantly upregulated ([107]Figure 4D), while the genes in the
   magenta module were significantly downregulated in HCM ([108]Figure
   4E).

Figure 4.

   [109]Figure 4
   [110]Open in a new tab

   Identification of key modules. (A) Module-trait relationships in the
   constructed network. (B and C) MM versus GS plot of the key modules. (D
   and E) Gene expression profiles of the key modules.

   Abbreviations: MM, module membership; GS, gene significance; N, normal
   myocardial tissue; H, hypertrophic cardiomyopathy myocardial tissue.

Functional Enrichment Analysis of Key Modules

   Analysis of BP in the blue module revealed that the genes were enriched
   in extracellular matrix organization, transmembrane receptor protein
   serine/threonine kinase signaling pathway, cellular response to
   transforming growth factor beta (TGF-β) stimulus, TGF-β receptor
   signaling pathway, collagen fibril organization, and keratan sulfate
   catabolic process. As for CC and MF, genes in the blue module were
   mainly enriched in extracellular matrix (CC), endoplasmic reticulum
   lumen (CC), extracellular matrix structural constituent (MF), and
   receptor regulator activity (MF) ([111]Figure 5A). The enriched
   pathways of the blue module were mainly associated with protein
   digestion and absorption and the thyroid hormone signaling pathway
   ([112]Figure 5B). Detailed information is listed in [113]Supplementary
   Tables 3 and [114]4. Since the genes of the magenta module were not
   enriched into any GO terms and KEGG pathway, we did not analyse the
   module further.

Figure 5.

   [115]Figure 5
   [116]Open in a new tab

   Functional enrichment analysis of the blue modules and identification
   of hub genes. (A) GO enrichment analysis of the blue module. (B) KEGG
   pathway enrichment analysis of the blue module. (C) Interaction network
   of the hub genes in the blue module; red nodes represent the hub genes;
   green nodes represent other related genes in the module.

   Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and
   Genomes; BP, biological process; CC, cellular component; MF, molecular
   function.

Identification of Hub Genes

   In the blue module, frizzled related protein (FRZB, IC = 27.6),
   collagen type XIV alpha 1 chain (COL14A1, IC = 24.7), cysteine rich
   secretory protein LCCL domain containing 1 (CRISPLD1, IC = 23.3),
   lumican (LUM, IC = 22.8), and secreted frizzled related protein 4
   (sFRP4, IC = 21.8) were identified as the hub genes ([117]Figure 5C).
   Simultaneously, the hub genes we identified based on high IC also have
   high GS and high MM, and the results are as follows: FRZB (MM
   cor=0.9378, GS 0.7264), COL14A1 (MM cor=0.8891, GS 0.7429), CRISPLD1
   (MM cor=0.9209, GS 0.6862), LUM (MM cor=0.9049, GS 0.6621), and sFRP4
   (MM cor=0.8462, GS 0.6997). These hub genes may play key roles in HCM.

ROC Curve Analysis of Hub Genes

   In order to verify the diagnostic values of hub genes, we performed ROC
   curve analysis. In the merged dataset, FRZB (AUC = 0.949), COL14A1 (AUC
   = 0.953), CRISPLD1 (AUC = 0.917), LUM (AUC = 0.902), sFRP4 (AUC =
   0.923), and the logistic regression model of hub genes (AUC = 0.971)
   showed an excellent predictive power ([118]Figure 6A). Through the
   verification of [119]GSE36961 datasets, we further proved that FRZB
   (AUC = 0.954), COL14A1 (AUC = 0.718), CRISPLD1 (AUC = 0.758), LUM (AUC
   = 0.888), sFRP4 (AUC=0.923), and the logistic regression model of hub
   genes (AUC = 0.889) had an excellent predictive power ([120]Figure 6B).

Figure 6.

   [121]Figure 6
   [122]Open in a new tab

   Analysis of the disease-predicting abilities of hub genes. (A) ROC
   curve analysis of hub genes in the merged dataset. (B) ROC curve
   analysis of hub genes in the verification dataset [123]GSE36961.

   Abbreviations: AUC, area under the curve; ROC, receiver operating
   characteristic.

Validation of Hub Gene Expression

   In the merged dataset, FRZB, COL14A1, CRISPLD1, LUM, and sFRP4 were
   significantly upregulated in HCM samples ([124]Figure 7A).
   Additionally, in the verification dataset [125]GSE36961, the expression
   levels of FRZB, COL14A1, CRISPLD1, LUM, and sFRP4 were significantly
   upregulated in HCM samples ([126]Figure 7B).

Figure 7.

   [127]Figure 7
   [128]Open in a new tab

   Validation of the expression of hub genes. (A) Expression of hub genes
   in the merged dataset. (B) Expression of hub genes in the verification
   dataset [129]GSE36961. (C) Masson’s trichrome staining of mouse hearts.
   (D) Expression of hub genes in TAC mice. Error bars indicate
   mean ± standard deviation. **P < 0.01; ***P < 0.001; ****P < 0.0001.

   Abbreviations: Normal, normal myocardial tissue; HCM, hypertrophic
   cardiomyopathy myocardial tissue; TAC, transverse aortic constriction.

   Myocardial fibrosis was observed in TAC mice ([130]Figure 7C), and
   qRT-PCR was performed on the sham and TAC myocardial samples. The
   results showed that ANP, FRZB, COL14A1, CRISPLD1, LUM, and sFRP4 were
   significantly upregulated in hypertrophic myocardium ([131]Figure 7D).

   Furthermore, we also found that protein levels of LUM and sFRP4
   increased in myocardium of TAC mice, which is consistent with changes
   of genetic transcription ([132]Figure 8). However, myocardial protein
   expression assay showed no difference of FRZB, Col14A1, and CRISPLD1
   between sham and TAC mice.

Figure 8.

   [133]Figure 8
   [134]Open in a new tab

   Western blotting tests representing protein levels of hub genes (FRZB,
   COL14A1, CRISPLD1, LUM, sFRP4) and remodeling gene RCAN in myocardium
   of sham and TAC mice.

   Abbreviations: FRZB, frizzled related protein; COL14A1, collagen type
   XIV alpha 1 chain; CRISPLD1, cysteine rich secretory protein LCCL
   domain containing 1; LUM, lumican; sFRP4, secreted frizzled related
   protein 4; RCAN, calcipressin regulatory protein.

Discussion

   Currently, bioinformatics has become an important tool to explore the
   potential pathogenesis of various diseases.[135]^25^,[136]^26 In this
   study, we combined multiple independent datasets for analysis and
   verification, and we identified the module and hub genes associated
   with HCM by using WGCNA. To the best of our knowledge, this is the
   first study to integrate combined analysis and WGCNA to explore the
   pathogenesis and potential hub genes of HCM.

   In the differential expression analysis, we identified a total of 455
   DEGs, including typical molecular markers of pathological cardiac
   hypertrophy, such as ANP and BNP.[137]^27 Enrichment analysis of DEGs
   revealed that the regulation of extracellular matrix organization,
   neurotransmitter levels, and renin-angiotensin system play an important
   role in HCM. Additionally, we performed WGCNA, and we identified the
   blue module having the highest correlation with HCM. Enrichment
   analysis of the blue module revealed that the extracellular matrix
   organization, collagen fibril organization, TGF-β receptor signaling
   pathway, transmembrane receptor protein serine/threonine kinase
   signaling pathway, and keratan sulfate catabolic process were highly
   related to HCM. Previous studies have shown that the occurrence and
   development of hypertrophic cardiomyopathy is significantly related to
   the changes in extracellular matrix and myocardial fibrosis, and the
   cardiac muscle-specific TGF-β signal plays an important role in
   fibrosis remodeling.[138]^28–30 Additionally, neurohormonal pathways,
   such as the renin-angiotensin-aldosterone system, are involved in the
   pathogenesis of cardiac fibrosis.[139]^31^,[140]^32 Our enrichment
   results are consistent with previous conclusions, which further
   emphasizes the important roles of extracellular matrix, fibrosis, and
   neurohormonal pathways in the pathogenesis of HCM.

   Based on the ICs, we identified five hub genes, namely FRZB, COL14A1,
   CRISPLD1, LUM, and sFRP4, in the blue module. In order to explore their
   potential functions, we analysed the expression levels of these genes.
   In the merged dataset, the expression of hub genes was significantly
   upregulated in HCM patients. Furthermore, in the validation dataset
   [141]GSE36961, the gene expression increased in HCM. qRT-PCR results
   revealed that the hub genes were significantly upregulated in TAC mice.
   These results indicate that hub genes are highly expressed in both
   human and mice hypertrophic myocardium. These genes may play an
   important role in the occurrence and development of HCM. FRZB belongs
   to the secreted glycoproteins family, and it modulates Wnt signaling
   activity.[142]^33 Previously, it has been reported that FRZB can
   inhibit Wnt-mediated cell proliferation in cardiac cushions and
   regulate mesenchymal cell proliferation.[143]^34 Additionally,
   congenital heart defect is associated with FRZB.[144]^35 However, the
   role of FRZB in HCM needs to be further verified. COL14A1 is a
   fibril-associated collagen that regulates fibrillogenesis.[145]^36 It
   has been reported that COL14A1 affects arterial remodeling, and it may
   be involved in the regulation of essential hypertension, acute aortic
   dissection, coronary heart disease, and idiopathic pulmonary
   hypertension.[146]^37–39 Importantly, COL14A1 plays a vital role in
   collagen fibril organization and ventricular morphogenesis.[147]^40
   Therefore, COL14A1 may play a regulatory role in fibrosis and vascular
   remodeling in HCM. CRISPLD1 is likely a protease that can target the
   extracellular matrix, and it may play an important role in cellular
   adhesion.[148]^41^,[149]^42 Remarkably, during the development of heart
   failure in patients with aortic stenosis, the expression of CRISPLD1
   increases progressively, and it may play a functional role by
   regulating Ca2+ cycling.[150]^43 LUM, a glycoprotein belonging to the
   family of small leucine-rich repeat proteoglycans, is highly expressed
   in connective tissues.[151]^44–46 It has been reported that the LUM
   levels in tissues may have a significant effect on the differentiation
   of monocytes into fibrocytes; thus, modulating LUM signaling may be
   useful to treat fibrosis.[152]^47 During proteomic analysis of the
   myocardium in HCM, LUM increases,[153]^48 which is consistent with our
   results. Additionally, LUM controls cardiomyocyte growth by regulating
   the pericellular extracellular matrix, and it may coordinate multiple
   factors of collagen assembly in the murine heart.[154]^49 A moderate
   lack of LUM attenuates cardiac fibrosis and improves diastolic
   dysfunction following pressure overload in mice.[155]^50 Meanwhile,
   mechanical and proinflammatory stimulations can induce fibroblasts to
   produce a large amount of LUM and participate in cardiac remodeling in
   the process of heart failure.[156]^51 These studies suggest that LUM
   may be an important profibrotic molecule in the heart that can serve as
   a potential therapeutic target. sFRP4, a member of the sFRP family,
   functions as a soluble modulator in Wnt signaling.[157]^52–55 Apart
   from binding with Wnt ligands, sFRP4 interacts with extracellular
   matrix receptors such as integrin receptors.[158]^56 Additionally,
   sFRP4 participates in multiple biological processes such as apoptosis,
   angiogenesis, and proliferation.[159]^57–59 Previously, it has been
   reported that sFRP4 inhibits the expression of cardiac genes and
   cardiac differentiation via the activated focal adhesion kinase
   pathway.[160]^60 Particularly, sFRP4 has been proposed in many
   cardiovascular diseases, such as dilated cardiomyopathy,
   atherosclerosis, myocardial ischemia, reperfusion injury, coronary
   artery disease, and heart failure.[161]^61–65 MYH7 is one of the
   sFRP4-targeted genes,[162]^61 which is closely related to HCM, as
   mentioned before. Therefore, sFRP4 may play a role in the development
   of HCM, and it may be a novel candidate for HCM.

   We also explored the diagnostic values of hub genes. The five hub genes
   showed a good predictive value in the merged and validation dataset,
   indicating the characteristic expressions of these genes, suggesting
   their potential usage as biomarkers. Although we combined multiple
   independent datasets for analysis, part of the information in the
   network may have been obscured due to the limitation of sample size.
   Generally speaking, hub genes may play an important role in the
   development of HCM. In our study, we found that the hub genes were
   significantly associated with HCM. However, further experimental
   validation is needed to verify whether these genes are the primary
   drivers of HCM. The protein levels of FRZB, Col14A1, and CRISPLD1
   showed no difference in myocardium between sham and TAC mice. This
   results may partially related to species difference and the limited
   sample size of TAC mice. Additionally, we only focus on the genes in
   the module with the highest correlation, which may cause some hub genes
   in other related modules to be neglected. This is another limitation
   inherent to our study.

Conclusion

   Our results revealed that FRZB, COL14A1, CRISPLD1, LUM, and sFRP4 are
   highly expressed in HCM, and these genes can be potentially used as
   biomarkers and even therapeutic targets for HCM. Additionally, the
   biological process obtained through an enrichment analysis of DEGs and
   key modules may provide a new direction of HCM research.

Funding Statement

   This work was supported by the National Natural Science Foundation of
   China (No. 81974025, NO.81941003); Natural Science Foundation of
   Zhejiang Province (No. LY19H020007); Medical and Health Science Program
   of Zhejiang Province (No. 2020RC016, No. 2019RC183).

Abbreviations

   BP, biological process; CC, cellular component; COL14A1, collagen type
   XIV alpha 1 chain; CRISPLD1, cysteine rich secretory protein LCCL
   domain containing 1; DEGs, differentially expressed genes; FRZB,
   frizzled related protein; GEO, Gene Expression Omnibus; HCM,
   hypertrophic cardiomyopathy; IC, intramodular connectivity; LUM,
   lumican; MEs, module eigengenes; MF, molecular function; MYH7, β-myosin
   heavy chain; PCA, principal component analysis; qRT-PCR, quantitative
   real-time polymerase chain reaction; RCAN, calcipressin regulatory
   protein; sFRP4, secreted frizzled related protein 4; TAC, transverse
   aortic constriction; TGF-β, transforming growth factor beta; WGCNA,
   weighted gene coexpression network analysis.

Data Sharing Statement

   The datasets analyzed during the current study are available in the GEO
   database ([163]http://www.ncbi.nlm.nih.gov/geo/).

Ethics Approval

   Experiments were performed under a project license (NO.: SRRSH20201215)
   granted by institutional ethics committee of Sir Run Run Shaw Hospital,
   in compliance with institutional guidelines for the care and use of
   animals.

Disclosure

   The authors declare that they have no conflict of interest.

References