Multi-omic profiling of primary mouse neutrophils predicts a pattern of sex- and age-related functional regulation

Abstract

Neutrophils are the most abundant human white blood cell and constitute a first line of defense in the innate immune response. Neutrophils are short-lived cells and thus the impact of organismal aging on neutrophil biology, especially as a function of biological sex, remains poorly understood. Here, we describe a multi-omic resource of mouse primary bone-marrow neutrophils from young and old female and male mice, at the transcriptomic, metabolomic and lipidomic levels. We identify widespread regulation of neutrophil ‘omics’ landscapes with organismal aging and biological sex. In addition, we leverage our resource to predict functional differences, including changes in neutrophil responses to activation signals. This dataset represents a large multi-omics resource for neutrophils across sex and ages and identifies neutrophil characteristics that could be targeted to improve immune responses as a function of sex and/or age.

Main

Neutrophils are the most abundant cells in human blood, representing 50–70% of leukocytes¹. These cells are continually produced in the bone marrow and released into circulation^2,3. Neutrophils are short-lived, with estimated cellular lifespans of ~10–18 h once released in the bloodstream^4,5, although they can survive longer^4,6. Throughout their cellular lifespan, neutrophils undergo ‘neutrophil aging’, a process distinct from organismal aging^3,7. Neutrophils perform many key functions, including production of antimicrobial granules and of ‘neutrophil extracellular traps’ (NETs)^3,8. Although neutrophils are essential against infections as the ‘first line of defense’, their aberrant activation can also aggravate inflammatory disease^2,3. Indeed, emerging evidence suggests that neutrophils play important roles in chronic inflammation⁹.

Organismal aging in mammals is characterized by systematic immune dysfunction and chronic inflammation, a phenomenon dubbed ‘inflamm-aging’¹⁰, which may be partially mediated by dysfunction of innate immune cells¹¹. Indeed, emerging evidence suggests that neutrophils from aged organisms are dysfunctional^{12,13,14,15,16}. Observed age-related dysfunctions in neutrophils include decreased NETosis in tumor necrosis factor (TNF)-α-primed conditions^13,14, reduced chemotaxis to sites of inflammation¹⁵ and secretion of intracellular granule proteases^15,17. Although gene expression changes throughout the lifespan have been reported across many cell types^18,19, how organismal aging (rather than ‘daily’ cell aging) affects neutrophil landscapes remains largely unknown.

Females and males present with many biological differences, which may underlie lifelong health disparities between sexes²⁰ and could result from differential ‘omic’ regulation^21,22. Although transcriptional differences between young female and male murine neutrophils have been profiled by ImmGen²², how these differences interplay with organismal aging and whether they are accompanied by other phenotypical differences remains largely unknown. Accumulated studies suggest that aspects of neutrophil biology are sex dimorphic, such as inflammatory mediator production⁹ or functional modulation by testosterone²³. However, the pathways underlying sex dimorphism in neutrophils, as well as the extent of sex dimorphism, are still unclear.

To gain insights into how neutrophils are regulated as a function of age and sex, we generated a multi-omic resource covering transcriptome, metabolome and lipidome profiling of primary bone-marrow mouse neutrophils. We identified widespread age-related and sex-dimorphic ‘omic’ regulation, including transcriptional regulation of chromatin-related pathways. Using the assay for transposase-accessible chromatin using sequencing (ATAC-seq), we showed that remodeling of chromatin-related pathways was associated with overall differences in the chromatin architecture of neutrophils from young versus old and female versus male mice. Consistently, we observed age- and sex-linked differences in NETosis inducibility. Machine learning showed that specific factors could predict age-related and sex-dimorphic gene regulation in neutrophils. Finally, we leveraged our resource and predicted sex differences in serum levels of neutrophil elastase in control and sepsis-like conditions.

Results

Multi-omics of bone-marrow neutrophils with age and sex

To understand how neutrophils are regulated as a function of age and sex, we obtained primary bone-marrow neutrophils from young (4–5-month-old) and old (20–21-month-old) C57BL/6Nia female and male mice (Fig. 1a). Primary neutrophils were isolated from bone marrow using magnetic-activated cell sorting (MACS) and profiling was then performed on purified neutrophils: (1) transcriptome profiling by RNA-sequencing (RNA-seq), (2) metabolomic profiling by untargeted liquid chromatography coupled with mass spectrometry (LC–MS) and (3) lipidomic profiling by targeted MS (Fig. 1a).

Fig. 1: A multi-omic analysis of primary mouse bone marrow neutrophils during aging and with respect to sex.

a, Experimental setup scheme. b–d, MDS analysis results of RNA expression by RNA-seq (b), untargeted metabolomics (c) or targeted lipidomics (d). e–g, Heat map of significant (FDR < 5%) sex-dimorphic genes (e), metabolic features (f) or lipid species (g). Significance of gene regulation by RNA-seq was calculated by DESeq2 and significance of metabolic features or lipid species regulation were calculated by limma.

As a first-level analysis to evaluate the similarity of our datasets, we utilized multidimensional scaling (MDS). MDS analysis for RNA-seq, metabolomics and lipidomics datasets showed clear separation of samples by animal sex, regardless of age (Fig. 1b–d). In contrast, although young and old samples separated within each sex, global separation by age regardless of sex was not clearly observed for each omics (Fig. 1b–d). To better understand the nature of differences between neutrophils from young versus old (with sex as a covariate) and female versus male animals (with age as a covariate), we identified transcriptional, metabolic and lipidomic features with significant age- or sex-related regulation at a false discovery rate (FDR) < 5% using multivariate linear modeling (Fig. 1e–g, Extended Data Fig. 1a–f and Supplementary Table 1a–f). We quality-checked our dataset for appropriate expression of sex-specific genes (Extended Data Fig. 1b,c). Finally, we confirmed differential gene expression trends between groups using a small replicate RNA-seq cohort (Extended Data Fig. 2a,b) and comparing our data with published datasets from female versus male mouse spleen neutrophils²² and human blood neutrophils²⁴ (Extended Data Fig. 2a). Thus, our analyses suggest that observed transcriptional differences with respect to organismal age and sex in neutrophils are robust.

Consistent with MDS, we identified genes, metabolic features and lipids with significant differential regulation with respect to sex (FDR < 5%; Fig. 1e–g, Extended Data Fig. 1a–d and Supplementary Table 1b,d,e). Differentially expressed genes with respect to sex were located throughout the genome, suggesting that sex-dimorphic gene expression was not just a byproduct of genomic location (Extended Data Fig. 1a,e). Notably, we observed no biases in the purity of MACS-purified neutrophils across groups (Extended Data Fig. 3a,b) and no difference in neutrophil abundance in bone marrow (Extended Data Fig. 3c) across groups. Thus, sex-dimorphic and age-related ‘omic’ phenotypes are unlikely to result from a systematic purification bias between groups.

Although aging led to clear changes in neutrophil gene expression profiles (Extended Data Fig. 1a,f and Supplementary Table 1), few to no metabolomic and lipidomic changes were observed (Extended Data Fig. 1a and Supplementary Table 1c,f). This is consistent with overall poorer age-related separation observed for metabolomics and lipidomics (Fig. 1c,d). Higher biological variability in metabolites and lipids might explain the lack of changes with age. Alternatively, it is possible that the paucity of age-related metabolic and lipidomic differences may stem from biological specificities of neutrophils (such as short cellular lifespan). When female and male neutrophils were analyzed separately, we found that the transcriptional impact of aging on neutrophils correlated strongly between sexes (Spearman Rho = 0.593 and P ~ 0 in significance of correlation test; Extended Data Fig. 1g). Thus, although neutrophils show clear sex differences throughout life, the trajectories of neutrophils with aging are highly similar regardless of sex. To note, we observed a greater number of transcriptional changes with organismal age in male versus female neutrophils (Extended Data Fig. 1h), suggesting that male neutrophils are more susceptible to aging. Together, our results indicate that neutrophils are highly sex dimorphic at the molecular level and that these sex differences persist (or may be amplified) with organismal aging.

Finally, as bone-marrow neutrophils can be heterogenous, we evaluated the composition of MACS-purified neutrophils across groups (Extended Data Fig. 4a–d and Supplementary Table 2a,b). We leveraged flow cytometry markers to estimate the proportions of neutrophil progenitors (pre- and pro-neutrophils), immature versus mature neutrophils and fresh versus ‘aged’ neutrophils in our population. Neutrophil progenitors were extremely rare among MACS-purified cells (pre- and pro-neutrophils represent <0.5% of cells; Extended Data Fig. 4a and Supplementary Table 2a) and are unlikely to drive large ‘omic’ differences. Notably, we did not observe significant differences in the presence of mature versus immature neutrophils across groups (Extended Data Fig. 4a,c and Supplementary Table 2). Finally, we observed a significant trend for the increased presence of senescent/aged neutrophils (CD62L⁻) among mature neutrophils (Extended Data Fig. 4a,d and Supplementary Table 2a). Our observation is consistent with a previous report of accumulation of senescent/aged neutrophils in old mice due to decreased clearance by efferocytosis²⁵. Based on previous studies of senescent/aged neutrophils, changes in the relative proportion of ‘fresh’ versus ‘aged’ neutrophils among bone-marrow neutrophils may have functional impacts on inflammation, chemotaxis and NETosis^26,27.

Age-associated changes in chromatin pathways in neutrophils

We first focused on the impact of organismal aging on neutrophils using pathway enrichment analysis (Fig. 2a–c, Extended Data Fig. 5a,b and Supplementary Table 3). To note, we only analyzed age-related functional remodeling at the transcriptomic level, as metabolomics and lipidomics showed few changes with aging (Extended Data Fig. 1a).

Fig. 2: Regulated pathways in bone-marrow neutrophils during aging reveals downregulation of chromatin-related pathways.

a, NetworkAnalyst predicted PPI network for significant age-regulated genes in neutrophils. Network edges are based on IMEx/InnateDB data, a knowledgebase specifically geared for analyses of innate immune networks¹¹⁹. Blue (red) shades are associated with decreased (increased) gene expression during aging. b,c, Top enriched gene sets from Reactome (b) and transcription factor targets (c) using GSEA for differential RNA expression. Only the top ten most-upregulated and top ten most-downregulated gene sets are plotted for readability. Full lists and statistics available in Supplementary Table 3. Also see Extended Data Fig. 5. Shown are pathways with FDR < 5%. NES, normalized enrichment score.

We assessed the coordination of age-related changes by performing a network analysis (Fig. 2a). The network was constructed using age-regulated genes and protein–protein interaction (PPI) information from IMEx/InnateDB. Our analysis revealed that significantly age-regulated genes form a clearly interconnected network (Fig. 2a). The strong interconnection of age-related genes is consistent with the notion that coordinated shifts occur in neutrophils with organismal aging, rather than stochastic changes and suggests that these genes are co-regulated or are involved in common processes.

We then asked which pathways were regulated with aging (Fig. 2b, Extended Data Fig. 5a,b and Supplementary Table 3). Immune-related pathways were significantly upregulated according to gene set enrichment analysis (GSEA) (Fig. 2b, Extended Data Fig. 5a,b and Supplementary Table 3a–c) and ingenuity pathway analysis (IPA) (Supplementary Table 3e). This suggests that neutrophil-associated immunity in older animals could have different characteristics, consistent with age-related immune dysfunction¹⁰.

