Novel Genes Involved in Hypertrophic Cardiomyopathy: Data of Transcriptome and Methylome Profiling

Hypertrophic cardiomyopathy (HCM) is the most common inherited heart disease; its pathogenesis is still being intensively studied to explain the reasons for the significant genetic and phenotypic heterogeneity of the disease. To search for new genes involved in HCM development, we analyzed gene expression profiles coupled with DNA methylation profiles in the hypertrophied myocardia of HCM patients. The transcriptome analysis identified significant differences in the levels of 193 genes, most of which were underexpressed in HCM. The methylome analysis revealed 1755 nominally significant differentially methylated positions (DMPs), mostly hypomethylated in HCM. Based on gene ontology enrichment analysis, the majority of biological processes, overrepresented by both differentially expressed genes (DEGs) and DMP-containing genes, are involved in the regulation of locomotion and muscle structure development. The intersection of 193 DEGs and 978 DMP-containing genes pinpointed eight common genes, the expressions of which correlated with the methylation levels of the neighboring DMPs. Half of these genes (AUTS2, BRSK2, PRRT1, and SLC17A7), regulated by the mechanism of DNA methylation, were underexpressed in HCM and were involved in neurogenesis and synapse functioning. Our data, suggesting the involvement of innervation-associated genes in HCM, provide additional insights into disease pathogenesis and expand the field of further research.


Introduction
Hypertrophic cardiomyopathy (HCM) is the most common inherited heart disease, with an estimated prevalence of about 1:200-1:500, and is characterized by an increase in the wall thickness of the left ventricle (LV) in the absence of obvious loading or metabolic causes for the observed magnitude of myocardial hypertrophy [1,2]. The clinical spectrum of HCM is very broad, from asymptomatic status during long life to progressive heart failure and sudden cardiac death at a young age [1].
HCM is believed to be a classic Mendelian monogenic disease, mostly caused by mutations in sarcomere genes. However, it still remains unrevealed as to what determines the high clinical heterogeneity of the disease, even among carriers of the same mutation and members of the same family, including monozygotic twins [3,4]. Moreover, the rate of identification of pathogenic/likely pathogenic (P/LP) variants in sarcomere or sarcomereassociated genes in HCM patients usually does not exceeds 40-60% [5], and it is still unclear as to what drives the development of HCM in other cases. Furthermore, many loci associated with HCM are also involved in the formation of dilated cardiomyopathy, which is characterized by opposite changes in the LV tissue [6]. These facts gave origin to a new opinion on HCM as an oligogenic disease with a complex inheritance, the phenotype of which is determined via intergenic interactions and the involvement of gene expression regulatory mechanisms at various levels, combined with environmental factors [7]. Recent genome-wide association studies (GWASs) support this hypothesis via the identification of non-sarcomeric HCM-associated loci, which are involved in the modulation of the disease phenotype [8,9]. Additional insights into HCM etiology and its pathogenesis may be provided through an investigation of transcriptome and epigenome changes in the hypertrophied myocardia of HCM patients.
To date, several transcriptomic studies of hypertrophied myocardium have been performed, most of them being in animal models; however, the correlation of the results in the context of HCM pathogenesis is complicated by described differences in gene expression between mouse and human HCM [10]. In humans, genome-wide expression profiles of hypertrophied myocardium were also evaluated [11][12][13][14][15][16]. An investigation into epigenetic modifications that cause long-lasting alterations of gene expression in myocardium can expand the understanding of HCM molecular mechanisms provided by transcriptomic studies, and thus help in elucidating additional genes that contribute to disease development. Such an attempt was performed only in [14], where histone acetylation was analyzed in parallel with the transcriptome. The DNA methylation of CpG dinucleotides in the C5 position of a cytosine ring is one of the most stable and well-studied epigenetic modifications that regulate the expression of the nearest genes through different mechanisms [17]. There are two studies of DNA methylation in the LV tissues of HCM patients that have been performed to date [18,19].
In the present study, we aimed to identify novel HCM-driving genes by employing a dual-omics approach, which includes RNA sequencing and genome-wide DNAmethylation analysis. As a comparison group, patients with aortic stenosis (AS) were used, in whom secondary LV hypertrophy develops compensatorily to maintain wall stress and cardiac function in the conditions of the increased afterload. By this, we managed to measurably level out the contribution of biological processes participating in the hypertrophy pathogenesis, regardless of its origin (such as cardiac energy balance, response to oxidative stress, apoptosis, fibrosis, etc.), and thus, to point towards new genes that may play an important role in regulating the pathological mechanisms underlying HCM.

