Sexual-dimorphism in human immune system aging

Abstract

Differences in immune function and responses contribute to health- and life-span disparities between sexes. However, the role of sex in immune system aging is not well understood. Here, we characterize peripheral blood mononuclear cells from 172 healthy adults 22-93 years of age using ATAC-seq, RNA-seq, and flow cytometry. These data reveal a shared epigenomic signature of aging including declining naïve T cell and increasing monocyte and cytotoxic cell functions. These changes are greater in magnitude in men and accompanied by a male-specific decline in B-cell specific loci. Age-related epigenomic changes first spike around late-thirties with similar timing and magnitude between sexes, whereas the second spike is earlier and stronger in men. Unexpectedly, genomic differences between sexes increase after age 65, with men having higher innate and pro-inflammatory activity and lower adaptive activity. Impact of age and sex on immune phenotypes can be visualized at https://immune-aging.jax.org to provide insights into future studies.

Conflict of interest statement

The authors declare no competing interests.

Figures

Fig. 1. Age is the main driver of variation in PBMC genomic data.

a Schematic summary of study design: PBMCs were isolated from blood samples of 172 healthy community-dwelling individuals (ages 22–93), from which ATAC-seq, RNA-seq and flow-cytometry data were generated. b Principal component 1 scores (PC1) were calculated for each individual from ATAC-seq (left) and RNA-seq (right) principal component analyses (PCA) results. Marker genes are selected using PC1 scores that are ≥25th percentile: top positive scores (ATAC-seq n = 9392 peaks, RNA-seq n = 1413 genes) and top negative scores (ATAC-seq n = 10,151 peaks, RNA-seq n = 1675 genes). PC1 scores from ATAC-seq and RNA-seq data increased with increasing age. Furthermore, we detected differences in PC1 scores of sexes in older subjects, where older men had higher scores than older women. c Functional enrichment of PC1-related genes based on immune modules. Note that myeloid/inflammation-related genes were associated with high and positive scores, whereas adaptive immunity/lymphoid related genes were associated with high and negative PC1 scores. These enrichments align well with an age-related increase in the myeloid lineage and inflammation, and an age-related decline in naive T cell activity. Hypergeometric test was used to calculate enrichment p-values. d Flow cytometry data from 128 individuals that reflect proportions of major cell types within PBMCs. Subjects were stratified based on age group and sex: Young: 22–40, Middle-aged: 41–64, Older: 65+ years old. Wilcoxon rank-sum test was used to compare data from female (n = 67) and male (n = 61) subjects. Note that T cell proportions decline with age in both sexes, whereas the decline in B cell proportions is specific to older men. Box plots represent median and IQR values, whiskers extend to 1.5 times the IQR; ***p < 0.001, **p < 0.01, *p < 0.05, n.s.: non-significant. Source data are provided as a Source Data file for (c).

Fig. 2. Shared and sex-specific epigenomic signatures of aging.

a Fold change (older-young) vs. read counts per million (logCPM) from ATAC-seq differential analyses in women (left) and men (right) based on 78,167 ATAC-seq peaks (5% FDR after Benjamini-Hochberg P-value adjustment, 58 subjects). Differentially opening/closing peaks are represented in orange/purple. b Functional annotations of differentially accessible (DA) peaks in men and women using chromHMM states in PBMCs and cell subsets. Closing peaks tend to be enhancers/promoters across all cell types, whereas opening peaks tend to be enhancers in monocytes and NK cells. c Distribution of DA peaks among cell-specific loci from chromHMM states. Innate cell-specific loci (monocytes, NK cells) are enriched in opening peaks, T cell-specific loci are enriched in closing peaks. B cell-specific loci are more likely to be opening peaks in women and closing peaks in men. d Left: Average chromatin accessibility profiles (grouped by age and sex) around the IL7R locus—is associated with chromatin closing with age in women (top, n = 32) and men (bottom, n = 26). Red bars indicate peaks and their significance value from differential analyses (*p < 0.05, **p < 0.01, ***p < 0.001, no symbols: non-significant). Right: Normalized chromatin accessibility (combined from three peaks shown in the figure) and gene expression levels for IL7R in young (n = 16) and older (n = 16) women (top) and young (n = 12) and older (n = 14) men (bottom). In both sexes, chromatin accessibility and gene expression levels decline with age albeit more significantly in men. e Enrichment of differentially expressed (DE) genes using cell-specific gene sets from scRNA-seq data. Enrichment p-values are calculated using hypergeometric test. Gene expression programs for innate cells (NK, monocytes) are activated, whereas programs for T cells are inactivated in both sexes. B-cell specific genes were downregulated specifically in men. f Average expression levels of T cell-specific genes grouped by age group and sex. Note the decline in both sexes. Box plots on the top summarize the data from all genes represented in the heatmap. All box plots represent median and IQR (interquartile range) values, with whiskers extending to 1.5 times the IQR. Source data are provided for (a, b, c, e).

