Proteogenomics reveals sex-biased aging genes and coordinated splicing in cardiac aging

Abstract

The risks of heart diseases are significantly modulated by age and sex, but how these factors influence baseline cardiac gene expression remains incompletely understood. Here, we used RNA sequencing and mass spectrometry to compare gene expression in female and male young adult (4 mo) and early aging (20 mo) mouse hearts, identifying thousands of age- and sex-dependent gene expression signatures. Sexually dimorphic cardiac genes are broadly distributed, functioning in mitochondrial metabolism, translation, and other processes. In parallel, we found over 800 genes with differential aging response between male and female, including genes in cAMP and PKA signaling. Analysis of the sex-adjusted aging cardiac transcriptome revealed a widespread remodeling of exon usage patterns that is largely independent from differential gene expression, concomitant with upstream changes in RNA-binding protein and splice factor transcripts. To evaluate the impact of the splicing events on cardiac proteoform composition, we applied an RNA-guided proteomics computational pipeline to analyze the mass spectrometry data and detected hundreds of putative splice variant proteins that have the potential to rewire the cardiac proteome. Taken together, the results here suggest that cardiac aging is associated with 1) widespread sex-biased aging genes and 2) a rewiring of RNA splicing programs, including sex- and age-dependent changes in exon usages and splice patterns that have the potential to influence cardiac protein structure and function. These changes contribute to the emerging evidence for considerable sexual dimorphism in the cardiac aging process that should be considered in the search for disease mechanisms.NEW & NOTEWORTHY Han et al. used proteogenomics to compare male and female mouse hearts at 4 and 20 mo. Sex-biased cardiac genes function in mitochondrial metabolism, translation, autophagy, and other processes. Hundreds of cardiac genes show sex-by-age interactions, that is, sex-biased aging genes. Cardiac aging is accompanied with a remodeling of exon usage in functionally coordinated genes, concomitant with differential expression of RNA-binding proteins and splice factors. These features represent an underinvestigated aspect of cardiac aging that may be relevant to the search for disease mechanisms.

Keywords: aging; alternative splicing; proteoforms; proteogenomics; sex difference.

Figure 1.

Study design. A: animal groups and schema of proteogenomics workflow. B: body mass, heart mass, and anterior wall thicknesses in each group (n = 5 biological replicates/group). P values: Student’s t test. C and D: echocardiographic measures associated with systolic (C) and diastolic (D) function.

Figure 2.

Age- and sex-associated cardiac genes. A: scatter plot of first two principal components of total transcriptome profiles showing linear separation of samples by age and sex groupings (n = 3 per group). B: point ranges of log2 fold changes and standard errors (y-axis) in male (M) vs. female (F) C57BL/6J hearts, showing the top 50 significant genes (x-axis). Point glyphs denote sex chromosome/autosome of the gene in mice. Color: log2 fold changes (FC). C: point ranges of log2 fold changes and standard errors (y-axis) in sex-adjusted 20-mo vs. 4-mo C57BL/6J hearts, showing the top 50 significant genes (x-axis). D: overrepresentation analysis against Reactome pathways (y-axis) among differentially expressed genes (10% FDR of |log2FC| ≥ 10%) in male (M) vs. female (F) hearts. Red dashed line, 10% FDR; point color, –log10-adjusted P; point size, number of quantified genes in a pathway. Red pathways are significant within 10% FDR. E: pathway diagram for Reactome DAG and IP₃ signaling pathway genes. Associated dot plots show the normalized read counts of sexually dimorphic genes. P values are from individual t test for reference; refer to Supplemental Data S1 for DESeq2 model P values. Genes that are more highly expressed in female are in blue; male, in red. F: overrepresentation against Reactome pathways (y-axis) among differentially expressed genes (10% FDR of |log2FC| ≥ 10%) in 20-mo vs. 4-mo hearts.

Figure 3.

Cardiac genes and pathways with significant age-by-sex interaction. A: heat maps showing the top 60 sexually dimorphic aging genes that show significant differences in aging response across sexes (DESeq2 age:sex 10% s-value). Colors denote standardized expression (Z score of normalized read counts). B: normalized read counts of four genes from panel a were visualized across groups individually. P values represent individual t tests in M vs. F comparison in 4-mo and 20-mo groups for reference; refer to Supplemental Data S3 for DESeq2 model P values. C: gene-concept network of 2 groups of significantly enriched Reactome terms among genes with significant age-by-sex interaction. Brown nodes are Reactome terms, linked through edges to their annotated genes in the foreground of the functional enrichment analysis. Gene node colors, log2 FC in male aged/adult vs. female aged/adult; Reactome term node size, gene count; edge color, Reactome terms. D: gene-concept network involving top enriched term in gene set enrichment analysis of the ordered gene list (by log2FC, male:aged vs. female:aged), linked to their annotated genes as in C.

Figure 4.