Clinical Characteristics of the Studied Individuals
A total of 13 HCM and 14 AS patients undergoing surgical interventions were enrolled in the study. The discovery sample consisted of eight HCM and five AS patients, while the overall sample was used for the RT-qPCR validation of the findings.
The main clinical characteristics of all of the studied patients are shown in Table 1 and  Table S1. The maximal LV wall thickness differed significantly between the studied groups (Mann-Whitney p-value < 0.05) in the overall and discovery samples. In the overall sample, differences in left atrium diameter, end-systolic volume index, Sokolow-Lyon index, and level of NT-proBNP between HCM and AS were significant.

Analysis of Gene Differential Expression Using RNA Sequencing
RNA sequencing (RNA-seq) was performed in LV myocardial samples of eight HCM and five AS patients (Table 1 and Table S1). Principal component analysis after variance stabilizing transformations showed that not only disease status, but also sex (at least in case of HCM) was among the main sources of variability between the samples ( Figure S1). Therefore, sex was included in the design formula in DESeq2 analysis as a covariate. Levels of 20887 transcripts were analyzed.
We found 193 genes, the expression levels of which significantly differed (padj < 0.05) between HCM and AS patients ( Figure 1A and Table S2). A total of 149 (77.2%) of these differentially expressed genes (DEGs) were shown to be downregulated in HCM. A heatmap of DEG expression levels in all participants of the study is shown in Figure 1B. A hierarchi-cal clustering of RNA samples showed that HCM and AS patients aggregated into separate clusters, with the only exception being the patient HCM-8, who had the smallest LV wall thickness among all of the HCM patients. Parameters are marked in bold if they significantly differ between HCM and AS patients in the sample (Mann-Whitney p-value < 0.05). ACE-angiotensin-converting enzyme, AS-aortic stenosis, ARB-angiotensin receptor blocker, BMI-body mass index, eGFR-estimated glomerular filtration rate, HCM-hypertrophic cardiomyopathy, LA-left atrium, LV-left ventricle, MRA-mineralocorticoid receptor antagonist, NT-proBNP-N-terminal pro-B-type natriuretic peptide.
When considering DEGs via the most pronounced changes in expression between HCM and AS patients (|Log2FC| > 1), 38 genes were characterized by lower levels in HCM patients (Log2FC from −1.00 to −2.10), and 14 genes-by higher levels (Log2FC from 1.05 to 1.82) ( Table 2).

Genome-Wide DNA Methylation Analysis
To search for the DNA methylation patterns specific for the myocardium of HCM patients, we performed genome-wide DNA methylation profiling of the same LV samples obtained from eight HCM and five AS patients (see Table 1 and Table S1). Methylation data from 729,969 individual CpG sites passing filtering were included in the analysis. SVD analysis showed that subjects' sex and age, as well as the batch effects of microarrays, do not significantly contribute to the observed variance ( Figure S2); thus, we did not include them in the further analysis as covariates. When considering DEGs via the most pronounced changes in expression between HCM and AS patients (|Log2FC| > 1), 38 genes were characterized by lower levels in HCM patients (Log2FC from −1.00 to −2.10), and 14 genes-by higher levels (Log2FC from 1.05 to 1.82) ( Table 2).  No differentially methylated positions (DMPs) crossed the genome-wide significance threshold p ≤ 5 × 10 −8 . Using a nominal significance threshold (p < 0.01, deltaBeta > 0.05), we detected 1755 DMPs ( Figure 2A and Table S3). A heatmap of DMP methylation levels is shown in Figure 2B. The hierarchical clustering of DNA samples on a heatmap demonstrated their well-defined aggregation into separate clusters corresponding to HCM and AS patients. In total, 1367 DMPs (77.9% of total) were hypomethylated in HCM; 1107 DMPs (63.1% of total) were located in/near 978 known genes.

