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. To date, this dataset
   represents the largest multi-omics resource for neutrophils across sex
   and ages. This resource identifies neutrophil characteristics which
   could be targeted to improve immune responses as a function of sex
   and/or age.

Introduction

   Neutrophils are the most abundant cells in human blood, representing
   50–70% of leukocytes^[38]1. These cells are continually produced in the
   bone marrow and released into circulation^[39]2, [40]3. Neutrophils are
   short-lived, with estimated cellular lifespans of ~10–18 hours once
   released in the bloodstream^[41]4, [42]5, although they can survive
   longer^[43]4, [44]6. Throughout their cellular lifespan, neutrophils
   undergo “neutrophil aging”, a process distinct from organismal
   aging^[45]3, [46]7. Neutrophils perform many key functions, including
   production of antimicrobial granules and of “Neutrophil Extracellular
   Traps” [NETs]^[47]3, [48]8. Although neutrophils are essential against
   infections as the “first line of defense”, their aberrant activation
   can also aggravate inflammatory disease^[49]2, [50]3. Indeed, emerging
   evidence suggests that neutrophils play important roles in chronic
   inflammation^[51]9.

   Organismal aging in mammals is characterized by systematic immune
   dysfunction and chronic inflammation, a phenomenon dubbed
   “inflamm-aging”^[52]10, which may be partially mediated by dysfunction
   of innate immune cells^[53]11. Indeed, emerging evidence suggests that
   neutrophils from aged organisms are dysfunctional^[54]12, [55]13,
   [56]14, [57]15, [58]16. Observed age-related dysfunctions in
   neutrophils include decreased NETosis in TNFα-primed conditions^[59]13,
   [60]14, reduced chemotaxis to sites of inflammation^[61]15 and
   secretion of intracellular granule proteases^[62]15, [63]17. Although
   gene expression changes throughout lifespan have been reported across
   many cell types^[64]18, [65]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^[66]20 and could
   result from differential “omic” regulation^[67]21, [68]22. Although
   transcriptional differences between young female and male murine
   neutrophils have been profiled by ImmGen^[69]22, 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,
   e.g. inflammatory mediators production^[70]9, or functional modulation
   by testosterone^[71]23. 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 ATAC-seq, we showed that
   remodeling of chromatin-related pathways was associated with overall
   differences in the chromatin architecture of neutrophils from young vs.
   old and female vs. 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
   months) and old (20–21 months) C57BL/6Nia female and male mice
   ([72]Figure 1a). Primary neutrophils were isolated from bone marrow
   using magnetic-activated cell sorting [MACS], and profiling was then
   performed on purified neutrophils: (i) transcriptome profiling by
   RNA-sequencing [RNA-seq], (ii) metabolomic profiling by untargeted
   liquid chromatography coupled with mass spectrometry [LC-MS], and (iii)
   lipidomic profiling by targeted MS ([73]Figure 1a).

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

   Figure 1:
   [74]Open in a new tab

   (a) Experimental setup scheme. (b-d) Multidimensional Scaling analysis
   results of RNA expression by RNA-seq (b), untargeted metabolomics (c),
   or targeted lipidomics (d). (e-g) Heatmap 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. MACS: Magnetic-Activated Cell Sorting. MDS:
   Multidimensional Scaling.

   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 ([75]Figure 1b–[76]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 ([77]Figure 1 b–[78]d). To better understand the nature of
   differences between neutrophils from young vs. old [with sex as a
   covariate], and female vs. 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 ([79]Figure 1e–[80]g and
   [81]Extended Data Figure 1a–[82]f; [83]Supplementary Table S1A–[84]F).
   We quality-checked our dataset for appropriate expression of
   sex-specific genes ([85]Extended Data Figure 1b,[86]c). Finally, we
   confirmed differential gene expression trends between groups using a
   small replicate RNA-seq cohort ([87]Extended Data Figure 2a,[88]b), and
   comparing our data with published datasets from female vs. males mouse
   spleen neutrophils^[89]22 and human blood neutrophils^[90]24
   ([91]Extended Data Figure 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%;
   [92]Figure 1e–[93]g and [94]Extended Data Figure 1a–[95]d;
   [96]Supplementary Table S1B,[97]D,[98]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 ([99]Extended Data Figure 1a,[100]e). Importantly,
   we observed no biases in the purity of MACS-purified neutrophils across
   groups ([101]Extended Data Figure 3a–[102]b), and no difference in
   neutrophil abundance in bone-marrow ([103]Extended Data Figure 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 ([104]Extended Data Figure 1a,[105]f; [106]Supplementary Table
   S1), few to no metabolomic and lipidomic changes were observed
   ([107]Extended Data Figure 1a; [108]Supplementary Table S1C,[109]F).
   This is consistent with overall poorer age-related separation observed
   for metabolomics and lipidomics ([110]Figure 1c–[111]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 (e.g. 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; [112]Extended Data Figure 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 vs. female neutrophils ([113]Extended Data
   Figure 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, since bone marrow neutrophils can be heterogenous, we
   evaluated the composition of MACS-purified neutrophils across groups
   ([114]Extended Data Figure 4a–[115]d; [116]Supplementary Table
   S2A,[117]B). We leveraged flow cytometry markers to estimate the
   proportions of neutrophil progenitors (i.e. pre- and pro-neutrophils),
   immature vs. mature neutrophils, and fresh vs. “aged” neutrophils in
   our population. Neutrophil progenitors were extremely rare among
   MACS-purified cells (i.e. pre- and pro-neutrophils represent < 0.5% of
   cells; [118]Extended Data Figure 4a; [119]Supplementary Table S2A), and
   are unlikely to drive large “omic” differences. Importantly, we did not
   observe significant differences in the presence of mature vs. immature
   neutrophils across groups ([120]Extended Data Figure 4a,[121]c;
   [122]Supplementary Table S2). Finally, we observed a significant trend
   for the increased presence of senescent/aged neutrophils (CD62L^−)
   among the mature neutrophils ([123]Extended Data Figure 4a,[124]d;
   [125]Supplementary Table S2A). Our observation is consistent with a
   previous report of accumulation of senescent/aged neutrophils in old
   mice due to decreased clearance by efferocytosis^[126]25. Based on
   previous studies of senescent/aged neutrophils, changes in the relative
   proportion of “fresh” vs. “aged” neutrophils among bone marrow
   neutrophils may have functional impacts on inflammation, chemotaxis,
   and NETosis^[127]26, [128]27.

Age-associated changes in chromatin pathways in neutrophils

   We first focused on the impact of organismal aging on neutrophils using
   pathway enrichment analysis ([129]Figure 2a–[130]c and [131]Extended
   Data Figure 5a–[132]b; [133]Supplementary Table S3). To note, we only
   analyzed age-related functional remodeling at the transcriptomic level,
   since metabolomics and lipidomics showed little changes with aging
   ([134]Extended Data Figure 1a).

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

   Figure 2:
   [135]Open in a new tab

   (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 ^[136]109. Blue (red) shades are associated to 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 10 most up- and top 10 most
   downregulated gene sets are plotted for readability. Full lists and
   statistics available in [137]Supplementary Table S3. Also see
   [138]Extended Data Figure 5. Shown pathways with FDR < 5%. NES:
   Normalized Enrichment Score (for GSEA analysis). FDR: False Discovery
   Rate.

   We assessed the coordination of age-related changes by performing a
   network analysis ([139]Figure 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 ([140]Figure
   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 ([141]Figure 2b
   and [142]Extended Data Figure 5a,[143]b; [144]Supplementary Table S3).
   Immune-related pathways were significantly upregulated according to
   Gene Set Enrichment Analysis [GSEA] ([145]Figure 2b and [146]Extended
   Data Figure 5a,[147]b; [148]Supplementary Table S3A–[149]C) and
   Ingenuity Pathway Analysis [IPA] ([150]Supplemental Table S3E). This
   suggests that neutrophil-associated immunity in older animals could
   have different characteristics, consistent with age-related immune
   dysfunction^[151]10.

   Intriguingly, the top downregulated pathways with organismal aging were
   overwhelmingly related to chromatin and cell cycle ([152]Figure 2b and
   [153]Extended Data Figure 5a,[154]b; [155]Supplementary Table
   S3A–[156]C,[157]E). Specifically, 18 histone-encoding genes were
   significantly downregulated with aging ([158]Supplemental Table S1A).
   Importantly, we did not observe significant differences in the
   cell-cycle phase distribution of MACS-isolated neutrophils across
   groups ([159]Extended Data Figure 6a,[160]b), suggesting that
   differential expression of cell-cycle pathways is not a byproduct of
   differential proliferation.

   Changes in the regulation of chromatin components are especially
   relevant to neutrophil biology. Indeed, neutrophils have a unique
   chromatin organization^[161]28, with increased chromatin compaction or
   “supercontraction”^[162]29. In addition, neutrophils can extrude their
   chromatin to produce NETs in response to pathogens^[163]30. The
   extrusion of chromatin directly participates in pathogen killing,
   notably thanks to antimicrobial properties of histones^[164]31.
   Chromatin relaxation can also be potentiated by nuclear-translocated
   granule-derived proteases (e.g. Elane), ultimately helping bacterial
   killing^[165]32. Finally, although neutrophils are post-mitotic,
   cell-cycle genes can control NET production^[166]33. 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 ([167]Figure 2b; [168]Supplementary
   Table S3B,[169]C,[170]E). This is consistent with loss of proteostasis
   being a “hallmark of aging”^[171]34. Control of autophagy is critical
   for neutrophil differentiation^[172]35, regulation of NET
   formation^[173]36 and of degranulation^[174]37.

   To identify candidate upstream regulators, we investigated whether
   target genes of specific transcription factors [TFs] were
   differentially regulated with aging ([175]Figure 2c; [176]Supplementary
   Table S3D). We obtained lists of TF target genes derived from ChEA,
   JASPAR or GEO through the Harmonizome (see [177]methods). Importantly,
   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 (e.g.
   chromatin-remodeling enzymes, 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^[178]38,
   [179]39 ([180]Figure 2c). Although E2f7 has not been studied in
   neutrophils, it is highly expressed in committed progenitors^[181]40.
   Targets of Foxo1 were also significantly downregulated with aging in
   neutrophils ([182]Figure 2c). Foxo TFs are the primary targets of the
   aging-regulating insulin/insulin-like growth factor 1
   signaling^[183]41, and Foxo1 regulates neutrophil-mediated bacterial
   immunity^[184]42. Finally, known targets of Padi4 were significantly
   downregulated with aging ([185]Figure 2c). Padi4 is a peptidylarginine
   deaminase, catalyzing histone citrullination, which plays a key role in
   NETosis in vivo^[186]43, [187]44. Transcript levels of Padi4 were
   significantly upregulated with aging (FDR = 6.9×10^−3; [188]Extended
   Data Figure 5c; [189]Supplementary Table S1B). Protein levels of Padi4
   were also significantly upregulated with aging regardless of sex (p =
   8.5×10^−3; [190]Extended Data Figure 5d,[191]e). This is consistent
   with reports of histone citrullination associating to transcriptional
   repression^[192]45, [193]46, [194]47, although the directionality of
   regulation may depend on cofactors and cell types^[195]48. 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 ([196]Figure
   3a–[197]f and [198]Extended Data Figure 7a–[199]c; [200]Supplementary
   Table S4). We first assessed the interconnection of sex-dimorphic genes
   by constructing putative PPI networks for female vs. male-biased gene
   expression ([201]Figure 3a,[202]b). Sex-biased gene programs showed
   clear interconnection, suggesting coherent differences with likely
   functional impact ([203]Figure 3a,[204]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^[205]49
   ([206]Figure 3a). Interestingly, Irf8 was the top connected node in the
   male-biased gene network ([207]Figure 3b). Although not specifically
   studied in neutrophils, Irf8 is an interferon signaling-related TF that
   regulates myeloid fate determination, usually suppressing neutrophil
   differentiation^[208]50, [209]51. In this context, Irf8 is likely to
   accomplish functions unrelated to myeloid differentiation (e.g.
   cytokine production)^[210]52, since no significant changes in relative
   bone marrow neutrophil abundances were observed ([211]Extended Data
   Figure 3c).

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

   Figure 3:
   [212]Open in a new tab

   (a-b) NetworkAnalyst predicted PPI networks for genes displaying
   significant bias in gene expression by RNA-seq towards female (a) or
   male (b) neutrophils. Network edges are based on IMEx/InnateDB data, a
   knowledgebase specifically geared for analyses of innate immune
   networks ^[213]109. (c-d) Top enriched gene sets from Reactome (c) and
   transcription factor targets (d) using GSEA for differential RNA
   expression. (e) Phenotype Set Enrichment Analysis (PSEA) for
   differential metabolite regulation. (f) Top Lipid Ontology (LION)
   functional enrichment analysis for differential lipid regulation. For
   (c-f), the top 10 most up- and top 10 most downregulated gene sets (if
   that many) are plotted for readability. Full lists and statistics
   available in [214]Supplementary Table S4. Also see [215]Extended Data
   Figure 7. Shown pathways with FDR < 5%. NES: Normalized Enrichment
   Score (for GSEA and PSEA analysis). FDR: False Discovery Rate.

   Next, we asked which pathways were regulated in a sex-dimorphic fashion
   in neutrophils ([216]Figure 3c and [217]Extended Data Figure
   7a,[218]b;` [219]Supplementary Table S4). Functional enrichment
   analysis revealed that significant female-biased pathways encompassed
   extra-cellular matrix [ECM] and cell surface-related pathways
   ([220]Figure 3c and [221]Extended Data Figure 7a,[222]b;
   [223]Supplementary Table S4A–[224]D). Indeed, collagen-encoding genes
   Col1a1 and Col1a2 were expressed at higher levels in female neutrophils
   ([225]Supplementary Table S1B). Collagen production, usually from
   fibroblasts, is a key event in fibrosis^[226]53, and neutrophils play
   an important role in the development of fibrotic lesions^[227]54.
   Alternatively, differential expression of cell surface-related genes
   may also impact the ability of neutrophils to migrate through
   endothelial barriers^[228]55. Interestingly, autophagy-related gene
   sets were female-biased in our dataset ([229]Supplementary Table
   S4A,[230]C,[231]E). As mentioned above, autophagy control is critical
   for neutrophil biology^[232]35, [233]36, including neutrophil
   degranulation^[234]37.

   Intriguingly, chromatin- and cell cycle-related pathways were
   overwhelmingly biased for higher expression in male neutrophils
   ([235]Figure 3c; [236]Extended Data Figure 7a,[237]b;
   [238]Supplementary Table S4A–[239]C). Consistently, 21 histone-encoding
   genes showed significant male-biased expression ([240]Supplementary
   Table S1B). 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 ([241]Figure 3d;
   [242]Supplementary Table S4D). Consistent with pathway analysis,
   targets of cell cycle related E2f7 and Mybl1^[243]39, [244]56 were more
   highly expressed in male neutrophils ([245]Figure 3d). Female
   neutrophils showed higher expression of Foxo4 targets ([246]Figure 3d).
   Interestingly, the Insulin/Igf1 signaling pathway, which regulates Foxo
   TF activity, has been linked to multiple sex-dimorphic
   phenotypes^[247]57, [248]58.

   We next examined sex-dimorphism based on metabolomics and identified a
   number of pathways with significant sex-bias ([249]Figures 3e;
   [250]Extended Data Figure 7c; [251]Supplementary Table S4F).
   Interestingly, metabolic pathways related to nucleotide metabolism were
   differentially regulated between sexes ([252]Figures 3e; [253]Extended
   Data Figure 7c; [254]Supplementary Table S4F). Indeed, a number of
   sex-dimorphic metabolic peaks were predicted to represent nucleotide
   species (e.g. AMP, ATP, GMP; see [255]Supplementary Table S1D).
   Sex-differences in nucleotide pools are consistent with transcriptomic
   observations of sex-dimorphism in (i) pathways directly related to
   nucleotide metabolism, which could lead to direct changes in nucleotide
   pools, and (ii) pathways responsive to nucleotide levels, such as
   signaling by Rho-GTPases, which can regulate neutrophil
   recruitment^[256]59 ([257]Supplementary Table S4A–[258]C,[259]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’^[260]60. Adenine
   nucleotides from neutrophils exert anti-inflammatory effects^[261]61,
   increase endothelial barrier function, attenuate neutrophil adhesion to
   endothelial cells, and modulate transendothelial migration^[262]62.
   Functional analysis of metabolomic data also suggested significant
   sex-differences in amino-acid metabolism, such as tryptophan and
   arginine/proline metabolism ([263]Figures 3e; [264]Extended Data Figure
   7c). Interestingly, arginine and tryptophan metabolism play important
   immunoregulatory roles^[265]63. 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 ([266]Figures 3f, [267]Extended
   Data Figure 1d, [268]Supplemental Table S4G). Interestingly, male
   neutrophils were strongly enriched for triacylglycerols [TAG], stored
   in lipid droplets, whereas female neutrophils were enriched for
   diacylglycerols [DAG] ([269]Figures 3f; [270]Supplemental Table S4G;
   [271]Extended Data Figure 1d). Consistently, genes associated to the GO
   term “negative regulation of lipid storage” were expressed at
   significantly higher levels in female neutrophils (FDR = 4.33×10^−3;
   [272]Supplemental Table S4C). Lipid droplets are crucial during
   neutrophil differentiation and as a source for inflammatory
   mediators^[273]35, [274]64. Adipose triglyceride lipase (Atgl), encoded
   by Pnpla2, regulates neutral lipid storage into lipid droplets in
   neutrophils^[275]65. Interestingly, Pnpla2 is expressed at higher
   levels in female neutrophils (FDR = 0.01; [276]Extended Data Figure 7d;
   [277]Supplemental Table S1B). Indeed, Atgl catalyzes the conversion of
   TAG to DAG, consistent with higher levels of DAG in female neutrophils
   and of TAG in male neutrophils ([278]Figure 3f, [279]Extended Data
   Figure 1d and [280]Extended Data Figure 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^[281]65. In addition, Atg7 is crucial
   for generation of free fatty acids [FFA] from lipid-droplets in
   neutrophils, and Atg7-deficient neutrophils show increased
   lipid-droplet storage^[282]35. Consistently, Atg7 is expressed at
   higher levels in female neutrophils (FDR = 0.03; [283]Extended Data
   Figure 7d; [284]Supplemental Table S1B). Finally, Lpin1 is crucial for
   synthesis of DAG from phosphatidic acid, and plays a role in the
   regulation of inflammation^[285]66. Consistent with increased DAG
   levels, Lpin1 is also expressed at higher levels in female neutrophils
   (FDR = 0.03; [286]Extended Data Figure 7d; [287]Supplemental Table
   S1B). Thus, both lipidomics and transcriptomics data suggest that male
   vs. female neutrophils have increased neutral lipid storage (or
   decreased usage) ([288]Extended Data Figure 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
   ([289]Supplementary Table S4H). We observed joint enrichment of
   male-biased molecules in pathways linked to cell cycle, DNA metabolism
   and chromatin-related pathways ([290]Supplementary Table S4H), 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 ([291]Supplementary Table
   S4H). 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^[292]28, [293]29. More than just a gene expression regulatory
   layer, neutrophil chromatin directly participates in antimicrobial
   responses through NETs^[294]30, [295]31. Thus, chromatin-metabolism
   differences could lead to profound changes in neutrophil biology. To
   directly evaluate neutrophil chromatin, we utilized the Assay for
   Transposon-Accessible Chromatin followed by sequencing
   [ATAC-seq]^[296]67. ATAC-seq takes advantage of adapter-loaded Tn5
   particles inserting into accessible chromatin, yielding information
   about local chromatin accessibility^[297]67 and higher-order
   organization^[298]68.

   To evaluate chromatin landscapes, we isolated neutrophils from an
   independent cohort of young (4–5 months) and old (20–21 months)
   C57BL/6Nia female and male mice and performed ATAC-seq ([299]Figure 4a;
   [300]Extended Data Figure 8a–[301]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 ([302]Extended Data Figure 8d).
   In addition, there were substantially more regions with age-related vs.
   sex-dimorphic differential accessibility ([303]Extended Data Figure
   8e–[304]h), consistent with numbers of differentially expressed genes
   with respect to age and sex ([305]Extended Data Figure 1a).
   Interestingly, enriched GO terms associated to differentially
   accessible regions with respect to organismal age (FDR < 5%) were
   mostly associated to regions with increased accessibility, and
   encompassed terms related to immune response ([306]Supplementary Table
   S3F). 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 ([307]Supplementary Table S4I).

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

   Figure 4:
   [308]Open in a new tab

   (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) Boxplot of
   nucleosome occupancy as calculated by NucleoATAC for high-confidence
   nucleosomes across groups. Also see [309]Extended Data Figure 8i. (e)
   Boxplot of nucleosome fuzziness as calculated by NucleoATAC for
   high-confidence nucleosomes across experimental groups. (f) Boxplot of
   nucleosome inter-dyad distance as calculated by NucleoATAC across
   experimental groups. The horizontal red line in panels (d-f) shows the
   median value in neutrophils from young females for ease of comparison.
   Significance in non-parametric two-sided Wilcoxon rank-sum tests is
   reported for panels (d-f). Black p-values represent differences between
   male vs. female neutrophils in each age group; pink (blue) p-values
   represent age-related differences in female (male) neutrophils. Also
   see [310]Extended Data Figure 8i,[311]j. (g) Workflow for the in vitro
   NETosis inducibility experiment. SYTOX Green was used to stain
   extracellular DNA, a proxy for NETosing cells. (h) Boxplot of NETosis
   induction in naïve neutrophils. Each point represents one animal,
   median of 4 technical replicates. Animals from 4 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) Boxplot of NETosis induction in neutrophils, primed with
   10ng/mL mouse TNFα for 10 min prior to NETosis induction. Each point
   represents one animal, median of 4 technical replicates. Animals from 4
   independent cohorts, n = 17 per group. Significance in non-parametric
   two-sided Wilcoxon rank-sum tests are reported for panels (h-i). Black
   p-values represent differences between male vs. female neutrophils in
   each age group; pink (blue) p-values represent age-related differences
   in female (male) neutrophils. For all boxplots in panels (d-f, h-i),
   the center line represents the sample median, the box limits consist of
   the 25^th and 75^th percentiles, the whiskers span 1.5x 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.

   Since enrichment analyses suggested that overall aspects of chromatin
   are differentially regulated with organismal age and sex, we leveraged
   NucleoATAC^[312]68 to analyze neutrophil chromatin architecture.
   Interestingly, nucleosome “v-plots” generated revealed a region at
   ~145–150bp which captured more reads in male vs. female neutrophils,
   regardless of animal age, suggesting global differences in nucleosome
   packaging ([313]Extended Data Figure 8a). Next, we leveraged indicators
   linked to chromatin compaction and organization around nucleosomes: (i)
   nucleosome occupancy, which indicates how frequently a position is
   occupied by a nucleosome across cells, where higher occupancy is
   associated to tighter chromatin structure, (ii) nucleosome fuzziness,
   which indicates how well-positioned the nucleosome is, where fuzzier
   nucleosomes correspond to looser chromatin, and (iii) inter-dyad
   distance, which measure DNA length between neighboring nucleosomes,
   with increased distance associated to looser chromatin.

   Occupancy metaplots around transcriptional start sites [TSS] of
   expressed genes revealed that, regardless of age, male neutrophils
   showed higher median nucleosomal occupancy compared to female
   neutrophils ([314]Figure 4b). Similarly, regardless of sex, aging was
   associated to increased median nucleosomal occupancy ([315]Figure 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 ([316]Figure 4d, [317]Extended Data Figure 8i). We observed
   decreased nucleosome fuzziness in male vs. female neutrophils, as well
   as in old vs. young neutrophils ([318]Figure 4e). Finally, we observed
   shorter inter-dyad distance in male vs. female, as well as old vs.
   young neutrophils ([319]Figure 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 vs. female and old vs. 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 ([320]Extended Data Figure
   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 ([321]Extended Data Figure
   8k–[322]l; [323]Supplementary Table S1G,[324]H). Eukaryotic genomes
   contain large proportions of repeats, including transposons, whose
   transcription is usually tightly repressed through compaction^[325]69.
   Consistent with a more compact chromatin architecture, male neutrophils
   showed lower transcription of repeats vs. females (350 female-biased
   elements vs. 12 male-biased; [326]Extended Data Figure 8l;
   [327]Supplementary Table S1H). Similarly, old neutrophils also showed
   lower transcription of repeats vs. young neutrophils (115 elements
   downregulated vs. 27 upregulated with age; [328]Extended Data Figure
   8k; [329]Supplementary Table S1G). Interestingly, chromatin
   decompaction is a limiting step of NETosis^[330]70, 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 vs. old, female vs. male mice ([331]Figure 4g–[332]i). We used
   flow cytometry to quantify cells with extracellularized DNA (i.e.
   NETosed or NETosing cells), after a 2h incubation in the presence of
   Phorbol 12-myristate 13-acetate [PMA], a known activator of NETosis,
   compared to vehicle ([333]Figure 4g). We first evaluated NETosis
   induction on unstimulated neutrophils ([334]Figure 4h). NETosis
   inducibility was significantly increased with organismal age in female
   cells (p = 4.5×10^−5), while higher animal-to-animal variability and no
   significant change were observed in male cells ([335]Figure 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;
   [336]Figure 4h).

   Increased NETosis in neutrophils from older female mice may seem
   contradictory with previous reports of decreased NETosis with
   organismal aging^[337]13, [338]14. To note, previous studies evaluated
   NETosis after in vitro TNFα priming^[339]13 or in conditions known to
   mimic TNFα priming^[340]14, [341]71, and did not evaluate sex as a
   biological variable. TNFα is a pro-inflammatory cytokine and a potent
   inducer of NETosis^[342]13. To evaluate how TNFα-priming may impact
   NETosis across groups, we also performed NETosis experiments on cells
   pre-treated with 10ng/mL TNFα ([343]Figure 4i). Although TNFα-priming
   yielded noisier data, NETosis inducibility was blunted with organismal
   age in female cells stimulated with TNFα ([344]Figure 4i). In contrast,
   TNFα-primed male neutrophils trended towards decreased inducibility of
   NETosis with organismal age (p = 0.092; [345]Figure 4i). In TNFα-primed
   conditions, young male vs. female neutrophils no longer showed
   differences in NETosis inducibility, but old female neutrophils still
   showed the highest overall NETosis inducibility ([346]Figure 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 ([347]Figure 5a–[348]c,
   [349]Figure 6a–[350]c, [351]Extended Data Figure 9a–[352]f).
   Machine-learning can help (i) to determine whether features contain
   information predictive of a state of interest (e.g. sex-dimorphic or
   age-regulated expression), and (ii) to identify predictive features of
   such states. We used seven algorithms to learn models discriminating
   between (i) up- and down-regulated expression with age, and (ii) male-
   or female-biased gene expression (i.e. conditional inference Trees
   [cTree], Linear Discriminant Analysis [LDA], neural networks [NNET],
   Logistic Regression [LogReg], random forests [RF], support vector
   machines [SVM], and gradient boosting [GBM]). The models were trained
   with features related to (i) genomic sequence (e.g. GC-richness), (ii)
   chromatin accessibility (e.g. changes in ATAC-seq promoter
   accessibility), and (iii) TF target genes (i.e. same as used with
   GSEA).

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

   Figure 5:
   [353]Open in a new tab

   (a) Scheme of our machine-learning pipeline for up vs. 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 [354]Extended
   Data Figure 9a–[355]c and [356]Supplementary Table S5A). (c) Top 20
   most important features from Random Forest models and Gradient Boosting
   Machine 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. cTree: conditional
   inference tree, LDA: linear discriminant analysis, SVM: support vector
   machine, NNET: neural network, LogReg: logistic regression, GBM:
   gradient boosting machine, RF: random forest.

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

   Figure 6:
   [357]Open in a new tab

   (a) Scheme of our machine-learning pipeline for male vs. 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 [358]Extended
   Data Figure 9d–[359]f and [360]Supplementary Table S5C). (c) Top 20
   most important features from Random Forest models and Gradient Boosting
   Machine 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. cTree: conditional
   inference tree, LDA: linear discriminant analysis, SVM: support vector
   machine, NNET: neural network, LogReg: logistic regression, GBM:
   gradient boosting machine, RF: random forest.

   We first evaluated our ability to discriminate between genes
   upregulated vs. downregulated during aging in neutrophils ([361]Figure
   5a). Our models assigned genes to the correct transcriptional change
   with age >64% of the time on testing data ([362]Figure 5b;
   [363]Supplementary Table S5A), significantly above random (50%). High
   accuracy suggests that clear biological features differentiate between
   genes that tend to be up- vs. 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 ([364]Figure 5c; [365]Supplementary
   Table S5B). Consistent with the tight link between chromatin and gene
   expression^[366]72, the average promoter ATAC-seq signal (i.e.
   describing how “open” the promoter is), as well as the fold difference
   in ATAC-seq signal between young and old neutrophils (i.e. age-related
   changes in openness) were among the top predictors for age-related
   changes ([367]Figure 5c; [368]Supplementary Table S5B). Key predictors
   also included CG/CpG content in the promoter and gene sequences
   ([369]Figure 5c), which is consistent with our previous work
   identifying promoter CpG content as a top predictor of age-related gene
   expression in mouse tissues^[370]18. Cytosines in CpG configuration are
   the primary targets of DNA methylation^[371]73. CpG DNA methylation
   [DNAme], catalyzed by DNA-methyltransferases (i.e. writers), can be
   removed by TET family proteins (i.e. erasers) and is interpreted
   through recognition by methyl-CpG binding proteins (i.e.
   readers)^[372]73. 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%; [373]Supplementary Table S1A). The
   presence of putative age-related changes in DNAme is consistent with
   the notion that DNAme patterns can serve as “aging clocks”^[374]74.

   Key predictors of age-related expression changes also included whether
   a gene was a target of specific TFs: Stat5a, Mtf2, Sirt6, Foxo1, etc.
   ([375]Figure 5c; [376]Supplementary Table S5B). 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 (e.g. Stat5a, Mtf2, Nod2; [377]Extended Data
   Figure 9g; [378]Supplementary Table S1A). Thus, it is possible that the
   activity of TF predictors without RNA-level regulation may occur
   through post-translational regulation. Interestingly, Stat5a mediates
   the effects of GM-CSF on granulocyte differentiation^[379]75, and is
   crucial for mature neutrophils survival^[380]76. Although the role of
   Mtf2 in neutrophils has not been studied, Mtf2 interacts with Jarid2
   (another top predictor) to promote repressive H3K27
   methylation^[381]77, which mediates recruitment to unmethylated
   CpGs^[382]78. Nod2 is a pattern recognition receptor that recognizes
   muramyl dipeptide, the minimal common motif of bacterial
   peptidoglycans^[383]79. Although Nod2 is not a TF, it has profound
   transcriptional effects on gene expression, with target genes including
   cytokines and chemokines^[384]80, [385]81. Indeed, Nod2 is an important
   contributor of neutrophil-mediated innate immunity^[386]81. Sirt6 is a
   histone deacetylase tightly linked to the regulation of mammalian aging
   and longevity^[387]82, [388]83. Sirt6 can limit inflammation by
   deacetylating the promoters of NF-kB targets^[389]84 and reducing
   cytokine production^[390]85. Consistently, Sirt6 knockout mice show
   liver inflammation with neutrophil infiltration^[391]86. 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
   ([392]Figure 6a). Importantly, as an additional feature, we also
   included the chromosomal location of a gene (i.e. sex-chromosomes vs.
   autosomes). Our models assigned genes to the correct sex-based
   expression bias >70% of the time on held-out testing data ([393]Figure
   6b; [394]Supplementary Table S5C), 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 ([395]Figure 6c; [396]Supplementary Table S5D). Importantly,
   being located on a sex chromosome (X or Y) was a poor predictor of
   sex-dimorphic gene expression (166^th combined importance rank;
   [397]Supplementary Table S5D), consistent with the fact that genetic
   context was not a crucial driver ([398]Extended Data Figure 1a).
   Similar to our age-related models, the average promoter ATAC-seq signal
   and the fold change between females vs. males in ATAC-seq signal were
   also top predictors for sex-biased gene expression ([399]Figure 6c;
   [400]Supplementary Table S5D). The GC/CpG content of promoter and genes
   sequences also showed strong predictive power, consistent with the
   sex-dimorphic expression of DNAme regulators (i.e. male-biased
   expression for Dnmt1 and Tet3; [401]Supplementary Table S1B).

   Finally, we asked which TFs were most predictive of sex-dimorphic gene
   expression in neutrophils ([402]Figure 6c; [403]Supplementary Table
   S5D). Being a known target of Foxm1, Sirt6, Mtf2 or Mybl2 was highly
   predictive of male vs. 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 (e.g. Foxm1, Irf8; [404]Extended Data Figure 9h;
   [405]Supplementary Table S1B). Although its role hasn’t been explored
   in neutrophils, Foxm1 regulates the expression of cell cycle-related
   genes^[406]87 and is highly expressed in late committed neutrophil
   precursors^[407]40. The transcription factor Mybl2 also regulates cell
   cycle-related genes^[408]88, and can modulate differentiation and cell
   fate decision of myeloid progenitors^[409]89. 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 ([410]Figure 3, [411]Extended Data Figure 7).
   Being a target of Sirt6 was also predictive of sex-dimorphic gene
   expression, in line with sex-dimorphic phenotypes of Sirt6 knock-out
   and overexpression mice^[412]82, [413]83. In addition, consistent with
   sex-dimorphism in lipid storage ([414]Figure 3f), Sirt6 is tightly
   linked to lipid metabolism^[415]90, and can limit lipid droplet
   accumulation in foam cells^[416]91. In addition, Sirt6 activity can be
   directly regulated by FFAs (e.g. myristic, oleic and linoleic
   acid)^[417]92. Consistently, female neutrophils have higher levels of
   FFAs, including FFA(14:0) and FFA(18:2) (i.e. myristic and linoleic
   acid), and may thus have higher basal Sirt6 activity levels
   ([418]Supplementary Table S1E,[419]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 vs. male neutrophils

   Intriguingly, when analyzing the top genes with the largest
   fold-difference between male vs. female neutrophils, we noticed that
   genes encoding primary neutrophil granule components showed top
   male-biased expression (e.g. Elane, Mpo, Prtn3 and Ctsg;
   [420]Supplementary Table S1B), suggesting potential sex differences in
   granule loading. Thus, we investigated potential differences in granule
   biogenesis ([421]Figure 7a–[422]b, [423]Extended Data Figure
   10a–[424]d). Neutrophils store a number of proteases and antimicrobial
   proteins in (i) primary/azurophilic granules, (ii) secondary/specific
   granules, and (iii) tertiary/gelatinase granules. To unbiasedly assess
   whether neutrophil granule biology is sex-dimorphic, we compiled a list
   of granule components from the literature ([425]Supplementary Table
   S6). Consistent with our first observation, GSEA revealed a robust
   trend for increased RNA expression of primary granule-related genes in
   male vs. female neutrophils, regardless of age ([426]Figure 7a–[427]b).
   In contrast, no significant enrichments were observed for genes
   encoding components of secondary or tertiary granules (FDR > 5%;
   [428]Extended Data Figure 10a–[429]d). Importantly, we confirmed that
   protein levels of neutrophil elastase Elane were higher in young males
   vs. young females by Western blot (p = 3.9×10^−3 in two-sided
   Wilcoxon’s test between young females and male samples), although this
   sex-difference disappeared in older animals ([430]Extended Data Figure
   5d,[431]e).

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

   [432]Figure 7:
   [433]Open in a new tab

   (a) Heatmap of normalized gene expression for primary (azurophilic)
   granule-related gene expression in our RNA-seq dataset. Also see
   [434]Extended Data Figure 10. (b) GSEA analysis 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. The data is derived from 3 independent cohorts of
   young male vs. female mice (n = 5 per sex and time point for 1, 3 and
   6h PBS injection; n = 8 per sex and time point for 1 and 3h LPS, and n
   = 9 per sex and time point for 6h LPS injection). For simplicity’s
   sake, all PBS-injected animals are reported as “0h of LPS exposure” in
   (d) and are replotted with time-based color-coding in (e). Significance
   in non-parametric two-sided Wilcoxon rank-sum test. For boxplots in
   panels (d-e), the center line represents the sample median, the box
   limits consist of the 25^th and 75^th percentiles, the whiskers span
   1.5x the interquartile range.

   Neutrophil degranulation is a key aspect of the response of neutrophils
   to pathogens^[435]93, and mice without a functional copy of Elane have
   increased mortality upon sepsis^[436]94. However, excessive levels of
   circulating neutrophil elastase can be deleterious, amplifying septic
   shock^[437]95. 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 LPS (a major component of
   gram-negative bacteria cell wall), for 1–6 hours, to mimic a septic
   state^[438]96 ([439]Figure 7c). Interestingly, males showed
   significantly higher levels of serum neutrophil elastase at 3 and 6
   hours post LPS exposure ([440]Figure 7d), and in the PBS-injected
   controls ([441]Figure 7d,[442]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^[443]97. Since 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
   (i.e. higher levels of DDX3Y associated with males). Consistent with
   our mouse RNA-seq and serum ELISA, we found significant correlation
   between DDX3Y levels (i.e. likelihood of the sample coming from a male
   donor) and expression of primary granule proteins ELANE, MPO and CTSG,
   although not of PRTN3 (Figure [444]Extended Data Figure 10e–[445]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. To our knowledge, this dataset is the largest multi-omic
   dataset for the study of neutrophils and is one of the rare cases to
   include both sexes and organismal aging, rather than focusing only on
   one sex or age group. Using this resource, we observed that (i)
   neutrophils are extremely sex-dimorphic, and that (ii) organismal age
   leads to large scale transcriptomic and epigenomic remodeling of
   neutrophils. Interestingly, when analyzed separately, males had >
   10-fold 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^[446]98. 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^[447]8. 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^[448]19,
   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^[449]10.

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 (i) genes induced vs. repressed during aging, and (ii) 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 (e.g. Foxo1, Sirt6). We also identified putative
   mediators of sex-biased gene expression of neutrophils (e.g. Foxm1,
   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^[450]99. This is
   especially relevant to aging, with large sex-differences in baseline
   lifespan and in response to pro-longevity interventions^[451]58.
   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^[452]21, [453]22.
   Consistently, accumulating evidence has shown that immune responses
   differ between sexes^[454]100, [455]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^[456]102. Interestingly, NETs may drive aspects of
   inflammatory disease, lead to tissue damage^[457]8, [458]30, and
   contribute to severe COVID-19^[459]103. In addition, neutrophil
   dysfunction has been linked to the pathogenesis of a number of chronic
   diseases (e.g. macular degeneration^[460]104, stroke^[461]105,
   obstructive pulmonary disease^[462]106, atherosclerosis^[463]107,
   cancer^[464]108, etc.). 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 protocols #20770, #20804 and #21004.

   For ‘omics’ analyses, male and female C57BL/6Nia mice (4–5 and 20–24
   months 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 months old animals) were obtained
   from Jackson Laboratories. Both Charles Rivers and Jackson Laboratories
   have specific-pathogen-free (SPF) facilities. Animals were acclimated
   at the SPF animal facility at USC for 2–4 weeks before any processing.
   The facility is on a 12-hour light/dark cycle, and animal housing rooms
   are maintained at 72°F and 30% humidity. All animals were euthanized
   between 8–11am for ‘omics’ experiments, flow cytometry and NETosis
   assays. For the sepsis model experiments, animals were injected in the
   morning, and euthanized between 4–6pm. 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[2] asphyxiation followed by cervical dislocation.

Isolation of primary neutrophils from mouse bone marrow

   The long bones of each mouse were harvested 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 the
   bone marrow from cleaned bones was collected into clean tubes^[465]110.
   Red blood cells from the marrow were removed using Red Blood Cell Lysis
   buffer (Miltenyi Biotec #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 #130–110-916) to retain only single cells for
   downstream processing. Neutrophils were isolated from other bone marrow
   cells using magnetic-assisted cell sorting (Miltenyi Biotec kit
   #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 #130–092-575), MACS-purified
   neutrophils were then stained using APC-Ly6G (Invitrogen #17–9668-80)
   and Vioblue-Cd11b (Miltenyi Biotec #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 V7
   software. Purity of MACS-isolated bone marrow neutrophils in young and
   aged male and female animals are reported in [466]Extended Data Figure
   3. Raw cytometry data was deposited on Figshare
   (doi.org/[467]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^[468]111, [469]112,
   [470]113 (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 vs. “aged” subsets of mature
   neutrophils were defined respectively as
   CD11b^+Ly6G^+c-Kit^−Cxcr4^−Cxcr2^+CD62L^+ vs.
   CD11b^+Ly6G^+c-Kit^−Cxcr4^−Cxcr2^+CD62L^− cells.

   Prior to analyzing compositional heterogeneity on MACS-purified
   neutrophils, to account for spill over from different lasers,
   compensation was performed using appropriate compensation beads
   (Miltenyi Biotec #130–104-693 for Miltenyi antibodies; ThermoFisher #
   01–3333-42 for others). The compensation file was saved and loaded
   prior to 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^5 cells were washed once by adding 1mL
   of PBS/EDTA 0.1% BSA (i.e. MACS resuspension buffer) in a 5mL
   polystyrene round bottom tube (Falcon #352054), then centrifuged 300g
   for 10 minutes at 4°C. After an Fc-blocking step (Miltenyi Biotec
   #130–092-575), MACS-purified neutrophils were then stained with
   antibodies for 20 minutes 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 was analyzed using FlowLogic V7. Gating was determined
   using fluorescence-minus-one [FMO] controls for each color used in the
   experiment to ensure that positive populations were solely associated
   with the antibody for that specific marker ([471]Extended Data Figure
   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 labelling by this
   antibody occurs normally on bone marrow cells prior to purification
   ([472]Extended Data Figure 4b). Both sets of raw cytometry data,
   including FMO controls and compensation files, were deposited to
   Figshare (doi.org/[473]10.6084/m9.figshare.14043929.v1,
   doi.org/[474]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 (i.e. G0/G1, S and G2/M). For each biological sample,
   5×10^5 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 minutes at room temperature in labelling
   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 was analyzed using FlowLogic
   V7. Raw cytometry data was deposited to Figshare
   (doi.org/[475]10.6084/m9.figshare.14043926.v1).

RNA purification and RNA-seq library preparation

   For RNA isolation, frozen cell pellets were resuspended in 1mL of
   Trizol reagent (Thermo-Fisher), 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 (RIN). 500ng 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 75bp reads were generated on the Illumina
   NextSeq500 platform at the SC^2 Core at CHLA (original cohort) or
   paired-end 150bp reads were generated on the Illumina HiSeq-Xten
   platform at the Novogene Corporation (USA) (replicate cohort).

RNA-seq analysis pipeline

   To best mimic the first RNA-seq cohort, paired end 150bp reads from the
   replicate cohort of neutrophil RNA-seq were hard-trimmed to 75bp using
   Fastx Trimmer. Both sets of paired-end reads were processed using
   Trimgalore 0.4.4 ([476]github.com/FelixKrueger/TrimGalore) (i) to
   retain only high-quality bases with phred score > 15, and (ii) to
   eliminate biases due to priming by hard clipping the first 6 bases of
   each read. Only pairs with both reads retaining a length of > 45bp
   after trimming were retained for further processing. Trimmed reads were
   mapped to the mm10 genome reference using STAR 2.5.0a^[477]114. Read
   counts were assigned to genes from the UCSC mm10 reference using
   subread1.5.3^[478]115 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/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^[479]116. R package ‘sva’ v3.34 was used to estimate surrogate
   variables, and the removeBatchEffect function from ‘limma’ v3.42.2 was
   used to regress out the effects of surrogate variables and
   RNA-integrity differences (RIN scores) from raw read counts. The
   ‘DESeq2’ R package (DESeq2 1.26.0) was used for further processing of
   the RNA-seq data in R^[480]117. Importantly, sex was encoded as a
   categorical variable (female vs. male), and age was encoded as a
   continuous numerical variable. Genes with a false discovery rate < 5%
   were considered statistically significant and are reported in
   [481]Supplementary Table S1A–[482]B.

Dimensionality reduction for exploratory data analysis

   To perform Multidimensional Scaling (MDS) analysis^[483]118, 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 (i) a small replicate
   RNA-seq cohort of primary bone marrow neutrophils with all four
   conditions, (ii) a dataset from mouse spleen neutrophils from young
   adult female vs. males from the ImmGen consortium (6 weeks-old animals;
   GEO Datasets [484]GSE124829)^[485]22 and (iii) a dataset from human
   blood neutrophils (ages unknown; GEO Datasets [486]GSE145231)^[487]24
   ([488]Extended Data Figure 2a,[489]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’ 2.42.1 package (accessed
   2020–12-18), and DESeq2 normalized fold-changes were then compared.

Putative Protein-Protein Interaction [PPI] network analysis

   Genes with significant regulation according to our DEseq2 analysis (FDR
   < 5%) were used as input for network analysis with NetworkAnalyst
   3.0^[490]119. For age-related gene regulation, significant genes and
   DEseq2 calculated log[2](fold-change) were used as input for a single
   network analysis. For sex-dimorphic gene regulation, the 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^[491]109, 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 Gene Set Enrichment Analysis (GSEA) paradigm^[492]120
   through the ‘phenotest’ 1.34.0 R package. Gene Ontology term annotation
   were obtained from ENSEMBL through Biomart (Ensembl 99; downloaded on
   2020–04-10), 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 v7.0 (e.g. Reactome,
   KEGG)^[493]120, [494]121 and the Harmonizome Database (e.g. ChEA,
   JASPAR, ENCODE, GEO TF targets; accessed 2020–03-25)^[495]122. To
   improve target coverage, we summarized putative transcription factor
   [TF] targets, including FOXO TF targets compiled in our previous
   study^[496]18, so as to have a reduced, unique, and non-redundant 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 10 most significant pathways with negative NES and top 10 most
   significant pathways with positive NES are shown in figures if more
   than that pass the FDR < 5% significance threshold. All significant
   gene sets are reported in [497]Supplementary Tables S3 (aging) and
   [498]S4 (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, 1x Halt Protease and
   Phosphatase Inhibitor Cocktail [Life Technologies # 78442]) and
   incubated on ice for 30 minutes. Samples were sonicated (Fisher
   Scientific # FB120; 60% power, 30 sec ON/ 120 sec OFF, 4 cycles) and
   centrifuged at 16,000g for 20 minutes at 4°C. Supernatant was collected
   and boiled in 4x Laemmli Sample Buffer (Biorad # 1610747). Proteins
   were separated on 10% SDS-PAGE gels and wet-transferred onto PVDF
   membranes (GE Healthcare # 10600021). To assess loading and transfer
   efficiency, membranes were stained using the Ponceau S solution (Sigma
   P7170) for 2 minutes at room temperature. Membranes were blocked using
   5% milk (Carnation) in 1xTBST buffer (TBS [Alfa Aesar J62938],
   supplemented with 0.05% Tween-20 [Bio-rad #161–0781]) for 1 hour 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 hours at 4°C.
   After three-5 minute washes using 1xTBST 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 hour at room
   temperature. For visualization, membranes were treated with
   chemiluminescence solution according to manufacturer’s manual (Biorad #
   1705062 and Thermo Scientific # 34580) and images were taken using
   Azure Biosystems c280.

   For band intensity quantification, ImageJ (version 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
   ([499]https://doi.org/10.6084/m9.figshare.14154665.v1) and as a source
   data file for [500]Extended Data Figure 5.

Chemicals for LC-MS

   LC-MS-grade solvents and mobile phase modifiers were obtained from
   Fisher Scientific (water, acetonitrile, methanol, methyl tert-butyl
   ether) 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, methyl tert-butyl
   ether (MTBE) and water^[501]123. Briefly, 300μL of methanol spiked-in
   with 54 deuterated internal standards provided with the Lipidyzer
   platform (SCIEX, cat #5040156, LPISTDKIT-101) was added to the cell
   pellet, samples were vigorously vortexed for 20 seconds and sonicated
   in a water bath 3 times for 30 seconds on ice. Lipids were solubilized
   by adding 1000μL of MTBE and incubated under agitation for 1h 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 4 times of ice-cold 1:1:1 isopropanol:acetonitrile:water
   spiked in with 17 labeled internal standards and incubating for 2 hours
   at −20°C. The supernatant was dried down to completion under nitrogen
   and re-suspended in 100μL of 1:1 MeOH: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^[502]124. 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 10mM 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 (i)
   injecting 6 and 12-pool samples to equilibrate the LC-MS system prior
   to run the sequence for RPLC and HILIC, respectively, (ii) sample
   randomization for metabolite extraction and data acquisition, and (iii)
   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 v2.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^[503]125.

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 HMDB ([504]http://www.hmdb.ca/), MoNA
   ([505]http://mona.fiehnlab.ucdavis.edu/) and MassBank
   ([506]http://www.massbank.jp/) using the R package ‘MetID’
   (v0.2.0)^[507]126 (Level 2). Briefly, metabolic feature tables from
   Progenesis QI were matched to fragmentation spectra with a m/z and a
   retention time window of ±15 ppm and ±30 s (HILIC) and ± 20 s (RPLC),
   respectively. When multiple MS/MS spectra match 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 above 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 (DMS) interface (SCIEX) and a high flow LC-30AD
   solvent delivery unit (Shimazdu). A full description of the method is
   available elsewhere^[508]123. Briefly, the lipid molecular species were
   identified and quantified using multiple reaction monitoring (MRM) and
   positive/negative switching. Two acquisition methods were employed
   covering 10 lipid classes; method 1 had SelexION voltages turned on,
   while method 2 had SelexION voltages turned off. Lipidyzer data were
   reported by the Lipidomics Workflow Manager (LWM) v1.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 (i) tuning the DMS
   compensation voltages using a set of lipid standards (SCIEX #5040141)
   after each cleaning, more than 24 hours of idling or 3 days of
   consecutive use, (ii) performing a quick system suitability test (QSST)
   (SCIEX #50407) before each batch to ensure acceptable limit of
   detection for each lipid class, iii) sample randomization for lipid
   extraction and data acquisition, and iv) triplicate injection of lipids
   extracted from a reference plasma sample (SCIEX #4386703) at the
   beginning of the batch.

Data pre-processing.

   Lipids detected in less 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^[509]125.

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 #23225) to account
   for differential starting material quantity. Then, Variance Stabilizing
   Normalization was applied to the data using ‘limma’ 3.42.2, as
   recommended by previous studies^[510]127, [511]128. Differential
   analysis for metabolomic or lipidomic features was performed using
   ‘limma’ in R. Features with a false discovery rate (FDR) < 5% were
   considered statistically significant.

   For analysis of the lipidomics data by lipid class ([512]Extended Data
   Figure 1; [513]Supplementary Table S1F), 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 a false discovery rate (FDR) < 5% were considered
   statistically significant.

Functional enrichment analysis for metabolomics

   Since only 1 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’ 2.0.4^[514]129 to perform Phenotype Set Enrichment
   analysis (PSEA)^[515]130 of KEGG metabolic pathways. Using the
   ‘mummichog’ method^[516]131, 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 [517]Supplementary Table S4F.

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^[518]132. We used the IMPaLA v12 web interface
   ([519]http://impala.molgen.mpg.de/), which evaluates overrepresentation
   of genes and metabolites across 5055 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 (server access on 06–30-2020). All
   pathways with an overall combined FDR < 5% are reported in
   [520]Supplementary Table S4H.

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^[521]133. Lipid features with
   FDR < 5% with a male or female bias were uploaded to the server (access
   on 03–24-2020). For ease of reading, only the 10 most significant
   pathways with male and with female bias are shown in [522]Figure 3f.
   All significant LION terms are reported in [523]Supplementary Table
   S4G.

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^[524]134. 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 150bp reads at the
   Novogene Corporation (USA).

ATAC-seq preprocessing pipeline

   Paired-end ATAC-seq reads were adapter-trimmed using
   NGmerge0.2^[525]135, 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 (2.2.6)^[526]136. After alignment, PCR
   duplicates were removed using the ‘rmdup’ function of samtools (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 (2.20.3) or
   samtools (1.5). We used the peak finding function of HOMER (4.11) to
   identify ATAC-seq accessible regions^[527]137.

Differential accessibility analysis with ATAC-seq

   We extracted a normalized read count matrix from our the downsampled
   bam files over merged HOMER regions using R package ‘Diffbind’ 2.14
   ^[528]138. We used this extracted matrix for downstream analyses.
   Because of sample-to-sample variability, we used surrogate variable
   analysis ^[529]116 to estimate and correct for experimental noise,
   similar to the RNA-seq analysis. R package ‘sva’ 3.34 was used to
   estimate surrogate variables, and the ‘removeBatchEffect’ function from
   ‘limma’ 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 (1.26.0) was used for further processing of ATAC-seq data.
   Regions with a false discovery rate < 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’ 4.0.4
   tool^[530]139, a tool 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 to default parameters. Results were exported as ‘tsv’ and
   processed in R to filter significant regions at FDR < 5%. Filtered
   results are reported in [531]Supplementary Tables S3F and [532]S4I.

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 v0.3.4 software^[533]68. Since 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 the authors of the software ^[534]68. We extracted
   similar metrics from the NucleoATAC output to what was described in
   previous studies using this software to investigate chromatin
   architecture^[535]140. All comparisons between groups were tested for
   statistical significance using two-tailed Wilcoxon rank-sum test
   (non-parametric 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^[536]137. To allow for
   whole library normalization, reads mapping over transcript were also
   counted using the same HOMER script. Read counts were imported into R
   to estimate differential repeat expression using the ‘DESeq2’ (1.26.0)
   R package^[537]117.

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^[538]141 and a 96-well plate plate-reader protocol^[539]142,
   both with neutrophils in suspension. Briefly, MACS-purified neutrophils
   were resuspended in a concentration of 2×10^6 cells/mL in neutrophil
   culture medium (RPMI 1640 without phenol red [Hyclone], supplemented
   with 1% Pennicilin/Streptomycin [Corning] and 0.1% BSA [Akron
   Biotechnology]). We then used 1×10^6 cells, aliquoted into sterile
   microcentrifuge tubes, one per starting biological sample. SYTOX green
   (ThermoFisher Scientific # S7020) was added to each sample to a final
   concentration of 200nM. Neutrophil media supplemented with DMSO
   [Vehicle] or 50nM Phorbol 12-myristate 13-acetate [PMA] (Sigma # P1565)
   was added to each tube. Tubes were slowly inverted three times to mix.
   Then, 2×10^5 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[2] at 37°C for 2 hours. The fraction of
   cells positive for SYTOX Green in each well was quantified using the
   MACSQuant10 and analyzed using Flowlogic V7. To account for differences
   in basal levels of NETosis across samples, NETosis was expressed as
   induction of NETosis: (median of PMA technical quadruplicates) /
   (median of DMSO technical quadruplicates). The raw cytometry dataset
   was deposited to Figshare
   (doi.org/[540]10.6084/m9.figshare.14043923.v1).

TNFα-primed NETosis analysis

   To compare our results with previously published results on NETosis in
   aged neutrophils^[541]13, [542]14 [which were obtained with direct TNFα
   priming^[543]13 or in conditions known to mimic TNFα priming^[544]14
   (i.e. thioglycolate elicitation of peritoneal neutrophils^[545]71)], we
   also performed NETosis analysis in primed conditions, using a slightly
   modified protocol from above to include a priming step. Specifically,
   the MACS-purified neutrophils cells were resuspended in neutrophil
   culture medium and first pre-warmed after isolation in a humidified
   incubator with 5% CO2 at 37°C for 15 minutes. Then, 2×10^6 cells were
   aliquoted into a sterile microcentrifuge tube, and primed with 10ng/mL
   mouse TNFα (PeProTech # 315–01A) in a humidified incubator with 5% CO2
   at 37°C for 15 minutes. Cells were then further aliquoted, supplemented
   with SYTOX green, DMSO/PMA, incubated and analyzed as above. The raw
   cytometry dataset was deposited to Figshare
   (doi.org/[546]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 to 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 (i) the average promoter
   accessibility by ATAC-seq in FPKM across all 4 conditions (promoters
   were defined as −500bp;+150bp with regards to annotated transcriptional
   start sites in mm10 build according to HOMER), as well as (ii) the
   average log[2](Fold Change) in accessibility between old and young
   neutrophils (age-regulation models) or the average log[2](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^[547]122
   and transcriptional targets of FOXO transcription factors from GEO
   experiments compiled in our previous study^[548]18. To reduce
   extraneous features, we engineered a TF target feature following these
   steps: (i) TF targets were summarized from all putative sources, (ii)
   only TFs with evidence of expression in primary neutrophils according
   to RNA-seq were retained as features, and (iii) 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 to 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 the 3,653 genes with significant age
   regulation, and 1,636 of the 1,734 genes with significant sex-dimorphic
   regulation. Finally, to determine whether location of the gene on
   autosomes vs. 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 2 questions: (i)
   models for age-regulated gene expression, and (ii) models for
   sex-dimorphic gene expression. For (i), age-regulation machine-learning
   models were built to learn whether up- or down-regulated genes could be
   discriminated using genomic features and predict potential master
   regulators. For (ii), 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 7 classification
   algorithms as implemented through R package ‘caret’ (caret 6.0–86).
   Auxiliary R packages were used with ‘caret’ to implement neural
   networks (NNET; ‘nnet’ 7.3–13), random forests (RF; ‘randomForest’
   4.6–14), gradient boosting (GBM; ‘gbm’ 2.1.5), radial kernel support
   vector machines (SVM; kernlab 0.9–29), sparse Linear Discriminant
   Analysis (LDA; ‘sparseLDA’ 0.1–9), conditional inference Trees (cTree;
   ‘party’ 1.3–4), Regularized Logistic Regression (LogReg; ‘LiblineaR’
   2.10–8). ‘Caret’ was allowed to optimize final model parameters on the
   training data using 10-fold cross validation. Accuracies,
   sensitivities, specificities and AUC (using package ‘pROC’ 1.16.2) for
   all trained classifiers were estimated using a test set of randomly
   held out 1/3 of the data (not used in the training phase) obtained
   using the ‘createDataPartition’ function. Feature importance estimation
   was only done using tree-based RF and GBM methods, since 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 DIA proteomics results from the supplemental material of a
   recent study^[549]97. We used normalized expression values from the
   article’s supplemental material and used the reported “donor-specific
   Spectronaut protein expression values” for further analysis. To exclude
   potential confounds linked to disease, we only used data from the 68
   healthy donors [HD] and excluded the data from diseased patients. Since
   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 (i.e. ELANE, MPO, PRTN3 and
   CTSG). Significance of the Spearman rank correlation is reported.

Mouse sepsis model through intra-peritoneal LPS injection

   Young adult (3–4 months) C57BL/6J mice were exposed to LPS, a
   pathogen-associated molecular pattern (PAMP) found in the cell-wall of
   Gram-negative bacteria, to elicit a sepsis-like response^[550]96.
   Briefly, mice were injected intra-peritoneally with 2μg of LPS per gram
   of body weight, using a sterile LPS stock (Sigma #L5293) diluted in
   PBS, and control animals were injected with sterile PBS^[551]96.
   Animals were monitored hourly for up to 6 hours after injection for the
   occurrence and severity of endotoxemia^[552]96. We did not observe
   gross differences in the state of female vs. 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
   euthanasia by CO[2] asphyxiation and cervical dislocation at 1, 3 and 6
   hours post injection, blood was collected from the heart. A total of 80
   animals across 3 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 hour at room
   temperature in MiniCollect Serum SeparatorTubes (Greiner Bio-One). The
   serum was then separated from the clot using centrifugation at 2,000g
   for 10 minutes, 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. 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 version 3.6.0–3.6.3.
   For all ‘omic’ analyses, results were corrected for multiple testing by
   the use of False Discovery Rate [FDR] significance correction. All
   statistical analyses were performed as specified in the figure legends
   and in the corresponding methods section. We used non-parametric
   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^[553]143. As one exception, RNA-seq of
   bone marrow neutrophils was independently replicated on a smaller
   cohort, since neutrophils are known to be RNA-poor, which could add
   more than usual biological variability ([554]Extended Data Figure 2).
   With the exception of a test of the heterogeneity flow cytometry panel
   on bone marrow cells (used to validate c-kit staining; [555]Extended
   Data Figure 4b; [556]Supplementary Table S2B), all other experiments
   were performed on animals from at least 2 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 is possible since groups are biologically determined. For
   LPS injection experiments, animals were randomly allocated to the PBS
   or LPS group for a euthanasia time point, on a cage-level basis (i.e.
   animals from the same cage randomly attributed to LPS vs. PBS
   injections; for each cohort, co-housed animals were used for the same
   time point; cages attributed to a specific euthanasia time point
   randomly). This process was performed independently for each of the 3
   cohorts for LPS-induced endotoxemia. Blinding was not relevant to the
   studies conducted here, as the data are collected by automated systems
   (i.e. sequencing, mass-spectrometry, flow cytometry) or in numeric form
   (i.e. ImageJ quantification of Western Blot bands, ELISA absorbance
   values), which cannot be influenced by subjectivity from the
   experimenter.

Data Availability

   Sequencing data has been submitted to the Sequence Read Archive (SRA)
   accessible through BioProject PRJNA630663. Untargeted metabolomics data
   was uploaded to metabolomics workbench DataTrack ID 2089. The
   lipidomics data is available as [557]Supplemental Table S7 and was
   deposited in Figshare
   ([558]https://doi.org/10.6084/m9.figshare.14524278). Raw flow cytometry
   data was deposited to Figshare
   (doi.org/[559]10.6084/m9.figshare.14043938.v1,
   doi.org/[560]10.6084/m9.figshare.14043932.v1,
   doi.org/[561]10.6084/m9.figshare.14043929.v1,
   doi.org/[562]10.6084/m9.figshare.14043926.v1, and
   doi.org/[563]10.6084/m9.figshare.14043923.v1). Western Blot raw and
   cropped images have been deposited to Figshare
   ([564]https://doi.org/10.6084/m9.figshare.14154665).

   We reanalyzed publicly available neutrophil RNA-seq data from GEO
   Datasets GSE1248296 (6 weeks-old mouse spleen neutrophil samples), and
   [565]GSE145231 (human blood neutrophil samples), and human neutrophil
   DIA proteomics data
   ([566]https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6442368/bin/141289_1
   _supp_251888_pjqfps.xlsx). Functional gene annotations were obtained
   from ENSEMBL Biomart ([567]https://www.ensembl.org/biomart/martview/),
   the Molecular Signature Database
   ([568]http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) and the
   Harmonizome ([569]https://maayanlab.cloud/Harmonizome/download).
   Mass-spectrometry peaks were identified by matching experimental m/z
   and fragmentation spectra to publicly available databases including
   HMDB ([570]http://www.hmdb.ca/), MoNA
   ([571]http://mona.fiehnlab.ucdavis.edu/) and MassBank
   ([572]http://www.massbank.jp/).

Code Availability

   The analytical code is available on the Benayoun lab Github repository
   ([573]https://github.com/BenayounLaboratory/Neutrophil_Omics_2020). All
   R code was run using R version 3.6.0–3.6.3.

Extended Data

Extended Data Figure 1: A multi-omic analysis of primary mouse bone marrow
neutrophils during aging and with respect to sex (continued).

   Extended Data Figure 1:
   [574]Open in a new tab

   (a) Table of significant “omic” features as a function of organismal
   age and sex in our datasets based on DEseq2 (RNA-seq) or limma
   (Metabolomics and lipidomics), at FDR < 5%. For the analysis of
   sex-dimorphism in gene expression, the number of significant genes
   located on autosomes (i.e. not on chromosomes X or Y) is also reported.
   (b-c) Barplot of DESeq2-normalized log2 counts for Xist (b) and Ddx3y
   (c), showing the expected pattern between male and female samples,
   which acts as a quality check for our dataset. (d) Heatmap of lipidomic
   changes summarized by lipid classes. Significance of difference as a
   function of age or sex was evaluated using a linear model, and the
   significance of the sex coefficient (as reported by the R ‘lm’
   function) is reported on the line of the heatmap. *: p < 0.05; ** : p <
   0.01; n.s. : p ≥ 0.05. See also [575]Supplementary Table S1F. (e)
   Circular genome plot of the positions of genes with significant
   sex-biased gene expression in neutrophils (FDR < 5%). (f) Heatmap of
   significant age-regulated genes (DESeq2 FDR < 5%). (g-h) Correlation
   plot of age-related gene expression change according to DESeq2 in
   female vs. male neutrophils from RNA-seq, showing genes with (g)
   significant and concordant age-regulation in both sexes at FDR 5%, or
   (h) genes with divergent age-regulation between sexes (FDR < 5% in one
   sex, and FDR >15% in the other). Spearman Rank correlation (Rho), and
   significance of this correlation are reported in (g). Importantly, age
   was inputted into the DEseq2 model as a continuous numerical variable
   (expressed in months), which yields fold change values per month of
   life.

Extended Data Figure 2: Comparison to other neutrophil RNA-seq datasets.

   Extended Data Figure 2:
   [576]Open in a new tab

   (a) DESeq2 normalized log[2] fold changes as a function of biological
   sex for differentially expressed mouse genes from our original cohort
   (FDR <5%) (and their human orthologs) across datasets [see
   [577]Methods]. The original mouse bone marrow neutrophil cohort data
   are plotted for comparison as the leftmost panel. On the right are
   plotted corresponding log[2] fold change values from our own smaller
   replication cohort of mouse bone marrow neutrophils (n = 3 per group),
   from mouse spleen neutrophils from ImmGen^[578]23 and a human blood
   neutrophil cohort^[579]25. P-values were calculated using a two-sided
   Wilcoxon test between female-biased and male-biased genes from our
   original cohort across new datasets, to test the robustness of such
   differences between our original cohort and new datasets. (b) DESeq2
   normalized log[2] fold changes per month during aging for
   differentially age-regulated mouse genes from our original cohort (FDR
   <5%). The original mouse bone marrow neutrophil cohort data are plotted
   for comparison as the leftmost panel. On the right are plotted
   corresponding log[2] fold change values from our own smaller
   replication cohort of mouse bone marrow neutrophils (n = 3 per group).
   P-values were calculated using a two-sided Wilcoxon test between
   upregulated vs. downregulated genes, to test the robustness of such
   differences between our original cohort log[2] fold changes and a new
   dataset. The center line represents the sample median, the box limits
   consist of the 25th and 75th percentiles, the whiskers span 1.5x the
   interquartile range, and each RNA-seq sample is represented by a point
   on the graph.

Extended Data Figure 3: Bone marrow neutrophil purity is not impacted by
animal sex or organismal age.

   Extended Data Figure 3:
   [580]Open in a new tab

   (a) Representative flow cytometry gating strategy of bone-marrow
   neutrophils purified using MACS from a young female mouse. Neutrophils
   are expected to be double positive for Cd11b and Ly6G. (b) Neutrophil
   purity from 2 independent cohorts of aging male vs. female mice (n = 8
   for old females, and n = 10 per group for all other groups), determined
   as the Cd11b+ Ly6G+ population. (c) Proportion of neutrophils among
   bone marrow nucleated cells according to the Hemavet 950FS from 2
   independent cohorts of aging male vs. female mice (n = 9 for old
   females, and n = 10 per group for all other groups). Statistical
   analysis for (b) and (c) derived from a linear modeling analysis with
   sex and age as covariates, similar to our “omic” analysis models, with
   significance of coefficients as reported by the R ‘lm’ function. The
   center line represents the sample median, the box limits consist of the
   25th and 75th percentiles, the whiskers span 1.5x the interquartile
   range, and each RNA-seq sample is represented by a point on the graph.

Extended Data Figure 4: Analysis of MACS-purified bone marrow neutrophil
heterogeneity.

   Extended Data Figure 4:
   [581]Open in a new tab

   (a) Representative flow cytometry gating strategy of bone-marrow
   neutrophils purified using MACS from a young female mouse. Populations
   were defined according to previously published markers^[582]136,
   [583]137, [584]138. Specifically, among live cell singlets,
   pro-neutrophils were defined as c-Kit^+Ly6g^−Cd81^+, with subpopulation
   of pro-neutrophil 1 defined as c-Kit^+Ly6G^−Cd81^+CD11b^−CD106^− and
   pro-neutrophils 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 vs. “aged”
   subsets of mature neutrophils were defined respectively as
   CD11b^+Ly6G^+c-Kit^−Cxcr4^−Cxcr2^+CD62L^+ vs.
   CD11b^+Ly6G^+c-Kit^−Cxcr4^−Cxcr2^+CD62L^− cells. The gate from which
   cells are plotted is indicated above each flow cytometry plot. (b)
   Representative flow cytometry gating strategy of RBC-lysed bone marrow
   cells from a young female mouse. The gate from which cells are plotted
   is indicated above each flow cytometry plot. This analysis was run as a
   control for the presence/abundance of cKit^+ cells, as MACS-purified
   cells had extremely low numbers of these cells. (c) Amounts of immature
   (left) vs. mature (right) neutrophil among all MACS-purified cells from
   2 independent cohorts of aging male vs. female mice (n = 10 per group).
   (d) Proportion of “fresh” vs. “aged” neutrophils among mature
   neutrophils from 2 independent cohorts of aging male vs. female mice (n
   = 10 per group). Statistical analysis for (c) and (d) derived from a
   linear modeling analysis with sex and age as covariates, similar to our
   “omic” analysis models, with significance of coefficients as reported
   by the R ‘lm’ function. The center line represents the sample median,
   the box limits consist of the 25th and 75th percentiles, the whiskers
   span 1.5x the interquartile range, and each RNA-seq sample is
   represented by a point on the graph. Also see [585]Supplementary Table
   S2 for quantifications of all populations.

Extended Data Figure 5: Regulated pathways in bone marrow neutrophils during
aging reveals downregulation of chromatin-related pathways (continued).

   Extended Data Figure 5:
   [586]Open in a new tab

   (a-b) Top enriched gene sets from Gene Ontology (a) and KEGG (b) using
   GSEA for differential RNA expression. Only the top 10 most up- and top
   10 most downregulated gene sets are plotted for readability. Full lists
   and statistics available in [587]Supplementary Table S2. Shown pathways
   with FDR < 5%. (c) Boxplot of Padi4 transcriptional levels from
   RNA-seq. Significance: DESeq2 FDR values for regulation as a function
   of aging or biological sex. The center line represents the sample
   median, the box limits consist of the 25^th and 75^th percentiles, the
   whiskers span 1.5x the interquartile range, and each RNA-seq sample is
   represented by a point on the graph (n = 4 RNA-seq samples per group).
   (d) Representative Western Blot images for primary bone marrow
   neutrophils for Elane, Padi4 and β-actin (loading control), as well as
   Ponceau S staining of PVDF blotting membrane. The membrane was cut
   along predicted molecular weights, and each strip was probed for each
   specific protein in that range. (e) Boxplot of quantification of
   Western Blot for Padi4 (left) and Elane (right) from primary bone
   marrow neutrophils. Data from 2 independent cohorts of aging male vs.
   female mice (n = 9 per group). Two concentrations of the same sample
   were loaded side by side for each biological sample. Statistical
   analysis derived from a linear modeling analysis with sex and age as
   covariates, similar to our “omic” analysis models, with significance of
   coefficients as reported by the R ‘lm’ function. The center line
   represents the sample median, the box limits consist of the 25th and
   75th percentiles, the whiskers span 1.5x the interquartile range, and
   each RNA-seq sample is represented by a point on the graph. All
   uncropped western blot images and quantification data is been made
   available on Figshare
   ([588]https://doi.org/10.6084/m9.figshare.14154665.v1). NES: Normalized
   Enrichment Score (for GSEA analysis). FDR: False Discovery Rate.

Extended Data Figure 6: The distribution of cell cycle phases of
MACS-purified bone marrow neutrophil is not impacted by animal sex or
organismal age.

   Extended Data Figure 6:
   [589]Open in a new tab

   (a) Representative flow cytometry gating strategy of bone-marrow
   neutrophils purified using MACS from a young male mouse stained using
   Propidium Iodide to assess cell-cycle phase based on DNA content.
   Neutrophils are expected to be largely post-mitotic. (b) Boxplots of
   MACS-purified bone marrow neutrophil cell cycle distribution (G0/G1, S
   or G2/M) from 2 independent cohorts of aging male vs. female mice (n =
   9 for old males, and n = 10 per group for all other groups), determined
   by DNA content. The overwhelming majority of samples have > 90% cells
   in G0/G1, and there are no trends associated to organismal age or
   biological sex. Statistical analysis for (b) derived from a linear
   modeling analysis with sex and age as covariates, similar to our “omic”
   analysis models, with significance of coefficients as reported by the R
   ‘lm’ function. The center line represents the sample median, the box
   limits consist of the 25th and 75th percentiles, the whiskers span 1.5x
   the interquartile range, and each RNA-seq sample is represented by a
   point on the graph.

Extended Data Figure 7: Sex-dimorphic pathways in bone marrow neutrophils
reveal differential regulation of chromatin-related pathways (continued).

   [590]Extended Data Figure 7:
   [591]Open in a new tab

   (a-b) Top enriched gene sets from Gene Ontology (a) and KEGG (b) using
   GSEA for differential RNA expression as a function of sex. Only the top
   10 most up- and top 10 most downregulated gene sets are plotted for
   readability. Full lists and statistics available in [592]Supplementary
   Table S2. All shown pathways and genes such as FDR < 5%. (c)
   MetaboAnalyst PSEA scatterplot for KEGG pathways from metabolomics.
   NES: Normalized Enrichment Score (for GSEA analysis). FDR: False
   Discovery Rate. (d) Boxplot of transcriptional levels from our RNA-seq
   for Pnpla2/Atgl, Atg7 and Lpin1, which are sex-dimorphic lipid
   metabolism-related genes. The significance is derived from DESeq2, and
   we report the corresponding FDR values for regulation as a function of
   aging or biological sex. The center line represents the sample median,
   the box limits consist of the 25th and 75th percentiles, the whiskers
   span 1.5x the interquartile range, and each RNA-seq sample is
   represented by a point on the graph (n = 4 RNA-seq samples per group).
   (e) Summary scheme of discussed lipid usage in female vs. male
   neutrophils based on lipidomic and transcriptomic data. Catalyzed
   reactions are derived from Wikipathway WP3901 (Lipid droplet
   metabolism) and ^[593]37.

Extended Data Figure 8: ATAC-seq analysis reveals age- and sex-related
differences in the chromatin architecture of bone marrow neutrophils
(continued).

   [594]Extended Data Figure 8:
   [595]Open in a new tab

   (a) NucleoATAC v-plots. The light green box overlay reveals a v-plot
   region where males have higher signal than females, regardless of age,
   suggesting differences in nucleosomal architecture.

   (b-c) Barplot of DESeq2-normalized log2 counts at ATAC-seq peaks
   associated to Xist (b) or situated on the Y chromosome (c), showing the
   expected pattern between male and female samples. (d) Multidimensional
   Scaling analysis results of chromatin accessibility from ATAC-seq. (e)
   Table of significantly differentially accessible ATAC-seq peaks as a
   function of organismal age and sex. The number of significant peaks
   located on autosomes (i.e. not on chromosomes X or Y) is also reported.
   (f-g) Heatmap of accessibility at peaks with significant age-regulated
   (f) or sex-dimorphic (g) change in accessibility (FDR < 5%). (h)
   Circular genome plot of the genomic positions of peaks with significant
   sex-bias in neutrophil ATAC-seq (FDR < 5%). (i) Boxplot of nucleosome
   occupancy as calculated by NucleoATAC for all called nucleosomes across
   groups. Also see [596]Figure 4d. (j) Boxplot of median ATAC-seq
   fragment length at ATAC-seq peaks across experimental groups. The
   horizontal red line in panels (i,j) shows the median value in
   neutrophils from young females for ease of comparison. Significance in
   non-parametric two-sided Wilcoxon rank-sum tests are reported for
   panels (i,j). Black p-values represent differences between male vs.
   female neutrophils in each age group; pink (blue) p-values represent
   age-related differences in female (male) neutrophils. For boxplots in
   panels (i,j), the center line represents the sample median, the box
   limits consist of the 25th and 75th percentiles, the whiskers span 1.5x
   the interquartile range, 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. (k-l) Heatmap of
   repetitive elements with significant age-regulated (k) or sex-dimorphic
   (l) change in transcription from RNA-seq (FDR < 5%).

Extended Data Figure 9: Machine-learning analysis reveals that age-regulated
and sex-dimorphic gene expression can be predicted by genomic and epigenomic
features.

   Extended Data Figure 9:
   [597]Open in a new tab

   (a-c) Age-related machine-learning model performance metrics: balanced
   classification accuracy over the 10 cross-validation folds during model
   training (a), balanced classification accuracy on held-out testing data
   (b), and Area Under the Curve [AUC] on held-out testing data (c). Also
   see [598]Supplementary Table S5A. (d-f) Sex-dimorphism machine-learning
   model performance metrics: balanced classification accuracy over the 10
   cross-validation folds during model training (d), balanced
   classification accuracy on held-out testing data (e), and Area Under
   the Curve [AUC] on held-out testing data (f). Also see
   [599]Supplementary Table S5C. (g-h) Heatmap of RNA-seq expression level
   for top TF predictors in age-regulated (g) or sex-dimorphic (h) gene
   expression models with FDR < 10%. DESeq2 FDR is reported on each line,
   with # : FDR < 0.1; * : FDR < 0.05; ** : FDR < 0.01; *** : FDR < 0.001.
   FDR: False Discovery Rate. For boxplots in panels (a,d), the center
   line represents the sample median, the box limits consist of the 25^th
   and 75^th percentiles, the whiskers span 1.5x the interquartile range.

Extended Data Figure 10: Male neutrophils express higher levels of primary
granule genes but not secondary and tertiary granule genes.

   Extended Data Figure 10:
   [600]Open in a new tab

   (a-c) Heatmap of normalized gene expression for primary (a), secondary
   (b) and tertiary (c) granule-related gene expression in our RNA-seq
   dataset. The estimated False Discovery Rate [FDR] using GSEA for each
   gene set is reported. (d) Heatmap of normalized gene expression for
   X-linked (Kdm5c, Kdm6a, Xist) and Y-linked (Uty, Ddx3y, Eif2s3y) genes
   in our RNA-seq dataset. Panel (a) is identical to [601]Figure 7a to
   help direct comparison between sex-bias expression patterns of
   different granule types as well as canonical sex-biased genes. (e-h)
   Scatterplots of protein levels of DDX3Y (as a proxy for likely of a
   sample being derived from a male donor), compared to protein levels for
   primary granule-related proteins ELANE (e), MPO (f), PRTN3 (g) and CTSG
   (h) in neutrophils from healthy human donors. Normalized data was
   obtained from the supplementary material of ^[602]117. Estimates of
   Spearman Rank correlation (Rho) and significance of correlation are
   reported for each correlation on the scatterplot.

Supplementary Material

   Supplementary information inventory
   [603]NIHMS1725952-supplement-Supplementary_information_inventory.docx^
   (18.7KB, docx)
   Combined supplementary table file
   [604]NIHMS1725952-supplement-Combined_supplementary_table_file.xlsx^
   (1.2MB, xlsx)

Acknowledgements