Comparisons of proteomics and transcriptomics data. A: running enrichment score (ES) plot of significant subcellular localizations in gene set enrichment analysis (GSEA) of the quantitative mass spectrometry comparisons of 20-mo vs. 4-mo hearts. Proteins were ordered from highest to lowest log2 fold change (log2FC). Left to right: age-associated proteins are enriched in extracellular (trending upregulated), mitochondrial (downregulated), nuclear speckle (downregulated), and actin cytoskeletal (upregulated) localizations. B: scatter plot of log2FC among differentially regulated proteins (limma P ≤ 0.05) and the corresponding log2FC of their corresponding mRNA. Data are divided into 4 quadrants based on whether they show concordant or discordant protein and mRNA changes as labeled. Point colors, data density; bars, mRNA/protein counts in each quadrant. C: bar charts showing significantly enriched Reactome terms (at 10% FDR) among proteins in each quadrant; x-axis, –log10 FDR-adjusted P value, Fisher’s exact test; bar color, log2 fold enrichment against all quantified proteins.

Figure 5.

mRNA-processing and splicing-related pathways in sex-adjusted aging hearts. A: volcano plot showing the limma –log10 P value (y-axis) vs. log2 fold change (x-axis) of gene set variance analysis (GSVA)-extracted gene sets. Gene sets related to mRNA splicing are labeled and colored in red and labeled with letters corresponding to the annotated pathways printed. B: volcano plot showing the limma –log10 P value (y-axis) vs. log2 fold change (x-axis) of PLIER-learned latent variables. Latent variables related to mRNA splicing are labeled and colored in red and labeled with letters corresponding to the annotated pathways printed. C: heat maps showing a subset of genes with highest feature scores in the mRNA splicing associated latent variables that are differentially regulated in aged hearts, i.e., genes belonging to mRNA processing and splicing pathways that are represented in latent variables 437 and 54 in B. The left matrix denotes the membership of the gene in the PLIER prior knowledge matrix (annotated pathways); colors in the right heat maps denote relative read counts in 20-mo vs. 4-mo female (F) and male (M) hearts. D: individual plots of normalized gene counts across each group among selected significant transcripts encoding RNA-binding proteins and splice factors. Numbers represent DESeq2 model s-value in 20-mo and 4-mo heart comparisons after adjusting for sex. See Supplemental Data S2 for details.

Figure 6.

Differential exon usages in aging hearts. A–C: exon usages in sex-adjusted comparisons between 20-mo (blue) and 4-mo (red) mouse hearts. Y-axis shows normalized read counts across disjoint exonic parts in DEXSeq after adjusting for gene-level differential expression. The disjoint exonic parts in DEXSeq are mapped to annotated exons in the Ensembl gene model below with the genomic coordinates labeled. Exonic parts with significant differential usage in 20-mo vs. 4-mo mouse hearts are colored in purple. D: intersection between genes with differential exon usage (DEU) at 20 mo vs. 4 mo (DEU_AGE) at 10% FDR with genes with differential gene expression (DGE) across age and sex comparisons. E: dot plot showing the enrichment ratio and adjusted P values of overrepresented Reactome pathways among genes with age-associated exon usage difference. Exon usage changes are enriched along RNA metabolism, ribosomal, and mitochondrial metabolic pathways. F: gene concept network showing top 15-enriched Reactome terms and their associated genes with significant differences in exon usage in 20-mo vs. 4-mo hearts.

Figure 7.

Protein-level consequences of exon usage changes. A: Sankey diagram illustrating the fate of all modeled splice junctions in the RNA-guided protein database methods. Node sets from left to right denote the modeled splice type of the exon junctions, whether the skipped junction count is above a threshold modeled from total junction read counts, translation status (see methods), and identification status. Inset: modeled read count cutoff. B, top: peptide spectrum match (PSM) of a noncanonical peptide belonging to titin from a database search using the RNA-guided protein database in A. Matched peptide fragment b- and y-ions are labeled and colored in red. PEP, Percolator PSM posterior error probability. B, bottom: best-hit candidate sequence matched to the identical spectrum (top) when queried against UniProt SwissProt canonical and isoform entries. C: differential exon usage of Acsl1 in 20-mo vs. 4-mo mouse heart. A disjoint exonic part E030 (green) shows higher usage in aged hearts and maps to an exon absent in the canonical transcript (Acsl1-201). D: peptide-spectrum match of a noncanonical peptide corresponding to a cassette exon matching to E030 identified in the quantitative mass spectrometry data using the RNA-guided protein database. E: predicted protein structure of the predicted full-length alternative ACSL1 containing the translated alternative region from which the noncanonical peptide was identified, joined to the full-length UniProt canonical sequence via a 10-amino acid joint in the upstream and downstream exons. E, left: predicted protein structure with the alternative region highlighted. E, right: predicted structure of the alternative region overlaid on the predicted structure of the canonical protein.