Correlation of Gene Expression and DNA Methylation Data
In order to search for the possible effects of DNA methylation levels on the expression of neighboring genes, we matched the results obtained by RNA-seq and methylation analyses. Intersecting the lists of 193 DEGs and 978 DMP-containing genes identified 14 genes that were both differentially expressed in HCM and contained DMPs. Spearman's correlation analysis ( Figure 3) identified significant correlations of gene expression and DNA methylation levels for 8 out of 14 genes, namely ANKRD33B, AUTS2, BRSK2, CBFA2T3, GIGYF1, HIVEP3, PRRT1, and SLC17A7. The expression levels of BRSK2 and HIVEP3 correlated with the methylation levels of two DMPs. The expression of HIVEP3 correlated positively with methylation level of cg15644324 (r = 0.59, p = 0.036) and negatively-with methylation level of cg20075528 (r = −0.63, p = 0.025). The expression of BRSK2 positively correlated with methylation levels of both cg14064268 and cg10590925 (r = 0.74 and 0.68, p = 0.0046 and 0.018, respectively). Other positive correlations were observed between the expression of AUTS2, GIGYF1, PRRT1, and SLC17A7; and the methylation levels of cg20713355, cg17047222, cg16040614, and cg16909495, respectively (r from 0.70 to 0.59, p from 0.0093 to 0.038). Negative correlations were identified between ANKRD33B and CBFA2T3, and cg24959663 and cg03704006, respectively (r = −0.71 and −0.60, and p = 0.0086 and 0.032, respectively). Table 2. Differentially expressed genes identified in the LV myocardial samples of patients with HCM when compared to AS (padj < 0.05, |Log2FC| > 1).

No.
Gene

RT-qPCR Validation
In order to confirm high-throughput sequencing data via independent methods, we evaluated the expression levels of several genes using RT-qPCR in the overall samples of 13 HCM patients and 14 AS patients (see Table 1 and Table S1). NOTCH3 and GADD45G were randomly selected among the top five genes that were down-and upregulated in HCM, respectively, and CBFA2T3, GIGYF1, and SLC17A7 were selected among genes with correlated expression and DNA methylation levels. Expression changes were confirmed for all genes, with the only exception being for GIGYF1, the level of which did not differ significantly between the studied groups ( Figure 4).

Correlation of Gene Expression and DNA Methylation Data
In order to search for the possible effects of DNA methylation levels on the expression of neighboring genes, we matched the results obtained by RNA-seq and methylation analyses. Intersecting the lists of 193 DEGs and 978 DMP-containing genes identified 14 genes that were both differentially expressed in HCM and contained DMPs.

Gene Ontology Enrichment Analysis
In order to evaluate functional patterns affected by 193 DEGs and 978 DMP-containing genes, we performed gene ontology (GO) enrichment analysis in the "biological process" category. DEGs were significantly enriched in 83 GO terms (Table S4), and DMP-containing genes in 285 GO terms (Table S5). Sixty-four biological processes were overrepresented by genes from both sets, accounting for 77.1% of all DEG-enriched GO terms and 22.5% of GO terms enriched by DMP-containing genes (Table S6).
The clustering of identified biological processes using REVIGO software ( Figure 5) showed that the majority of GO terms, overrepresented by both DEGs and DMP-containing genes, were associated with the regulation of locomotion and muscle structure development. Biological processes enriched only with DEGs are involved in the regulation of neuron migration, and in cytoskeleton organization and functioning, while processes enriched only by DMP-containing genes are involved in the regulation of binding and cell junction organization, and in responses to endogenous stimuli and neurotransmitter transport.

