The landscape of gene regulatory elements active in human hearts
Heart CAGE dataset contained 103 samples across 21 healthy (10 males and 11 females) and 10 failing (7 males and 3 females) human hearts, with 20 of the healthy samples (3 males and 4 females) supported with RNA-seq (Fig. 1A, Supplementary Table 1). This study focused on promoters and enhancers with alternative TREs and the resulting gene expression (Fig. 1B).
Comparing the distance to previously identified TSS for male and female unidirectional TREs, more female unidirectional TREs coincide with the earlier predicted TSS: 800 of the female unidirectional TREs are within 10 bp from the previously identified TSS, and 724 of the male unidirectional TREs are within 10 base pairs from the TSS. The mean distance to previously identified TSS is -143 for females and + 1560 for males.
The principal component analysis (PCA) across all TREs showed the most significant separation between atria and ventricles, explaining 47.76% of the total variance (Fig. 1C). The second largest contributor to gene expression variance, 9.47%, was health status (Fig. 1D). This variation could be associated with the left ventricular ejection fraction (Fig. 1E). The third principal component separated samples by sex, explaining 4.77% of the variance (Fig. 1F). Kmeans clustering of the samples in the space of the first and the second principal components for K = 2 has partitioned the samples into clearly separated groups. Kmeans clustering of the samples in the space of the third and the fourth principal components for K = 2 has partitioned the samples into groups that were more difficult to separate: group 1 contained 46 male samples and no female samples; group 2 included 50 female samples and 7 male samples. The 7 samples came from four individuals: all chambers of a 51-year-old failing male, the left atrium of the 48 years of failing male, the left ventricle of the 61 years of a healthy male, and the right atrium of the 57 years old healthy male hearts.
Pearson’s correlation analysis between all samples demonstrated two high-correlation groups – atrial and ventricular except for 9 out of 15 failing atria samples (Supplementary Fig. 1). All failing female atria (n = 6) clustered with ventricles, while only 4 out of 9 failing male atria shifted to ventricular profile. In Deviatiiarov et al.3, we have shown that failing atrial reprogramming to ventricular phenotype was primarily due to downregulation of atria-specific genes, which also confirmed findings by other researchers4, who termed it “ventricularization” of atrial myocytes. Interestingly, female atria showed higher susceptibility to this reprogramming or ventricularization, potentially indicating higher plasticity of atrial myocytes in failing female hearts.
Differential expression analysis reveals distinct sex-specific regulatory elements and genes in healthy atria and ventricles.
Since the atria and ventricles showed the most significant degree of variance, we analyzed those regions separately. We identified 364 and 558 differentially expressed TREs between male and female healthy atria and ventricles, respectively (Fig. 2A). TREs located in sex chromosomes represented 15.7% (57/364) of all differentially expressed peaks in healthy atria (Supplementary Fig. 2A) and 11.1% (62/558) in ventricles (Supplementary Fig. 2B). Thus, most of the male-female TRE differences are in autosomal chromosomes. Supplementing our findings with RNA-seq data, we identified 96 differentially expressed genes in the atria and 191 in the ventricles (Fig. 2A). Comparing CAGE data with RNA-seq suggests a sex-specific regulatory role without a large difference in gene expression. It is important to note that there is a lag between activation and transcription. Thus, a perfect correlation should not be expected. The top overexpressed TREs across both regions in females were associated with XIST, a non-coding RNA on the X chromosome which regulates the X-inactivation process (Supplementary Table 2). The top male TREs were associated with the Y chromosome and protein synthesis function. Interestingly, mostly TREs located in sex chromosomes overlapped between atria and ventricles, while overall overlap indicated unique differences in each cardiac region (Fig. 2B).
To understand the function of sex-specific regulatory elements, we performed gene ontology (GO) analysis separately, considering the effect of each sex on each region. Among the top sex-specific pathways passing FDR < 0.05, we identified immune system response, muscle cell-related pathways, and responses to stimuli (Fig. 2C). Compared to males, healthy female atria had a significantly more active immune system, while healthy female ventricles had enrichment in the cellular response to stimulus, specifically to cytokines. Both healthy male atria and ventricles showed enrichment in pathways related to cardiac muscle function, such as muscle cell development and differentiation, which was only partially enriched in female ventricles. These differences indicate unique functional differences in the heart chambers between males and females.
To further understand the differential expression of TREs, we identified genes differentially expressed in CAGE and RNA-seq data. One of them, HOPX (HOP homeobox gene, atypical homeodomain protein), promotes cardiomyogenesis by activating BMP to inhibit WNT5. We found that both, the main and alternative promoters, were upregulated (FDR < 0.05) in healthy male atria and ventricles compared to females. RNA-seq data supported this finding (Fig. 2D).
Differential expression analysis reveals distinct sex-specific regulatory elements in failing atria and ventricles
After identifying numerous sex differences in healthy hearts, we investigated whether these differences persist during heart failure. We identified 636 differentially expressed TREs in the atria and 258 in the ventricles (Fig. 3A). Failing male atria had more differentially expressed TREs than females in both regions. Again, differential expression in sex chromosomes represented only a minor portion of all statistically different TREs, 69 in failing atria and 49 in failing ventricles (Supplementary Fig. 2C-D). Similar to healthy samples, there was a small overlap of differentially expressed TREs between atria and ventricles of the same sex (Fig. 3B).
Interestingly, only a small number of TREs had consistent sex differences in healthy and failing samples in both atria and ventricles (Supplementary Fig. 3A-D), suggesting that heart failure pathogenesis manifests differently in males and females. We noticed the opposite enrichment of GO pathways compared to healthy hearts. This time, failing male ventricles showed enrichment in cytokine response compared to female ventricles, while failing female ventricles had enrichment in cardiac muscle processes compared to male ventricles (Fig. 3C). Failing male atria enriched pathways belonged to several homeostatic processes and cell communication and signaling while failing female atria showed enrichment in developmental processes. These differentiation and developmental processes can potentially be explained by the fact that all failing female atria exhibited reprogramming or ventricularization.
One of the most striking examples of sex dimorphism during heart failure was PCSK1N, a proprotein convertase subtilisin/kexin type 1 inhibitor that functions as a prohormone convertase 1 inhibitor which regulates the proteolytic cleavage of neuroendocrine peptide precursors. Its deficiency was shown to lead to obesity6. Additionally, PCSK1N interacts with PCSK1 which genetic variants are associated with the risk of coronary artery disease in type 2 diabetic patients of Han Chinese decent7. In healthy hearts of both sexes, PCSK1N was equally expressed in atria, but absent in ventricles (Fig. 2D). During heart failure, the PCSK1N expression in male atria increased, while in female atria, the expression almost completely vanished. These findings indicate PCSK1N chamber specificity and opposite expression changes in each sex during heart failure.
Sex-specific promoter usage in heart
CAGE allows for identifying differential promoter usage between male and female hearts. In total, sex-specific promoter usage passed 0.05 FDR for 720 genes (Fig. 4A, Supplementary Fig. 4, Supplementary Table 3), from 45 (failing atria) to 282 (healthy atria) sex-specific cases of differential promoter usage for sample subsets corresponding to individual heart chambers and health states. Notably, 40 genes demonstrated promoter switching events, where the dominant promoter contributing the most to the total expression of a gene was different between the two sexes.
The promoter-switching events within the same gene varied between health states and chambers. We highlighted 35 genes located at the autosomes and demonstrated sex-specific promoter switching (Fig. 4A and Supplementary Fig. 4). There were several functionally distinct groups exhibiting sex-specific promoter switch: cardiomyocyte function and remodeling (MYOT8, ABR9, PRRX110,11, and CLCN312,13), metabolic processes (CRAMD1B (cholesterol)14, GRB14 (insulin)15, THAP4 (detoxification of ROS)16, USP217 (protein metabolism), and PRKAG218 (energy sensor)), immune system function (CD3619,20, CABIN1, and NFATC121), and developmental formation of the heart (PACSIN222, RASGRP323,24, ABL125, LTBP126, RGS327, and TBX528).
One important example is TBX5, encoding a critically important transcription factor, which along with strong chamber-specific promoter usage, demonstrates sex-specific differential promoter usage in failing atria, with the female samples mostly expressing short TBX5-203 isoform and no full-length TBX5-201 and TBX5-202 (CAGE differential promoter usage FDR < 10− 5) (Fig. 4B). Thus, unlike male failing atria, female atria start to resemble ventricles under failing conditions. Interestingly, TBX5 has an associated CAGE peak located ~ 40kbp upstream without any gene models in between, and therefore it was previously assigned to TBX5 as a proximal gene3. This regulatory region may function as a promoter as well as an enhancer29. Previously, it has been classified either as a promoter30 or enhancer exclusively31,32. This TBX5 finding could provide the mechanistic basis for sex-specific susceptibility and outcome of atrial fibrillation in the settings of heart failure.
Moreover, the CAGE signal of this upstream element is located in the same topologically associated domain with the bidirectional promoter of TBX5 and TBX5-AS1 and has a significant correlation with its expression (FDR < 0.05). We found that expression of this upstream regulatory region is atria-specific in healthy males and females. However, in the failing state, its activity is significantly downregulated in female atria compared to male samples, highlighting the high sensitivity of this region to heart failure specifically in females, but not males (Fig. 4C).
We performed RNA-seq differential gene and isoform expression analyses to validate the sex-specific promoter switching, which detected 443 and 167 genes, respectively, with sex-specific expression at FDR < 0.05 (Fig. 4A). Differential gene and isoform expression results from RNA-seq significantly overlapped with the CAGE-based differential gene and promoter expression analyses (Fisher’s exact test odds ratio = 6.8, p < 10− 15 for sex-specific differential gene expression in healthy samples, odds ratio = 7.6, p < 10− 15 for genes with a sex-specific expression of promoters and RNA isoforms).
The expression of CRE modulator isoforms in the heart is sex-specific
As a particularly important case study, we selected CREM, CAMP Responsive Element Modulator, that showed significant and consistent sex-specificity in multiple tests (FDR < 10− 4 comparing all male and female samples, see CAGE tags profiles and isoform abundance estimates in Fig. 4D-E). CREM encodes a transcription factor that binds cAMP-responsive elements in many promoter regions and can either activate or repress transcription depending on a particular CREM protein isoform33. CREM encodes about 50 isoforms according to Ensembl34. However, the majority are involved in spermatogenesis and exclusively expressed in testis35. The rest of the commonly expressed CREM protein-coding isoforms36 encode repressors that block the functioning of both CREB and CREM 'activator' isoforms37–39.
All 9 protein-coding CREM isoforms out of the top 10 expressed in the heart encode repressors lack either trans-activating or DNA-binding domain (Fig. 4D). The first promoter (rDPI_5954 of TRE) corresponds to CREM-239 (ENST00000489321), CREM-229 (ENST00000474362), CREM-218 (ENST00000460270), and CREM-211 (ENST00000374726) and is active in all samples regardless of sex. The other two promoters (rDPI_5961 and rDPI_5963) are male-specific and produce CREM-237 (ENST00000488328), CREM-228 (ENST00000473940), CREM-207 (ENST00000356917), CREM-230 (ENST00000474931), and CREM-242 (ENST00000490511). This expression pattern is confirmed (Fig. 4E) by RNA-seq abundance estimates, where CREM-237 is significantly overexpressed in male compared to female hearts (log2 fold change = -1.3, FDR < 0.05).
Sex-specific TRE changes associated with age
Women tend to live longer than men, yet it remains poorly understood whether cardiac aging presents differently in males and females40,41. Thus, we aimed to examine the relationship between genome regulatory network and age. The age range of healthy male hearts is 26–61, and of healthy female hearts is 26–71. We analyzed contributions of transcription by comparing age-related changes detected in TREs for healthy male and female atria and ventricles separately (Supplementary Table 4). We found that male atria have the fewest changes in TREs expression due to age (Fig. 5A), while female atria have the most (Fig. 5B). TREs downregulated with age were related to metabolic function, while TREs upregulated with age were involved in signaling and developmental processes (Fig. 5C).
Next, we computed correlation coefficients between TREs and the age of healthy samples stratified by sex. Figure 5D shows the histogram of correlations between expression and age for all TREs and Fig. 5E for unidirectional TREs. To identify the sex-specific age-related trends, we ranked the TREs by the difference in correlation coefficients between female and male samples. We selected the top 1,000 ‘Female’ TREs (correlation of expression with age is larger in females than in males) and top 1,000 ‘Male’ TREs (correlation of expression with age is larger in males than in females) and investigated the features of these age-related TREs.
In Deviatiiarov et al., we compared TREs with known epigenetic data3. Here, we analyzed whether there was an age-related change in epigenetic modification specific to one sex (Supplementary Table 5A). We found that the most significant histone differences were associated with the DNA packaging protein H3 (H3K36me3): in the 'Female' group, 984/1,000 TREs showed no enrichment, while in the 'Male' group, 643/1,000 TREs showed no enrichment. This histone modification is responsible for maintaining gene expression stability and is vital throughout aging, impacting longevity. Genes that change their expression during aging have much lower levels of H3K36me3 in their gene bodies42, meaning aging females have potentially more genes with changing expression due to lower enrichment of H3K36me3 compared to males. Comparing the KEGG assignment for 'Female' and 'Male' TREs, we found more 'Male' TREs belonged to the PPAR signaling and Fatty acid degradation pathway43 (Supplementary Table 5B). In contrast, more 'Female' TREs belonged to the Insulin resistance pathway.
Next, we looked at individual genes with a significant correlation between age and promoter expression that also had significant differences in expression trends between healthy females and healthy males (Supplementary Results). In total, we identified 2,482 promoters passing adjusted p-value threshold of 0.05, of which 1,632 promoters had a positive expression correlation with age in females and negative in males, and only 81 promoters had a reverse trend (Supplementary Table 5C). The top 10 genes with opposite correlation (positive in females, negative in males) were crucial cardiac developmental transcription factors GATA444 and TBX2045 along with RERE46, which help to regulate genes in early development (Supplementary Fig. 5). These findings indicate the sex dimorphism of gene regulation that progresses with aging.
Differential promoter usage highlights sex-specific aging of human hearts
Similar to the sex-specific differential promoter usage (DPU) analysis but with age as the numerical variable, we identified age-dependent promoter usage in female and male atria and ventricles. The number of genes passing FDR < 0.05 was 1,100 in total and varied between the test groups ranging from 232 (for male atria) to 700 (for female ventricles) with 103 common genes, revealing both general and sex-specific effects of aging of atria and ventricles (Fig. 6A and Supplementary Table 6).
To highlight and validate the aging differences between male and female hearts, we performed the differential isoform usage (DIU) analysis using RNA-seq. Overall, we detected 162 and 205 genes with age-dependent isoform changes (FDR < 0.05 for female and male hearts, respectively) that significantly overlapped with the results of the CAGE differential promoter usage analysis (Fig. 6A, Fisher's exact test odds ratio = 3, p < 10− 3 for females, and odds ratio = 2.5, p = 0.007 for males).
To reveal the biological processes affected by heart aging in a sex-specific way, we performed gene set enrichment analysis (GSEA). Overall, 113 GO-BP terms passed FDR < 0.05 and NES (normalized enrichment score) > 1 (Fig. 6B). Differential promoter usage in female ventricles and atria was associated with membrane protein translation, muscle contraction, mitochondrial metabolism, and ribonucleoside metabolic processes. In contrast, the male age-dependent differential promoter usage was enriched with such terms as cell response to external stimulus, differentiation, developmental processes, and anatomical structure morphogenesis. GSEA results of RNA-seq differential isoform usage did not provide individual terms passing FDR < 0.05 but were generally similar (Fig. 6C). Particularly, genes with female age-dependent differential isoform usage were associated with protein localization to membrane raft (related to membrane protein translation found for differential promoter usage), muscle contraction, mitochondrial calcium transport (related to mitochondrial metabolism), and cellular response to ATP (related to ribonucleoside metabolic processes). Genes showing age-dependent differential isoform usage in males were mostly enriched by tissue morphogenesis and GPCR signaling related to developmental processes and cell response to an external stimulus detected for differential promoter usage, respectively (Fig. 6C).
We illustrate age-dependent differential promoter usage in Fig. 6D-F and Supplementary Fig. 6A-C with UCKL1 (Fig. 6D) and HAND2 (Fig. 6F), which show female- and male-specific age-dependent differential promoter usage, respectively, and CDH2 (Fig. 6E) with sex-independent aging-associated changes in the promoter usage. UCKL1 encodes a uridine kinase that catalyzes the phosphorylation of uridine to uridine monophosphate47. The relative expression of its protein-coding isoform UCKL1-203 (ENST00000369908, Fig. 6D, Supplementary Fig. 6A) is reduced with age in females: CAGE differential promoter usage FDR < 10− 5 for ventricles and < 0.01 for atria, RNA-seq differential isoform usage p-value < 0.04 (compared to male DPU FDR = 0.4 for ventricles, 0.4 for atria and DIU p-value = 0.3). HAND2 encodes a transcription factor essential for cardiac morphogenesis. It demonstrates a significant male-specific age-dependent increase in the relative abundance of its short isoform HAND2-202 (ENST00000503024, Fig. 6F, Supplementary Fig. 6C) that is predicted to go through nonsense-mediated decay: CAGE differential promoter usage FDR < 0.2 for ventricles and < 0.007 for atria, RNA-seq differential isoform usage p-value < 0.03 (compared to female DPU FDR = 0.5 for ventricles, 0.9 for atria and DIU p-value = 0.9)48. CDH2 encodes cell adhesion protein cadherin that mediates adherens junctions connecting actomyosin networks of cardiomyocytes49. In contrast to UCKL1 and HAND2, this gene demonstrates no sex-specificity, but the fraction of short isoform CDH2-211 increases with age in both males and females (ENST00000676445, Fig. 6E, Supplementary Fig. 6B, CAGE differential promoter usage FDR < 10− 4 for female and < 10− 6 for male ventricles as well as < 10− 5 for female and < 5*10− 4 for male atria, RNA-seq differential isoform usage p-value < 0.09 for females and < 0.07 for males).