Notably, the top downregulated pathways with organismal aging were overwhelmingly related to chromatin and cell cycle (Fig. 2b, Extended Data Fig. 5a,b and Supplementary Table 3a–c,e). Specifically, 18 histone-encoding genes were significantly downregulated with aging (Supplementary Table 1a). We did not observe significant differences in the cell-cycle phase distribution of MACS-isolated neutrophils across groups (Extended Data Fig. 6a,b), suggesting that differential expression of cell-cycle pathways is not a byproduct of differential proliferation.

Changes in regulation of chromatin components are especially relevant to neutrophil biology. Indeed, neutrophils have a unique chromatin organization²⁸, with increased chromatin compaction or ‘supercontraction’²⁹. In addition, neutrophils can extrude their chromatin to produce NETs in response to pathogens³⁰. The extrusion of chromatin directly participates in pathogen killing, notably thanks to antimicrobial properties of histones³¹. Chromatin relaxation can also be potentiated by nuclear-translocated granule-derived proteases (such as Elane), ultimately helping bacterial killing³². Finally, although neutrophils are post-mitotic, cell-cycle genes can control NET production³³. Thus, our analysis suggests that neutrophils may experience changes in chromatin organization with organismal aging, a feature that may impact their ability to respond to infection.

We also observed that several autophagy-related pathways were upregulated in old neutrophils (Fig. 2b and Supplementary Table 3b,c,e). This is consistent with loss of proteostasis being a ‘hallmark of aging’³⁴. Control of autophagy is critical for neutrophil differentiation³⁵, regulation of NET formation³⁶ and of degranulation³⁷.

To identify candidate upstream regulators, we investigated whether target genes of specific transcription factors (TFs) were differentially regulated with aging (Fig. 2c and Supplementary Table 3d). We obtained lists of TF target genes derived from ChEA, JASPAR or Gene Expression Omnibus (GEO) through the Harmonizome (Methods). Notably, we restricted our testing to TFs with detectable expression in neutrophils according to RNA-seq. To note, the GEO ‘TF’ set from Harmonizome includes several non-TF regulators (such as chromatin-remodeling enzymes and signaling receptors), which are referred to hereafter as TFs for simplicity. Consistent with functional pathway enrichments, we observed significant age-related downregulation of E2f7 targets, which are known to regulate cycle-related expression^38,39 (Fig. 2c). Although E2f7 has not been studied in neutrophils, it is highly expressed in committed progenitors⁴⁰. Targets of Foxo1 were also significantly downregulated with aging in neutrophils (Fig. 2c). Foxo TFs are the primary targets of aging-regulating insulin/insulin-like growth factor 1 signaling⁴¹ and Foxo1 regulates neutrophil-mediated bacterial immunity⁴². Finally, known targets of Padi4 were significantly downregulated with aging (Fig. 2c). Padi4 is a peptidylarginine deaminase, catalyzing histone citrullination, which plays a key role in NETosis in vivo^43,44. Transcript levels of Padi4 were significantly upregulated with aging (FDR = 6.9 × 10⁻³; Extended Data Fig. 5c and Supplementary Table 1b). Protein levels of Padi4 were also significantly upregulated with aging regardless of sex (P = 8.5 × 10⁻³; Extended Data Fig. 5d,e). This is consistent with reports of histone citrullination associating to transcriptional repression^45,46,47, although the directionality of regulation may depend on cofactors and cell types⁴⁸. Together, our functional pathway analysis suggests that major aspects of chromatin metabolism are shifting in aging bone-marrow neutrophils.

Sex dimorphism in chromatin pathways in neutrophils

We next analyzed the functional impact of neutrophil sex dimorphism at the transcriptomic, metabolomic and lipidomic levels (Fig. 3a–f, Extended Data Fig. 7a–c and Supplementary Table 4). We first assessed the interconnection of sex-dimorphic genes by constructing putative PPI networks for female versus male-biased gene expression (Fig. 3a,b). Sex-biased gene programs showed clear interconnection, suggesting coherent differences with likely functional impact (Fig. 3a,b). The top-connected node in the female-biased gene network, Iqcb1, encodes a primary cilia component and has been linked to severe kidney disease⁴⁹ (Fig. 3a). Notably, Irf8 was the top connected node in the male-biased gene network (Fig. 3b). Although not specifically studied in neutrophils, Irf8 is an interferon signaling-related TF that regulates myeloid fate determination, usually suppressing neutrophil differentiation^50,51. In this context, Irf8 is likely to accomplish functions unrelated to myeloid differentiation (such as cytokine production)⁵², as no significant changes in relative bone-marrow neutrophil abundances were observed (Extended Data Fig. 3c).

Fig. 3: Sex-dimorphic pathways in bone-marrow neutrophils reveal differential regulation of chromatin-related pathways.

a,b, NetworkAnalyst predicted PPI networks for genes displaying significant bias in gene expression by RNA-seq toward female (a) or male (b) neutrophils. Network edges are based on IMEx/InnateDB data, a knowledgebase specifically geared for analyses of innate immune networks¹¹⁹. c,d, Top enriched gene sets from Reactome (c) and transcription factor (TF) targets (d) using GSEA for differential RNA expression. DSCAM, down syndrome cell adhesion molecule; MET, MET proto-oncogene, receptor tyrosine kinase; VEGF, vascular endothelial growth factor. e, Phenotype set enrichment analysis (PSEA) for differential metabolite regulation. f, Top LION functional enrichment analysis for differential lipid regulation. For c–f, the top ten most-upregulated and top ten most-downregulated gene sets (if that many) are plotted for readability. Full lists and statistics available in Supplementary Table 4. Also see Extended Data Fig. 7. Shown are pathways with FDR < 5%.

Next, we asked which pathways were regulated in a sex-dimorphic fashion in neutrophils (Fig. 3c, Extended Data Fig. 7a,b and Supplementary Table 4). Functional enrichment analysis revealed that significant female-biased pathways encompassed extracellular matrix (ECM) and cell surface-related pathways (Fig. 3c, Extended Data Fig. 7a,b and Supplementary Table 4a–d). Indeed, collagen-encoding genes Col1a1 and Col1a2 were expressed at higher levels in female neutrophils (Supplementary Table 1b). Collagen production, usually from fibroblasts, is a key event in fibrosis⁵³ and neutrophils play an important role in the development of fibrotic lesions⁵⁴. Alternatively, differential expression of cell-surface-related genes may also impact the ability of neutrophils to migrate through endothelial barriers⁵⁵. Notably, autophagy-related gene sets were female-biased in our dataset (Supplementary Table 4a,c,e). As mentioned above, autophagy control is critical for neutrophil biology^35,36, including neutrophil degranulation³⁷.

Chromatin- and cell-cycle-related pathways were overwhelmingly biased for higher expression in male neutrophils (Fig. 3c, Extended Data Fig. 7a,b and Supplementary Table 4a–c). Consistently, 21 histone-encoding genes showed significant male-biased expression (Supplementary Table 1b). The sex-dimorphic regulation of chromatin-related pathways suggests that sex may lead to differences in neutrophil chromatin metabolism.

We next investigated whether sex-dimorphic expression was correlated with predicted activity of TFs in neutrophils (Fig. 3d and Supplementary Table 4d). Consistent with pathway analysis, targets of cell cycle-related E2f7 and Mybl1^39,56 were more highly expressed in male neutrophils (Fig. 3d). Female neutrophils showed higher expression of Foxo4 targets (Fig. 3d). Notably, the Insulin/Igf1 signaling pathway, which regulates Foxo TF activity, has been linked to multiple sex-dimorphic phenotypes^57,58.

We next examined sex-dimorphism based on metabolomics and identified a number of pathways with significant sex-bias (Fig. 3e, Extended Data Fig. 7c and Supplementary Table 4f). Notably, metabolic pathways related to nucleotide metabolism were differentially regulated between sexes (Fig. 3e, Extended Data Fig. 7c and Supplementary Table 4f). Indeed, a number of sex-dimorphic metabolic peaks were predicted to represent nucleotide species (such as AMP, ATP, GMP; Supplementary Table 1d). Sex differences in nucleotide pools are consistent with transcriptomic observations of sex-dimorphism in (1) pathways directly related to nucleotide metabolism, which could lead to direct changes in nucleotide pools and (2) pathways responsive to nucleotide levels, such as signaling by Rho-GTPases, which can regulate neutrophil recruitment⁵⁹ (Supplementary Table 4a–c,e). Among nucleotides, the impact of adenine derivatives in neutrophils has been studied most extensively. ATP/ADP levels are regulated in neutrophils as a function of activation and ‘cellular age’⁶⁰. Adenine nucleotides from neutrophils exert anti-inflammatory effects⁶¹, increase endothelial barrier function, attenuate neutrophil adhesion to endothelial cells and modulate transendothelial migration⁶². Functional analysis of metabolomic data also suggested significant sex differences in amino-acid metabolism, such as tryptophan and arginine/proline metabolism (Fig. 3e and Extended Data Fig. 7c). Of note, arginine and tryptophan metabolism play important immunoregulatory roles⁶³; however, the final effect of differential metabolite levels from these pathways between sexes will need further investigation.

We then used the Lipid Ontology (LION) framework to analyze sex dimorphism in lipidomic profiles (Fig. 3f, Extended Data Fig. 1d and Supplementary Table 4g). Notably, male neutrophils were strongly enriched for triacylglycerols (TAGs) stored in lipid droplets, whereas female neutrophils were enriched for diacylglycerols (DAGs) (Fig. 3f, Supplementary Table 4g and Extended Data Fig. 1d). Consistently, genes associated with the Gene Ontology (GO) term ‘negative regulation of lipid storage’ were expressed at significantly higher levels in female neutrophils (FDR = 4.33 × 10⁻³; Supplementary Table 4c). Lipid droplets are crucial during neutrophil differentiation and as a source for inflammatory mediators^35,64. Adipose triglyceride lipase (Atgl), encoded by Pnpla2, regulates neutral lipid storage into lipid droplets in neutrophils⁶⁵. Of note, Pnpla2 is expressed at higher levels in female neutrophils (FDR = 0.01; Extended Data Fig. 7d and Supplementary Table 1b). Indeed, Atgl catalyzes the conversion of TAGs to DAGs, consistent with higher levels of DAGs in female neutrophils and of TAGs in male neutrophils (Fig. 3f, Extended Data Fig. 1d and Extended Data Fig. 7e). Lower levels of Pnpla2 (mimicking a ‘masculinized’ state) can lead to abnormal neutral lipid accumulation in neutrophils, increased chemotactic ability and reduced release of proinflammatory lipids⁶⁵. In addition, Atg7 is crucial for generation of free fatty acids (FFAs) from lipid droplets in neutrophils and Atg7-deficient neutrophils show increased lipid-droplet storage³⁵. Consistently, Atg7 is expressed at higher levels in female neutrophils (FDR = 0.03; Extended Data Fig. 7d and Supplementary Table 1b). Finally, Lpin1 is crucial for synthesis of DAG from phosphatidic acid and plays a role in the regulation of inflammation⁶⁶. Consistent with increased DAG levels, Lpin1 is also expressed at higher levels in female neutrophils (FDR = 0.03; Extended Data Fig. 7d and Supplementary Table 1b). Thus, both lipidomics and transcriptomics data suggest that male versus female neutrophils have increased neutral lipid storage (or decreased usage) (Extended Data Fig. 7e). Because of the functional importance of lipid droplets, these observations are likely to have deep functional consequences on neutrophil function.