RT-qPCR Validation
In order to confirm high-throughput sequencing data via independent methods, we evaluated the expression levels of several genes using RT-qPCR in the overall samples of 13 HCM patients and 14 AS patients (see Table 1 and Table S1). NOTCH3 and GADD45G were randomly selected among the top five genes that were down-and upregulated in HCM, respectively, and CBFA2T3, GIGYF1, and SLC17A7 were selected among genes with correlated expression and DNA methylation levels. Expression changes were confirmed for all genes, with the only exception being for GIGYF1, the level of which did not differ significantly between the studied groups ( Figure 4).  The data are presented using a scatter plot (median with interquartile ranges). The differences in RNA levels between HCM patients and AS patients were considered to be significant if the p-value in the Mann-Whitney U test was less than 0.05.

Gene Ontology Enrichment Analysis
In order to evaluate functional patterns affected by 193 DEGs and 978 DMP-containing genes, we performed gene ontology (GO) enrichment analysis in the "biological process" category. DEGs were significantly enriched in 83 GO terms (Table S4), and DMP- Figure 4. Expression levels of NOTCH3, GADD45G, GIGYF1, SLC17A7, and CBFA2T3 in 13 HCM patients and 14 AS patients, as determined via RT-qPCR. RNA levels were calculated relative to GAPDH. The data are presented using a scatter plot (median with interquartile ranges). The differences in RNA levels between HCM patients and AS patients were considered to be significant if the p-value in the Mann-Whitney U test was less than 0.05.

Figure 5. Visualization of GO terms overrepresented by DEGs and DMP-containing gen
Euler-Venn diagram schematically represents sets of identified GO terms enriched with taining genes (orange) or with DEGs (blue) and their overlap (green). TreeMap diagrams IGO clasterization of the identified GO term sets by similarity. Each colored rectangle o represents one GO term; the size of the rectangle indicates the 'uniqueness' of the GO negative of its average similarity to all other input terms. Groups of similar processes fo cluster are marked with the same color.

Discussion
Although HCM was described several decades ago, the understanding of it and pathogenesis still evolves. New data obtained using high-throughput genom ods refine the picture on the molecular mechanisms underlying this common dis significant social consequences. The results of the present study expand a limi оf transcriptome and DNA methylome investigations performed in the heart HCM patients.
In order to advance the understanding of what processes are involved i mation of primary LV hypertrophy, we compared the transcription profiles o with HCM and patients with AS. Although AS patients also developed LV hyp as a result of a long-lasting overloading condition, it was much less pronoun than in HCM patients (the maximal LV wall thickness was 15.0 ± 2.0 and 23.5 respectively; p = 0.0039). We identified significant differences in myocardium tion profiles between HCM and AS patients: 193 DEGs passed the threshold fo corrections, and most of them were underexpressed in HCM. When considerin DEGs with the most pronounced changes in expression between HCM and AS (| > 1), C4B was shown to be the most significantly downregulated gene in HCM; Figure 5. Visualization of GO terms overrepresented by DEGs and DMP-containing genes. Top left Euler-Venn diagram schematically represents sets of identified GO terms enriched with DMPcontaining genes (orange) or with DEGs (blue) and their overlap (green). TreeMap diagrams show REVIGO clasterization of the identified GO term sets by similarity. Each colored rectangle on the map represents one GO term; the size of the rectangle indicates the 'uniqueness' of the GO term-the negative of its average similarity to all other input terms. Groups of similar processes forming the cluster are marked with the same color.