Fig. 3. Sex-dimorphic changes in monocyte- and B cell-associated loci.

Correlation of age-related ATAC-seq (a) and RNA-seq (b) remodeling between women and men PBMCs. Note the overall large and positive Pearson correlation coefficients. Genes are associated to ATAC-seq peaks based on nearest TSS, and are color coded in both plots according to their association to immune modules (purple: T cells, green: B cells, pink: NK cells, yellow: monocytes). Only regulatory (TSS/enhancer) peaks are included and peaks-genes are matched between both plots (n = 10,707 loci). Blue-red gradient on data points represents their relative local density. (c) Correlation between sexes for age-related ATAC-seq remodeling stratified by cell-specific loci from chromHMM annotations. Note that the highest correlation is observed in naive T cells (n = 833 peaks), which is associated with negative fold changes (i.e., loss in chromatin accessibility with age) in both sexes. Chromatin remodeling correlates the least between sexes for B cell- (n = 1645 peaks) and monocyte-specific loci (n = 6861 peaks), with NK cells (n = 3,008 peaks) showing a positive trend in both sexes. d Correlation between sexes for age-related RNA-seq remodeling stratified by cell-specific genes from single-cell RNA-seq data. Note that the highest correlation is observed in naive T cells (n = 70 genes), which is associated with negative fold changes (i.e., downregulation with age) in both sexes. On the other hand, NK cells (n = 403 genes) are highly correlated between sexes and associated with increased expression with age in both sexes. In agreement with ATAC-seq data, gene expression remodeling correlates the least between sexes for B cell- (n = 144) and monocyte-specific genes (n = 748). Values at the top left corner represent the Pearson correlation coefficient between sexes only for the genes/loci associated to that cell type. e Normalized expression and accessibility levels for important molecules associated to T, cytotoxic, and B cells. LogFC = mean log2 of fold changes between old and young samples; MA = middle-aged. Box plots represent median and IQR values, with whiskers extending to 1.5 times the IQR. Source data are provided as a Source Data file for (a–d).

Fig. 4. Epigenome and transcriptome changes over human adult lifespan.

a Heatmap of ATAC-seq data (fitted values from ARIMA models) with a significant chronological trend in women (left, n = 13,297) and men (right, n = 13,295), as a function of age in years. Values represent z-score normalized accessibility values relative to the row (i.e., peak) mean. K-means clustering was used to group these peaks into three clusters in men and women (M1/F1, M2/F2, M3/F3). Color bar on the top represents discrete age groupings as defined in this study (young, middle-aged, older). Rows are annotated according to their position relative to the nearest TSS: proximal if <1 kbp distance, distal otherwise. b ChromHMM state annotations of temporal peak clusters in women (F1–F3) and men (M1–M3). Colors represent hypergeometric enrichment test p-values; light gray cells indicate insufficient number of genes to run an enrichment test. c Heatmap of ARIMA-fitted expression values for genes with a significant chronological trend in women (left, n = 1068) and men (right, n = 1471), as a function of age in years. Three clusters per sex are identified using k-means clustering: F1–3 and M1–3. Values represent z-score normalized expression values relative to the row (i.e., gene) mean. d Annotation of temporal gene clusters using cell-specific gene sets derived from single-cell RNA-seq data. Colors represent hypergeometric enrichment test p-values; light gray cells in plot indicate that there were insufficient genes in the gene sets to run an enrichment test. e Inverse log p-value distributions from breakpoint analysis for each cluster for women (top) and men (bottom), where curve height indicates magnitude of differences between preceding and succeeding age windows. Note that there are two age brackets where epigenomic changes take place abruptly both in men and women. Points and vertical lines mark median age estimates for a breakpoint integrated over multiple scales, and C1–C3 correspond to temporal clusters M1–M3 or F1–F3 from each sex (see the “Methods” section and Fig. S5c for details). Source data are provided as a Source Data file for (a–e).