Finally, we performed an integrated analysis using IMPaLA, combining information from transcriptomics, metabolomics and lipidomics (Supplementary Table 1h). We observed joint enrichment of male-biased molecules in pathways linked to cell cycle, DNA metabolism and chromatin-related pathways (Supplementary Table 4h), further supporting the notion that neutrophil chromatin architecture may be regulated in a sex-dimorphic manner. In contrast, female-biased pathways were enriched for lipid metabolism (Supplementary Table 4h). These joint observations provide an integrated confirmation of our observations for individual ‘omic’ layers. Together, our analyses suggest that major aspects of neutrophils are regulated in a sex-dimorphic fashion and may ultimately underlie sex differences in immune responses.

Neutrophil chromatin organization changes with sex and age

Our enrichment analyses revealed that chromatin-related pathways were significantly modulated by sex and organismal age in neutrophils, which is expected to result in profound changes in chromatin architecture. Neutrophils hold their chromatin in a compacted polylobular nuclear architecture, which earned them the name of ‘polymorphonuclear’ cells^28,29. More than just a gene expression regulatory layer, neutrophil chromatin directly participates in antimicrobial responses through NETs^30,31. Thus, chromatin-metabolism differences could lead to profound changes in neutrophil biology. To directly evaluate neutrophil chromatin, we utilized ATAC-seq⁶⁷, which takes advantage of adaptor-loaded Tn5 particles inserting into accessible chromatin, yielding information about local chromatin accessibility⁶⁷ and higher-order organization⁶⁸.

To evaluate chromatin landscapes, we isolated neutrophils from an independent cohort of young (4–5-month-old) and old (20–21-month-old) C57BL/6Nia female and male mice and performed ATAC-seq (Fig. 4a and Extended Data Fig. 8a–h). In contrast with transcriptomic, metabolomic and lipidomic data, MDS analysis on ATAC-seq showed that aging led to better sample separation than sex when examining local differences in chromatin accessibility (Extended Data Fig. 8d). In addition, there were substantially more regions with age-related versus sex-dimorphic differential accessibility (Extended Data Fig. 8e–h), consistent with numbers of differentially expressed genes with respect to age and sex (Extended Data Fig. 1a). Notably, enriched GO terms associated with differentially accessible regions with respect to organismal age (FDR < 5%) were mostly associated with regions with increased accessibility and encompassed terms related to immune response (Supplementary Table 3f). Significantly enriched GO terms for differential regions with respect to sex (FDR < 5%) were all female-biased and featured terms related to MAPK signaling and chromatin (Supplementary Table 4i).

Fig. 4: ATAC-seq analysis reveals age- and sex-related differences in the chromatin architecture of bone-marrow neutrophils.

a, Setup scheme for ATAC-seq experiment. b,c, Metaplot analysis of median nucleosome occupancy (as calculated by NucleoATAC) around the TSS of neutrophil-expressed genes, to analyze sex- (b) or age-related (c) differences in neutrophil nucleosome occupancy. Note that increased occupancy is observed in male compared to female neutrophils, as well as in old compared to young neutrophils. Significance of the difference between median occupancy profiles at TSSs assessed using a two-sided Kolmogorov–Smirnov goodness-of-fit test for b,c. d, Box plot of nucleosome occupancy as calculated by NucleoATAC for high-confidence nucleosomes across groups. Also see Extended Data Fig. 8i. e, Box plot of nucleosome fuzziness as calculated by NucleoATAC for high-confidence nucleosomes across experimental groups. f, Box plot of nucleosome inter-dyad distance as calculated by NucleoATAC across experimental groups. The horizontal red line in d–f shows the median value in neutrophils from young females for ease of comparison. Significance in nonparametric two-sided Wilcoxon rank-sum tests is reported for d–f. Grey P values represent differences between male versus female neutrophils in each age group; pink (blue) P values represent age-related differences in female (male) neutrophils. Also see Extended Data Fig. 8i,j. g, Workflow for in vitro NETosis inducibility experiment. SYTOX Green was used to stain extracellular DNA, a proxy for NETosing cells. h, Box plot of NETosis induction in naive neutrophils. Each point represents one animal, median of four technical replicates. Animals from four independent cohorts, n = 16–19 per group (variation due to animal death prior to experiment and technical failures of the flow cytometer on some samples; n = 17 for young males and females; n = 16 for old females; n = 19 for old males). i, Box plot of NETosis induction in neutrophils, primed with 10 ng ml⁻¹ mouse TNF-α for 10 min before NETosis induction. Each point represents one animal, median of four technical replicates. Animals from four independent cohorts, n = 17 per group. Significance in nonparametric two-sided Wilcoxon rank-sum tests are reported for h,i. Grey P values represent differences between male versus female neutrophils in each age group; pink (blue) P values represent age-related differences in female (male) neutrophils. For all box plots in d–f,h,i, the center line represents the sample median, the box limits consist of the 25th and 75th percentiles, the whiskers span 1.5× the interquartile range. For d–f, for readability and consistent with practices in genomics, outliers are not shown on the graph although they are always taken into account for the statistics.

As enrichment analyses suggested that overall aspects of chromatin are differentially regulated with organismal age and sex, we leveraged NucleoATAC⁶⁸ to analyze neutrophil chromatin architecture. Notably, nucleosome ‘v-plots’ generated revealed a region at ~145–150 bp that captured more reads in male versus female neutrophils, regardless of animal age, suggesting global differences in nucleosome packaging (Extended Data Fig. 8a). Next, we leveraged indicators linked to chromatin compaction and organization around nucleosomes: (1) nucleosome occupancy, which indicates how frequently a position is occupied by a nucleosome across cells, where higher occupancy is associated to tighter chromatin structure, (2) nucleosome fuzziness, which indicates how well positioned the nucleosome is, where fuzzier nucleosomes correspond to looser chromatin and (3) inter-dyad distance, which measure DNA length between neighboring nucleosomes, with increased distance associated to looser chromatin.

Occupancy metaplots around transcriptional start sites (TSSs) of expressed genes revealed that, regardless of age, male neutrophils showed higher median nucleosomal occupancy compared to female neutrophils (Fig. 4b). Similarly, regardless of sex, aging was associated to increased median nucleosomal occupancy (Fig. 4c). More generally, male neutrophils had a slight, but significant, increase in nucleosomal occupancy compared to female neutrophils and median nucleosomal occupancy also slightly increased with organismal age (Fig. 4d and Extended Data Fig. 8i). We observed decreased nucleosome fuzziness in male versus female neutrophils, as well as in old versus young neutrophils (Fig. 4e). Finally, we observed shorter inter-dyad distance in male versus female, as well as old versus young neutrophils (Fig. 4f). Together, small (but consistent) increases in occupancy, decreases in fuzziness and decreases in inter-dyad distance, are indicative of an overall more compacted chromatin architecture in male versus female and old versus young neutrophils. Finally, when analyzing ATAC-seq fragments at accessible regions (including subnucleosomal fragments), we observed lower median fragment length in male and old neutrophils (Extended Data Fig. 8j). While this observation may seem contradictory, overall shorter ATAC-seq fragments in the libraries may result from relatively ‘more accessible’ nucleosome-free regions in the context of overall tighter chromatin.

To independently evaluate chromatin compaction, we leveraged our RNA-seq to examine transcription of repeats (Extended Data Fig. 8k,l and Supplementary Table 1g,h). Eukaryotic genomes contain large proportions of repeats, including transposons, whose transcription is usually tightly repressed through compaction⁶⁹. Consistent with a more compact chromatin architecture, male neutrophils showed lower transcription of repeats versus females (350 female-biased elements versus 12 male-biased; Extended Data Fig. 8l and Supplementary Table 1h). Similarly, old neutrophils also showed lower transcription of repeats versus young neutrophils (115 elements downregulated versus 27 upregulated with age; Extended Data Fig. 8k and Supplementary Table 1g). Interestingly, chromatin decompaction is a limiting step of NETosis⁷⁰, suggesting that baseline differences in chromatin compaction in neutrophils of different ages and sex may impact the dynamics and timeline of NETosis. However, chromatin profiling alone cannot predict the direction in which NETosis could be affected.

Based on transcriptome and ATAC-seq results, we predicted that NETosis inducibility should be regulated as a function of sex and organismal age. We evaluated NETosis inducibility in neutrophils isolated from young versus old, female versus male mice (Fig. 4g–i). We used flow cytometry to quantify cells with extracellularized DNA (NETosed or NETosing cells), after a 2 h incubation in the presence of phorbol 12-myristate 13-acetate (PMA), a known activator of NETosis, compared to vehicle (Fig. 4g). We first evaluated NETosis induction on unstimulated neutrophils (Fig. 4h). NETosis inducibility was significantly increased with organismal age in female cells (P = 4.5 × 10⁻⁵), while higher animal-to-animal variability and no significant change were observed in male cells (Fig. 4h). Young male neutrophils showed higher NETosis capacity than young female neutrophils (P = 0.057), a trend reversed with organismal aging, with old female neutrophils showing highest NETosis inducibility (P = 0.029; Fig. 4h).

Increased NETosis in neutrophils from older female mice may seem contradictory with previous reports of decreased NETosis with organismal aging^13,14. To note, previous studies evaluated NETosis after in vitro TNF-α priming¹³ or in conditions known to mimic TNF-α priming^14,71 and did not evaluate sex as a biological variable. TNF-α is a proinflammatory cytokine and a potent inducer of NETosis¹³. To evaluate how TNF-α priming may impact NETosis across groups, we also performed NETosis experiments on cells pretreated with 10 ng ml⁻¹ TNF-α (Fig. 4i). Although TNF-α priming yielded noisier data, NETosis inducibility was blunted with organismal age in female cells stimulated with TNF-α (Fig. 4i). In contrast, TNF-α-primed male neutrophils trended toward decreased inducibility of NETosis with organismal age (P = 0.092; Fig. 4i). In TNF-α-primed conditions, young male versus female neutrophils no longer showed differences in NETosis inducibility, but old female neutrophils still showed the highest overall NETosis inducibility (Fig. 4i). Together, we find that NETosis is differentially modulated as a function of organismal age and sex, consistent with observed differences in chromatin architecture.

Machine learning predicts features of neutrophils regulation

To understand whether age-related and sex-dimorphic gene expression can be predicted from genomic sequence, chromatin and/or regulation by TFs, we took advantage of machine learning (Figs. 5a–c and 6a–c and Extended Data Fig. 9a–f). Machine learning can help (1) to determine whether features contain information predictive of a state of interest (such as sex-dimorphic or age-regulated expression) and (2) to identify predictive features of such states. We used seven algorithms to learn models discriminating between (2) up- and downregulated expression with age and (3) male- or female-biased gene expression (conditional inference trees, linear discriminant analysis, neural networks, logistic regression, random forests (RFs), support vector machines and gradient-boosting machines (GBMs)). The models were trained with features related to (1) genomic sequence (such as GC richness), (2) chromatin accessibility (such as changes in ATAC-seq promoter accessibility) and (3) TF target genes (the same as used with GSEA).

Fig. 5: Machine-learning analysis reveals that age-regulated gene expression can be predicted by genomic and epigenomic features.