Discussion
Although HCM was described several decades ago, the understanding of its etiology and pathogenesis still evolves. New data obtained using high-throughput genomic methods refine the picture on the molecular mechanisms underlying this common disease with significant social consequences. The results of the present study expand a limited series of transcriptome and DNA methylome investigations performed in the heart tissues of HCM patients.
In order to advance the understanding of what processes are involved in the formation of primary LV hypertrophy, we compared the transcription profiles of patients with HCM and patients with AS. Although AS patients also developed LV hypertrophy as a result of a long-lasting overloading condition, it was much less pronounced in AS than in HCM patients (the maximal LV wall thickness was 15.0 ± 2.0 and 23.5 ± 5.6 mm, respectively; p = 0.0039). We identified significant differences in myocardium transcription profiles between HCM and AS patients: 193 DEGs passed the threshold for multiple corrections, and most of them were underexpressed in HCM. When considering only 52 DEGs with the most pronounced changes in expression between HCM and AS (|Log2FC| > 1), C4B was shown to be the most significantly downregulated gene in HCM; it encodes complement 4b, which is abundantly expressed in the healthy heart muscle [20] and was previously shown to be downregulated in the peripheral blood plasma of HCM patients [21]. NOTCH3, a top-two downregulated gene in HCM, also attracts attention since it is well known to play an essential role in the protection of heart tissue from apoptosis and fibrotic changes [22][23][24]. Too little is currently known about the most upregulated genes, EIF4EBP3 and CTXND1; their role in cardiac tissue still remains to be uncovered. The expression level of the topthree upregulated gene, GADD45G, negatively correlates with LV reverse remodeling [25]. In general, the most significant DEGs participate in the processes of the inflammation and regeneration of the heart; the levels of NOTCH3 and GADD45G in the HCM group, in comparison with AS, may indicate decreased compensatory abilities of the primary hypertrophied myocardium.
An important benefit of our study was in the fact that the DNA methylation profiling was performed in addition to transcriptome profiling in the same myocardium samples from HCM and AS patients. We showed that the profile of DMPs clearly distinguished HCM patients from AS patients, most likely indicating the involvement of DNA methylation in the formation of primary LV hypertrophy.
By comparing the lists of DEGs and genes containing DMPs, we identified a significant correlation of the expression of eight genes with methylation levels of the neighboring DMPs; these genes can be involved in the formation of the HCM clinical phenotype. Half of these genes-AUTS2, BRSK2, PRRT1, and SLC17A7-are implicated in neurodevelopment and synapse formation, suggesting differences in heart innervation between HCM and AS patients. AUTS2 is highly expressed in the brain and mostly plays a role in transcription regulation in neurons during development [26]. BRSK2 encodes serine/threonine kinase SAD-A, which intermediates external signals, polarizing neurons and promoting axonogenesis [27]. The PRRT1 gene product slows the deactivation and desensitization of the synaptic AMPA-regulated glutamate receptors, which are necessary for synapse development and function [28]. SLC17A7 (also known as VGLUT1) encodes a multifunctional transporter mainly studied in the neuron-rich regions of the brain, where it transports L-glutamate from the cytoplasm into synaptic vesicles at presynaptic nerve terminals [29]. The decreased expression of the above-mentioned genes in HCM was positively correlated with the hypomethylation of CpG sites in gene bodies. This phenomenon is not surprising, since gene body methylation is known to be mostly associated with gene expression upregulation [30]. The remaining identified genes play roles in tyrosine kinase receptor signaling (GIGYF1) and transcription regulation (CBFA2T3 and HIVEP3), or they have unknown functions (ANKRD33B).
Gene Ontology enrichment analysis identified a broad variety of biological processes significantly enriched by DEGs (83 pathways) and by DMP-containing genes (285 pathways). The majority of GO terms enriched by genes from both gene sets are associated with the regulation of locomotion and muscle structure development (see Figure 4), which points toward the differences in the molecular mechanisms of hypertrophy formation in HCM and AS, both at the transcriptomic and epigenomic levels.
Interestingly, a REVIGO analysis highlighted quite a large cluster of GO terms such as "regulation of neuron migration" enriched by DEGs and a cluster "neurotransmitter transport" enriched by DMP-containing genes. Although similar innervation-associated clusters were not observed among shared GO terms, pathway GO:0022008 "neurogenesis" was identified during a manual search (see Table S6). Thus, GO terms associated with neuron functioning were enriched by both DEGs-and DMP-containing genes. This fact supports the assumption on the differences in heart innervation between HCM and AS patients, which was previously made based on the functions of genes identified in a correlation analysis. As expected, two of these genes, AUTS2 and BRSK2, participate in the pathway GO:0022008 "neurogenesis", mentioned above. Several other genes associated with cardiomyopathies were earlier shown to be involved in a variety of neuromuscular disorders [31]; this points additionally to the link between innervation-related mechanisms and heart pathology. Sympathetic and parasympathetic innervation of the heart is known to play a crucial role in heart maturation during ontogenesis, through regulation of the proliferation, physiological hypertrophy, metabolism, and several electrophysiologic properties of cardiomyocytes [32]. In adult hearts, sympathetic/parasympathetic imbalance is associated with pathological remodeling and heart failure [33,34]. An innervation imaging may have a prognostic value in patients with cardiomyopathy [35]. Overall, these data suggest the involvement of innervation-associated mechanisms in HCM pathogenesis.
In contrast to previous transcriptomic studies searching for the differences between HCM myocardium and healthy heart tissue, we used samples from AS patients as a control. However, it seems that the pattern of biological processes affected by identified DEGs is quite similar to previous studies. Indeed, the main biological pathways previously shown to be deregulated in HCM were associated with cell motility, muscle cell development, cytoskeleton reorganization and adhesion, angiogenesis, response to wounding, etc. [11,13,14]; similar GO terms were highlighted in the present work (see Table S4). The processes involved in neurogenesis and neuromuscular disorders were also reported [11,14], however, they did not receive much attention. At the same time, in contrast to published studies [14,16], we did not observe the dysregulation of pathways affecting ATP synthesis and mitochondrial functions in HCM. Most likely, these processes are involved not only in the pathogenesis of HCM, but also in AS-associated secondary LV hypertrophy. Thus, by using AS the comparison group, we were unable to identify any and all aspects of HCM pathogenesis, but we revealed the processes differently involved in genetically driven myocardial hypertrophy in HCM, and the more common secondary hypertrophy instead.
Our study included HCM patients with identified P/LP variants in the ALPK3 gene (the patient HCM-2) and the MYBPC3 gene (the patient HCM-5), along with six genotype-negative patients. However, the DEG profiles of HCM patients (except for the patient HCM-8), as well as their DMP profiles (without exception) clustered together (see Figure 1 and Figure 2). The published data also show that gene expression changes are mostly similar among HCM patients with P/LP variants in different causative genes [12]. Altogether, it suggests the existence of the universal molecular mechanisms of HCM pathogenesis, regardless of the differences in genetic background.
By choosing AS patients as a comparison group for our study, we were able to exclude the confounding effects of universal LV hypertrophy pathogenic factors, such as arterial hypertension, which was equally prevalent in HCM and AS patients (84.6% and 85.7% of patients, respectively). However, adding the data on transcriptome and methylome profiling in healthy heart samples would certainly expand the informativity of the research. In this connection, we have to say that the inaccessibility of healthy LV tissues for subsequent analysis is the discernible limitation of our study. Another limitation is the small number of participants, which is determined by the number of surgical interventions performed; it entails a lack of statistical power to discover genome-wide DNA methylation changes. Nevertheless, the prominent similarity of the biological processes affected by DEGs and DMPs in our study and the processes reported in previous works, in our opinion, testify to the validity of our data.
Generally, we observed clear differences in the gene expression and DNA methylation profiles of LV myocardial samples obtained from HCM and AS patients, which indicate the differences in the molecular mechanisms that define primary and secondary myocardial hypertrophy. Such a multilevel analysis pinpointed eight genes involved in HCM, which might be the subjects of further study. We believe that our data, suggesting the possible involvement of innervation-associated genes in HCM, can provide additional insights into disease pathogenesis and expand the field of research.