Fig. 5. Sex-specific patterns in temporal peaks/genes.

a Heatmap of ATAC-seq peaks (fitted values from ARIMA models) with a chronological trend in both women and men (n = 3197), as a function of age in years. Values represent z-score normalized accessibility values relative to the row (i.e., peak) mean. K-means clustering was used to group these peaks into three clusters (C1–C3) using concatenated data from men and women. Color bar on the top represents discrete age groupings as defined in this study (young, middle-aged, older). Rows are annotated according to their position relative to the nearest TSS: proximal if <1 kbp distance, distal otherwise. b, c Annotations of shared temporal peaks using chromHMM states (b) and gene sets from DICE database (c). Colors represent hypergeometric enrichment test p-values; light gray cells in plot indicate that there were insufficient genes in the gene sets to run an enrichment test. The enrichment pattern strongly associates cluster C1 to T cells, suggesting a delayed loss of accessibility in women relative to men; C2 to CD19+ cells, suggesting the presence of CD19+ specific loci with opposing temporal behavior in men and women; and C3 to monocytes and NK cells. d Transcription factor (TF) motif enrichment results for each temporal cluster (C1–C3), relative to the other two clusters. Motif enrichment analyses carried out on 1388 PWMs, grouped into families based on the sequence similarity, and most significant p-value for each motif family is represented here. Tests were done using HOMER. e, f Expression levels of TFs associated to cluster 1 (C1) and cluster 3 (C3) grouped by age group and sex whose expression follows the same pattern as the peak temporal clusters where they are enriched. Cluster 2 (C2) is omitted since all TFs in this group show a significant increase with age in females or both sexes. Box plots represent median and IQR values, with whiskers extending to 1.5 times the IQR. Wilcoxon rank-sum test used to compare expression levels between sexes (significance value below boxes) and age groups (above boxes): *p < 0.05, **p < 0.01, ***p < 0.001, ns: non-significant. Sample sizes for young individuals n = 11 F, 6 M, middle-aged n = 10 F, 20 M, older n = 13 F, 14 M. Source data are provided as a Source Data file for (a, b, d).

Fig. 6. Genomic differences between sexes at different age groups.

MA plots representing mean log2 fold change (male-female) vs. average log2 read counts per million reads (logCPM) for ATAC-seq (a) and RNA-seq (b) data at three age groups using 78,167 peaks (12,199 genes) at autosomal chromosomes. Both epigenomic and transcriptomic differences between sexes increase with age. Peaks and genes significantly upregulated in women (men) relative to men (women) are represented in red (blue). Fold changes obtained via GLM modeling of the data sets; statistical significance assessed at a 5% FDR threshold based on Benjamini-Hochberg p-value adjustment. Tests based on n = 100 (ATAC-seq) and n = 74 (RNA-seq) independent samples. (c) Enrichment of significantly sex-biased peaks/genes in older individuals using cell-specific gene sets obtained from single-cell RNA-seq data. Note the male bias toward increased accessibility/expression for monocytes and DCs and the female bias toward increased accessibility/expression for T and B cells. P-values are based on hypergeometric enrichment tests. Numbers on bars represent the number of differential genes overlapping each gene set. d Left: Selected pathways/module enrichments for male-biased and female-biased genes/loci. The complete list of enrichments is presented in Tables S12, S13. * represents annotations obtained from ATAC-seq data analyses. Right: Arrows indicating whether the same pathway/module has been significantly associated with age-related changes. Significant enrichments (FDR 5%) and their directions are represented with arrows. Red arrows for women, blue arrows for men. NS = not significant. The complete list of enrichments is presented in Tables S6, S8. e ELISA data from 500 Functional Genomics (500FG) project to measure serum levels of pro-inflammatory cytokines (n = 267). Protein expression levels are compared between sexes and between young and older individuals. P-values were calculated using Wilcoxon rank sum (two sided). Box plots represent median and IQR values, with whiskers extending to 1.5 times the IQR. Source data are provided as a Source Data file for (a–c).