a, Scheme of our machine-learning pipeline for upregulated versus downregulated genes in neutrophils with aging. b, Balanced classification accuracy over different machine-learning model algorithms. The accuracy of the models is measured using held-out testing data. ‘Random’ accuracy illustrates the accuracy of a meaningless model (~50%) and ‘perfect’ that of a model with no mistakes (100%). All tests were more accurate than random (see other measures of model performance in Extended Data Fig. 9a–c and Supplementary Table 5a). CTree, conditional inference tree; LDA, linear discriminant analysis; SVM, support vector machine; NNET, neural network; LogReg, logistic regression. c, Top 20 most-important features from RF models and GBM models. A rank product approach was used to determine overall top predictive features from both models. High values for variable importance indicate most influential predictors.

Fig. 6: Machine-learning analysis reveals that sex-dimorphic gene expression can be predicted by genomic and epigenomic features.

a, Scheme of our machine-learning pipeline for male versus female-biased gene expression in neutrophils. b, Balanced classification accuracy over different machine-learning model algorithms. The accuracy of the models is measured using held-out testing data. ‘Random’ accuracy illustrates the accuracy of a meaningless model (~50%) and ‘perfect’ that of a model with no mistakes (100%). All tests were more accurate than random (see other measures of model performance in Extended Data Fig. 9d–f and Supplementary Table 5c). c, Top 20 most-important features from RF models and GBM models. A rank product approach was used to determine overall top predictive features from both models. High values for variable importance indicate most influential predictors.

We first evaluated our ability to discriminate between genes upregulated versus downregulated during aging in neutrophils (Fig. 5a). Our models assigned genes to the correct transcriptional change with age >64% of the time on testing data (Fig. 5b and Supplementary Table 5a), significantly above random (50%). High accuracy suggests that clear biological features differentiate between genes that tend to be upregulated versus downregulated with aging. To assess which features were most predictive of age-related changes, we examined predictor contribution from RF and GBM models, which provide native evaluation of feature importance (Fig. 5c and Supplementary Table 5b). Consistent with the tight link between chromatin and gene expression⁷², the average promoter ATAC-seq signal (describing how ‘open’ the promoter is), as well as the fold difference in ATAC-seq signal between young and old neutrophils (age-related changes in openness) were among the top predictors for age-related changes (Fig. 5c and Supplementary Table 5b). Key predictors also included CG/CpG content in the promoter and gene sequences (Fig. 5c), which is consistent with our previous work identifying promoter CpG content as a top predictor of age-related gene expression in mouse tissues¹⁸. Cytosines in CpG configuration are the primary targets of DNA methylation⁷³. CpG DNA methylation (DNAme), catalyzed by DNA-methyltransferases (writers), can be removed by TET family proteins (erasers) and is interpreted through recognition by methyl-CpG binding proteins (readers)⁷³. Consistent with the notion that the predictive power of CpG content is related to changes in CpG methylation regulation, RNA-seq revealed significant downregulation with organismal aging of DNA-methyltransferase encoding Dnmt1, demethylation-catalyzing TET encoding Tet1 and Tet2, as well as methyl-CpG binding domain proteins encoding Mbd4 and Mbd5 (FDR < 5%; Supplementary Table 1a). The presence of putative age-related changes in DNAme is consistent with the notion that DNAme patterns can serve as ‘aging clocks’⁷⁴.

Key predictors of age-related expression changes also included whether a gene was a target of specific TFs: Stat5a, Mtf2, Sirt6 and Foxo1 (Fig. 5c and Supplementary Table 5b). Although causality cannot be inferred from machine learning, predictors provide a list of candidate drivers/mediators of programs related to organismal aging. Only a subset of top TF predictors were themselves differentially regulated with aging (such as Stat5a, Mtf2 and Nod2; Extended Data Fig. 9g and Supplementary Table 1a). Thus, it is possible that the activity of TF predictors without RNA-level regulation may occur through post-translational regulation. Notably, Stat5a mediates the effects of granulocyte–macrophage colony-stimulating factor (GM–CSF) on granulocyte differentiation⁷⁵ and is crucial for mature neutrophil survival⁷⁶. Although the role of Mtf2 in neutrophils has not been studied, Mtf2 interacts with Jarid2 (another top predictor) to promote repressive H3K27 methylation⁷⁷, which mediates recruitment to unmethylated CpGs⁷⁸. Nod2 is a pattern-recognition receptor that recognizes muramyl dipeptide, the minimal common motif of bacterial peptidoglycans⁷⁹. Although Nod2 is not a TF, it has profound transcriptional effects on gene expression, with target genes including cytokines and chemokines^80,81. Indeed, Nod2 is an important contributor of neutrophil-mediated innate immunity⁸¹. Sirt6 is a histone deacetylase tightly linked to the regulation of mammalian aging and longevity^82,83. Sirt6 can limit inflammation by deacetylating the promoters of nuclear factor (NF)-κB targets⁸⁴ and reducing cytokine production⁸⁵. Consistently, Sirt6 knockout mice show liver inflammation with neutrophil infiltration⁸⁶. Thus, changes in activity for top predictive TFs may drive aspects of age-related transcriptional remodeling.

We next asked whether sex-dimorphic gene expression could be predicted from genomic information, local chromatin states and/or TF regulation (Fig. 6a). As an additional feature, we also included the chromosomal location of a gene (sex chromosomes versus autosomes). Our models assigned genes to the correct sex-based expression bias >70% of the time on held-out testing data (Fig. 6b and Supplementary Table 5c), consistent with the notion that sex-biased genes share a common regulatory signature. To understand which features were most predictive of male- or female-biased expression, we again examined predictor contribution from RF and GBM models (Fig. 6c and Supplementary Table 5d). Notably, being located on a sex chromosome (X or Y) was a poor predictor of sex-dimorphic gene expression (166th combined importance rank; Supplementary Table 5d), consistent with the fact that genetic context was not a crucial driver (Extended Data Fig. 1a). Similar to our age-related models, the average promoter ATAC-seq signal and the fold change between females versus males in ATAC-seq signal were also top predictors for sex-biased gene expression (Fig. 6c and Supplementary Table 5d). The GC/CpG content of promoter and genes sequences also showed strong predictive power, consistent with the sex-dimorphic expression of DNAme regulators (male-biased expression for Dnmt1 and Tet3; Supplementary Table 1b).

Finally, we asked which TFs were most predictive of sex-dimorphic gene expression in neutrophils (Fig. 6c and Supplementary Table 5d). Being a known target of Foxm1, Sirt6, Mtf2 or Mybl2 was highly predictive of male versus female neutrophil-biased gene expression, suggesting that these may drive aspects of the sex-dimorphic neutrophil biology. Similar to above, only a small subset of top predictors from our models were themselves significantly regulated at the transcriptional level (for example Foxm1, Irf8; Extended Data Fig. 9h and Supplementary Table 1b). Although its role hasn’t been explored in neutrophils, Foxm1 regulates the expression of cell cycle-related genes⁸⁷ and is highly expressed in late committed neutrophil precursors⁴⁰. The transcription factor Mybl2 also regulates cell cycle-related genes⁸⁸ and can modulate differentiation and cell fate decision of myeloid progenitors⁸⁹. Thus, top predictors of sex-dimorphic neutrophil gene expression Foxm1 and Mybl2 might be linked to observed sex-dimorphic expression of cell cycle and chromatin-related genes (Fig. 3 and Extended Data Fig. 7). Being a target of Sirt6 was also predictive of sex-dimorphic gene expression, in line with sex-dimorphic phenotypes of Sirt6 knockout and overexpression mice^82,83. In addition, consistent with sex dimorphism in lipid storage (Fig. 3f), Sirt6 is tightly linked to lipid metabolism⁹⁰ and can limit lipid-droplet accumulation in foam cells⁹¹. In addition, Sirt6 activity can be directly regulated by FFAs (such as myristic, oleic and linoleic acid)⁹². Consistently, female neutrophils have higher levels of FFAs, including FFA(14:0) and FFA(18:2) (myristic and linoleic acid) and may thus have higher basal Sirt6 activity levels (Supplementary Table 1e,f). Thus, our machine-learning analysis reveals candidate regulators that may drive neutrophil sex dimorphism and will deserve further investigation.

Primary granule differences in female versus male neutrophils

Notably, when analyzing the top genes with the largest fold difference between male versus female neutrophils, we noticed that genes encoding primary neutrophil granule components showed top male-biased expression (such as Elane, Mpo, Prtn3 and Ctsg; Supplementary Table 1b), suggesting potential sex differences in granule loading. Thus, we investigated potential differences in granule biogenesis (Fig. 7a,b and Extended Data Fig. 10a–d). Neutrophils store a number of proteases and antimicrobial proteins in (1) primary/azurophilic granules, (2) secondary/specific granules and (3) tertiary/gelatinase granules. To unbiasedly assess whether neutrophil granule biology is sex dimorphic, we compiled a list of granule components from the literature (Supplementary Table 6). Consistent with our first observation, GSEA revealed a robust trend for increased RNA expression of primary granule-related genes in male versus female neutrophils, regardless of age (Fig. 7a,b). In contrast, no significant enrichments were observed for genes encoding components of secondary or tertiary granules (FDR > 5%; Extended Data Fig. 10a–d). Notably, we confirmed that protein levels of neutrophil elastase Elane were higher in young males versus young females by western blot (P = 3.9 × 10⁻³ in two-sided Wilcoxon’s test between young females and male samples), although this sex difference disappeared in older animals (Extended Data Fig. 5d,e).

Fig. 7: Male neutrophils express higher levels of primary granule genes, correlating with increased serum ELANE levels in control and septic mice.

a, Heat map of normalized gene expression for primary (azurophilic) granule-related gene expression in our RNA-seq dataset. Also see Extended Data Fig. 10. b, GSEA of primary (azurophilic) granule-related gene expression reveals biased expression to male neutrophils. c, Setup scheme for serum ELANE measurement in control and sepsis-like mice. d,e, Analysis of ELANE protein levels in mouse serum by ELISA. Data are derived from three independent cohorts of young male versus female mice (n = 5 per sex and time point for 1, 3 and 6 h PBS injection; n = 8 per sex and time point for 1 and 3 h LPS and n = 9 per sex and time point for 6 h LPS injection). For simplicity, all PBS-injected animals are reported as ‘0 h of LPS exposure’ in d and are replotted with time-based color coding in e. Significance in nonparametric two-sided Wilcoxon rank-sum test. For box plots in d,e, the center line represents the sample median, the box limits consist of the 25th and 75th percentiles and the whiskers span 1.5× the interquartile range.

Neutrophil degranulation is a key aspect of the response of neutrophils to pathogens⁹³ and mice without a functional copy of Elane have increased mortality upon sepsis⁹⁴. However, excessive levels of circulating neutrophil elastase can be deleterious, amplifying septic shock⁹⁵. Based on RNA-seq, we hypothesized that higher Elane expression in neutrophils could lead to increased Elane circulating levels, both in basal conditions and upon an immune challenge. To test this hypothesis, we obtained female and male young adult C57BL/6J mice and injected them with sterile PBS or lipopolysaccharide (LPS; major component of Gram-negative bacteria cell wall), for 1–6 h, to mimic a septic state⁹⁶ (Fig. 7c). Notably, males showed significantly higher levels of serum neutrophil elastase at 3 and 6 h after LPS exposure (Fig. 7d) and in the PBS-injected controls (Fig. 7d,e). Thus, young males had overall higher serum protein levels of Elane.