Patients and Controls
Thirteen patients with severe obstructive HCM (mean age 56.5 ± 12.0 years; women 46.2%) were recruited into the study. HCM was diagnosed based on the 2014 recommendations of the European Society of Cardiology [36]: increased LV wall thickness ≥ 15 mm that was not solely explained by abnormal loading conditions. Patients suspicious of phenocopies of HCM were excluded from the study. Other exclusion criteria were systemic autoimmune, viral, bacterial, and oncological diseases. All HCM patients were previously studied by a targeted massively parallel sequencing method using a panel that included robust validated HCM-associated genes and the genes of its main phenocopies (ACTC1, DES, FHL1, FHOD3, GLA, LAMP2, MYBPC3, MYH7, MYL2, MYL3, PRKAG2, PTPN11, TNNC1, TNNI3, TNNT2, TPM1, TRIM63, and TTR) in certified genetic laboratories. One HCM patient carried a variant in the MYBPC3 gene, and one patient in the ALPK3 gene; in other patients, P/LP variants were not identified. The control group consisted of 14 patients with severe AS (mean age 60.4 ± 9.0 years; women 35.7%).

Sample Processing
Myocardial samples were obtained from HCM and AS patients undergoing surgical interventions. In HCM patients, the sampling area was located under the right aortic coronary leaflet and corresponded to the zone of Morrow septal myectomy. In AS patients, myocardial biopsies sized 5 × 7 × 10 mm were taken from the left-side interventricular septum during open-heart surgery before aortic valve replacement. All interventions were performed in the Federal Scientific and Clinical Center of Specialized Medical Care and Medical Technologies (Moscow, Russia) by the same surgeon. After extraction, tissue samples were immediately cut into fragments up to ≤ 0.5 cm in any single dimension, submerged in 5 volumes of RNAlater solution (Thermo Fisher Scientific, Waltham, MA, USA), and stored at -20 • C for further processing. The disruption of tissue samples was performed using the TissueLyser LT bead mill (Qiagen, Germantown, MD, USA). Total RNA was extracted from disrupted samples using the RNeasy Mini Kit (Qiagen, USA), and genomic DNA using the DNA Mini Kit (Qiagen, USA). RNA and DNA purity and quantity were assessed using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA). RNA integrity was assessed by the QIAxcel Advanced System (Qiagen, USA); samples with RNA integrity numbers of above eight were included in subsequent experiments.