Finally, we then took advantage of a published proteomics dataset of human blood neutrophils from 68 healthy donors⁹⁷. As the biological sex of human donors was not specified, we used the reported protein levels of DDX3Y, a protein encoded on the Y chromosome, as a proxy for the likelihood that the sample was derived from a male donor (higher levels of DDX3Y associated with males). Consistent with our mouse RNA-seq and serum ELISA, we found significant correlation between DDX3Y levels (likelihood of the sample coming from a male donor) and expression of primary granule proteins ELANE, MPO and CTSG, although not of PRTN3 (Extended Data Fig. 10e–g). Thus, our results suggest that neutrophils derived from human males have higher protein expression of key primary granule components. These observations suggest that, at least for these components, transcriptomic trends are predictive of protein-level trends and that the male bias in primary granule components observed in our mouse data may be conserved in humans.

Discussion

A resource for the study of neutrophils across sex and aging

We have generated transcriptomic, metabolomic, lipidomic and epigenomic datasets using bone-marrow neutrophils from young and old, female and male animals. This dataset is a large multi-omic dataset for the study of neutrophils and includes both sexes and organismal aging, rather than focusing only on one sex or age group. Using this resource, we observed that (1) neutrophils are extremely sex dimorphic and that (2) organismal age leads to large scale transcriptomic and epigenomic remodeling of neutrophils. Of note, when analyzed separately, males had over tenfold more significantly regulated genes with age than females, suggesting that the molecular rate of neutrophil aging differs between sexes. Observed changes may likely be downstream of signals that drive aging itself, but these changes will impact immune responses in aged animals, participating in immune dysfunction.

We then predicted increased release of neutrophil elastase in males in control and sepsis-like conditions. Excess release of neutrophil elastase is known to exacerbate inflammation and cause tissue damage⁹⁸. Thus, this resource should help open new avenues of research and identify candidate mediators that underlie sex differences in lifelong immunity. Future studies should investigate how the differences in bone marrow neutrophils are maintained, erased or amplified in circulating blood neutrophils.

Sex and aging influence neutrophil chromatin metabolism

We observed that regulation of chromatin metabolism is a hallmark of neutrophil aging and a key aspect with sex-dimorphic regulation. Far from being a mere regulatory layer for gene expression, chromatin organization has profound implications on NETosis⁸. Our analyses suggest that female neutrophils have an overall ‘looser’ chromatin, associated to the transcriptional upregulation of a number of transposable elements (TEs). Conversely (and regardless of sex), neutrophils from old individuals show increased chromatin condensation accompanied by reduced TE transcription. This contrasts with observation of age-related TE derepression in other contexts¹⁹ and may reflect the unique neutrophil biology. Thus, it will be important to elucidate the mechanism driving differences in NETosis between sexes and throughout lifespan, with implications in aberrant chronic age-related inflammation¹⁰.

Machine learning as a powerful candidate-identification tool

By using machine learning, we showed that both age regulation and sex dimorphism in gene expression can be predicted accurately. The high accuracy of our models supports the notion of coordinated differences between (1) genes induced vs. repressed during aging and (2) genes expressed in a sex-dimorphic manner. Among important predictors of transcriptional states, we identified regulators whose activity change during organismal aging may underlie ‘omic’ profile remodeling throughout life (such as Foxo1 and Sirt6). We also identified putative mediators of sex-biased gene expression of neutrophils (such as Foxm1 and Mybl2). Although machine learning does not provide information about causality, predictors represent candidate mediators that could help shape neutrophil landscapes with aging or as a function of sex. Thus, it will be important to elucidate the potential role of these predictors to understand the emergence of age-related and sex-dimorphic phenotypes.

Neutrophils as mediators for sex differences in immune aging

Accumulating evidence has shown that, even outside of reproduction, mammalian biology is extremely sex dimorphic⁹⁹. This is especially relevant to aging, with large sex differences in baseline lifespan and in response to pro-longevity interventions⁵⁸. Although this has not been broadly investigated with aging, pioneering studies have started to compare male and female cells, revealing profound sex dimorphism in gene regulation^21,22. Consistently, accumulating evidence has shown that immune responses differ between sexes^100,101.

Sex dimorphism in immunity is exemplified in the current health crisis, with men representing the majority of severe cases and deaths from COVID-19 (ref. ¹⁰²). Notably, NETs may drive aspects of inflammatory disease, lead to tissue damage^8,30 and contribute to severe COVID-19 (ref. ¹⁰³). In addition, neutrophil dysfunction has been linked to the pathogenesis of a number of chronic diseases (such as macular degeneration¹⁰⁴, stroke¹⁰⁵, obstructive pulmonary disease¹⁰⁶, atherosclerosis¹⁰⁷ and cancer¹⁰⁸). Thus, sex differences in neutrophils could constitute targets to optimize immune responses throughout life and help tailor therapeutics to men and women.

Methods

Mouse husbandry

All animals were treated and housed in accordance with the Guide for Care and Use of Laboratory Animals. All experimental procedures were approved by the University of Southern California’s Institutional Animal Care and Use Committee (IACUC) and are in accordance with institutional and national guidelines. Samples were derived from animals on approved IACUC protocol nos. 20770, 20804 and 21004.

For ‘omics’ analyses, male and female C57BL/6Nia mice (4–5- and 20–24-month-old animals) were obtained from the National Institute on Aging (NIA) colony at Charles Rivers. For the sepsis assays in young animals, male and female C57BL/6J mice (3–4-month-old animals) were obtained from Jackson Laboratories. Both Charles Rivers and Jackson Laboratories have specific-pathogen-free facilities. Animals were acclimated at the specific-pathogen-free animal facility at University of Southern California (USC) for 2–4 weeks before any processing. The facility is on a 12-h light/dark cycle and animal housing rooms are maintained at 72 °F and 30% humidity. All animals were euthanized between 8–11 am for ‘omics’ experiments, flow cytometry and NETosis assays. For the sepsis model experiments, animals were injected in the morning and euthanized between 4–6 pm. In all cases, animals were euthanized using a ‘snaking order’ across all groups to minimize batch-processing confounds due to circadian processes (including potential confounds due to neutrophil cellular ‘age’). All animals were euthanized by CO₂ asphyxiation followed by cervical dislocation.

Isolation of primary neutrophils from mouse bone marrow

The long bones of each mouse were collected and kept on ice in D-PBS (Corning) supplemented with 1% penicillin/streptomycin (Corning) until further processing. Muscle tissue was removed from the bones and bone marrow from cleaned bones was collected into clean tubes¹⁰⁹. Red blood cells from the marrow were removed using Red Blood Cell Lysis buffer (Miltenyi Biotec, no. 130-094-183), according to the manufacturer’s instructions, albeit with no vortexing step to avoid unscheduled neutrophil activation. The suspension was filtered on 70-μm mesh filters (Miltenyi Biotec, no. 130-110-916) to retain only single cells for downstream processing. Neutrophils were isolated from other bone-marrow cells using MACS (Miltenyi Biotec, no. 130-097-658). Viability and yield were assessed using trypan blue exclusion and an automated COUNTESS cell counter (Thermo Fisher Scientific). Purified cells were pelleted at 300g and snap-frozen in liquid nitrogen until processing for RNA, lipid or metabolite isolation.

Neutrophil purity estimates by flow cytometry