RNA-Seq Analysis
RNA libraries were constructed from 1 µg RNA using MGIEasy RNA Library Prep Set following manufacturer instructions. Global transcriptome profiling was performed via high throughput RNA-seq on the MGISEQ-200 platform.
The quality control of raw sequencing reads was assessed with FASTQC. Low quality nucleotides were trimmed from the raw fastq files using Trimmomatic v. 0.39. Samples were mapped on Genome primary assembly GRCh38 with the Gencode annotation using STAR v. 2.7.6a [37]. Reads per gene were counted using the GeneCounts option of the STAR software. Genes with a small number of reads (<100 per sample) were excluded from the analysis.
DESeq2 R package v. 1.32.0 was used for normalization, principal component analysis, and differential gene expression analysis (design = ∼ Sex + Condition) [38]. The amplitude of expression changes was represented in the log2FC format. The False Discovery Rate (FDR) procedure of Benjamini-Hochberg was used to adjust for multiple testing corrections [39].

Genome-Wide DNA Methylation Analysis
Bisulfite conversion of the genomic DNA was performed using the EZ DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA). DNA methylation levels were analyzed on the Infinium MethylationEPIC BeadChip using the iScan array scanner (Illumina, USA) at the SB RAS Genomics Core Facility (Institute of Chemical Biology and Fundamental Medicine SB RAS, Novosibirsk).
Methylation data processing and normalization were performed using the ChAMP R package v. 3.14 [40]. Probes with detection p-values above 0.01, multi-hit probes, probes close to SNPs, and non-autosomal probes were filtered out from the subsequent analy-sis using the champ.filter() function with default parameters. A beta-mixture quantile normalization method (BMIQ) [41] provided by the champ.norm() function was used to normalize methylation data for probe type bias (quality control plots for non-normalized and normalized methylation data are provided in the Figure S3). The examination of the methylation data for possible covariates, such as sex and age, as well as the batch effects of different microarray slides, was performed using singular value decomposition (SVD) analysis implemented in the champ.SVD() function. The evaluation of methylation level for each probe was performed by calculating its beta-value-the ratio of the methylation signal intensity and the overall intensity (the sum of intensities of the methylated and unmethylated probes). Beta-values ranged for each probe from 0 (fully unmethylated probes) to 1 (fully methylated probes). DMPs were identified using the general linear model. Along with the standard p-value, FDR correction was applied.