We used aging C57BL/6Nia mice from two independent cohorts to estimate potential differences in neutrophil purity with respect to sex and age. After an Fc-blocking step (Miltenyi Biotec, no. 130-092-575), MACS-purified neutrophils were then stained using APC-Ly6G (Invitrogen, no. 17-9668-80) and Vioblue-Cd11b (Miltenyi Biotec, no. 130-113-238) at a 1:50 dilution according to the manufacturer’s instructions. Stained cells were then analyzed by flow cytometry on a MACS Quant Analyzer 10 (Miltenyi Biotec). Flow cytometry results were processed using the FlowLogic v.7 software. Purity of MACS-isolated bone-marrow neutrophils in young and aged male and female animals are reported in Extended Data Fig. 3. Raw cytometry data were deposited on Figshare (https://doi.org/10.6084/m9.figshare.14043932.v1).

Purified neutrophil heterogeneity estimate by flow cytometry

We used flow cytometry analysis to evaluate the heterogeneity/maturity of our isolated neutrophils from two independent cohorts of aging C57BL/6Nia mice. We used a panel of antibodies designed to evaluate population heterogeneity of cells purified by MACS based on defined subpopulations of neutrophils in the literature^110,111,112 (see below). Specifically, among live cell singlets, pro-neutrophils were defined as c-Kit⁺Ly6g^-Cd81⁺, with subpopulations of pro-neutrophil 1 defined as c-Kit⁺Ly6G⁻Cd81⁺CD11b⁻CD106⁻ and pro-neutrophil 2 as c-Kit⁺Ly6G⁻Cd81⁺CD11b⁺CD106⁺. Pre-neutrophils were defined as Cd11b⁺Ly6G⁺c-Kit⁺Cxcr4⁺ cells. Immature neutrophils were defined as CD11b⁺Ly6G⁺c-Kit⁻Cxcr4⁻Cxcr2⁻ cells and mature neutrophils as CD11b⁺Ly6G⁺c-Kit⁻Cxcr4⁻Cxcr2⁺. Fresh versus ‘aged’ subsets of mature neutrophils were defined respectively as CD11b⁺Ly6G⁺c-Kit⁻Cxcr4⁻Cxcr2⁺CD62L⁺ versus CD11b⁺Ly6G⁺c-Kit⁻Cxcr4⁻Cxcr2⁺CD62L⁻ cells.

Before analyzing compositional heterogeneity on MACS-purified neutrophils, to account for spill-over from different lasers, compensation was performed using appropriate compensation beads (Miltenyi Biotec, no. 130-104-693 for Miltenyi Biotec antibodies; Thermo Fisher Scientific, no. 01-3333-42 for others). The compensation file was saved and loaded before analyzing neutrophils. Compensation was rerun anytime at least one fresh vial of antibody had to be used for the staining experiment.

Antibodies used for this heterogeneity panel were:

• CD184 (CXCR4)–PE-Vio770 (clone REA107; Miltenyi Biotec, 130-102-914) at 1:10 dilution

• CD81–PE (clone EAT2; Miltenyi Biotec, 130-102-632) at 1:10 dilution

• CD11b–Vioblue (clone REA592; Miltenyi Biotec, 130-113-810) at 1:50 dilution

• Ly6G–PerCP–Vio700 (clone REA526; Miltenyi Biotec, 130-117-500) at 1:50 dilution

• CD182 (CXCR2)–APC–Vio770 (clone REA942; Miltenyi Biotec, 130-115-637) at 1:50 dilution

• CD62L–FITC (clone REA828; Miltenyi Biotec, 130-112-835) at 1:50 dilution

• CD106 (VCAM-1)–APC (clone REA971; Miltenyi Biotec, 130-116-324) at 1:50 dilution

• CD117 (c-kit)–Brilliant Violet 510 (clone ACK2; BioLegend, 135119) at 1:20 dilution

After isolation by MACS, 2.5 × 10⁵ cells were washed once by adding 1 ml of PBS/EDTA 0.1% BSA (MACS resuspension buffer) in a 5-ml polystyrene round-bottom tube (Falcon, no. 352054), then centrifuged 300g for 10 min at 4 °C. After an Fc-blocking step (Miltenyi Biotec, no. 130-092-575), MACS-purified neutrophils were then stained with antibodies for 20 min at 4 °C and excess antibody was washed away using resuspension buffer. Stained cells were then analyzed by flow cytometry on a MACS Quant Analyzer 10 (Miltenyi Biotec) and flow cytometry data were analyzed using FlowLogic v.7. Gating was determined using fluorescence-minus-one controls for each color used in the experiment to ensure that positive populations were solely associated with the antibody for that specific marker (Extended Data Fig. 4a). Due to the low amount of c-kit⁺ events in MACS-purified neutrophils, we also ran whole bone marrow (pooled by group in one cohort) with the same scheme to confirm that positive labeling by this antibody occurs normally on bone-marrow cells before purification (Extended Data Fig. 4b). Both sets of raw cytometry data, including fluorescence-minus-one controls and compensation files, were deposited to Figshare (https://doi.org/10.6084/m9.figshare.14043929.v1, https://doi.org/10.6084/m9.figshare.14043938.v1).

Estimate of neutrophil proportion in bone marrow in young female and male mice

We used aging C57BL/6Nia mice from two independent cohorts to estimate potential differences in neutrophil proportions within nucleated bone-marrow cells with respect to sex and age. An aliquot of bone-marrow cell suspension was obtained after bone-marrow extraction, red blood cell lysis and filtration on 70-µm mesh filters (see above). Cell composition analysis was obtained using the Hemavet 950FS at the USC Leonard Davis School of Gerontology mouse phenotyping core. Percentage of cells was used instead of absolute cell numbers to avoid confounding results due to animal size.

Evaluation of cell-cycle patterns of MACS-purified neutrophils by flow cytometry

We used aging C57BL/6Nia mice from two independent cohorts to estimate cell cycle proportions of MACS-purified neutrophils with respect to sex and age. We utilized a standard method using DNA content, as measured by propidium iodide staining of fixed cells, to estimate the phases of the cell cycle (G0/G1, S and G2/M). For each biological sample, 5 × 10⁵ cells were aliquoted and fixed with 70% ethanol at −20 °C overnight. The next day, cells were washed with DPBS (Corning) and nuclear DNA was stained for 30 min at room temperature in labeling buffer (50 µg ml⁻¹ propidium iodide (Alfa Aesar), 100 µg ml⁻¹ RNaseA (Invitrogen), 0.05% Triton X-100 (Sigma) in DPBS (Corning)). Stained cells were then analyzed by flow cytometry on a MACS Quant Analyzer 10 (Miltenyi Biotec) and flow cytometry data were analyzed using FlowLogic v.7. Raw cytometry data were deposited to Figshare (https://doi.org/10.6084/m9.figshare.14043926.v1).

RNA purification and RNA-seq library preparation

For RNA isolation, frozen cell pellets were resuspended in 1 ml of Trizol reagent (Thermo Fisher Scientific) and total RNA was purified following the manufacturer’s instructions. RNA quality was assessed using the Agilent Bioanalyzer platform at the USC Genome Core using the RNA integrity number. Then, 500 ng of total RNA was subjected to ribosomal-RNA depletion using the NEBNext rRNA Depletion kit (New England Biolabs), according to the manufacturer’s protocol. Strand specific RNA-seq libraries were then constructed using the SMARTer Stranded RNA-seq kit (Clontech), according to the manufacturer’s protocol. Libraries were quality controlled on the Agilent Bioanalyzer 2100 platform at the USC Genome Core before multiplexing the libraries for sequencing. Paired-end 75-bp reads were generated on the Illumina NextSeq500 platform at the SC2 Core at CHLA (original cohort) or paired-end 150-bp reads were generated on the Illumina HiSeq-Xten platform at the Novogene Corporation (replicate cohort).

RNA-seq analysis pipeline

To best mimic the first RNA-seq cohort, paired-end 150-bp reads from the replicate cohort of neutrophil RNA-seq were hard-trimmed to 75 bp using Fastx Trimmer. Both sets of paired-end reads were processed using Trimgalore v.0.4.4 (github.com/FelixKrueger/TrimGalore) (1) to retain only high-quality bases with phred score >15 and (2) to eliminate biases due to priming by hard clipping the first six bases of each read. Only pairs with both reads retaining a length of >45 bp after trimming were retained for further processing. Trimmed reads were mapped to the mm10 genome reference using STAR v.2.5.0a¹¹³. Read counts were assigned to genes from the UCSC mm10 reference using subread v.1.5.3 (ref. ¹¹⁴) and were imported into R to perform differential gene expression analysis. Based on general RNA-seq processing guidelines, only genes with mapped reads in at least 6 out of 16 RNA-seq libraries were considered to be expressed and retained for downstream analysis. Due to high sample-to-sample variability, we used surrogate variable analysis to estimate and correct for unwanted experimental noise¹¹⁵. R package ‘sva’ v.3.34 was used to estimate surrogate variables and the removeBatchEffect function from ‘limma’ v.3.42.2 was used to regress out the effects of surrogate variables and RNA integrity differences (RNA integrity number scores) from raw read counts. The ‘DESeq2’ R package (DESeq2 v.1.26.0) was used for further processing of the RNA-seq data in R¹¹⁶. Notably, sex was encoded as a categorical variable (female versus male) and age was encoded as a continuous numerical variable. Genes with FDR < 5% were considered statistically significant and are reported in Supplementary Table 1a,b.

Dimensionality reduction for exploratory data analysis

To perform MDS analysis¹¹⁷, we used a distance metric between samples based on the Spearman’s rank correlation value (1-Rho), which was then provided to the core R command ‘cmdscale’. Dimensionality reduction was applied to pre-normalized data, as described in relevant sections.

Replicate and public neutrophil RNA-seq comparison

To assess the robustness of transcriptional changes linked to biological sex and/or organismal age, we used (1) a small replicate RNA-seq cohort of primary bone-marrow neutrophils with all four conditions, (2) a dataset from mouse spleen neutrophils from young adult female versus males from the ImmGen consortium (6-week-old animals; GEO datasets GSE124829)²² and (3) a dataset from human blood neutrophils (ages unknown; GEO Datasets GSE145231)²⁴ (Extended Data Fig. 2a,b). Reads were mapped to mm10 and hg38, respectively using STAR and reads were summarized to genes using subread as before. DESeq2-normalized fold changes were then used to estimate differential gene expression as a function of age and of sex. For mouse datasets, the comparison directly assessed the DESeq2-normalized fold changes for genes with FDR < 5% in the original cohort RNA-seq. For the human dataset, orthologs of significant mouse genes in humans were identified using the R ‘biomaRt’ v.2.42.1 package and DESeq2-normalized fold changes were then compared.

Putative protein–protein interaction network analysis

Genes with significant regulation according to our DEseq2 analysis (FDR < 5%) were used as input for network analysis with NetworkAnalyst v.3.0 (ref. ¹¹⁸). For age-related gene regulation, significant genes and DEseq2 calculated log₂(fold change) were used as input for a single network analysis. For sex-dimorphic gene regulation, female- and male-biased gene lists were used as separate inputs for the analysis. To avoid network clusters due to sex-chromosome-encoded genes, only significant autosomal genes were included in the sex-dimorphism gene expression networks. In both cases, NetworkAnalyst was set up to use PPI information derived from IMEx/InnateDB data, a knowledgebase specifically geared for analyses of innate immune networks¹¹⁹, to construct putative PPI networks. In each case, the largest subnetwork determined by NetworkAnalyst was used in figures and analyses.

Functional enrichment analysis (transcriptomics)

We used the GSEA paradigm¹²⁰ through the ‘phenotest’ v.1.34.0 R package. GO term annotations were obtained from ENSEMBL through Biomart (Ensembl 99) and gene-term association with only author statement support (GO evidence codes ‘TAS’ and ‘NAS’) or unclear evidence (GO evidence code ‘ND’) were filtered out. Other annotations were obtained from the Molecular Signature Database v.7.0 (Reactome and KEGG)^120,121 and the Harmonizome database (ChEA, JASPAR, ENCODE, GEO TF targets)¹²². To improve target coverage, we summarized putative TF targets, including FOXO TF targets compiled in our previous study¹⁸, to have a reduced, unique and nonredundant list of TF targets summarized from all these sources. The DEseq2 t-statistic was used to generate the ranked list of genes for functional enrichment analysis, both for the sex and aging effects. For ease of reading, only the top ten most-significant pathways with negative NES and top ten most-significant pathways with positive NES are shown in figures if more than that passed the FDR < 5% significance threshold. All significant gene sets are reported in Supplementary Table 3 (aging) and Supplementary Table 4 (sex). In parallel, functional enrichment was also independently assessed using the ingenuity pathway analysis portal, using genes with a significant male or female bias (FDR < 5%).

Western blot analysis of Padi4 and Elane protein levels

To prepare neutrophil lysates, ~5 million cells were collected from aging C57BL/6Nia mice from two independent cohorts (n = 9 animals per group). Flash-frozen cell pellets were resuspended in 200 µl ice-cold lysis buffer (10 mM Tris-HCl (pH 8.0), 1% SDS, 1× Halt Protease and Phosphatase Inhibitor Cocktail (Life Technologies, no. 78442)) and incubated on ice for 30 min. Samples were sonicated (Thermo Fisher Scientific, no. FB120; 60% power, 30 s ON/120 s OFF, four cycles) and centrifuged at 16,000g for 20 min at 4 °C. Supernatant was collected and boiled in 4× Laemmli Sample Buffer (Bio-Rad, no. 1610747). Proteins were separated on 10% SDS–PAGE gels and wet-transferred onto PVDF membranes (GE Healthcare, no. 10600021). To assess loading and transfer efficiency, membranes were stained using the Ponceau S solution (Sigma, P7170) for 2 min at room temperature. Membranes were blocked using 5% milk (Carnation) in 1× TBST buffer (TBS (Alfa Aesar J62938), supplemented with 0.05% Tween-20 (Bio-Rad no. 161-0781)) for 1 h at room temperature. Membranes were incubated with primary antibodies diluted in 5% milk in TBST (anti-β-actin: Cell Signaling D6A8 at 1:5,000 dilution; anti-ELANE: Abcam ab68672 at 1:1,000 dilution; and anti-PADI4: Abcam, ab96758 at 1:3,000 dilution) for 16 h at 4 °C. After three 5-min washes using 1× TBST buffer, membranes were incubated with an HRP-conjugated secondary antibody (goat anti-rabbit IgG H&L, Abcam, ab205718 at 1:10,000 dilution) for 1 h at room temperature. For visualization, membranes were treated with chemiluminescence solution according to manufacturer’s manual (Bio-Rad, no. 1705062 and Thermo Fisher Scientific, no. 34580) and images were taken using Azure Biosystems c280.

For band intensity quantification, ImageJ (v.1.53) was used. Grayscale converted images were imported to ImageJ and intensity was measured from all bands. Background intensity was subtracted from each measurement. Anti-β-actin-relative values from each membrane were normalized by the median value of the specific membrane to mitigate membrane-specific variations.

All raw original images and cropped equivalents, as well as ImageJ quantification, have been made available on Figshare (https://doi.org/10.6084/m9.figshare.14154665.v1) and as a source data file for Extended Data Fig. 5.

Chemicals for LC–MS

LC–MS-grade solvents and mobile-phase modifiers were obtained from Thermo Fisher Scientific (water, acetonitrile, methanol, methyl tert-butyl ether (MTBE)) and Sigma-Aldrich (ammonium acetate).

Metabolite and lipid extraction from neutrophils

Metabolites and lipids were extracted from neutrophil cell pellets and analyzed in a randomized order. Extraction was performed using a biphasic separation protocol with ice-cold methanol, MTBE and water¹²³. Briefly, 300 μl of methanol spiked-in with 54 deuterated internal standards provided with the Lipidyzer platform (SCIEX, cat. no. 5040156, LPISTDKIT-101) was added to the cell pellet, samples were vigorously vortexed for 20 s and sonicated in a water bath three times for 30 s on ice. Lipids were solubilized by adding 1,000 μl of MTBE and incubated under agitation for 1 h at 4 °C. After addition of 250 μl of ice-cold water, the samples were vortexed for 1 min and centrifuged at 14,000g for 5 min at 20 °C. The upper phase containing the lipids was then collected and dried down under nitrogen. The dry lipid extracts were reconstituted with 300 μl of 10 mM ammonium acetate in 9:1 methanol:toluene for analysis. The lower phase containing metabolites was subjected to further protein precipitation by adding four times of ice-cold 1:1:1 isopropanol:acetonitrile:water spiked in with 17 labeled internal standards and incubating for 2 h at −20 °C. The supernatant was dried down to completion under nitrogen and resuspended in 100 μl of 1:1 methanol:water for analysis.

Untargeted LC–MS metabolomics

Data acquisition

Metabolic extracts were analyzed four times using hydrophilic liquid chromatography (HILIC) and reverse-phase liquid chromatography (RPLC) separation in both positive and negative ionization modes as previously described¹²⁴. Data were acquired on a Thermo Q Exactive plus mass spectrometer equipped with a HESI-II probe and operated in full MS scan mode. MS/MS data were acquired on pool samples consisting of an equimolar mixture of all the samples in the study. HILIC experiments were performed using a ZIC-HILIC column 2.1 × 100 mm, 3.5 μm, 200 Å (Merck Millipore) and mobile-phase solvents consisting of 10 mM ammonium acetate in 50:50 acetonitrile:water (A) and 10 mM ammonium acetate in 95:5 acetonitrile:water (B). RPLC experiments were performed using a Zorbax SBaq column 2.1 × 50 mm, 1.7 μm, 100 Å (Agilent Technologies) and mobile-phase solvents consisting of 0.06% acetic acid in water (A) and 0.06% acetic acid in methanol (B). Data quality was ensured by (1) injecting 6- and 12-pool samples to equilibrate the LC–MS system before running the sequence for RPLC and HILIC, respectively, (2) sample randomization for metabolite extraction and data acquisition and (3) checking mass accuracy, retention time and peak shape of internal standards in every samples.

Data processing

Data from each mode were independently analyzed using Progenesis QI software v.2.3 (Nonlinear Dynamics). Metabolic features from blanks and that didn’t show sufficient linearity upon dilution were discarded. Only metabolic features present in >33% of the samples in each group were kept for further analysis and missing values were imputed by drawing from a random distribution of small values in the corresponding sample¹²⁵.

Metabolic feature annotation

Annotation confidence levels for each metabolite were provided following the Metabolomics Standards Initiative (MSI) confidence scheme. Peak annotation was first performed by matching experimental m/z, retention time and MS/MS spectra to an in-house library of analytical-grade standards (level 1). Remaining peaks were identified by matching experimental m/z and fragmentation spectra to publicly available databases including the Human Metabolome Database (HMDB) (http://www.hmdb.ca/), MassBank of North America (MoNA) (http://mona.fiehnlab.ucdavis.edu/) and MassBank (http://www.massbank.jp/) using the R package ‘MetID’ (v.0.2.0)¹²⁶ (level 2). Briefly, metabolic feature tables from Progenesis QI were matched to fragmentation spectra with an m/z and a retention time window of ±15 ppm and ± 30 s (HILIC) and ± 20 s (RPLC), respectively. When multiple MS/MS spectra matched a single metabolic feature, all matched MS/MS spectra were used for the identification. Next, MS1 and MS2 pairs were searched against public databases and a similarity score was calculated using the forward dot–product algorithm, which takes into account both fragments and intensities. Metabolites were reported if the similarity score was >0.4. Level 3 corresponds to unknown metabolites.

Targeted lipidomics with the Lipidyzer platform

Data acquisition

Lipid extracts were analyzed using the Lipidyzer platform that comprises a 5500 QTRAP System equipped with a SelexION differential mobility spectrometry interface (SCIEX) and a high flow LC-30AD solvent delivery unit (Shimazdu). A full description of the method is available elsewhere¹²³. Briefly, the lipid molecular species were identified and quantified using multiple reaction monitoring (MRM) and positive/negative switching. Two acquisition methods were employed covering ten lipid classes; method 1 had SelexION voltages turned on, whereas method 2 had SelexION voltages turned off. Lipidyzer data were reported by the Lipidomics Workflow Manager (LWM) v.1.0.5.0 software, which calculates concentration in nmol g⁻¹ for each detected lipid as average intensity of the analyte MRM/average intensity of the most structurally similar internal standard MRM multiplied by its concentration. Data quality was ensured by (1) tuning the differential mobility spectrometry compensation voltages using a set of lipid standards (SCIEX, no. 5040141) after each cleaning, more than 24 h of idling or 3 d of consecutive use, (2) performing a quick system suitability test (QSST) (SCIEX, no. 50407) before each batch to ensure acceptable limit of detection for each lipid class, (3) sample randomization for lipid extraction and data acquisition and (4) triplicate injection of lipids extracted from a reference plasma sample (SCIEX, no. 4386703) at the beginning of the batch.

Data preprocessing

Lipids detected in fewer than 66% of the samples in each group were discarded and missing values were imputed in each class by drawing from a random distribution of small values in the corresponding sample¹²⁵.

Differential analysis of metabolomics and lipidomics data

Metabolomics and lipidomics datasets were first normalized to the total protein content as determined by BCA assay (Pierce, no. 23225) to account for differential starting material quantity. Then, variance stabilizing normalization was applied to the data using ‘limma’ v.3.42.2, as recommended by previous studies^127,128. Differential analysis for metabolomic or lipidomic features was performed using ‘limma’ in R. Features with an FDR < 5% were considered to be statistically significant.

For analysis of the lipidomics data by lipid class (Supplementary Fig. 1d and Supplementary Table 1f), lipids were summarized by class after BCA and VSN corrections. To determine regulation at the class level, because of the small number of analyzed classes, a basic linear model approach was used through base R ‘lm’ function and FDR correction was applied using the base R ‘p.adjust’ function. Lipid classes with an FDR < 5% were considered statistically significant.

Functional enrichment analysis for metabolomics

As only one metabolite was significantly regulated with aging, we focused on analyzing differentially regulated metabolite sets with respect to sex. To analyze the directional regulation of metabolic pathways from untargeted metabolomics, we used R package ‘MetaboAnalystR’ v.2.0.4 (ref. ¹²⁹) to perform PSEA¹³⁰ of KEGG metabolic pathways. Using the ‘mummichog’ method¹³¹, we provided an input table of metabolic peaks represented by mass over charge ratios (m/z) and retention time, limma-derived P value and t-scores and analysis mode (negative or positive ion) to have maximum sensitivity for the functional enrichment analysis of untargeted metabolomic data. All significant gene sets (FDR < 5%) are reported in Supplementary Table 4f.

Integrated functional enrichment analysis for RNA-seq and metabolomics using IMPaLA

Based on the differential analyses results in RNA-seq and metabolomics, we focused on analyzing differentially regulated genes and metabolites with respect to sex. To provide an integrated view of RNA-seq and metabolomics results, we took advantage of the IMPaLA paradigm¹³². We used the IMPaLA v.12 web interface (http://impala.molgen.mpg.de/), which evaluates overrepresentation of genes and metabolites across 5,055 annotated pathways from 12 databases, with default parameters. For the joint analysis, we used genes with FDR < 5% with respect to biological sex, using gene symbols as input IDs. For metabolic features, we selected only the manually validated species (levels 1 and 2) with FDR < 5% with respect to biological sex. Finally, we also selected lipid species FDR < 5% with respect to biological sex. Importantly, we used HMDB IDs as input IDs for the metabolomic/lipidomic side of the analysis, thus analyzing only features with corresponding HMDB IDs. Joint analysis of significant features with a male bias, or, separately, those with female bias, were uploaded separately to the server. All pathways with an overall combined FDR < 5% are reported in Supplementary Table 4h.

Functional enrichment analysis for lipidomics

Since there was no significant difference in lipid composition with aging, we focused on analyzing differentially regulated lipid sets with respect to sex. The Lipid Ontology (LION) website was used to perform functional enrichment analysis of lipids¹³³. Lipid features with FDR < 5% with a male or female bias were uploaded to the server. For ease of reading, only the ten most-significant pathways with male and with female bias are shown in Fig. 3f. All significant LION terms are reported in Supplementary Table 4g.

Neutrophil ATAC-seq library generation

An independent cohort of C57BL/6Nia mice was used to assess age- and sex-related differences in chromatin architecture using ATAC-seq. To assay potential differences in the chromatin landscape of mouse bone marrow neutrophils across age and sex, we used the omni-ATAC protocol starting from 50,000 MACS-purified cells¹³⁴. Libraries were quality controlled on the Agilent Bioanalyzer 2100 platform at the USC Genome Core before pooling. Libraries were multiplexed and sequenced on the Illumina HiSeq-Xten platform as paired-end 150-bp reads at the Novogene Corporation.

ATAC-seq preprocessing pipeline

Paired-end ATAC-seq reads were adapter-trimmed using NGmerge0.2 (ref. ¹³⁵), which clips overhanging reads in a sequence-independent fashion and has been recommended for use in ATAC-seq preprocessing. Trimmed paired-end were mapped to the mm10 genome build using bowtie2 (v.2.2.6)¹³⁶. After alignment, PCR duplicates were removed using the ‘rmdup’ function of samtools (v.1.5). To minimize analytical artifacts from uneven sequencing depth between biological samples/libraries, libraries were randomly downsampled to the same depth for downstream analyses using PicardTools (v.2.20.3) or samtools (v.1.5). We used the peak finding function of HOMER (v.4.11) to identify ATAC-seq accessible regions¹³⁷.

Differential accessibility analysis with ATAC-seq

We extracted a normalized read count matrix from our downsampled bam files over merged HOMER regions using R package ‘Diffbind’ v.2.14 (ref. ¹³⁸). We used this extracted matrix for downstream analyses. Because of sample-to-sample variability, we used surrogate variable analysis¹¹⁵ to estimate and correct for experimental noise, similar to the RNA-seq analysis. R package ‘sva’ v.3.34 was used to estimate surrogate variables and the ‘removeBatchEffect’ function from ‘limma’ v.3.42.2 was used to regress out the effects of surrogate variables, PCR duplicate content and variation in reads mapping to peaks according to Diffbind (fraction of reads in peaks). The ‘DESeq2’ R package (v.1.26.0) was used for further processing of ATAC-seq data. Regions with FDR < 5% were considered statistically significant.

Functional enrichment analysis of differentially accessible regions from ATAC-seq

We used the significantly differentially accessible peaks identified by DESeq2 at FDR < 5% to analyze functional enrichments linked to these regions. For this purpose, we leveraged the ‘GREAT’ v.4.0.4 tool¹³⁹, which is specifically optimized to identify functional enrichment starting from genomic regions. We used all ATAC-seq accessible regions as the background for enrichment. All other options were left as default parameters. Results were exported as ‘tsv’ and processed in R to filter significant regions at FDR < 5%. Filtered results are reported in Supplementary Tables 3f and 4i.

Chromatin architecture analysis from ATAC-seq datasets

To analyze the underlying chromatin architecture differences between neutrophils isolated from different biological groups, we used the NucleoATAC v.0.3.4 software⁶⁸. As NucleoATAC requires high sequencing depth to reliably measure nucleosome profiles, we pooled depth-matched reads from each biological group to attain this depth as recommended by authors of the software⁶⁸. We extracted similar metrics from the NucleoATAC output to what was described in previous studies using this software to investigate chromatin architecture¹⁴⁰. All comparisons between groups were tested for statistical significance using a two-tailed Wilcoxon rank-sum test (nonparametric test).

Repetitive element transcription analysis

To evaluate potential differential regulation of transposable elements, we used the ‘analyzeRepeats.pl’ functionality of the HOMER software to count STAR-mapped reads over repetitive elements¹³⁷. To allow for whole-library normalization, reads mapping over transcripts were also counted using the same HOMER script. Read counts were imported into R to estimate differential repeat expression using the DESeq2 (v.1.26.0) R package¹¹⁶.

Neutrophil NETosis assay

We used a flow cytometry-based NETosis assay protocol using extracellularized DNA staining by SYTOX Green to quantify cells undergoing NETosis, adapted from a flow-cytometry-based protocol¹⁴¹ and a 96-well plate plate-reader protocol¹⁴², both with neutrophils in suspension. Briefly, MACS-purified neutrophils were resuspended in a concentration of 2 × 10⁶ cells ml⁻¹ in neutrophil culture medium (RPMI 1640 without phenol red (Hyclone), supplemented with 1% penicillin/streptomycin (Corning) and 0.1% BSA (Akron Biotechnology)). We then used 1 × 10⁶ cells, aliquoted into sterile microcentrifuge tubes, one per starting biological sample. SYTOX green (Thermo Fisher Scientific, no. S7020) was added to each sample to a final concentration of 200 nM. Neutrophil medium supplemented with dimethylsulfoxide (DMSO) (vehicle) or 50 nM PMA (Sigma, no. P1565) was added to each tube. Tubes were slowly inverted three times to mix. Then, 2 × 10⁵ cells were seeded in technical quadruplicates in wells of a sterile black 96-well plate (Greiner Bio-One) and incubated in a humidified incubator with 5% CO₂ at 37 °C for 2 h. The fraction of cells positive for SYTOX Green in each well was quantified using the MACSQuant10 and analyzed using Flowlogic v.7. To account for differences in basal levels of NETosis across samples, NETosis was expressed as induction of NETosis as (median of PMA technical quadruplicates) / (median of DMSO technical quadruplicates). The raw cytometry dataset was deposited to Figshare (https://doi.org/10.6084/m9.figshare.14043923.v1).

TNF-α-primed NETosis analysis

To compare our results with previously published results on NETosis in aged neutrophils^13,14, which were obtained with direct TNF-α priming¹³ or in conditions known to mimic TNF-α priming¹⁴ (thioglycolate elicitation of peritoneal neutrophils⁷¹), we also performed NETosis analysis in primed conditions, using a slightly modified protocol from above to include a priming step. Specifically, MACS-purified neutrophils cells were resuspended in neutrophil culture medium and first pre-warmed after isolation in a humidified incubator with 5% CO₂ at 37 °C for 15 min. Then, 2 × 10⁶ cells were aliquoted into a sterile microcentrifuge tube and primed with 10 ng ml⁻¹ mouse TNF-α (PeProTech, no. 315-01 A) in a humidified incubator with 5% CO₂ at 37 °C for 15 min. Cells were further aliquoted, supplemented with SYTOX green, DMSO/PMA, incubated and analyzed as above. The raw cytometry dataset was deposited to Figshare (https://doi.org/10.6084/m9.figshare.14043923.v1).

Feature extraction for machine-learning analysis

For each significantly age-regulated or sex-dimorphic gene in neutrophils according to DEseq2 (FDR < 5%), we extracted a number of features associated with that gene. First, we took advantage of our ATAC-seq data to evaluate the chromatin architecture of neutrophils (see above). For each gene, we extracted (1) the average promoter accessibility by ATAC-seq in fragments per kilobase of transcript per million mapped reads across all four conditions (promoters were defined as the region between −500 bp to +150 bp with regard to annotated transcriptional start sites in mm10 build according to HOMER), as well as (2) the average log₂(fold change) in accessibility between old and young neutrophils (age-regulation models) or the average log₂(fold change) in accessibility between female and male neutrophils (sex-dimorphism models). Second, we took advantage of gene sets coverage known of predicted targets of transcriptions factors: ChEA, ENCODE, JASPAR and GEO TF perturbation information from the Harmonizome database¹²² and transcriptional targets of FOXO transcription factors from GEO experiments compiled in our previous study¹⁸. To reduce extraneous features, we engineered a TF target feature following these steps: (1) TF targets were summarized from all putative sources, (2) only TFs with evidence of expression in primary neutrophils according to RNA-seq were retained as features and (3) only TFs with at least 25 targets across significant age-related genes (age-regulation models) or sex-dimorphic genes (sex-dimorphism models) were retained. This process yielded 357 (age-regulation models) and 249 (sex-dimorphism models) TF target sets used as features for model training. Finally, we included several DNA sequence features for each gene: the percentage of CG nucleotide and the percentage of CpG dinucleotide in promoters and exons, computed using HOMER. To note, HOMER was only able to provide information for 3,421 of 3,653 genes with significant age regulation and 1,636 of 1,734 genes with significant sex-dimorphic regulation. Finally, to determine whether location of the gene on autosomes versus chromosomes was a key predictive factor for sex-dimorphic gene expression, we also included, for the sex-dimorphic gene expression models, a feature encoding whether the gene is located on autosomes or sex chromosomes (X or Y).

Machine-learning analysis for age-regulated or sex-dimorphic genes

We trained machine-learning classification models for two questions: (1) models for age-regulated gene expression and (2) models for sex-dimorphic gene expression. For (1), age-regulation machine-learning models were built to learn whether up- or downregulated genes could be discriminated using genomic features and predict potential master regulators. For (2), sex-dimorphism machine-learning models were built to learn whether female and male-biased genes could be discriminated using genomic features and predict potential master regulators.

In both cases, we built classification models using seven classification algorithms as implemented through R package ‘caret’ (v.6.0-86). Auxiliary R packages were used with caret to implement neural networks (v.7.3-13), RFs (v.4.6-14), GBMs (v.2.1.5), radial kernel support vector machines (v.0.9-29), sparse linear discriminant analysis (v.0.1-9), conditional inference trees (‘party’, v.1.3-4) and regularized logistic regression (‘LiblineaR’, v.2.10-8). Caret was allowed to optimize final model parameters on training data using tenfold cross validation. Accuracies, sensitivities, specificities and area under the curve (using package ‘pROC’, v.1.16.2) for all trained classifiers were estimated using a test set of randomly held-out one-third of the data (not used in the training phase) obtained using the ‘createDataPartition’ function. Feature importance estimation was only conducted using tree-based RF and GBM methods, as other algorithms do not natively allow for it. The feature importance was computed and scaled by caret for the RF and GBM models.

Reanalysis of DIA proteomics data from human neutrophils

We obtained data-independent acquisition (DIA) proteomics results from the supplemental material of a recent study⁹⁷. We used normalized expression values from the article’s supplemental material and used reported ‘donor-specific Spectronaut protein expression values’ for further analysis. To exclude potential confounders linked to disease, we used only data from 68 healthy donors and excluded data from patients with disease. As the article did not report the biological sex of donors, we used detected expression levels of DDX3Y, a Y-chromosome-encoded protein, as a proxy for the likelihood of the sample belonging to a male donor. We then examined rank correlation statistics between protein expression of DDX3Y compared to that of human orthologs of top male-biased primary granules genes from our mouse RNA-seq data (ELANE, MPO, PRTN3 and CTSG). The significance of the Spearman rank correlation is reported.

Mouse sepsis model through intraperitoneal LPS injection

Young adult (3–4-month-old) C57BL/6J mice were exposed to LPS, a pathogen-associated molecular pattern found in the cell wall of Gram-negative bacteria, to elicit a sepsis-like response⁹⁶. Briefly, mice were injected intraperitoneally with 2 µg of LPS per gram of body weight, using a sterile LPS stock (Sigma, no. L5293) diluted in PBS and control animals were injected with sterile PBS⁹⁶. Animals were monitored hourly for up to 6 h after injection for the occurrence and severity of endotoxemia⁹⁶. We did not observe gross differences in the state of female versus male animals injected with LPS over the course of the experiment, although it is possible that such differences would have emerged at longer time points. After killing by CO₂ asphyxiation and cervical dislocation at 1, 3 and 6 h after injection, blood was collected from the heart. A total of 80 animals across three independent experimental days were used (n = 8 animals per sex, age and time group for LPS injections; n = 5 animals per sex, age and time group for control PBS injections). To obtain serum for downstream analysis, the blood was left to clot for ≥1 h at room temperature in MiniCollect Serum SeparatorTubes (Greiner Bio-One). Serum was then separated from the clot using centrifugation at 2,000g for 10 min and frozen at −80 °C until further use.

Quantification of serum neutrophil elastase (ELANE) by ELISA

Quantitative evaluation of circulating ELANE levels was performed from serum samples. ELISA quantification of serum ELANE levels was performed using Abcam’s Mouse Neutrophil Elastase SimpleStep ELISA kit (ab252356) according to the manufacturer’s instructions. Technical replicates from the same sample were averaged as one value before statistical analysis and plotting.

Statistics and reproducibility

All statistical analyses were conducted using R v.3.6.0–3.6.3. For all ‘omic’ analyses, results were corrected for multiple testing by the use of FDR significance correction. All statistical analyses were performed as specified in figure legends and in the corresponding methods section. We used nonparametric statistical tests whenever possible to avoid assuming normality of data distributions. No statistical methods were used to predetermine sample size, since effect sizes and variance were not known a priori. No data passing quality control checks were excluded from the analyses. All ‘omic’ samples were performed on independent biological samples, usually from one cohort of animals to minimize known issues linked to batch effects in such analyses¹⁴³. As one exception, RNA-seq of bone marrow neutrophils was independently replicated on a smaller cohort, as neutrophils are known to be RNA-poor, which could add more than usual biological variability (Extended Data Fig. 2). With the exception of a test of the heterogeneity flow cytometry panel on bone-marrow cells (used to validate c-kit staining; Extended Data Fig. 4b and Supplementary Table 2b), all other experiments were performed on animals from at least two independent cohorts.

For each experiment, animals were processed in an alternating/snaking order rather than in large homogeneous groups to minimize group, batch and circadian effects. For observational aging studies, no randomization was possible as groups were biologically determined. For LPS injection experiments, animals were randomly allocated to the PBS or LPS group for time point of euthanasia, on a cage-level basis (animals from the same cage were randomly attributed to LPS or PBS injections; for each cohort, co-housed animals were used for the same time point; cages were randomly attributed to a specific time point for euthanasia). This process was performed independently for each of the three cohorts for LPS-induced endotoxemia experiments. Blinding was not relevant to the studies conducted as data were collected by automated systems (sequencing, MS, flow cytometry) or were in numeric form (ImageJ quantification of western blot bands, ELISA absorbance values), which cannot be influenced by subjectivity from the experimenter.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.