RT-qPCR
Total RNA was reverse transcribed using the High Capacity RNA-to-cDNA Kit (Ther-moFisher Scientific, USA), and cDNA was then subjected to qPCR in triplicate using the TaqMan Gene Expression Assays (ThermoFisher Scientific, USA) according to the manufacturer's protocol. GAPDH was used as an endogenous control. The relative expression of each gene in the sample was calculated as 2-deltaCt, where deltaCt for the gene is equal to its average Ct minus the average Ct of GAPDH. The comparison of gene expression levels between groups was performed using the Mann-Whitney U test (p < 0.05 were considered significant).

Gene Set Enrichment Analysis and Data Visualization
GO term enrichment analysis was performed using the Panther classification system [42]. GO terms with FDR corrected p-values of less than 0.05 were considered as being significant. In order to shorten the list of GO terms and to find the most representative ones, we used the REVIGO service, which implements a simple clustering algorithm based on measures of semantic similarity and the following redundancy reduction; a reduction setting of "Large (0.9)" was used for the analysis [43].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms232315280/s1, Figure S1: Two-dimensional clustering of samples based on principal component analysis of 500 most variable genes. Percentage of variance explained by each principal component (PC) is indicated on the corresponding axis; Figure S2: Singular value decomposition analysis (SVD) plots for normalized methylation data from 729969 probes passing filtering; Figure S3: Quality control plots for non-normalized and normalized methylation data from 729969 probes passing filtering; Table S1: Characteristics of HCM and AS patients included in the study; Table S2: Differentially expressed genes identified in myocardia of HCM patients when compared to AS (padj < 0.05); Table S3: Differentially methylated CpG sites identified in myocardia of HCM patients when compared to AS (p < 0.01, delta β > 0.05); Table S4: GO Biological processes significantly enriched (adj. p < 0.05) by differentially expressed genes (DEGs), which were identified in the myocardia of HCM patients; Table S5: GO Biological processes significantly enriched (adj. p < 0.05) by genes containing differentially methylated positions (DMPs), which were identified in the myocardia of HCM patients; Table S6: GO Biological processes significantly enriched (adj. p < 0.05) by both differentially expressed genes (DEGs) and genes containing differentially methylated positions (DMPs).

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of City Clinical Hospital #17 (protocol #8, 25OCT2019).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: RNA-seq data are available in the Gene Expression Omnibus (GEO) database under accession number GSE206978, and genome-wide DNA methylation data-under number GSE207095.

Conflicts of Interest:
The authors declare that they have no conflicts of interest.