A molecular atlas of human granulopoiesis

A molecular atlas of human granulopoiesis | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Resource A molecular atlas of human granulopoiesis Christoph Klein, Sebastian Hesse, Jiadong Mao, Armin Hadziahmetovic, and 8 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6184761/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted You are reading this latest preprint version Abstract We developed a comprehensive molecular atlas of human granulopoiesis, encompassing the transcriptome (mRNA), proteome, and non-coding RNA (miRNA, lncRNA) of seven defined stages of neutrophil maturation. Two waves of mRNA transcription activity were observed. Whereas the first wave, extending from myeloblasts to metamyelocytes, showed correlated mRNA transcripts and protein translation, the second wave was characterized by uncoupling of transcription and translation. Integrated analysis revealed intricate dynamics of transcription and translation parallelization in granule subset proteins, ribosomes, and mitochondria, respectively. We identified the hsa-miR-106a-5p family and hsa-miR-125a-5p as potential repressors of myeloid differentiation in myeloblasts, and hsa-miR-193a-5p as a repressor of the SWI/SNF complex. LncRNA transcripts displayed remarkable stage specificity. Whereas NEAT1 and MAILR were expressed only from the band stage onwards, ITGB2-AS1 was expressed exclusively in metamyelocytes. Novel lncRNAs with high expression were found in the terminal S-stage and PMN (ENSG00000288819, ENSG00000203644). A comprehensive analysis of transcription factor signatures, including activity on target genes, revealed previously unrecognized factors as well as a yet underappreciated transcriptional role of lactotransferrin (LTF) in early maturation stages. Characterization of single-cell-clusters via transcriptional profiles of histologically defined maturation stages enabled stage-specific cell labeling. To ensure the availability of our data and the reproducibility of our analysis, our complete datasets, analysis codes and online tool for exploration are available on www.granulopoiesis.com. Biological sciences/Immunology/Innate immune cells/Granulocytes/Neutrophils Biological sciences/Immunology/Haematopoiesis/Myelopoiesis Neutrophil granulocytes Hematopoiesis Multi Omics data integration Transcription factors Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Objectives Creation of a stage-specific multiomics atlas of human granulopoiesis Integrative analysis of transcriptome (mRNA), proteome and non-coding RNAs (miRNA, lncRNA) Quantification of miRNA-mRNA interactions Identification of functional transcription factor signatures and global omics shift Key results A comprehensive molecular atlas of human granulopoiesis with 7 developmental stages at the mRNA, protein and non-coding RNA level, presented as FAIR-accessible data and code on www.granulopoiesis.com Revised understanding of granule protein synthesis, as transcription and translation are initiated in parallel for all granule subsets and terminated at the metamyelocyte stage Unraveling of hsa-miR-106a-5p and hsa-miR-125a-5p as potential inhibitors of myeloid differentiation in lineage committed neutrophil progenitors Identification of novel stage specific long non-coding RNAs (ENSG00000288819, ENSG00000203644) in terminal maturation stages Identification of LTF (lactotransferrin) as a central transcription factor in early granulopoiesis, a signature of 32 TFs active throughout each maturation step and an exclusive signature of 131 TFs active in the final progression from the bone marrow S-stage to peripheral blood neutrophils Stage-specific molecular markers on mRNA, miRNA, lncRNA and protein level Novelties Unprecedented resolution of the molecular landscapes of human granulopoiesis Integrative analysis workflows to define and interpret feature trajectories of cell-development Identification non-coding RNA governance factors on miRNA and lncRNA level Annotation of single cell profiles in accordance with histology based cell identities Introduction Neutrophil granulocytes are short-lived effector cells of the innate immune system that make up 50–70% of leukocytes in human peripheral blood. Their functions range from defense against pathogens to tissue remodeling and they play a role in cancer defense as well as carcinogenesis (Burn et al., 2021 ). Sophisticated attack mechanisms enable them to rapidly destroy microorganisms or tissue. Their eponymous granules are typically divided into primary (azurophilic), secondary (specific), tertiary, and ficolin-1-rich subsets that are thought to harbor mostly non-overlapping sets of proteins. Other effector functions include the generation of reactive oxygen species (ROS), the formation of neutrophil extracellular traps (NETs) and phagocytosis (Faurschou and Borregaard, 2003 ). Human diseases affecting neutrophils range from deficits in one of their specialized defense mechanisms, such as deficiency of specific granules or reactive oxygen species generation, to defects in granulopoiesis that include over 30 known monogenic causes of congenital neutropenia (Donadieu and Bellanné-Chantelot, 2022 ). Granulopoiesis, the differentiation of neutrophil granulocytes in the bone marrow, is thought to follow a linear developmental process that is divided into major cellular states of metamorphosis defined by specific states of transcriptome, proteome and non-coding RNA landscapes (Larsen et al., 2012 ). They differentiate from lineage committed stem-cell like myeloblasts (MB) to promyelocytes (PM), myelocytes (MC), metamyelocytes (MM), band (B) and segmented (S) nuclear forms. At the terminal maturation stage, cells emigrate from the bone marrow into the peripheral blood, referred to as polymorphonuclear neutrophils (PMN) (Elghetany et al., 2004 ). Granulopoiesis is a high-throughput differentiation process that transforms readily dividing stem cells into terminal effector cells with an enormous steady-state production rate of 10^11 cells per day (Strydom and Rankin, 2013 ). Previously, we and others showed that the highly specialized function of the neutrophil is mirrored in its proteome configuration, as less than 10 proteins make up over 50% of the total protein content (Grabowski et al., 2019 ) (Sollberger et al., 2024 ). Several groups have recently provided comprehensive multiomic analyses of neutrophil granulocytes subpopulations in diseases such as Covid (Schulte-Schrepping et al., 2020 ; Wu et al., 2021 ), or cancer (Ng et al., 2024 ). To better understand the molecular landscape of developing neutrophil granulocytes, we set out to characterize their stage-specific molecular signatures on the transcriptome, proteome and non-coding RNA level (Fig. 1a). Results Transcription-translation dynamics of granule-subtypes shows parallel and time-delayed synthesis Based on defined precursor stages of granulopoiesis (Kwok et al. 2020), we developed a panel for flow cytometric cell sorting to purify bone marrow cells from four healthy donors plus one biological replicate measurement (Fig. 1a). Using stringent quality criteria (sFig. 3II), we quantified 8,432 mRNA transcripts (transcriptome), 3,156 proteins (proteome) and 283 mature miRNA transcripts (miRNome). For initial analysis, we considered each omics level individually to examine its coherence with the hypothesized developmental stage association from our FACS and microscopy analysis. A high-level analysis was based on principal component analysis (PCA) for dimensionality reduction in combination with Ward-D based Euclidean clustering (Ward, 1963 ). For each omics level, we observed clear separation of stages in a developmentally coherent manner (Fig. 2I-III a). Clustering revealed near perfect separability of the individual stages at the level of the transcriptome (Fig. 2Ia) and miRNome (Fig. 2IIIa), while the proteome showed three main developmental clusters: MB, PM with MC and MM, B and PMN (Fig. 2IIa). Thus, the maturation process of granulocytes might be characterized by highly dynamic and stage-specific RNA landscapes, while changes in the proteome overlap considerably between stages, indicating a gradual developmental process. Since we observed a strong coherence of the samples with the assumed developmental stages, we next analyzed the changes between each successive developmental step with differential expression (DE) analysis. We detected a large number of differentially expressed features (1000–2000 combined) at the transcriptome level in each developmental stage, while significant changes at the proteome level were only detected during the initial differentiation (MB->PM->MC->MM). After the metamyelocyte stage, the proteome showed only minimal changes in protein abundance levels (0–9 DE proteins). At the miRNA level, the most pronounced shift was observed from the myeloblasts to the myelocytes, with a steady state during final maturation (Fig. 2IIIa). As clustering and DE indicated high stringency of the obtained omics signals for each maturation stage, we applied enrichment analysis to determine the biological themes of the up-regulated (Fig. 2I-III d) and down-regulated (Fig. 2I-III e) features. At the transcriptome level, all known granule subfractions (azurophilic = primary, specific = secondary, tertiary and ficolin-1 granules) were upregulated in myeloblasts and promyelocytes, followed by a gradual and stage-specific downregulation of granule subfraction features. At the proteome level, myeloblasts begin translation of the azurophil and ficolin-1 granules, while promyelocytes show translation of all granule fraction features. Myelocytes stop the translation of the primary granule proteins and simultaneously upregulate the translation of the specific, tertiary and ficolin-1-rich granules. In addition, we found an upregulation of ribosomal proteins from myeloblasts to promyelocytes, coinciding with granule protein synthesis, followed by a downregulation of ribosomal transcripts and loss of ribosomal proteins from the metamyelocyte stage onwards, suggesting termination of protein synthesis. To analyze how features of defined compartments (granule subfractions, plasma membrane, ribosomes and mitochondria) behave on mRNA and protein level throughout development, we plotted the summarized medians along the longitudinal development path (Fig. 2 VI). Granule proteins showed stage specific patterns that exhibited primarily coordinated transcriptome-proteome increases (primary and ficolin-1 granules from MB to PM) as well as translational time lags (secondary gr. in PM to MC) followed by transcriptional downregulation while protein levels remained stable. Ribosomes and mitochondria showed an initial increase of protein intensities with stable mRNA levels in promyelocytes, followed by continuous loss of features on transcriptome and proteome level. Thus we found evidence for developmental stage-specific signals in each of the three omics layers. Interestingly, transcription of all granule subtypes appeared to be initiated at the myeloblast stage, whereas translation initiation and transcription termination showed evidence for stage specificity. However, we found ample evidence for coordinated synthesis of granule subtype features. While significant proteome changes were only observed up to the metamyelocyte stage, transcriptome changes showed a second peak in the final stages of maturation while the proteome remained stable. Longitudinal analysis reveals developmental feature trajectories and allows unified pathway enrichment on multiomic plane Since differential expression indicated complex patterns, we conducted a longitudinal analysis of feature expression throughout development, here referred to as “trajectories”. To this end, we combined WGCNA with a multivariate time course analysis (described in Methods). We found that the majority of traits in each omics level showed a unidirectional expression pattern (only up or down). In the transcriptome, 1629 features ( ≙ 19% of the total, Fig. 3 Ia) showed an increasing trajectory and 1970 features (23%, Fig. 2Ib) showed a decreasing trajectory. At the proteome level, 483 (15%, Fig. 3 IIa) proteins increased and 504 (16%, Fig. 3 IIb) proteins decreased during human granulopoiesis. In the miRNome, 54 transcripts were identified that increased (19%, Fig. 3 IIIa) and 68 that decreased (24%, Fig. 3 IIIb). For a smaller subset of features, we observed bidirectional trajectories comprising a total of 674 transcripts (UpDown: 361, DownUp: 313, Fig. 2 Ic), 590 proteins (UpDown: 329, DownUp: 261, Fig. 3 IIc) and 33 miRNA transcripts (UpDown: 21, DownUp: 12, Fig. 3 IIIc). A minority of features showed multidirectional trajectories (sFig. 1 I-II). To derive the common biological themes of features in each trajectory, we filtered for overlapping terms after enriching for pathways in the Reactome database (Fig. 3 IV) (Milacic et al., 2024 ). Interestingly, we observed a negative correlation between miRNAs and features of transcriptome and proteome, consistent with the notion that miRNAs exert mainly down-regulatory effects (Shang et al., 2023). For example, pathways with up-regulated features at the transcriptome and proteome levels and down-regulated interacting miRNAs were VEGF, Ras, mTOR, MAPK and HIF1, while we also observed this pattern for features of cell-cell adhesion of leukocytes and extracellular trap formation (Fig. 3 IV, upper part). Conversely, we observed down-regulated mRNA and protein features with up-regulated interfering miRNAs for ribosomes, the translation apparatus and mitochondrial membrane proteins, as well as the signal recognition particle (SRP), nonsense-mediated decay (NMD) and the regulatory network for the expression of SLIT and ROBO proteins (Fig. 3 IV, lower part). Some of the pathways and cell compartments showed complex patterns of interacting miRNAs with bidirectionally regulated transcripts and bimodal protein trajectories (like ribosomes, translation, SRP complex). We compared the genes involved in acute myeloid leukemia (AML) and severe congenital neutropenia (SCN) with the features of the transcriptome and proteome (Fig. 3V), respectively. Intriguingly, we found that the majority of genes responsible for severe congenital neutropenia showed a downward trajectory of expression levels at both omics levels, demonstrating a reduction in expression during granulopoiesis. In summary, we found that the majority of features at all omics levels exhibit continuously increasing or decreasing trajectories with biologically coherent themes at the transcriptome, proteome and miRNA levels. Specific miRNA transcripts govern development stage associated processes in human granulopoiesis In an effort to determine hierarchical effects of miRNA transcripts on the granulopoiesis transcriptome, we probed the miRNAs of trajectory clusters for binding-site enrichment in differentially expressed RNA transcripts (using enrichMIR (Soutschek et al., 2022 ), see methods). Coupled with enrichment analysis of the mRNA transcripts, we could link miRNAs with the cellular processes they may control. For miRNAs with increasing trajectory (Fig. 3 IIIa), we selected DE mRNA transcripts between the two endpoints, myeloblasts (MB) and peripheral blood neutrophils (PMN). Here we found hsa-miR-193a-5p (Fig. 4 Ia) to interact with a total of 86 mRNA features (Fig. 4 Ib) that enriched for the ATPase complex and the SWI/SNF complex (Fig. 4I c,d). For the 68 miRNAs that decrease during granulopoiesis (Fig. 3 IIIb), we again enriched DE features that alternate between the two endpoints, but in reverse order, from PMN to MB to maintain the expected fold-change direction for the enrichMIR algorithm. We found two miRNA clusters, hsa-miR-19a/b-3p and hsa-miR-106a/17/20a/b-5p (Fig. 4 IIa), that enriched for a large group of mRNA transcripts in myeloblasts (Fig. 4 IIb). Interestingly, transcripts associated with the latter group (miR-106a) showed enrichment for myeloid cell differentiation, cell growth and histone modification (Fig. 4 II c,d). For miRNA features with up-down trajectory (Fig. 3 IIIc), we analyzed the features between band stage and myeloblasts (Fig. 4 III). We found that hsa-mir-454-3p binds a large group of mRNA transcripts associated with negative regulation of transcription by RNA polymerase II. Conversely, miRNAs with a down-up pattern were enriched for DE features between metamyelocytes and myeloblasts with terms similar to those of miRNAs with a decreasing trajectory, namely myeloid cell differentiation and cell growth (Fig. 4 IV). In the miRNA trajectory group showing a late up-down expression pattern (sFig. 1 II), we found that hsa-miR-1343-3p interacts with DE features between metamyelocytes and band stage granulocytes, that were enriched for translation and ribosomes (sFig. 1 III a-d). Transcription factor signatures in human granulopoiesis reveal a role of LTF in early differentiation steps In the final maturation step, from segmented nucleated neutrophils in the bone marrow to peripheral blood neutrophils, GO enrichment identified the downregulation of a group of 22 transcription factors. This finding encouraged us to further investigate transcription factor signatures in granulopoiesis. To this end we employed ChIP-X Enrichment Analysis 3 (Keenan et al., 2019 ), that enriches TF targets in DE transcriptome sets (details in Methods). First, we focused on the 10% most enriched TF signatures in DE transcripts of each successive development step (Fig. 5 Ia). In addition to previously described TFs orchestrating granulopoiesis (e.g. NFE2, CEBPs, MXD1), we identified LTF (lactotransferrin) most strongly induced and active in the two early stage progressions (MB to PM and PM to MC, Fig. 5 Ia). We extended the analysis with a comparison of all significantly enriched TF signatures between all developmental stages (Fig. 5 Ib). We identified a set of 32 TFs expressed and active at each step of granulopoiesis, in particular NFE2, MXD1, ZNF467, MXD3 and SP1. Another 40 TF signatures were shared between all initial development steps (MB to B), with prominent appearances of CEBPB and CEBPD. When we compared mature neutrophils in the bone marrow (S) with terminal neutrophils in peripheral blood (PMN), we identified an exclusive signal of 131 active TFs. Most of them, e.g. THAP11 (Thanatos-associated protein 11), associated with cell growth suppression (Zhu et al, 2009 ), have not previously been described in granulopoiesis. Stage specificity of long non-coding RNA expression lncRNAs emerge as critical factors in acute myeloid leukemia, driving proliferation, inhibiting differentiation and conferring treatment resistance (Connerty and Lock, 2023 ). We queried our comprehensive dataset to assess the abundance of lncRNA during defined stages of differentiation. Figure 5 IIa shows a heatmap of the top40 variably expressed lncRNA transcripts during granulopoiesis. They show highly specific expression patterns with SNHG (Small nucleolar RNA host genes) family members prominently expressed from MB to MC, NEAT1 and MAILR expression only from Band stage on, and various currently uncharacterized transcripts. Of special note, ENSG00000288836 was expressed specifically in myelocytes and metamyelocytes, and spans a genetic region annotated with variant associations that include neutrophil count ( The Association to Function Knowledge Portal , 2024), sFig. 2 IIa). Two poorly studied lncRNA were strongly expressed in segmented and polymorphonuclear stages, ENSG00000288819 and ENSG00000203644, respectively. While the targets of these ncRNA remain to be defined, ITGB2-AS1, a known antisense RNA transcript to Integrin subunit beta2, was found exclusively expressed from the metamyelocyte stage onwards. As our global transcriptome included all transcripts in unison, we separated protein-coding mRNA and long non-coding transcripts and plotted the number of differentially expressed features (Fig. 5 II b,c). As the protein coding mRNA constituted 88% of transcripts, the distribution of significantly up- and down-regulated features did not change in comparison to the total transcriptome (Fig. 1 Ic). However, the DE lncRNAs showed a markedly different pattern with an over-representation of continuously up-regulated transcripts from MB to PM, MC and MM and a transcriptionally silent phase from MM to B (mirroring the total transcriptome results) with few down-regulated features. A pronounced difference was observed from the B stage on, where 90 lncRNAs are significantly upregulated in transition to S. The final maturation, from S to PMN, was the only step where the majority of lncRNA transcripts were down-regulated. Neutrophil differentiation displays omics-shifts in cell compartments In previous studies, PMNs were found to have a uniform proteome configuration, with as few as 7 proteins accounting for 50% of the total protein content (Grabowski et al., 2019 ), a result that was recently confirmed (Sollberger et al., 2024 ). We hypothesized that the proteome remodeling might be detectable in our dataset as a stepwise process during differentiation. Our results were consistent with previous findings of a highly condensed proteome in peripheral blood neutrophils and a coherent proteome shift during granulopoiesis, with myeloblasts showing the most heterogeneous proteome (56 proteins with a cumulative abundance of 50% in MB) that concentrates to the observed homogenous state with each maturation step (13 proteins with an abundance of 50% in PMN, Fig. 5 IIIa). In the transcriptome, we found intriguing similarities and differences (Fig. 5 IIIb). In myeloblasts, the transcriptome configuration mirrored the proteome as it showed the most heterogeneous mRNA landscape (86 transcripts accounting for 50% of the cumulative abundance). Likewise, it showed reduced heterogeneity with progressing maturation until the band stage. Thereafter, a pronounced contrast to the proteome configuration was observed in the last two stages, segmented nuclei (S) and peripheral blood neutrophils (PMN). We noticed a reversal of transcriptome homogeneity with a sudden increase in heterogeneity of transcripts (over 50 transcripts accounting for 50% of the total transcriptome). As a follow up, we determined which fraction of the proteome (Fig. 5 IVa) and transcriptome (Fig. 5 IVb) is dedicated to distinct subcellular compartments. In the myeloblast proteome, we observed the highest summed up protein intensities annotated to nucleus, secretory vesicles and cytosol. As differentiation progressed, summed up protein intensities of granule proteins continuously increased to take up the majority of the proteome, while the intensities of ER and Golgi, mitochondria and nucleus continuously decreased. Interestingly, we found that primary granules made up the largest proportion in promyelocytes and decreased slowly but steadily as maturation progresses, while secretory vesicles represented the largest fraction from the metamyelocyte stage onwards. In the transcriptome, we again found similarities to the proteome during initial maturation but observed a striking difference in the last two stages, S and PMN. Mirroring the proteome, the transcriptome initially showed a continuous decrease in intensity of features of the ER and Golgi, cytosol, mitochondria and nucleus. The features of the granule subgroups increased with maturation and showed partial sub-stadial specificity by constituting the largest fractions: primary granule features in promyelocytes, ficolin-1 granules in myelocytes and secondary granule features in metamyelocytes. Up to the B stage, the shifts corresponded to those in the proteome. In the S and PMN stages, however, the transcriptome landscape reversed its granule focused specialization and returned to a heterogeneous state, similar to that of myeloblasts. Granulopoiesis transcriptome data enable refined clustering of neutrophil single-cell data We investigated whether our bulk-RNAseq data from histology defined maturation stages of granulopoiesis could lead to a more refined characterisation of single-cell (sc) neutrophil clusters. To this end, we applied a novel data integration method called “Φ-Space” (Mao, Deng and Lê Cao, 2024 ) to a dataset of single cell neutrophil sequences from COVID-19 patients ((Schulte-Schrepping et al., 2020 )). The authors identified distinctive clusters of neutrophils in severe COVID-19 patients and utilized marker genes from mouse granulopoiesis to categorize three crude maturation stages: pre-, pro- and mature neutrophils (Kwok et al., 2020). As an alternative, we used “Φ-Space” to integrate the sc based maturation-scores with our bulk-data of histology-defined maturation stages (analysis workflow in sFig. 2 IIIa). A heatmap based on these maturation scores (Fig. 5 Va) showed that the single cell clusters defined by Schulte-Schrepping et al. overwhelmingly matched the individual transcriptional profiles of our histologically defined maturation stages. Of note, we observed an overlapping profile match between the “ribosome high” cluster and our profiles of band-stage neutrophils and myeloblasts. We further visualized these maturation scores using Uniform Manifold Approximation and Projection (Fig. 5V b, UMAP, (McInnes et al., 2018 )). We re-clustered the single neutrophils based on their maturation scores which gave a better delineation of cell clusters compared to the original mapping (sFig. 2 IIIb). Discussion Here, we present a comprehensive molecular atlas of 7 defined developmental stages of human granulopoiesis on transcriptome, proteome, and miRNome levels. Our data characterize granulopoiesis as a finely coordinated process of extreme cellular differentiation from proliferating myeloblasts to postmitotic differentiated polymorphonuclear neutrophils. We observed a highly specialized proteome landscape that was reduced to a metabolic minimum with abundant effector proteins. While these changes occur in a stepwise manner, several neutrophil-specific proteins (primary and ficolin-1 granule components) are already detectable in myeloblasts, indicating that the developmental trajectory towards neutrophils may be initiated already upstream of myeloblasts (Kwok et al., 2020). Most changes in protein abundance occurred in early phases, whereas the proteome remained fairly stable at the MM, B, and PMN stage. By contrast, the granulopoietic transcriptome showed changes in each development step and was divided into two distinct phases. The first transcriptome phase, from MB to MM, coincides with a proteome landscape rich in ribosomes, and shows a strong correlation between mRNA and protein features. The second phase, after the MM stage, takes place in a static proteome with minimal ribosomal content and shows decoupling of mRNA and protein features. This dual wave pattern of transcripts in neutrophil precursor development appears to contradict the “targeting by time” hypothesis (Le Cabec et al., 1996), that proposed that granule synthesis may be a time-controlled, serial process for each granule subfraction individually. We find evidence of parallel processing from the promyelocyte stage onwards. Interestingly, there is a time-lag between transcription and translation for the secondary granule subfraction. The mechanisms explaining this apparent discordance remain to be discovered. In parallel with the up-regulation of granule proteins, differentiation is accompanied by loss of cell cycle machinery and DNA replication, with diminished cytoplasmic ribosomes and mitochondrial matrix. This reflects transition from replicative to non-replicative cells and coincides with morphological changes of neutrophils (Robinson et al., 1999). An unexpected finding was the loss of nonsense-mediated decay (NMD). These findings have implications to interpret proteomic imbalances in genetic diseases. For example, we previously found upregulation of members of the NMD machinery in the neutrophil proteome of patients with severe congenital neutropenia (Grabowski et al., 2019 ). Our current study highlights the need to distinguish between signatures of immaturity versus compensatory mechanisms related to disease states. Our analysis of non-coding RNAs (lncRNA and miRNA) highlights previously unrecognized features. Even though, due to technical reasons, we could not obtain representative miRNA samples from promyelocytes, the remaining stages revealed a characteristic pattern. Two independent strategies were used to integrate the miRNome with transcriptome and proteome. First, when we compared enrichment terms, we could infer pathways and processes that miRNAs may inhibit (Mathieu and Ruohola-Baker, 2013), including VEGF, Ras, mTOR, cell-cell adhesion, extracellular trap formation) or silence during maturation (SRP, NMD, ribosomes). Second, we determined miRNA binding site enrichment in DE transcripts and performed post-hoc enrichment analysis, allowing us to infer activity of specific miRNA transcripts on major biological themes. miRNAs with high expression in myeloblasts (miR-19a/b-3p, miR-106a-5p family, miR-125a-5p) showed binding-site enrichment for mRNAs involved in myeloid differentiation (NBEAL2, NEDD9, TET2, STAT3). As miR-19 was found to control NFkB expression (Gantier et al., 2012), we speculate that miR-19 may repress neutrophil differentiation in myeloblasts via inhibition of NFkB. miR-193a on the other hand, previously shown to increase during granulopoiesis and to down-regulate transcripts of the essential SWI/SNF-complex, may facilitate myeloid leukemogenesis in conjunction with LINC00152 and CDK9 (Witzel et al., 2017) (Zhang and Tao, 2019). Long non-coding (lnc) RNA biology may have specific effects on a wide variety of cellular processes, from cell structure to chromatin architecture and transcription (Mattick et al., 2023). We found that lncRNA transcripts were upregulated from early to late phases of neutrophil differentiation in the bone marrow. When mature neutrophils leave the bone marrow (transition from S to PMN stage) lncRNA species are down-regulated. One particular class of lncRNA are lncRNA-type small nucleolar RNA host genes (SNHG), known to play pivotal roles in cancer progression (Zheng et al., 2023). SNHGs are overexpressed in early granulopoiesis (MB to MC), in line with the notion that their prognostic value in AML may may stem from their activity in less differentiated cell states (Shi, Ding and Lu, 2020 ). Interestingly, a recent study found that the overexpression of ITGB2-TAS1 (ITGB2 antisense RNA 1) was associated with adverse prognosis in AML. In our study, ITGB2-AS1 was under-expressed in MB to MC, with increasing expression in the MM stage, an opposite expression pattern to the aforementioned SNHGs. Our dataset enabled us to conduct an in-depth analysis of transcription factor signatures in the developing granulopoiesis transcriptome. For example, Lactotransferrin (LTF) is a prominent constituent of secondary granules in PMN with antimicrobial and immunomodulatory functions (Legrand, 2016). Its isoform, delta-lactoferrin (ΔLf), acts as a transcription factor with antiproliferative properties (Mariller et al., 2012). A recent single-cell analysis in mice found a signature of LTF activity in early granulopoiesis (Grieshaber-Bouyer et al., 2021). Our study finds LTF to be transcriptionally active and highly expressed from myeloblasts to metamyelocytes, peaking in myeloblasts and promyelocytes. Of note, we also discovered a number of transcription factors in mature stages of neutrophil development that were not previously linked to neutrophil biology. The last granulocytic stage differentiation, from S to PMN, showed a large and exclusive transcription factor signature. The two most enriched, THAP11 and THAP5, belong to the thanatos-associated protein family and have not been described yet in granulopoiesis. THAP5 is described to cause cell cyle arrest (Balakrishnan et al., 2009 ) while THAP11 (formerly RONIN) is a epigenetic regulator and essential for pluripotency in embryonic stem cells (Dejosez et al., 2008 ) that downregulates cMYC (Zhu et al., 2009 ). Our study points to a general limitation in the interpretation of proteomic data using the ‘proteomic ruler’ approach to infer protein copy numbers from histone intensities (Grabowski et al., 2019 ; Sollberger et al., 2024 ). Conceptually, this is based on the assumption that the ratio of histones to DNA remains stable and can serve as a constant that allows the inference of protein copy numbers in proteome samples with known cell numbers (Wiśniewski et al., 2014 ). However, we found that histone expression at the proteome level was strongly induced during granulopoiesis at both the mRNA and protein levels, reaching a maximum in PMN. Accordingly, the ‘proteomic ruler' approach cannot easily be applied since it would introduce a bias in the interpretation of protein abundances. In summary, we here describe intricate processes of cellular specialization from proliferating myeloblasts to postmitotic differentiated polymorphonuclear neutrophils. Whereas previous studies have focused on either transcriptome (Grassi et al., 2018), proteome (Hoogendijk et al., 2019), or miRNome (Larsen et al., 2012 ) individually, joint analysis of multiomic data from matched samples remains challenging. Our website, www.granulopoiesis.com enables full access to raw and processed data, the complete analysis scripts and interactive data analysis for further use. Declarations Funding SH: Förderprogramm für Forschung und Lehre (FöFoLe, Project number 1052) RZ: SFB1123 Atherosclerosis CK: German Center for Child and Adolescent Health (DZKJ) Munich Site, Sonderforschungsbereich TRR332 "Neutrophils: origin, fate & function" References Allaire, J. (2012) ‘RStudio: integrated development environment for R’, Boston, MA , 770(394), pp. 165–171. Available at: https://www.academia.edu/download/49999670/IntR_-_Interactive_GUI_for_R20161031-11791-s0u02o.pdf#page=14. Aparicio-Puerta, E. et al. (2023) ‘miEAA 2023: updates, new functional microRNA sets and improved enrichment visualizations’, Nucleic acids research , 51(W1), pp. W319–W325. Available at: https://doi.org/10.1093/nar/gkad392. Balakrishnan, M.P. et al. (2009) ‘THAP5 is a human cardiac-specific inhibitor of cell cycle that is cleaved by the proapoptotic Omi/HtrA2 protease during cell death’, American journal of physiology. Heart and circulatory physiology , 297(2), pp. H643–53. Available at: https://doi.org/10.1152/ajpheart.00234.2009. Barrett, T. et al. (2013) ‘NCBI GEO: archive for functional genomics data sets--update’, Nucleic acids research , 41(Database issue), pp. D991–5. Available at: https://doi.org/10.1093/nar/gks1193. Bodein, A. et al. (2022) ‘timeOmics: an R package for longitudinal multi-omics data integration’, Bioinformatics , 38(2), pp. 577–579. Available at: https://doi.org/10.1093/bioinformatics/btab664. Brunson, J.C. (2020) ‘ggalluvial: Layered Grammar for Alluvial Plots’, Journal of open source software , 5(49). Available at: https://doi.org/10.21105/joss.02017. Burn, G.L. et al. (2021) ‘The neutrophil’, Immunity , 54(7), pp. 1377–1391. Available at: https://doi.org/10.1016/j.immuni.2021.06.006. Choi, J. et al. (2019) ‘Stemformatics: visualize and download curated stem cell data’, Nucleic acids research , 47(D1), pp. D841–D846. Available at: https://doi.org/10.1093/nar/gky1064. Connerty, P. and Lock, R.B. (2023) ‘The tip of the iceberg-The roles of long noncoding RNAs in acute myeloid leukemia’, Wiley interdisciplinary reviews. RNA , 14(6), p. e1796. Available at: https://doi.org/10.1002/wrna.1796. Cox, J. and Mann, M. (2008) ‘MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification’, Nature biotechnology , 26(12), pp. 1367–1372. Available at: https://doi.org/10.1038/nbt.1511. Dejosez, M. et al. (2008) ‘Ronin is essential for embryogenesis and the pluripotency of mouse embryonic stem cells’, Cell , 133(7), pp. 1162–1174. Available at: https://doi.org/10.1016/j.cell.2008.05.047. Deutsch, E.W. et al. (2023) ‘The ProteomeXchange consortium at 10 years: 2023 update’, Nucleic acids research , 51(D1), pp. D1539–D1548. Available at: https://doi.org/10.1093/nar/gkac1040. Dobin, A. et al. (2013) ‘STAR: ultrafast universal RNA-seq aligner’, Bioinformatics , 29(1), pp. 15–21. Available at: https://doi.org/10.1093/bioinformatics/bts635. Donadieu, J. and Bellanné-Chantelot, C. (2022) ‘Genetics of severe congenital neutropenia as a gateway to personalized therapy’, Hematology , 2022(1), pp. 658–665. Available at: https://doi.org/10.1182/hematology.2022000392. Elghetany, M.T. et al. (2004) ‘Flow cytometric study of neutrophilic granulopoiesis in normal bone marrow using an expanded panel of antibodies: correlation with morphologic assessments’, Journal of clinical laboratory analysis , 18(1), pp. 36–41. Available at: https://doi.org/10.1002/jcla.20001. Faurschou, M. and Borregaard, N. (2003) ‘Neutrophil granules and secretory vesicles in inflammation’, Microbes and infection / Institut Pasteur , 5(14), pp. 1317–1327. Available at: https://doi.org/10.1016/j.micinf.2003.09.008. Gene Ontology Consortium et al. (2023) ‘The Gene Ontology knowledgebase in 2023’, Genetics , 224(1). Available at: https://doi.org/10.1093/genetics/iyad031. Grabowski, P. et al. (2019) ‘Proteome analysis of human neutrophil granulocytes from patients with monogenic disease using data-independent acquisition’, Molecular & cellular proteomics: MCP [Preprint]. Available at: https://doi.org/10.1074/mcp.RA118.001141. Gu, Z. (2022) ‘Complex heatmap visualization’, iMeta , 1(3). Available at: https://doi.org/10.1002/imt2.43. Huber, W. et al. (2002) ‘Variance stabilization applied to microarray data calibration and to the quantification of differential expression’, Bioinformatics , 18 Suppl 1, pp. S96–104. Available at: https://doi.org/10.1093/bioinformatics/18.suppl_1.s96. Keenan, A.B. et al. (2019) ‘ChEA3: transcription factor enrichment analysis by orthogonal omics integration’, Nucleic acids research , 47(W1), pp. W212–W224. Available at: https://doi.org/10.1093/nar/gkz446. Kern, F. et al. (2021) ‘miRTargetLink 2.0-interactive miRNA target gene and target pathway networks’, Nucleic acids research , 49(W1), pp. W409–W416. Available at: https://doi.org/10.1093/nar/gkab297. Langmead, B. and Salzberg, S.L. (2012) ‘Fast gapped-read alignment with Bowtie 2’, Nature methods , 9(4), pp. 357–359. Available at: https://doi.org/10.1038/nmeth.1923. Larsen, M.T. et al. (2012) ‘miRNA Profiling in Human Granulopoiesis and in Activated Neutrophils From Skin Window’, Blood , 120(21), pp. 2136–2136. Available at: https://doi.org/10.1182/blood.V120.21.2136.2136. Makkar, H. et al. (2023) ‘Acute myeloid leukemia: novel mutations and their clinical implications’, American journal of blood research , 13(1), pp. 12–27. Available at: https://www.ncbi.nlm.nih.gov/pubmed/36937458. Mao, J., Deng, Y. and Lê Cao, K.-A. (2024) ‘Φ-Space: Continuous phenotyping of single-cell multi-omics data’, bioRxiv . Available at: https://doi.org/10.1101/2024.06.19.599787. McInnes, L. et al. (2018) ‘UMAP: Uniform Manifold Approximation and Projection’, Journal of open source software , 3(29), p. 861. Available at: https://doi.org/10.21105/joss.00861. McInnes, L., Healy, J. and Melville, J. (2018) ‘UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction’, arXiv [stat.ML] . Available at: http://arxiv.org/abs/1802.03426 (Accessed: 18 September 2024). Meissner, F. et al. (2013) ‘Direct proteomic quantification of the secretome of activated immune cells’, Science , 340(6131), pp. 475–478. Available at: https://doi.org/10.1126/science.1232578. Milacic, M. et al. (2024) ‘The Reactome Pathway Knowledgebase 2024’, Nucleic acids research , 52(D1), pp. D672–D678. Available at: https://doi.org/10.1093/nar/gkad1025. Ng, M.S.F. et al. (2024) ‘Deterministic reprogramming of neutrophils within tumors’, Science , 383(6679), p. eadf6493. Available at: https://doi.org/10.1126/science.adf6493. R Foundation for Statistical Computing, R. (2018) ‘R: A language and environment for statistical computing’, RA Lang Environ Stat Comput [Preprint]. Ritchie, M.E. et al. (2015) ‘limma powers differential expression analyses for RNA-sequencing and microarray studies’, Nucleic acids research , 43(7), p. e47. Available at: https://doi.org/10.1093/nar/gkv007. Roehr, J.T., Dieterich, C. and Reinert, K. (2017) ‘Flexbar 3.0 - SIMD and multicore parallelization’, Bioinformatics , 33(18), pp. 2941–2942. Available at: https://doi.org/10.1093/bioinformatics/btx330. Rohart, F. et al. (2017) ‘mixOmics: An R package for ’omics feature selection and multiple data integration’, PLoS computational biology , 13(11), p. e1005752. Available at: https://doi.org/10.1371/journal.pcbi.1005752. Rørvig, S. et al. (2013) ‘Proteome profiling of human neutrophil granule subsets, secretory vesicles, and cell membrane: correlation with transcriptome profiling of neutrophil precursors’, Journal of leukocyte biology , 94(4), pp. 711–721. Available at: https://doi.org/10.1189/jlb.1212619. Schulte-Schrepping, J. et al. (2020) ‘Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell Compartment’, Cell , 182(6), pp. 1419–1440.e23. Available at: https://doi.org/10.1016/j.cell.2020.08.001. Shi, J., Ding, W. and Lu, H. (2020) ‘Identification of long non-coding RNA SNHG family as promising prognostic biomarkers in acute myeloid leukemia’, OncoTargets and therapy , 13, pp. 8441–8450. Available at: https://doi.org/10.2147/OTT.S265853. Sollberger, G. et al. (2024) ‘Quantitative proteomics reveals tissue-specific, infection-induced and species-specific neutrophil protein signatures’, Scientific reports , 14(1), p. 5966. Available at: https://doi.org/10.1038/s41598-024-56163-6. Soutschek, M. et al. (2022) ‘enrichMiR predicts functionally relevant microRNAs based on target collections’, Nucleic acids research , 50(W1), pp. W280–W289. Available at: https://doi.org/10.1093/nar/gkac395. Strydom, N. and Rankin, S.M. (2013) ‘Regulation of circulating neutrophil numbers under homeostasis and in disease’, Journal of innate immunity , 5(4), pp. 304–314. Available at: https://doi.org/10.1159/000350282. The Association to Function Knowledge Portal (2024). Available at: a2fkp.org. Thul, P.J. et al. (2017) ‘A subcellular map of the human proteome’, Science , 356(6340). Available at: https://doi.org/10.1126/science.aal3321. UniProt Consortium (2023) ‘UniProt: the Universal Protein Knowledgebase in 2023’, Nucleic acids research , 51(D1), pp. D523–D531. Available at: https://doi.org/10.1093/nar/gkac1052. Wang, L., Wang, S. and Li, W. (2012) ‘RSeQC: quality control of RNA-seq experiments’, Bioinformatics , 28(16), pp. 2184–2185. Available at: https://doi.org/10.1093/bioinformatics/bts356. Ward, J.H., Jr (1963) ‘Hierarchical grouping to optimize an objective function’, Journal of the American Statistical Association , 58(301), pp. 236–244. Available at: https://doi.org/10.1080/01621459.1963.10500845. WGCNA: an R package for weighted correlation network analysis (no date) BioMed Central . Available at: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-559 (Accessed: 3 May 2024). Wiśniewski, J.R. et al. (2014) ‘A “proteomic ruler” for protein copy number and concentration estimation without spike-in standards’, Molecular & cellular proteomics: MCP , 13(12), pp. 3497–3506. Available at: https://doi.org/10.1074/mcp.M113.037309. Wu, P. et al. (2021) ‘The trans-omics landscape of COVID-19’, Nature communications , 12(1), p. 4543. Available at: https://doi.org/10.1038/s41467-021-24482-1. Zhu, C.-Y. et al. (2009) ‘Cell growth suppression by thanatos-associated protein 11(THAP11) is mediated by transcriptional downregulation of c-Myc’, Cell death and differentiation , 16(3), pp. 395–405. Available at: https://doi.org/10.1038/cdd.2008.160. Wong, A. C., Wong, J. J., Rasko, J. E., & Schmitz, U. (2024). SpliceWiz: interactive analysis and visualization of alternative splicing in R. Briefings in Bioinformatics , 25 (1). Materials and methods Sample collection The study protocol was approved by the Ethics Committee of the Medical Faculty of the Ludwig Maximilian University of Munich (Reg. No. 19-685 2). Four healthy German volunteers over the age of 25, who had not previously been diagnosed with cancer or a serious illness, were included in the study. All healthy volunteers had signed a written consent form prior to the collection of 300 ml of bone marrow from the iliac crest. The bone marrow was anticoagulated in EDTA (Thermo Scientific 005500, 0.1 ml EDTA / 50 ml bone marrow) and subjected to negative selection to remove the non-granulocytic cell fractions using Miltenyi's neutrophil isolation kit (Miltenyi Biotec, 130-104-434) according to the manufacturer's protocol. To isolate defined precursors of granulopoiesis, we created a novel flow cytometry panel (sFig. 3 Ia) to detect 6 substages of human granulopoiesis in human bone marrow (sFig. 3 Ib) based on published surface antigens (Elghetany et al, 2004; Fajtova and Babusikova, 2010; Gorczyca et al, 2011; Gustafson et al, 2015). The isolated cells were dispersed in 1 ml PBS (Gibco, 14190-169) containing 10 % FCS (Value FBS South-American, ThermoFisher A5256701) and stained with a ten-color antibody panel for sorting using a BD FACSAria III cell sorter with blue, violet, red and yellow/green laser. After exclusion of dead cells and monocytes with BD Horizon Fixable Viability Stain (Red 780, BD 565388) and CD14 (APC/Fire750, Biolegend 301854), the major maturation stages were differentiated using CD117+ (PE, BB700, BD 566548) and CD49d+ (PE, BD 555503). Double positive cells represented the earliest stages with myeloblasts (MB) and promyelocytes (PM), further differentiated by CD34+ (BV421, BD 562577) for MB and CD34- for PM. Intermediate maturation stages were selected via the CD117- and CD49d+ population with myelocytes (MC) CD11b+ (FITC, Biolegend 301330) and metamyelocytes CD101+ (BV605, BD 747548). The terminal maturation stages (CD117- / CD49d-) were differentiated as band-stage (B) via CD24+ (PE-Dazzle594, Biolegend 311134) / CD16 dim (PE-Cy7, BD557744) / CD10- (CD10 BV711 BD 740770), and as a segmented stage via CD24+/CD16+/CD10+. To represent the terminal stage from peripheral blood, neutrophil granulocytes were isolated from peripheral venous blood using the neutrophil isolation kit in combination with an erythrocyte depletion cocktail (Miltenyi Biotec, 130-098-196). Light microscopic preparations from each collected fraction showed >95% pure cell populations with the respective stage-typical morphologies (sFig. 3 Ib). Sequencing For sequencing, we aimed for a cell count of 1x10^6 cells per sample and prioritized the collection of transcriptome and proteome samples over miRNA samples. A complete overview of the samples and cell numbers can be found in sTable 1. Cells were collected in aliquots of 0.1x10^6 - 1x10^6 cells, washed once in PBS, dissolved in storage reagents optimized for downstream analysis and frozen at -80 C. For proteome analysis, samples were stored in 5ul 25x protease inhibitor (Roche #04693159001). Transcriptome samples were stored in RLT-plus buffer (Qiagen #1053393) with the addition of 10 ul/1ml beta-mercaptoethanol (Thermo Scientific #21985023). To preserve the microRNA, the cells were dissolved in 750 ul TRIZOL-LS (Invitrogen #15596026). The sequencing of transcriptomes (mRNA) and microRNA was performed by a commercial partner, Azenta life sciences (Im Leuschnerpark 1B, 64347 Griesheim, Germany). RNA was extracted using the Qiagen RNeasy Mini Kit according to the manufacturer's protocol and sequenced on the Illumina NovaSeq 6000 with a 2x150 pair-end configuration. Raw BCL data were converted to FASTQ files and demultiplexed using Illumina bcl2fastq 2.20 software. Detailed protocols can be found in the extended Methods. Proteome sequencing was done with a nanoflow liquid chromatography (Thermo Scientific, EASY-nLC) coupled to a Orbitrap Exploris 480 instrument with a nano-electrospray ion source (Thermo Fisher Scientific). Protein sample preparation was done as described previously (Meissner et al. , 2013), with Trypsin (Promega #VA9000) and LysC (FUJIFILM Wako #121-05063) digestion. For a detailed protocol, please refer to the extended Methods. After analyzing the raw sequencing data, described in depth in the extended methods, we applied stringent quality criteria for feature inclusion that resulted in 8,432 transcriptome features, 3,156 proteins, and 283 microRNAs that were reliably quantified in at least two samples per stage. The complete list of identified features and their expression values (stage median expression values) can be found in the downloadable supplement tables sT2 (transcriptome), sT3 (proteome) and sT4 (miRNnome). Bioinformatic analysis FAIR data commitment We strive to ensure the highest standards of FAIR principles for data sharing (Findability, Accessibility, Interoperability, and Reuse) and appreciate that a major contribution of this work is to make the datasets and analysis workflows available to the community. To this end, we published our data in various forms on different platforms and create a dedicated website for data and analsis sharing. Raw data are made available via ProteomeXchange (Deutsch et al. , 2023) and Gene Expression Omnibus (Barrett et al. , 2013). To enable comparisons with other published transcriptome data, we integrated our data with stemformatics.org (Choi et al. , 2019), an interactive portal of high quality stem cell datasets. Since our results included a wealth of data beyond the transcriptome, we created a dedicated website, www.granulopoiesis.com, to allow readers to download all intermediate and final data along with the complete, reproducible code used in all steps of analysis, from raw counts to analysis. We have also created several ways to query our data directly on granulopoiesis.com. The expression of single and grouped features of transcriptome, proteome and miRNome can be visualized directly on the website. Raw data processing Here we describe the raw data processing only in short, while a detailed description can be found in the Extended Materials and Methods. Proteome raw data were processed with MaxQuant (version 1.6.2.1, (Cox and Mann, 2008)) with enabled ‘match between runs' option to receive a Uniprot annotated a maxLFQ based intensity table. Transcriptome (mRNA) raw sequence reads were aligned using STAR (v2.7.9a_2021-06-25, (Dobin et al. , 2013)) with RSeQC for quality control (Wang, Wang and Li, 2012). Data were processed as counts for differential expression analysis and tpm (transcripts per million) for all other forms of analysis (clustering, trajectories). MicroRNA raw data were processed with Flexbar (v3.5.0, (Roehr, Dieterich and Reinert, 2017)) and Bowtie2 (v2.4.4, (Langmead and Salzberg, 2012)) and exported as a counts table. After raw data processing, all further data analysis was performed using R (4.3.2, 2023-10-31 (R Foundation for Statistical Computing, 2018)) and R-Studio (2023.12.0+369 "Ocean Storm" (Allaire, 2012)). Feature and sample selection In order to ensure consideration of highly stringent features only, we applied strict inclusion criteria at each omics level. For the transcriptome, from the original 61,552 features, we selected only those that had valid expression values in at least 3 samples for each stage, with the exception of the proteome, where at least 2 samples were required to accommodate the fact that of the PMN population, only two samples were successfully sequenced. Valid expression values were defined as more than 10 counts or a mean tpm value of 0.28 adjusted for library size (average 37 million) (Law et al., 2016), resulting in a total of 8,432 transcriptome features. From the total 3503 quantified proteins, we only selected features that were quantified with at least 2 peptides per protein. Furthermore, we only included proteins that were quantified in at least 2 samples per stage and excluded known contaminants (MaxQuant contaminant list and all keratins), resulting in 3,156 proteins. For microRNA, the raw data contained 2950 features. Filtering for features with at least 10 counts in at least 2 samples per stage resulted in 283 stringent features for downstream analysis. All individual data sets were analyzed for outliers using the signal intensity and principal component analysis (PCA) before unification of technical repeat samples via median expression values (sFig. 4-6). No batch effects were detected in any of the datasets and no correction methods were applied. For normalization, we applied VSN (variant stabilization normalization) to each dataset (Huber et al. , 2002). Sample analysis A detailed description of the analysis steps and packages used can be found in the Extended Materials and Methods as well as the provided markdown scripts. Here we provide only a short summary of the applied methods and packages. Unless specified, base-R methods were used. For clustering analysis of each omics dataset, we tested euclidean-based wardD and k-means and plotted the results on the first two principal components (Fig. 2 I-IIIa). Stage specific features were selected using PLS-DA of micOmics (Rohart et al. , 2017) and plotted with the ComplexHeatmap package (Gu, 2022). Differential expression analysis for each dataset was done with Limma (Ritchie et al. , 2015). As significance thresholds, we applied a Benjamini Hochberg corrected p.value < 0.05 and a fold change equal or below/above +/-1. For biological-term enrichment, we employed ClusterProfiler to probe for GO terms (Gene Ontology Consortium et al. , 2023) and the pathway annotations of Reactome (Milacic et al. , 2024). For miRNA enrichment, we used miEAA (Aparicio-Puerta et al. , 2023) and miRTargetLink 2.0 (Kern et al. , 2021). A significance threshold of adj.P =1 (upwards) or <=-1 (down) was applied for each of the three omics levels. Trajectory analysis To derive feature trajectories over the granulopoietic developmental course, we combined weighted gene coexpression network analysis (WGCNA, ( WGCNA: an R package for weighted correlation network analysis , no date)) with time course analysis (timeOmics, (Bodein et al. , 2022)). WGCNA enables the identification of features with correlated expression intensities in a series of groups. This approach has been successfully used to identify large clusters of features with similar expression trajectories in time series, but a post-hoc analysis of individual feature medians showed large variation of the individual patterns (data not shown). To reduce the variability of identified features, we employed timeOmics' multivariate analysis to delineate highly stringent expression clusters across time and omics layers. We found an advantageous combination of both approaches by first determining the number of correlated features showing similar trajectories using WGCNA and using the respective numbers of features per module to define the feature numbers of timeOmics components. Next, we conducted an enrichment analysis of the features of the derived components with clusterProfiler and miRTargetLink 2.0 (Kern et al. , 2021). The resulting terms per omics layer were compared for overlaps between all 3 strata and plotted in combination to visualize the multi omics based term enrichment (Fig. 3 IV). Next, we probed features of trajectories on transcriptome and proteome with granulopoiesis associated disease genes (severe congenital neutropenia (Donadieu and Bellanné-Chantelot, 2022) and acute myeloid leukemia (Makkar et al. , 2023)) (Fig. 2 IVb). Functional microRNA enrichment For microRNA target enrichment, we employed the online based tool EnrichMIR (Soutschek et al. , 2022). It enables integration of a set of microRNA transcripts against a set of differentially expressed features to identify enrichment of targets and the motif-specific cumulative fold-change distribution of targets and non-targets. We interrogated the miRNA features of each detected trajectory (from Figure 2 III a-d) with the appropriate differentially expressed transcripts (see results text for details). First, we plotted all enriched miRNAs with their enrichment factors (x-axis), p-values (y-axis) and the fold change of the miRNA in the respective comparison group (Fig. 4 I-IV a). Next, we plotted the dose-response curves for features of the main enriched miRNA to visualize differences of expression intensity for features with (colored by binding motive) and without (black line) binding sites for the respective miRNA transcript (Fig. 4 I-IV b). Lastly, we applied GO enrichment to the miRNA targeted features to probe for their respective biological themes (Fig. 4 I-IV c) and plotted selected terms with their connected features (Fig. 4 I-IV d). As miRNAs generally act as mRNA transcript suppressors, we modeled the input data in a way that maintained this relationship. This required inversing the fold changes of differential expression in cases where the miRNAs decrease with maturation. Transcription factor enrichment To derive transcription factor (TF) signatures, we used the R-package of ChEA3 (Keenan et al. , 2019) to investigate the differentially expressed features between differentiation steps. The results were first plotted as the top 10% of TF signatures per step (Fig. 5 Ia), while we secondly conducted a comparative analysis between all signatures to derive common and stage specific TF signatures using an Upset analysis (Fig. 5 Ib). We employed the integrated top-rank results of ChIP-X Enrichment Analysis 3 ( ChEA3 ), a recently updated tool that enriches user-defined gene sets with libraries based on RNAseq, ChIPseq and published data (Keenan et al. , 2019). Omics level shifts Based on our previous results of the unique proteome rank distribution in neutrophil granulocytes (Grabowski et al. , 2019), we plotted the intensity based cumulative protein abundance (y-axis) per rank abundance (x-axis) of the different maturation stages (Figure 5 III a,c) on proteome and transcriptome level. Of note, we found that neutrophil granulocytes actively synthesize histones during maturation, and concluded that previously derived molecular numbers by the histone ruler (Wiśniewski et al. , 2014) are not applicable, as this method assumes a stable ratio between DNA and histones. Next, we summarized the subcellular locations of all expressed features. To this end, we combined the location data of the human protein atlas (Thul et al. , 2017), UniProt (UniProt Consortium, 2023), the GO consortium (Gene Ontology Consortium et al. , 2023) and the neutrophil granule and membrane specific publication of Rorvig at al (Rørvig et al. , 2013). For a detailed description of the process, please refer to the extended methods and the annotation preparation script. We summarized the total intensities per omics layer and subcellular location for each stage and used an alluvium plot for the location specific intensity as a percentage of the stages total (Brunson, 2020) (Figure 5 IV a,b). Integration with single cell transcriptome data To showcase the utility of our neutrophil atlas to enhance our understanding of scRNA-seq data, we applied Φ-Space (Mao et al., 2024) to map the neutrophil subset of the peripheral blood mononuclear cells (PMMC) scRNA-seq data from (Schulte-Schrepping et al. , 2020)) to our bulk atlas. Φ-Space outputs scores measuring the similarity between each single neutrophil and each of our bulk cell stages. We first visualized the similarity scores as a heatmap. Next, we computed PC1 and PC2 of the similarity scores and visualized the PCs. Lastly, we visualized the similarity scores using uniform manifold approximation and projection (UMAP,(McInnes, Healy and Melville, 2018)). Additional Declarations There is NO Competing Interest. Supplementary Files GA.png Graphical Abstract Extendedmethods.docx Supplementlegeds.docx SupplementfiguresandTable.pdf Cite Share Download PDF Status: Under Review Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6184761","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Resource","associatedPublications":[],"authors":[{"id":431279608,"identity":"2c6a6f94-d414-4424-a0da-1677461d48ca","order_by":0,"name":"Christoph Klein","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABW0lEQVRIie2RP0vDQBTAXwm0y4WsJxX9BMIF4VQIfg7HhMB16T/orgEhLi2uWfwIQl1KxwsHzZIPkEEwLpkq1EVSKOhdqpBUxdUhv+Hx7s+P994dQE3NP8QIALhKCCCTy3hQ5AB4e97wVLTLCk6UYheKrZTjPxVy8bklFa7KOd5Wgd+V9uQ5XOdW/wR0j6+Gj537m+upps9P4ShwF+l6bvWh5fKysh8RgWw2OvMMCAOS9WbxYqjpMQaasI45idkIUFYuQzADAbZwptwAgYjozZIu0XQfX8mEthu+cDzcJTtKmBeK7okNER36tCwUWWXwJpV3qQxWlfEZcFQoiAsgwqYJ+lK6TalwVaX6yLIxxNhITRSOiTBnMRuGd0qJM7o38V3HR1m5MSNg2mtuWX0Sjc0034hDGomH9MW/BBq5GV77585ty03hOzZopRWvnDV/uL/7WTU1NTU1FT4AmAiHoWefnO0AAAAASUVORK5CYII=","orcid":"https://orcid.org/0000-0003-0956-0445","institution":"Dr. von Hauner Children's Hospital, LMU, Munich","correspondingAuthor":true,"prefix":"","firstName":"Christoph","middleName":"","lastName":"Klein","suffix":""},{"id":431279609,"identity":"f12cf365-eff2-4b70-b4ea-a5b5fcb61790","order_by":1,"name":"Sebastian Hesse","email":"","orcid":"https://orcid.org/0000-0002-7148-3068","institution":"Ludwig-Maximilians-Universität München","correspondingAuthor":false,"prefix":"","firstName":"Sebastian","middleName":"","lastName":"Hesse","suffix":""},{"id":431279610,"identity":"1b02fcab-a75d-415e-9a85-f24df838fc87","order_by":2,"name":"Jiadong Mao","email":"","orcid":"https://orcid.org/0000-0002-3818-1981","institution":"Melbourne Integrative Genomics and School of Mathematics and Statistics, University of Melbourne, Melbourne, Victoria, Australia","correspondingAuthor":false,"prefix":"","firstName":"Jiadong","middleName":"","lastName":"Mao","suffix":""},{"id":431279611,"identity":"90f2dac3-095e-42d5-8b27-7e260d425c9f","order_by":3,"name":"Armin Hadziahmetovic","email":"","orcid":"https://orcid.org/0000-0002-3765-5460","institution":"Institute of Bioinformatics and Practical Informatics, Department of Informatics, Ludwig-Maximilians-Universität München (LMU), Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Armin","middleName":"","lastName":"Hadziahmetovic","suffix":""},{"id":431279612,"identity":"d8f144ea-3d86-453b-9cd2-d366d0330fe7","order_by":4,"name":"Felix Offensperger","email":"","orcid":"","institution":"Institute of Bioinformatics and Practical Informatics, Department of Informatics, Ludwig-Maximilians-Universität München (LMU), Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Felix","middleName":"","lastName":"Offensperger","suffix":""},{"id":431279613,"identity":"bd52d21a-946d-4802-9ada-f46f835b74a0","order_by":5,"name":"Francisca Rojas Ringeling","email":"","orcid":"","institution":"Gene Center, Ludwig-Maximilians-Universität München, Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Francisca","middleName":"Rojas","lastName":"Ringeling","suffix":""},{"id":431279614,"identity":"82506745-49b7-40da-aabc-154f87b98530","order_by":6,"name":"Pablo Monteagudo-Mesas","email":"","orcid":"","institution":"Gene Center, Ludwig-Maximilians-Universität München, Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Pablo","middleName":"","lastName":"Monteagudo-Mesas","suffix":""},{"id":431279615,"identity":"00436d7c-c553-4fdf-b602-41a276e6ae5c","order_by":7,"name":"Annika Frauenstein","email":"","orcid":"","institution":"Systems Immunology Laboratory, Max Planck Institute of Biochemistry, Martinsried/Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Annika","middleName":"","lastName":"Frauenstein","suffix":""},{"id":431279616,"identity":"a2e9d07a-ebf5-477f-aa3b-4f677e207b8d","order_by":8,"name":"Felix Meissner","email":"","orcid":"https://orcid.org/0000-0003-1000-7989","institution":"Institute of Innate Immunity, University of Bonn","correspondingAuthor":false,"prefix":"","firstName":"Felix","middleName":"","lastName":"Meissner","suffix":""},{"id":431279617,"identity":"108fdc55-79c6-49cb-beab-e6e0432b40d9","order_by":9,"name":"Stefan Canzar","email":"","orcid":"","institution":"Gene Center, Ludwig-Maximilians-Universität München, Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Stefan","middleName":"","lastName":"Canzar","suffix":""},{"id":431279618,"identity":"00701094-4ba7-4d43-8f1a-870469846673","order_by":10,"name":"Ralf Zimmer","email":"","orcid":"","institution":"Institute of Bioinformatics and Practical Informatics, Department of Informatics, Ludwig-Maximilians-Universität München (LMU), Munich, Germany","correspondingAuthor":false,"prefix":"","firstName":"Ralf","middleName":"","lastName":"Zimmer","suffix":""},{"id":431279619,"identity":"25a2965e-d332-4fb9-aef2-3e9fbd188060","order_by":11,"name":"Kim-Anh le Cao","email":"","orcid":"","institution":"University of Melbourne","correspondingAuthor":false,"prefix":"","firstName":"Kim-Anh","middleName":"le","lastName":"Cao","suffix":""}],"badges":[],"createdAt":"2025-03-08 15:00:10","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6184761/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6184761/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":79819393,"identity":"1909e4e8-b533-45aa-9e97-ef8601af0eca","added_by":"auto","created_at":"2025-04-03 08:26:32","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":834104,"visible":true,"origin":"","legend":"\u003cp\u003eWorkflow\u003c/p\u003e\n\u003cp\u003ea) Workflow to create a multiomics based molecular atlas of granulopoiesis. Bone marrow of 4 healthy donors (one with a biological replicate) was subjected to flow cytometry sorting to collect 6 histology-based maturation stages, while PMN were isolated from venous blood. We prepared transcriptome (mRNA), microRNA and proteome sequences of each stage from each donor. A dedicated website www.granulopoiesis.com provide the comprehensive bioinformatic analysis, data and analysis scripts.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/323a5be507bbf608074ea698.png"},{"id":79818273,"identity":"a2ab06c8-24a4-464c-a845-25e56e3e88c5","added_by":"auto","created_at":"2025-04-03 08:18:32","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2861220,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTranscriptome (I), proteome (II) and miRNAnome (III) of human granulopoiesis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eMyeloblasts (MB), promyelocytes (PM), myelocytes (MC), metamyelocytes (MM), band- (B) and segmented- (S) nucleated forms and peripheral blood neutrophils (PMN).\u003c/p\u003e\n\u003cp\u003eI-III a) Principal component analysis (PC1/2) of transcriptome, proteome and microRNA samples with encirclement of euclidean clusters and coloring according to maturation stages. While clustering of transcriptome and miRNome shows stringent association with maturation stages, the proteome is clustered in 3 groups (MB, PM with MC, as well as B, S and PMN).\u003c/p\u003e\n\u003cp\u003eI-III b) Barplots of the number of significantly differentially expressed features between development steps with adj.pval \u0026lt; 0.05 and abs(logFC) \u0026gt; 1. While each development step shows large changes on transcriptome level, proteome changes are detectable only from MB to MM with minimal changes observed afterwards. Similarly, microRNA differences were observed mainly during initial development (MB to MC) while later stages show more nuanced changes.\u003c/p\u003e\n\u003cp\u003eI-III c, d) Dotplots of comparative enrichment analysis in up (c) and downregulated (d) features. Dot size represents fold enrichment, number label indicates total feature counts per term. Shown are interomics-overlapping terms with the lowest p-value in each group (combined results of GO biological process, cellular component and molecular function). For microRNA, enrichment by subcellular localisations is displayed.\u003c/p\u003e\n\u003cp\u003eVI) Mean Z-values of transcriptome and proteome expression in neutrophil subcellular compartments.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/26c9397e16511a66eb83a2df.png"},{"id":79819395,"identity":"02a4bd65-6b72-42ef-90e1-31afed2efb5f","added_by":"auto","created_at":"2025-04-03 08:26:32","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2125388,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eLongitudinal multiomic trajectories reveal molecular pathways of granulopoiesis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eMolecular trajectories of transcriptome (I), proteome (II), and miRNome (III) in human granulopoiesis. Graphical illustration longitudinal trajectory curves along maturation stages on the x-axis and relative expression intensities on the y-axis. Colors indicate trajectory cluster association of features (yellow = high, blue = low) based on the distance to the median expression curve (red).\u003c/p\u003e\n\u003cp\u003eI-II a,b) Features with unidirectional course, c) Features with bi-directional course, d) Table of features per trajectory as number and percentage of the total features per omics plane.\u003c/p\u003e\n\u003cp\u003eIV a) Reactome enrichment terms (y-axis) for trajectory features from all three omics planes along the x-axis. Shapes indicating omics levels with circle = transcriptome (mrn), filled circle = proteome (pro), diamond = miRNome (mic), colors indicate trajectory association with orange = up, purple = down, blue = up down.\u003c/p\u003e\n\u003cp\u003eVI b) Upset plot demonstrating overlaps in two sets of neutrophil associated genes, SCN (severe congenital neutropenia) and AML (acute myeloid leukemia), and trajectory features of transcriptome (mrn) and proteome (pro).\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/1ad22d30f17e4e85d55cdb3e.png"},{"id":79818275,"identity":"0acfabff-ecd6-4abe-8ab5-8c4aa5b24cab","added_by":"auto","created_at":"2025-04-03 08:18:32","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":2482295,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003emiRNA governance of human granulopoiesis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea) Dotplots of enrichment results with enrichment factor on x-axis and p-value on y-axis. Dot color indicates log fold change for miRNA in the probed development step (see below). Bold labeled dots represent miRNA with FDR \u0026lt; 0.05 enrichment, roman labeled dots miRNAs with p-values \u0026lt; 0.05.\u003c/p\u003e\n\u003cp\u003eb) Dose-response curves of miRNA targets (colored by binding motif) vs non-targets (black). Fold-change values of features in the probed development step along the x-axis, cumulative proportion of features on the y-axis.\u003c/p\u003e\n\u003cp\u003eA left-shift of the colored curves indicates a lower expression of features displaying miRNA binding sites.\u003c/p\u003e\n\u003cp\u003ec) Dotplots of GO enrichment results for features regulated by the respective miRNA transcript. GeneRatio (total in term / identified in term) on the x-axis, enriched terms on the y-axis, dot color represents adjusted p-values (Benjamini Hochberg corrected) for enrichment.\u003c/p\u003e\n\u003cp\u003ed) c-net plot of selected enrichment terms and the identified transcriptome features. Dot color represents log fold change of the feature in the probed development step\u003c/p\u003e\n\u003cp\u003eI) Enrichment of miRNAs with upward trajectory (Figure 3 IIIa) for differentially expressed (DE) transcriptome features between MB and PMN.\u003c/p\u003e\n\u003cp\u003eII) Enrichment of miRNAs with downward trajectory (Figure 3 IIIb) for DE transcriptome features between PMN and MB (resulting in inverted fold change directions to maintain the left shift of the response curves).\u003c/p\u003e\n\u003cp\u003eIII) Enrichment of miRNAs with up-down trajectory (Figure 3 IIIc upper) for DE transcriptome features between B and MB.\u003c/p\u003e\n\u003cp\u003eIV) Enrichment of miRNAs with down-up trajectory (Figure 3 IIIc lower) for DE transcriptome features between MM and MB.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/682975c5035941c69042065b.png"},{"id":79818276,"identity":"f2258abb-e541-46a0-b71a-41118cda0878","added_by":"auto","created_at":"2025-04-03 08:18:32","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":4236687,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTranscription factor activity, long non-coding RNAs and omics-shifts in human granulopoiesis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eI a) Top 10% of enriched transcription factor signatures per development step identified in upregulated transcriptome features. X-axis represents the enrichment factor (= 1-adj.pval), y-axis indicates gene symbols of transcription factors, color represents log fold change in development step.\u003c/p\u003e\n\u003cp\u003eI b) Upset plot to compare overlaps in transcription factor signatures in all development steps. Y-axis represents set names and total set size, bars at the top represent sizes of intersection sets. Single dots represent unique sets, connected dots represent overlaps between the sets. It becomes evident that the last development step (StoPMN) shows the largest signal of unique TF signatures, while signatures of 32 TFs are identified in every development step.\u003c/p\u003e\n\u003cp\u003eII a) Heatmap of the top40 variably expressed long non-coding RNAs. Expression normalized to z-values. Gene symbols or, if unavailable ENSG IDs at the right edge, samples colored by maturation stage at the lower edge\u003c/p\u003e\n\u003cp\u003eII b) Bar plots of the total number of differentially expressed long non-coding RNAs (up- and down-regulated) between the consecutive maturation stages (adj.pval \u0026lt; 0.05 and abs(logFC) \u0026gt; 1).\u003cbr\u003e\nII c) Bar plots of the total number of differentially expressed protein coding RNAs (up- and down-regulated) between the consecutive maturation stages (adj.pval \u0026lt; 0.05 and abs(logFC) \u0026gt; 1).\u003c/p\u003e\n\u003cp\u003eIII a,b) Cumulative abundance ranks in proteome (a) and transcriptome (b) for the different development stages (colored). Inset table shows total features of the 50th abundance percentile in the different stages.\u003c/p\u003e\n\u003cp\u003eIV a,b) Flow-barplots representing the summed up protein intensities and sequencing reads per subcellular localisation for proteome (c) and transcriptome, respectively (d). X-axis represents cell stages in physiological order, y-axis and colors represents subcellular locations. Percentages indicate the intensity fraction of the total signal assigned to proteins annotated to the indicated subcellular location in the respective maturation stage. \u003cbr\u003e\nV a) Heatmap of Φ-Space scores, measuring the similarity between clusters of single neutrophils from severe COVID-19 patients with our morphology-based maturation stages: the higher the Φ-Space score, the stronger the similarity. Original clusters from Schulte-Schrepping at the left edge, maturation-stage based clustering at the lower edge. \u003cbr\u003e\nV b) UMAP of re-clustered single cell data showing improved separability of cell types according to histological maturation stages.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/132b20f2e6a6a6e7a61b3be6.png"},{"id":79819765,"identity":"2bfcc756-c2b1-41bc-b16e-608d00875aaa","added_by":"auto","created_at":"2025-04-03 08:34:39","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":16419847,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/4e5156bd-1807-48be-a617-0d61eaa44dca.pdf"},{"id":79818266,"identity":"5cfb47b2-1a61-4087-a60f-93ffdcd91569","added_by":"auto","created_at":"2025-04-03 08:18:32","extension":"png","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":930000,"visible":true,"origin":"","legend":"\u003cp\u003eGraphical Abstract\u003c/p\u003e","description":"","filename":"GA.png","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/8cc37be7fb2310aa5fd3cd5b.png"},{"id":79818269,"identity":"173c9339-535a-457f-90ab-bcb18e179cd5","added_by":"auto","created_at":"2025-04-03 08:18:32","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":1949242,"visible":true,"origin":"","legend":"","description":"","filename":"Extendedmethods.docx","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/d3ca4d68d7a5500db0b4f812.docx"},{"id":79819394,"identity":"6597896b-3b2f-422f-961c-263fe3ce32e7","added_by":"auto","created_at":"2025-04-03 08:26:32","extension":"docx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":1686147,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementlegeds.docx","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/b2ce4a682b8b3feb912b38bb.docx"},{"id":79818278,"identity":"503c4246-f085-49ac-bec9-26778cd5957f","added_by":"auto","created_at":"2025-04-03 08:18:33","extension":"pdf","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":6307308,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementfiguresandTable.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6184761/v1/a28bf16cff79e1f3eb7cbb39.pdf"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"A molecular atlas of human granulopoiesis","fulltext":[{"header":"Objectives","content":"\u003cul\u003e\n \u003cli\u003eCreation of a stage-specific multiomics atlas of human granulopoiesis\u003c/li\u003e\n \u003cli\u003eIntegrative analysis of transcriptome (mRNA), proteome and non-coding RNAs (miRNA, lncRNA)\u003c/li\u003e\n \u003cli\u003eQuantification of miRNA-mRNA interactions\u003c/li\u003e\n \u003cli\u003eIdentification of functional transcription factor signatures and global omics shift\u003c/li\u003e\n\u003c/ul\u003e\n"},{"header":"Key results","content":"\u003cul\u003e\n \u003cli\u003eA comprehensive molecular atlas of human granulopoiesis with 7 developmental stages at the mRNA, protein and non-coding RNA level, presented as FAIR-accessible data and code on www.granulopoiesis.com\u003c/li\u003e\n \u003cli\u003eRevised understanding of granule protein synthesis, as transcription and translation are initiated in parallel for all granule subsets and terminated at the metamyelocyte stage\u003c/li\u003e\n \u003cli\u003eUnraveling of hsa-miR-106a-5p and hsa-miR-125a-5p as potential inhibitors of myeloid differentiation in lineage committed neutrophil progenitors\u003c/li\u003e\n \u003cli\u003eIdentification of novel stage specific long non-coding RNAs (ENSG00000288819, ENSG00000203644) in terminal maturation stages\u003c/li\u003e\n \u003cli\u003eIdentification of LTF (lactotransferrin) as a central transcription factor in early granulopoiesis, a signature of 32 TFs active throughout each maturation step and an exclusive signature of 131 TFs active in the final progression from the bone marrow S-stage to peripheral blood neutrophils\u003c/li\u003e\n \u003cli\u003eStage-specific molecular markers on mRNA, miRNA, lncRNA and protein level\u003c/li\u003e\n\u003c/ul\u003e"},{"header":"Novelties","content":"\u003cul\u003e\n \u003cli\u003eUnprecedented resolution of the molecular landscapes of human granulopoiesis\u0026nbsp;\u003c/li\u003e\n \u003cli\u003eIntegrative analysis workflows to define and interpret feature trajectories of cell-development\u003c/li\u003e\n \u003cli\u003eIdentification non-coding RNA governance factors on miRNA and lncRNA level \u0026nbsp;\u0026nbsp;\u003c/li\u003e\n \u003cli\u003eAnnotation of single cell profiles in accordance with histology based cell identities\u003c/li\u003e\n\u003c/ul\u003e"},{"header":"Introduction","content":"\u003cp\u003eNeutrophil granulocytes are short-lived effector cells of the innate immune system that make up 50\u0026ndash;70% of leukocytes in human peripheral blood. Their functions range from defense against pathogens to tissue remodeling and they play a role in cancer defense as well as carcinogenesis (Burn et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Sophisticated attack mechanisms enable them to rapidly destroy microorganisms or tissue. Their eponymous granules are typically divided into primary (azurophilic), secondary (specific), tertiary, and ficolin-1-rich subsets that are thought to harbor mostly non-overlapping sets of proteins. Other effector functions include the generation of reactive oxygen species (ROS), the formation of neutrophil extracellular traps (NETs) and phagocytosis (Faurschou and Borregaard, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2003\u003c/span\u003e). Human diseases affecting neutrophils range from deficits in one of their specialized defense mechanisms, such as deficiency of specific granules or reactive oxygen species generation, to defects in granulopoiesis that include over 30 known monogenic causes of congenital neutropenia (Donadieu and Bellann\u0026eacute;-Chantelot, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2022\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eGranulopoiesis, the differentiation of neutrophil granulocytes in the bone marrow, is thought to follow a linear developmental process that is divided into major cellular states of metamorphosis defined by specific states of transcriptome, proteome and non-coding RNA landscapes (Larsen et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). They differentiate from lineage committed stem-cell like myeloblasts (MB) to promyelocytes (PM), myelocytes (MC), metamyelocytes (MM), band (B) and segmented (S) nuclear forms. At the terminal maturation stage, cells emigrate from the bone marrow into the peripheral blood, referred to as polymorphonuclear neutrophils (PMN) (Elghetany et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Granulopoiesis is a high-throughput differentiation process that transforms readily dividing stem cells into terminal effector cells with an enormous steady-state production rate of 10^11 cells per day (Strydom and Rankin, \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2013\u003c/span\u003e).\u003c/p\u003e \u003cp\u003ePreviously, we and others showed that the highly specialized function of the neutrophil is mirrored in its proteome configuration, as less than 10 proteins make up over 50% of the total protein content (Grabowski et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2019\u003c/span\u003e) (Sollberger et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Several groups have recently provided comprehensive multiomic analyses of neutrophil granulocytes subpopulations in diseases such as Covid (Schulte-Schrepping et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Wu et al., \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), or cancer (Ng et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). To better understand the molecular landscape of developing neutrophil granulocytes, we set out to characterize their stage-specific molecular signatures on the transcriptome, proteome and non-coding RNA level (Fig.\u0026nbsp;1a).\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003eTranscription-translation dynamics of granule-subtypes shows parallel and time-delayed synthesis\u003c/p\u003e \u003cp\u003eBased on defined precursor stages of granulopoiesis (Kwok et al. 2020), we developed a panel for flow cytometric cell sorting to purify bone marrow cells from four healthy donors plus one biological replicate measurement (Fig.\u0026nbsp;1a). Using stringent quality criteria (sFig. 3II), we quantified 8,432 mRNA transcripts (transcriptome), 3,156 proteins (proteome) and 283 mature miRNA transcripts (miRNome).\u003c/p\u003e \u003cp\u003eFor initial analysis, we considered each omics level individually to examine its coherence with the hypothesized developmental stage association from our FACS and microscopy analysis. A high-level analysis was based on principal component analysis (PCA) for dimensionality reduction in combination with Ward-D based Euclidean clustering (Ward, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e1963\u003c/span\u003e). For each omics level, we observed clear separation of stages in a developmentally coherent manner (Fig.\u0026nbsp;2I-III a). Clustering revealed near perfect separability of the individual stages at the level of the transcriptome (Fig.\u0026nbsp;2Ia) and miRNome (Fig.\u0026nbsp;2IIIa), while the proteome showed three main developmental clusters: MB, PM with MC and MM, B and PMN (Fig.\u0026nbsp;2IIa). Thus, the maturation process of granulocytes might be characterized by highly dynamic and stage-specific RNA landscapes, while changes in the proteome overlap considerably between stages, indicating a gradual developmental process. Since we observed a strong coherence of the samples with the assumed developmental stages, we next analyzed the changes between each successive developmental step with differential expression (DE) analysis. We detected a large number of differentially expressed features (1000\u0026ndash;2000 combined) at the transcriptome level in each developmental stage, while significant changes at the proteome level were only detected during the initial differentiation (MB-\u0026gt;PM-\u0026gt;MC-\u0026gt;MM). After the metamyelocyte stage, the proteome showed only minimal changes in protein abundance levels (0\u0026ndash;9 DE proteins). At the miRNA level, the most pronounced shift was observed from the myeloblasts to the myelocytes, with a steady state during final maturation (Fig.\u0026nbsp;2IIIa).\u003c/p\u003e \u003cp\u003eAs clustering and DE indicated high stringency of the obtained omics signals for each maturation stage, we applied enrichment analysis to determine the biological themes of the up-regulated (Fig.\u0026nbsp;2I-III d) and down-regulated (Fig.\u0026nbsp;2I-III e) features. At the transcriptome level, all known granule subfractions (azurophilic\u0026thinsp;=\u0026thinsp;primary, specific\u0026thinsp;=\u0026thinsp;secondary, tertiary and ficolin-1 granules) were upregulated in myeloblasts and promyelocytes, followed by a gradual and stage-specific downregulation of granule subfraction features. At the proteome level, myeloblasts begin translation of the azurophil and ficolin-1 granules, while promyelocytes show translation of all granule fraction features. Myelocytes stop the translation of the primary granule proteins and simultaneously upregulate the translation of the specific, tertiary and ficolin-1-rich granules. In addition, we found an upregulation of ribosomal proteins from myeloblasts to promyelocytes, coinciding with granule protein synthesis, followed by a downregulation of ribosomal transcripts and loss of ribosomal proteins from the metamyelocyte stage onwards, suggesting termination of protein synthesis.\u003c/p\u003e \u003cp\u003eTo analyze how features of defined compartments (granule subfractions, plasma membrane, ribosomes and mitochondria) behave on mRNA and protein level throughout development, we plotted the summarized medians along the longitudinal development path (Fig.\u0026nbsp;2 VI). Granule proteins showed stage specific patterns that exhibited primarily coordinated transcriptome-proteome increases (primary and ficolin-1 granules from MB to PM) as well as translational time lags (secondary gr. in PM to MC) followed by transcriptional downregulation while protein levels remained stable. Ribosomes and mitochondria showed an initial increase of protein intensities with stable mRNA levels in promyelocytes, followed by continuous loss of features on transcriptome and proteome level.\u003c/p\u003e \u003cp\u003eThus we found evidence for developmental stage-specific signals in each of the three omics layers. Interestingly, transcription of all granule subtypes appeared to be initiated at the myeloblast stage, whereas translation initiation and transcription termination showed evidence for stage specificity. However, we found ample evidence for coordinated synthesis of granule subtype features. While significant proteome changes were only observed up to the metamyelocyte stage, transcriptome changes showed a second peak in the final stages of maturation while the proteome remained stable.\u003c/p\u003e \u003cp\u003eLongitudinal analysis reveals developmental feature trajectories and allows unified pathway enrichment on multiomic plane\u003c/p\u003e \u003cp\u003eSince differential expression indicated complex patterns, we conducted a longitudinal analysis of feature expression throughout development, here referred to as \u0026ldquo;trajectories\u0026rdquo;. To this end, we combined WGCNA with a multivariate time course analysis (described in Methods). We found that the majority of traits in each omics level showed a unidirectional expression pattern (only up or down). In the transcriptome, 1629 features (\u0026thinsp;≙\u0026thinsp;19% of the total, Fig.\u0026nbsp;3 Ia) showed an increasing trajectory and 1970 features (23%, Fig.\u0026nbsp;2Ib) showed a decreasing trajectory. At the proteome level, 483 (15%, Fig.\u0026nbsp;3 IIa) proteins increased and 504 (16%, Fig.\u0026nbsp;3 IIb) proteins decreased during human granulopoiesis. In the miRNome, 54 transcripts were identified that increased (19%, Fig.\u0026nbsp;3 IIIa) and 68 that decreased (24%, Fig.\u0026nbsp;3 IIIb). For a smaller subset of features, we observed bidirectional trajectories comprising a total of 674 transcripts (UpDown: 361, DownUp: 313, Fig.\u0026nbsp;2 Ic), 590 proteins (UpDown: 329, DownUp: 261, Fig.\u0026nbsp;3 IIc) and 33 miRNA transcripts (UpDown: 21, DownUp: 12, Fig.\u0026nbsp;3 IIIc). A minority of features showed multidirectional trajectories (sFig. 1 I-II).\u003c/p\u003e \u003cp\u003eTo derive the common biological themes of features in each trajectory, we filtered for overlapping terms after enriching for pathways in the Reactome database (Fig.\u0026nbsp;3 IV) (Milacic et al., \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Interestingly, we observed a negative correlation between miRNAs and features of transcriptome and proteome, consistent with the notion that miRNAs exert mainly down-regulatory effects (Shang et al., 2023). For example, pathways with up-regulated features at the transcriptome and proteome levels and down-regulated interacting miRNAs were VEGF, Ras, mTOR, MAPK and HIF1, while we also observed this pattern for features of cell-cell adhesion of leukocytes and extracellular trap formation (Fig.\u0026nbsp;3 IV, upper part). Conversely, we observed down-regulated mRNA and protein features with up-regulated interfering miRNAs for ribosomes, the translation apparatus and mitochondrial membrane proteins, as well as the signal recognition particle (SRP), nonsense-mediated decay (NMD) and the regulatory network for the expression of SLIT and ROBO proteins (Fig.\u0026nbsp;3 IV, lower part). Some of the pathways and cell compartments showed complex patterns of interacting miRNAs with bidirectionally regulated transcripts and bimodal protein trajectories (like ribosomes, translation, SRP complex).\u003c/p\u003e \u003cp\u003eWe compared the genes involved in acute myeloid leukemia (AML) and severe congenital neutropenia (SCN) with the features of the transcriptome and proteome (Fig.\u0026nbsp;3V), respectively. Intriguingly, we found that the majority of genes responsible for severe congenital neutropenia showed a downward trajectory of expression levels at both omics levels, demonstrating a reduction in expression during granulopoiesis.\u003c/p\u003e \u003cp\u003eIn summary, we found that the majority of features at all omics levels exhibit continuously increasing or decreasing trajectories with biologically coherent themes at the transcriptome, proteome and miRNA levels.\u003c/p\u003e \u003cp\u003eSpecific miRNA transcripts govern development stage associated processes in human granulopoiesis\u003c/p\u003e \u003cp\u003eIn an effort to determine hierarchical effects of miRNA transcripts on the granulopoiesis transcriptome, we probed the miRNAs of trajectory clusters for binding-site enrichment in differentially expressed RNA transcripts (using enrichMIR (Soutschek et al., \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), see methods). Coupled with enrichment analysis of the mRNA transcripts, we could link miRNAs with the cellular processes they may control. For miRNAs with increasing trajectory (Fig.\u0026nbsp;3 IIIa), we selected DE mRNA transcripts between the two endpoints, myeloblasts (MB) and peripheral blood neutrophils (PMN). Here we found hsa-miR-193a-5p (Fig.\u0026nbsp;4 Ia) to interact with a total of 86 mRNA features (Fig.\u0026nbsp;4 Ib) that enriched for the ATPase complex and the SWI/SNF complex (Fig.\u0026nbsp;4I c,d). For the 68 miRNAs that decrease during granulopoiesis (Fig.\u0026nbsp;3 IIIb), we again enriched DE features that alternate between the two endpoints, but in reverse order, from PMN to MB to maintain the expected fold-change direction for the enrichMIR algorithm. We found two miRNA clusters, hsa-miR-19a/b-3p and hsa-miR-106a/17/20a/b-5p (Fig.\u0026nbsp;4 IIa), that enriched for a large group of mRNA transcripts in myeloblasts (Fig.\u0026nbsp;4 IIb). Interestingly, transcripts associated with the latter group (miR-106a) showed enrichment for myeloid cell differentiation, cell growth and histone modification (Fig.\u0026nbsp;4 II c,d). For miRNA features with up-down trajectory (Fig.\u0026nbsp;3 IIIc), we analyzed the features between band stage and myeloblasts (Fig.\u0026nbsp;4 III). We found that hsa-mir-454-3p binds a large group of mRNA transcripts associated with negative regulation of transcription by RNA polymerase II. Conversely, miRNAs with a down-up pattern were enriched for DE features between metamyelocytes and myeloblasts with terms similar to those of miRNAs with a decreasing trajectory, namely myeloid cell differentiation and cell growth (Fig.\u0026nbsp;4 IV). In the miRNA trajectory group showing a late up-down expression pattern (sFig. 1 II), we found that hsa-miR-1343-3p interacts with DE features between metamyelocytes and band stage granulocytes, that were enriched for translation and ribosomes (sFig. 1 III a-d).\u003c/p\u003e \u003cp\u003eTranscription factor signatures in human granulopoiesis reveal a role of LTF in early differentiation steps\u003c/p\u003e \u003cp\u003eIn the final maturation step, from segmented nucleated neutrophils in the bone marrow to peripheral blood neutrophils, GO enrichment identified the downregulation of a group of 22 transcription factors. This finding encouraged us to further investigate transcription factor signatures in granulopoiesis. To this end we employed ChIP-X Enrichment Analysis 3 (Keenan et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), that enriches TF targets in DE transcriptome sets (details in Methods).\u003c/p\u003e \u003cp\u003eFirst, we focused on the 10% most enriched TF signatures in DE transcripts of each successive development step (Fig.\u0026nbsp;5 Ia). In addition to previously described TFs orchestrating granulopoiesis (e.g. NFE2, CEBPs, MXD1), we identified LTF (lactotransferrin) most strongly induced and active in the two early stage progressions (MB to PM and PM to MC, Fig.\u0026nbsp;5 Ia).\u003c/p\u003e \u003cp\u003eWe extended the analysis with a comparison of all significantly enriched TF signatures between all developmental stages (Fig.\u0026nbsp;5 Ib). We identified a set of 32 TFs expressed and active at each step of granulopoiesis, in particular NFE2, MXD1, ZNF467, MXD3 and SP1. Another 40 TF signatures were shared between all initial development steps (MB to B), with prominent appearances of CEBPB and CEBPD. When we compared mature neutrophils in the bone marrow (S) with terminal neutrophils in peripheral blood (PMN), we identified an exclusive signal of 131 active TFs. Most of them, e.g. THAP11 (Thanatos-associated protein 11), associated with cell growth suppression (Zhu et al, \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2009\u003c/span\u003e), have not previously been described in granulopoiesis.\u003c/p\u003e \u003cp\u003eStage specificity of long non-coding RNA expression\u003c/p\u003e \u003cp\u003elncRNAs emerge as critical factors in acute myeloid leukemia, driving proliferation, inhibiting differentiation and conferring treatment resistance (Connerty and Lock, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). We queried our comprehensive dataset to assess the abundance of lncRNA during defined stages of differentiation.\u003c/p\u003e \u003cp\u003eFigure\u0026nbsp;5 IIa shows a heatmap of the top40 variably expressed lncRNA transcripts during granulopoiesis. They show highly specific expression patterns with SNHG (Small nucleolar RNA host genes) family members prominently expressed from MB to MC, NEAT1 and MAILR expression only from Band stage on, and various currently uncharacterized transcripts. Of special note, ENSG00000288836 was expressed specifically in myelocytes and metamyelocytes, and spans a genetic region annotated with variant associations that include neutrophil count (\u003cem\u003eThe Association to Function Knowledge Portal\u003c/em\u003e, 2024), sFig. 2 IIa).\u003c/p\u003e \u003cp\u003eTwo poorly studied lncRNA were strongly expressed in segmented and polymorphonuclear stages, ENSG00000288819 and ENSG00000203644, respectively. While the targets of these ncRNA remain to be defined, ITGB2-AS1, a known antisense RNA transcript to Integrin subunit beta2, was found exclusively expressed from the metamyelocyte stage onwards.\u003c/p\u003e \u003cp\u003eAs our global transcriptome included all transcripts in unison, we separated protein-coding mRNA and long non-coding transcripts and plotted the number of differentially expressed features (Fig.\u0026nbsp;5 II b,c). As the protein coding mRNA constituted 88% of transcripts, the distribution of significantly up- and down-regulated features did not change in comparison to the total transcriptome (Fig.\u0026nbsp;1 Ic). However, the DE lncRNAs showed a markedly different pattern with an over-representation of continuously up-regulated transcripts from MB to PM, MC and MM and a transcriptionally silent phase from MM to B (mirroring the total transcriptome results) with few down-regulated features. A pronounced difference was observed from the B stage on, where 90 lncRNAs are significantly upregulated in transition to S. The final maturation, from S to PMN, was the only step where the majority of lncRNA transcripts were down-regulated.\u003c/p\u003e \u003cp\u003eNeutrophil differentiation displays omics-shifts in cell compartments\u003c/p\u003e \u003cp\u003eIn previous studies, PMNs were found to have a uniform proteome configuration, with as few as 7 proteins accounting for 50% of the total protein content (Grabowski et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), a result that was recently confirmed (Sollberger et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). We hypothesized that the proteome remodeling might be detectable in our dataset as a stepwise process during differentiation.\u003c/p\u003e \u003cp\u003eOur results were consistent with previous findings of a highly condensed proteome in peripheral blood neutrophils and a coherent proteome shift during granulopoiesis, with myeloblasts showing the most heterogeneous proteome (56 proteins with a cumulative abundance of 50% in MB) that concentrates to the observed homogenous state with each maturation step (13 proteins with an abundance of 50% in PMN, Fig.\u0026nbsp;5 IIIa). In the transcriptome, we found intriguing similarities and differences (Fig.\u0026nbsp;5 IIIb). In myeloblasts, the transcriptome configuration mirrored the proteome as it showed the most heterogeneous mRNA landscape (86 transcripts accounting for 50% of the cumulative abundance). Likewise, it showed reduced heterogeneity with progressing maturation until the band stage. Thereafter, a pronounced contrast to the proteome configuration was observed in the last two stages, segmented nuclei (S) and peripheral blood neutrophils (PMN). We noticed a reversal of transcriptome homogeneity with a sudden increase in heterogeneity of transcripts (over 50 transcripts accounting for 50% of the total transcriptome).\u003c/p\u003e \u003cp\u003eAs a follow up, we determined which fraction of the proteome (Fig.\u0026nbsp;5 IVa) and transcriptome (Fig.\u0026nbsp;5 IVb) is dedicated to distinct subcellular compartments. In the myeloblast proteome, we observed the highest summed up protein intensities annotated to nucleus, secretory vesicles and cytosol. As differentiation progressed, summed up protein intensities of granule proteins continuously increased to take up the majority of the proteome, while the intensities of ER and Golgi, mitochondria and nucleus continuously decreased. Interestingly, we found that primary granules made up the largest proportion in promyelocytes and decreased slowly but steadily as maturation progresses, while secretory vesicles represented the largest fraction from the metamyelocyte stage onwards.\u003c/p\u003e \u003cp\u003eIn the transcriptome, we again found similarities to the proteome during initial maturation but observed a striking difference in the last two stages, S and PMN. Mirroring the proteome, the transcriptome initially showed a continuous decrease in intensity of features of the ER and Golgi, cytosol, mitochondria and nucleus. The features of the granule subgroups increased with maturation and showed partial sub-stadial specificity by constituting the largest fractions: primary granule features in promyelocytes, ficolin-1 granules in myelocytes and secondary granule features in metamyelocytes. Up to the B stage, the shifts corresponded to those in the proteome. In the S and PMN stages, however, the transcriptome landscape reversed its granule focused specialization and returned to a heterogeneous state, similar to that of myeloblasts.\u003c/p\u003e \u003cp\u003eGranulopoiesis transcriptome data enable refined clustering of neutrophil single-cell data\u003c/p\u003e \u003cp\u003eWe investigated whether our bulk-RNAseq data from histology defined maturation stages of granulopoiesis could lead to a more refined characterisation of single-cell (sc) neutrophil clusters. To this end, we applied a novel data integration method called \u0026ldquo;Φ-Space\u0026rdquo; (Mao, Deng and L\u0026ecirc; Cao, \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2024\u003c/span\u003e) to a dataset of single cell neutrophil sequences from COVID-19 patients ((Schulte-Schrepping et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2020\u003c/span\u003e)). The authors identified distinctive clusters of neutrophils in severe COVID-19 patients and utilized marker genes from mouse granulopoiesis to categorize three crude maturation stages: pre-, pro- and mature neutrophils (Kwok et al., 2020). As an alternative, we used \u0026ldquo;Φ-Space\u0026rdquo; to integrate the sc based maturation-scores with our bulk-data of histology-defined maturation stages (analysis workflow in sFig. 2 IIIa). A heatmap based on these maturation scores (Fig.\u0026nbsp;5 Va) showed that the single cell clusters defined by Schulte-Schrepping et al. overwhelmingly matched the individual transcriptional profiles of our histologically defined maturation stages. Of note, we observed an overlapping profile match between the \u0026ldquo;ribosome high\u0026rdquo; cluster and our profiles of band-stage neutrophils and myeloblasts. We further visualized these maturation scores using Uniform Manifold Approximation and Projection (Fig.\u0026nbsp;5V b, UMAP, (McInnes et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2018\u003c/span\u003e)). We re-clustered the single neutrophils based on their maturation scores which gave a better delineation of cell clusters compared to the original mapping (sFig. 2 IIIb).\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eHere, we present a comprehensive molecular atlas of 7 defined developmental stages of human granulopoiesis on transcriptome, proteome, and miRNome levels.\u003c/p\u003e \u003cp\u003eOur data characterize granulopoiesis as a finely coordinated process of extreme cellular differentiation from proliferating myeloblasts to postmitotic differentiated polymorphonuclear neutrophils. We observed a highly specialized proteome landscape that was reduced to a metabolic minimum with abundant effector proteins. While these changes occur in a stepwise manner, several neutrophil-specific proteins (primary and ficolin-1 granule components) are already detectable in myeloblasts, indicating that the developmental trajectory towards neutrophils may be initiated already upstream of myeloblasts (Kwok et al., 2020). Most changes in protein abundance occurred in early phases, whereas the proteome remained fairly stable at the MM, B, and PMN stage. By contrast, the granulopoietic transcriptome showed changes in each development step and was divided into two distinct phases. The first transcriptome phase, from MB to MM, coincides with a proteome landscape rich in ribosomes, and shows a strong correlation between mRNA and protein features. The second phase, after the MM stage, takes place in a static proteome with minimal ribosomal content and shows decoupling of mRNA and protein features.\u003c/p\u003e \u003cp\u003eThis dual wave pattern of transcripts in neutrophil precursor development appears to contradict the \u0026ldquo;targeting by time\u0026rdquo; hypothesis (Le Cabec et al., 1996), that proposed that granule synthesis may be a time-controlled, serial process for each granule subfraction individually. We find evidence of parallel processing from the promyelocyte stage onwards. Interestingly, there is a time-lag between transcription and translation for the secondary granule subfraction. The mechanisms explaining this apparent discordance remain to be discovered.\u003c/p\u003e \u003cp\u003eIn parallel with the up-regulation of granule proteins, differentiation is accompanied by loss of cell cycle machinery and DNA replication, with diminished cytoplasmic ribosomes and mitochondrial matrix. This reflects transition from replicative to non-replicative cells and coincides with morphological changes of neutrophils (Robinson et al., 1999).\u003c/p\u003e \u003cp\u003eAn unexpected finding was the loss of nonsense-mediated decay (NMD). These findings have implications to interpret proteomic imbalances in genetic diseases. For example, we previously found upregulation of members of the NMD machinery in the neutrophil proteome of patients with severe congenital neutropenia (Grabowski et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Our current study highlights the need to distinguish between signatures of immaturity versus compensatory mechanisms related to disease states.\u003c/p\u003e \u003cp\u003eOur analysis of non-coding RNAs (lncRNA and miRNA) highlights previously unrecognized features. Even though, due to technical reasons, we could not obtain representative miRNA samples from promyelocytes, the remaining stages revealed a characteristic pattern. Two independent strategies were used to integrate the miRNome with transcriptome and proteome.\u003c/p\u003e \u003cp\u003eFirst, when we compared enrichment terms, we could infer pathways and processes that miRNAs may inhibit (Mathieu and Ruohola-Baker, 2013), including VEGF, Ras, mTOR, cell-cell adhesion, extracellular trap formation) or silence during maturation (SRP, NMD, ribosomes). Second, we determined miRNA binding site enrichment in DE transcripts and performed post-hoc enrichment analysis, allowing us to infer activity of specific miRNA transcripts on major biological themes. miRNAs with high expression in myeloblasts (miR-19a/b-3p, miR-106a-5p family, miR-125a-5p) showed binding-site enrichment for mRNAs involved in myeloid differentiation (NBEAL2, NEDD9, TET2, STAT3). As miR-19 was found to control NFkB expression (Gantier et al., 2012), we speculate that miR-19 may repress neutrophil differentiation in myeloblasts via inhibition of NFkB.\u003c/p\u003e \u003cp\u003emiR-193a on the other hand, previously shown to increase during granulopoiesis and to down-regulate transcripts of the essential SWI/SNF-complex, may facilitate myeloid leukemogenesis in conjunction with LINC00152 and CDK9 (Witzel et al., 2017) (Zhang and Tao, 2019).\u003c/p\u003e \u003cp\u003eLong non-coding (lnc) RNA biology may have specific effects on a wide variety of cellular processes, from cell structure to chromatin architecture and transcription (Mattick et al., 2023). We found that lncRNA transcripts were upregulated from early to late phases of neutrophil differentiation in the bone marrow. When mature neutrophils leave the bone marrow (transition from S to PMN stage) lncRNA species are down-regulated. One particular class of lncRNA are lncRNA-type small nucleolar RNA host genes (SNHG), known to play pivotal roles in cancer progression (Zheng et al., 2023).\u003c/p\u003e \u003cp\u003eSNHGs are overexpressed in early granulopoiesis (MB to MC), in line with the notion that their prognostic value in AML may may stem from their activity in less differentiated cell states (Shi, Ding and Lu, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Interestingly, a recent study found that the overexpression of ITGB2-TAS1 (ITGB2 antisense RNA 1) was associated with adverse prognosis in AML. In our study, ITGB2-AS1 was under-expressed in MB to MC, with increasing expression in the MM stage, an opposite expression pattern to the aforementioned SNHGs.\u003c/p\u003e \u003cp\u003eOur dataset enabled us to conduct an in-depth analysis of transcription factor signatures in the developing granulopoiesis transcriptome. For example, Lactotransferrin (LTF) is a prominent constituent of secondary granules in PMN with antimicrobial and immunomodulatory functions (Legrand, 2016). Its isoform, delta-lactoferrin (ΔLf), acts as a transcription factor with antiproliferative properties (Mariller et al., 2012). A recent single-cell analysis in mice found a signature of LTF activity in early granulopoiesis (Grieshaber-Bouyer et al., 2021). Our study finds LTF to be transcriptionally active and highly expressed from myeloblasts to metamyelocytes, peaking in myeloblasts and promyelocytes. Of note, we also discovered a number of transcription factors in mature stages of neutrophil development that were not previously linked to neutrophil biology. The last granulocytic stage differentiation, from S to PMN, showed a large and exclusive transcription factor signature. The two most enriched, THAP11 and THAP5, belong to the thanatos-associated protein family and have not been described yet in granulopoiesis. THAP5 is described to cause cell cyle arrest (Balakrishnan et al., \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2009\u003c/span\u003e) while THAP11 (formerly RONIN) is a epigenetic regulator and essential for pluripotency in embryonic stem cells (Dejosez et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2008\u003c/span\u003e) that downregulates cMYC (Zhu et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2009\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eOur study points to a general limitation in the interpretation of proteomic data using the \u0026lsquo;proteomic ruler\u0026rsquo; approach to infer protein copy numbers from histone intensities (Grabowski et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Sollberger et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Conceptually, this is based on the assumption that the ratio of histones to DNA remains stable and can serve as a constant that allows the inference of protein copy numbers in proteome samples with known cell numbers (Wiśniewski et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). However, we found that histone expression at the proteome level was strongly induced during granulopoiesis at both the mRNA and protein levels, reaching a maximum in PMN. Accordingly, the \u0026lsquo;proteomic ruler' approach cannot easily be applied since it would introduce a bias in the interpretation of protein abundances.\u003c/p\u003e\u003cp\u003eIn summary, we here describe intricate processes of cellular specialization from proliferating myeloblasts to postmitotic differentiated polymorphonuclear neutrophils. Whereas previous studies have focused on either transcriptome (Grassi et al., 2018), proteome (Hoogendijk et al., 2019), or miRNome (Larsen et al., \u003cspan class=\"CitationRef\"\u003e2012\u003c/span\u003e) individually, joint analysis of multiomic data from matched samples remains challenging. Our website, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ewww.granulopoiesis.com\u003c/span\u003e\u003c/span\u003e enables full access to raw and processed data, the complete analysis scripts and interactive data analysis for further use.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch3\u003eFunding\u003c/h3\u003e\n\u003cp\u003eSH: F\u0026ouml;rderprogramm f\u0026uuml;r Forschung und Lehre (F\u0026ouml;FoLe, Project number 1052)\u003cbr\u003e\u0026nbsp;RZ: SFB1123 Atherosclerosis\u003c/p\u003e\n\u003cp\u003eCK: German Center for Child and Adolescent Health (DZKJ) Munich Site, Sonderforschungsbereich TRR332 \u0026quot;Neutrophils: origin, fate \u0026amp; function\u0026quot;\u003c/p\u003e\n"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAllaire, J. (2012) \u0026lsquo;RStudio: integrated development environment for R\u0026rsquo;, \u003cem\u003eBoston, MA\u003c/em\u003e, 770(394), pp. 165\u0026ndash;171. Available at: https://www.academia.edu/download/49999670/IntR_-_Interactive_GUI_for_R20161031-11791-s0u02o.pdf#page=14.\u003c/li\u003e\n\u003cli\u003eAparicio-Puerta, E. \u003cem\u003eet al.\u003c/em\u003e (2023) \u0026lsquo;miEAA 2023: updates, new functional microRNA sets and improved enrichment visualizations\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 51(W1), pp. W319\u0026ndash;W325. Available at: https://doi.org/10.1093/nar/gkad392.\u003c/li\u003e\n\u003cli\u003eBalakrishnan, M.P. \u003cem\u003eet al.\u003c/em\u003e (2009) \u0026lsquo;THAP5 is a human cardiac-specific inhibitor of cell cycle that is cleaved by the proapoptotic Omi/HtrA2 protease during cell death\u0026rsquo;, \u003cem\u003eAmerican journal of physiology. Heart and circulatory physiology\u003c/em\u003e, 297(2), pp. H643\u0026ndash;53. Available at: https://doi.org/10.1152/ajpheart.00234.2009.\u003c/li\u003e\n\u003cli\u003eBarrett, T. \u003cem\u003eet al.\u003c/em\u003e (2013) \u0026lsquo;NCBI GEO: archive for functional genomics data sets--update\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 41(Database issue), pp. D991\u0026ndash;5. Available at: https://doi.org/10.1093/nar/gks1193.\u003c/li\u003e\n\u003cli\u003eBodein, A. \u003cem\u003eet al.\u003c/em\u003e (2022) \u0026lsquo;timeOmics: an R package for longitudinal multi-omics data integration\u0026rsquo;, \u003cem\u003eBioinformatics \u003c/em\u003e, 38(2), pp. 577\u0026ndash;579. Available at: https://doi.org/10.1093/bioinformatics/btab664.\u003c/li\u003e\n\u003cli\u003eBrunson, J.C. (2020) \u0026lsquo;ggalluvial: Layered Grammar for Alluvial Plots\u0026rsquo;, \u003cem\u003eJournal of open source software\u003c/em\u003e, 5(49). Available at: https://doi.org/10.21105/joss.02017.\u003c/li\u003e\n\u003cli\u003eBurn, G.L. \u003cem\u003eet al.\u003c/em\u003e (2021) \u0026lsquo;The neutrophil\u0026rsquo;, \u003cem\u003eImmunity\u003c/em\u003e, 54(7), pp. 1377\u0026ndash;1391. Available at: https://doi.org/10.1016/j.immuni.2021.06.006.\u003c/li\u003e\n\u003cli\u003eChoi, J. \u003cem\u003eet al.\u003c/em\u003e (2019) \u0026lsquo;Stemformatics: visualize and download curated stem cell data\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 47(D1), pp. D841\u0026ndash;D846. Available at: https://doi.org/10.1093/nar/gky1064.\u003c/li\u003e\n\u003cli\u003eConnerty, P. and Lock, R.B. (2023) \u0026lsquo;The tip of the iceberg-The roles of long noncoding RNAs in acute myeloid leukemia\u0026rsquo;, \u003cem\u003eWiley interdisciplinary reviews. RNA\u003c/em\u003e, 14(6), p. e1796. Available at: https://doi.org/10.1002/wrna.1796.\u003c/li\u003e\n\u003cli\u003eCox, J. and Mann, M. (2008) \u0026lsquo;MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification\u0026rsquo;, \u003cem\u003eNature biotechnology\u003c/em\u003e, 26(12), pp. 1367\u0026ndash;1372. Available at: https://doi.org/10.1038/nbt.1511.\u003c/li\u003e\n\u003cli\u003eDejosez, M. \u003cem\u003eet al.\u003c/em\u003e (2008) \u0026lsquo;Ronin is essential for embryogenesis and the pluripotency of mouse embryonic stem cells\u0026rsquo;, \u003cem\u003eCell\u003c/em\u003e, 133(7), pp. 1162\u0026ndash;1174. Available at: https://doi.org/10.1016/j.cell.2008.05.047.\u003c/li\u003e\n\u003cli\u003eDeutsch, E.W. \u003cem\u003eet al.\u003c/em\u003e (2023) \u0026lsquo;The ProteomeXchange consortium at 10 years: 2023 update\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 51(D1), pp. D1539\u0026ndash;D1548. Available at: https://doi.org/10.1093/nar/gkac1040.\u003c/li\u003e\n\u003cli\u003eDobin, A. \u003cem\u003eet al.\u003c/em\u003e (2013) \u0026lsquo;STAR: ultrafast universal RNA-seq aligner\u0026rsquo;, \u003cem\u003eBioinformatics \u003c/em\u003e, 29(1), pp. 15\u0026ndash;21. Available at: https://doi.org/10.1093/bioinformatics/bts635.\u003c/li\u003e\n\u003cli\u003eDonadieu, J. and Bellann\u0026eacute;-Chantelot, C. (2022) \u0026lsquo;Genetics of severe congenital neutropenia as a gateway to personalized therapy\u0026rsquo;, \u003cem\u003eHematology\u003c/em\u003e, 2022(1), pp. 658\u0026ndash;665. Available at: https://doi.org/10.1182/hematology.2022000392.\u003c/li\u003e\n\u003cli\u003eElghetany, M.T. \u003cem\u003eet al.\u003c/em\u003e (2004) \u0026lsquo;Flow cytometric study of neutrophilic granulopoiesis in normal bone marrow using an expanded panel of antibodies: correlation with morphologic assessments\u0026rsquo;, \u003cem\u003eJournal of clinical laboratory analysis\u003c/em\u003e, 18(1), pp. 36\u0026ndash;41. Available at: https://doi.org/10.1002/jcla.20001.\u003c/li\u003e\n\u003cli\u003eFaurschou, M. and Borregaard, N. (2003) \u0026lsquo;Neutrophil granules and secretory vesicles in inflammation\u0026rsquo;, \u003cem\u003eMicrobes and infection / Institut Pasteur\u003c/em\u003e, 5(14), pp. 1317\u0026ndash;1327. Available at: https://doi.org/10.1016/j.micinf.2003.09.008.\u003c/li\u003e\n\u003cli\u003eGene Ontology Consortium \u003cem\u003eet al.\u003c/em\u003e (2023) \u0026lsquo;The Gene Ontology knowledgebase in 2023\u0026rsquo;, \u003cem\u003eGenetics\u003c/em\u003e, 224(1). Available at: https://doi.org/10.1093/genetics/iyad031.\u003c/li\u003e\n\u003cli\u003eGrabowski, P. \u003cem\u003eet al.\u003c/em\u003e (2019) \u0026lsquo;Proteome analysis of human neutrophil granulocytes from patients with monogenic disease using data-independent acquisition\u0026rsquo;, \u003cem\u003eMolecular \u0026amp; cellular proteomics: MCP\u003c/em\u003e [Preprint]. Available at: https://doi.org/10.1074/mcp.RA118.001141.\u003c/li\u003e\n\u003cli\u003eGu, Z. (2022) \u0026lsquo;Complex heatmap visualization\u0026rsquo;, \u003cem\u003eiMeta\u003c/em\u003e, 1(3). Available at: https://doi.org/10.1002/imt2.43.\u003c/li\u003e\n\u003cli\u003eHuber, W. \u003cem\u003eet al.\u003c/em\u003e (2002) \u0026lsquo;Variance stabilization applied to microarray data calibration and to the quantification of differential expression\u0026rsquo;, \u003cem\u003eBioinformatics \u003c/em\u003e, 18 Suppl 1, pp. S96\u0026ndash;104. Available at: https://doi.org/10.1093/bioinformatics/18.suppl_1.s96.\u003c/li\u003e\n\u003cli\u003eKeenan, A.B. \u003cem\u003eet al.\u003c/em\u003e (2019) \u0026lsquo;ChEA3: transcription factor enrichment analysis by orthogonal omics integration\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 47(W1), pp. W212\u0026ndash;W224. Available at: https://doi.org/10.1093/nar/gkz446.\u003c/li\u003e\n\u003cli\u003eKern, F. \u003cem\u003eet al.\u003c/em\u003e (2021) \u0026lsquo;miRTargetLink 2.0-interactive miRNA target gene and target pathway networks\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 49(W1), pp. W409\u0026ndash;W416. Available at: https://doi.org/10.1093/nar/gkab297.\u003c/li\u003e\n\u003cli\u003eLangmead, B. and Salzberg, S.L. (2012) \u0026lsquo;Fast gapped-read alignment with Bowtie 2\u0026rsquo;, \u003cem\u003eNature methods\u003c/em\u003e, 9(4), pp. 357\u0026ndash;359. Available at: https://doi.org/10.1038/nmeth.1923.\u003c/li\u003e\n\u003cli\u003eLarsen, M.T. \u003cem\u003eet al.\u003c/em\u003e (2012) \u0026lsquo;miRNA Profiling in Human Granulopoiesis and in Activated Neutrophils From Skin Window\u0026rsquo;, \u003cem\u003eBlood\u003c/em\u003e, 120(21), pp. 2136\u0026ndash;2136. Available at: https://doi.org/10.1182/blood.V120.21.2136.2136.\u003c/li\u003e\n\u003cli\u003eMakkar, H. \u003cem\u003eet al.\u003c/em\u003e (2023) \u0026lsquo;Acute myeloid leukemia: novel mutations and their clinical implications\u0026rsquo;, \u003cem\u003eAmerican journal of blood research\u003c/em\u003e, 13(1), pp. 12\u0026ndash;27. Available at: https://www.ncbi.nlm.nih.gov/pubmed/36937458.\u003c/li\u003e\n\u003cli\u003eMao, J., Deng, Y. and L\u0026ecirc; Cao, K.-A. (2024) \u0026lsquo;\u0026Phi;-Space: Continuous phenotyping of single-cell multi-omics data\u0026rsquo;, \u003cem\u003ebioRxiv\u003c/em\u003e. Available at: https://doi.org/10.1101/2024.06.19.599787.\u003c/li\u003e\n\u003cli\u003eMcInnes, L. \u003cem\u003eet al.\u003c/em\u003e (2018) \u0026lsquo;UMAP: Uniform Manifold Approximation and Projection\u0026rsquo;, \u003cem\u003eJournal of open source software\u003c/em\u003e, 3(29), p. 861. Available at: https://doi.org/10.21105/joss.00861.\u003c/li\u003e\n\u003cli\u003eMcInnes, L., Healy, J. and Melville, J. (2018) \u0026lsquo;UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction\u0026rsquo;, \u003cem\u003earXiv [stat.ML]\u003c/em\u003e. Available at: http://arxiv.org/abs/1802.03426 (Accessed: 18 September 2024).\u003c/li\u003e\n\u003cli\u003eMeissner, F. \u003cem\u003eet al.\u003c/em\u003e (2013) \u0026lsquo;Direct proteomic quantification of the secretome of activated immune cells\u0026rsquo;, \u003cem\u003eScience\u003c/em\u003e, 340(6131), pp. 475\u0026ndash;478. Available at: https://doi.org/10.1126/science.1232578.\u003c/li\u003e\n\u003cli\u003eMilacic, M. \u003cem\u003eet al.\u003c/em\u003e (2024) \u0026lsquo;The Reactome Pathway Knowledgebase 2024\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 52(D1), pp. D672\u0026ndash;D678. Available at: https://doi.org/10.1093/nar/gkad1025.\u003c/li\u003e\n\u003cli\u003eNg, M.S.F. \u003cem\u003eet al.\u003c/em\u003e (2024) \u0026lsquo;Deterministic reprogramming of neutrophils within tumors\u0026rsquo;, \u003cem\u003eScience\u003c/em\u003e, 383(6679), p. eadf6493. Available at: https://doi.org/10.1126/science.adf6493.\u003c/li\u003e\n\u003cli\u003eR Foundation for Statistical Computing, R. (2018) \u0026lsquo;R: A language and environment for statistical computing\u0026rsquo;, \u003cem\u003eRA Lang Environ Stat Comput\u003c/em\u003e [Preprint].\u003c/li\u003e\n\u003cli\u003eRitchie, M.E. \u003cem\u003eet al.\u003c/em\u003e (2015) \u0026lsquo;limma powers differential expression analyses for RNA-sequencing and microarray studies\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 43(7), p. e47. Available at: https://doi.org/10.1093/nar/gkv007.\u003c/li\u003e\n\u003cli\u003eRoehr, J.T., Dieterich, C. and Reinert, K. (2017) \u0026lsquo;Flexbar 3.0 - SIMD and multicore parallelization\u0026rsquo;, \u003cem\u003eBioinformatics \u003c/em\u003e, 33(18), pp. 2941\u0026ndash;2942. Available at: https://doi.org/10.1093/bioinformatics/btx330.\u003c/li\u003e\n\u003cli\u003eRohart, F. \u003cem\u003eet al.\u003c/em\u003e (2017) \u0026lsquo;mixOmics: An R package for \u0026rsquo;omics feature selection and multiple data integration\u0026rsquo;, \u003cem\u003ePLoS computational biology\u003c/em\u003e, 13(11), p. e1005752. Available at: https://doi.org/10.1371/journal.pcbi.1005752.\u003c/li\u003e\n\u003cli\u003eR\u0026oslash;rvig, S. \u003cem\u003eet al.\u003c/em\u003e (2013) \u0026lsquo;Proteome profiling of human neutrophil granule subsets, secretory vesicles, and cell membrane: correlation with transcriptome profiling of neutrophil precursors\u0026rsquo;, \u003cem\u003eJournal of leukocyte biology\u003c/em\u003e, 94(4), pp. 711\u0026ndash;721. Available at: https://doi.org/10.1189/jlb.1212619.\u003c/li\u003e\n\u003cli\u003eSchulte-Schrepping, J. \u003cem\u003eet al.\u003c/em\u003e (2020) \u0026lsquo;Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell Compartment\u0026rsquo;, \u003cem\u003eCell\u003c/em\u003e, 182(6), pp. 1419\u0026ndash;1440.e23. Available at: https://doi.org/10.1016/j.cell.2020.08.001.\u003c/li\u003e\n\u003cli\u003eShi, J., Ding, W. and Lu, H. (2020) \u0026lsquo;Identification of long non-coding RNA SNHG family as promising prognostic biomarkers in acute myeloid leukemia\u0026rsquo;, \u003cem\u003eOncoTargets and therapy\u003c/em\u003e, 13, pp. 8441\u0026ndash;8450. Available at: https://doi.org/10.2147/OTT.S265853.\u003c/li\u003e\n\u003cli\u003eSollberger, G. \u003cem\u003eet al.\u003c/em\u003e (2024) \u0026lsquo;Quantitative proteomics reveals tissue-specific, infection-induced and species-specific neutrophil protein signatures\u0026rsquo;, \u003cem\u003eScientific reports\u003c/em\u003e, 14(1), p. 5966. Available at: https://doi.org/10.1038/s41598-024-56163-6.\u003c/li\u003e\n\u003cli\u003eSoutschek, M. \u003cem\u003eet al.\u003c/em\u003e (2022) \u0026lsquo;enrichMiR predicts functionally relevant microRNAs based on target collections\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 50(W1), pp. W280\u0026ndash;W289. Available at: https://doi.org/10.1093/nar/gkac395.\u003c/li\u003e\n\u003cli\u003eStrydom, N. and Rankin, S.M. (2013) \u0026lsquo;Regulation of circulating neutrophil numbers under homeostasis and in disease\u0026rsquo;, \u003cem\u003eJournal of innate immunity\u003c/em\u003e, 5(4), pp. 304\u0026ndash;314. Available at: https://doi.org/10.1159/000350282.\u003c/li\u003e\n\u003cli\u003e\u003cem\u003eThe Association to Function Knowledge Portal\u003c/em\u003e (2024). Available at: a2fkp.org.\u003c/li\u003e\n\u003cli\u003eThul, P.J. \u003cem\u003eet al.\u003c/em\u003e (2017) \u0026lsquo;A subcellular map of the human proteome\u0026rsquo;, \u003cem\u003eScience\u003c/em\u003e, 356(6340). Available at: https://doi.org/10.1126/science.aal3321.\u003c/li\u003e\n\u003cli\u003eUniProt Consortium (2023) \u0026lsquo;UniProt: the Universal Protein Knowledgebase in 2023\u0026rsquo;, \u003cem\u003eNucleic acids research\u003c/em\u003e, 51(D1), pp. D523\u0026ndash;D531. Available at: https://doi.org/10.1093/nar/gkac1052.\u003c/li\u003e\n\u003cli\u003eWang, L., Wang, S. and Li, W. (2012) \u0026lsquo;RSeQC: quality control of RNA-seq experiments\u0026rsquo;, \u003cem\u003eBioinformatics \u003c/em\u003e, 28(16), pp. 2184\u0026ndash;2185. Available at: https://doi.org/10.1093/bioinformatics/bts356.\u003c/li\u003e\n\u003cli\u003eWard, J.H., Jr (1963) \u0026lsquo;Hierarchical grouping to optimize an objective function\u0026rsquo;, \u003cem\u003eJournal of the American Statistical Association\u003c/em\u003e, 58(301), pp. 236\u0026ndash;244. Available at: https://doi.org/10.1080/01621459.1963.10500845.\u003c/li\u003e\n\u003cli\u003e\u003cem\u003eWGCNA: an R package for weighted correlation network analysis\u003c/em\u003e (no date) \u003cem\u003eBioMed Central\u003c/em\u003e. Available at: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-559 (Accessed: 3 May 2024).\u003c/li\u003e\n\u003cli\u003eWiśniewski, J.R. \u003cem\u003eet al.\u003c/em\u003e (2014) \u0026lsquo;A \u0026ldquo;proteomic ruler\u0026rdquo; for protein copy number and concentration estimation without spike-in standards\u0026rsquo;, \u003cem\u003eMolecular \u0026amp; cellular proteomics: MCP\u003c/em\u003e, 13(12), pp. 3497\u0026ndash;3506. Available at: https://doi.org/10.1074/mcp.M113.037309.\u003c/li\u003e\n\u003cli\u003eWu, P. \u003cem\u003eet al.\u003c/em\u003e (2021) \u0026lsquo;The trans-omics landscape of COVID-19\u0026rsquo;, \u003cem\u003eNature communications\u003c/em\u003e, 12(1), p. 4543. Available at: https://doi.org/10.1038/s41467-021-24482-1.\u003c/li\u003e\n\u003cli\u003eZhu, C.-Y. \u003cem\u003eet al.\u003c/em\u003e (2009) \u0026lsquo;Cell growth suppression by thanatos-associated protein 11(THAP11) is mediated by transcriptional downregulation of c-Myc\u0026rsquo;, \u003cem\u003eCell death and differentiation\u003c/em\u003e, 16(3), pp. 395\u0026ndash;405. Available at: https://doi.org/10.1038/cdd.2008.160.\u003c/li\u003e\n\u003cli\u003eWong, A. C., Wong, J. J., Rasko, J. E., \u0026amp; Schmitz, U. (2024). SpliceWiz: interactive analysis and visualization of alternative splicing in R. \u003cem\u003eBriefings in Bioinformatics\u003c/em\u003e, \u003cem\u003e25\u003c/em\u003e(1).\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Materials and methods","content":"\u003ch2\u003eSample collection\u003c/h2\u003e\n\u003cp\u003eThe study protocol was approved by the Ethics Committee of the Medical Faculty of the Ludwig Maximilian University of Munich (Reg. No. 19-685 2). Four healthy German volunteers over the age of 25, who had not previously been diagnosed with cancer or a serious illness, were included in the study. All healthy volunteers had signed a written consent form prior to the collection of 300 ml of bone marrow from the iliac crest. The bone marrow was anticoagulated in EDTA (Thermo Scientific 005500, 0.1 ml EDTA / 50 ml bone marrow) and subjected to negative selection to remove the non-granulocytic cell fractions using Miltenyi's neutrophil isolation kit (Miltenyi Biotec, 130-104-434) according to the manufacturer's protocol.\u003c/p\u003e\n\u003cp\u003eTo isolate defined precursors of granulopoiesis, we created a novel flow cytometry panel (sFig. 3 Ia) to detect 6 substages of human granulopoiesis in human bone marrow (sFig. 3 Ib) based on published surface antigens (Elghetany et al, 2004; Fajtova and Babusikova, 2010; Gorczyca et al, 2011; Gustafson et al, 2015).\u003c/p\u003e\n\u003cp\u003eThe isolated cells were dispersed in 1 ml PBS (Gibco, 14190-169) containing 10 % FCS (Value FBS South-American, ThermoFisher A5256701) and stained with a ten-color antibody panel for sorting using a BD FACSAria III cell sorter with blue, violet, red and yellow/green laser. After exclusion of dead cells and monocytes with BD Horizon Fixable Viability Stain (Red 780, BD 565388) and CD14 (APC/Fire750, Biolegend 301854), the major maturation stages were differentiated using CD117+ (PE, BB700, BD 566548) and CD49d+ (PE, BD 555503). Double positive cells represented the earliest stages with myeloblasts (MB) and promyelocytes (PM), further differentiated by CD34+ (BV421, BD 562577) for MB and CD34- for PM. Intermediate maturation stages were selected via the CD117- and CD49d+ population with myelocytes (MC) CD11b+ (FITC, Biolegend 301330) and metamyelocytes CD101+ (BV605, BD 747548). The terminal maturation stages (CD117- / CD49d-) were differentiated as band-stage (B) via CD24+ (PE-Dazzle594, Biolegend 311134) / CD16 dim (PE-Cy7, BD557744) / CD10- (CD10 BV711 BD 740770), and as a segmented stage via CD24+/CD16+/CD10+. To represent the terminal stage from peripheral blood, neutrophil granulocytes were isolated from peripheral venous blood using the neutrophil isolation kit in combination with an erythrocyte depletion cocktail (Miltenyi Biotec, 130-098-196).\u003cbr\u003e Light microscopic preparations from each collected fraction showed \u0026gt;95% pure cell populations with the respective stage-typical morphologies (sFig. 3 Ib).\u003c/p\u003e\n\u003cp\u003eSequencing\u003c/p\u003e\n\u003cp\u003eFor sequencing, we aimed for a cell count of 1x10^6 cells per sample and prioritized the collection of transcriptome and proteome samples over miRNA samples. A complete overview of the samples and cell numbers can be found in sTable 1. Cells were collected in aliquots of 0.1x10^6 - 1x10^6 cells, washed once in PBS, dissolved in storage reagents optimized for downstream analysis and frozen at -80 C. For proteome analysis, samples were stored in 5ul 25x protease inhibitor (Roche #04693159001). Transcriptome samples were stored in RLT-plus buffer (Qiagen #1053393) with the addition of 10 ul/1ml beta-mercaptoethanol (Thermo Scientific #21985023). To preserve the microRNA, the cells were dissolved in 750 ul TRIZOL-LS (Invitrogen #15596026).\u003c/p\u003e\n\u003cp\u003eThe sequencing of transcriptomes (mRNA) and microRNA was performed by a commercial partner, Azenta life sciences (Im Leuschnerpark 1B, 64347 Griesheim, Germany). RNA was extracted using the Qiagen RNeasy Mini Kit according to the manufacturer's protocol and sequenced on the Illumina NovaSeq 6000 with a 2x150 pair-end configuration. Raw BCL data were converted to FASTQ files and demultiplexed using Illumina bcl2fastq 2.20 software. Detailed protocols can be found in the extended Methods.\u003c/p\u003e\n\u003cp\u003eProteome sequencing was done with a nanoflow liquid chromatography (Thermo Scientific, EASY-nLC) coupled to a Orbitrap Exploris 480 instrument with a nano-electrospray ion source (Thermo Fisher Scientific). Protein sample preparation was done as described previously (Meissner \u003cem\u003eet al.\u003c/em\u003e, 2013), with Trypsin (Promega #VA9000) and LysC (FUJIFILM Wako #121-05063) digestion. For a detailed protocol, please refer to the extended Methods.\u003c/p\u003e\n\u003cp\u003eAfter analyzing the raw sequencing data, described in depth in the extended methods, we applied stringent quality criteria for feature inclusion that resulted in 8,432 transcriptome features, 3,156 proteins, and 283 microRNAs that were reliably quantified in at least two samples per stage. The complete list of identified features and their expression values (stage median expression values) can be found in the downloadable supplement tables sT2 (transcriptome), sT3 (proteome) and sT4 (miRNnome).\u003c/p\u003e\n\u003ch2\u003eBioinformatic analysis \u003c/h2\u003e\n\u003ch3\u003eFAIR data commitment \u003c/h3\u003e\n\u003cp\u003eWe strive to ensure the highest standards of FAIR principles for data sharing (Findability, Accessibility, Interoperability, and Reuse) and appreciate that a major contribution of this work is to make the datasets and analysis workflows available to the community. To this end, we published our data in various forms on different platforms and create a dedicated website for data and analsis sharing.\u003c/p\u003e\n\u003cp\u003eRaw data are made available via ProteomeXchange (Deutsch \u003cem\u003eet al.\u003c/em\u003e, 2023) and Gene Expression Omnibus (Barrett \u003cem\u003eet al.\u003c/em\u003e, 2013). To enable comparisons with other published transcriptome data, we integrated our data with stemformatics.org (Choi \u003cem\u003eet al.\u003c/em\u003e, 2019), an interactive portal of high quality stem cell datasets. \u003c/p\u003e\n\u003cp\u003eSince our results included a wealth of data beyond the transcriptome, we created a dedicated website, www.granulopoiesis.com, to allow readers to download all intermediate and final data along with the complete, reproducible code used in all steps of analysis, from raw counts to analysis. We have also created several ways to query our data directly on granulopoiesis.com. The expression of single and grouped features of transcriptome, proteome and miRNome can be visualized directly on the website.\u003c/p\u003e\n\u003ch3\u003eRaw data processing \u003c/h3\u003e\n\u003cp\u003eHere we describe the raw data processing only in short, while a detailed description can be found in the Extended Materials and Methods. Proteome raw data were processed with MaxQuant (version 1.6.2.1, (Cox and Mann, 2008)) with enabled ‘match between runs' option to receive a Uniprot annotated a maxLFQ based intensity table. Transcriptome (mRNA) raw sequence reads were aligned using STAR (v2.7.9a_2021-06-25, (Dobin \u003cem\u003eet al.\u003c/em\u003e, 2013)) with RSeQC for quality control (Wang, Wang and Li, 2012). Data were processed as counts for differential expression analysis and tpm (transcripts per million) for all other forms of analysis (clustering, trajectories). MicroRNA raw data were processed with Flexbar (v3.5.0, (Roehr, Dieterich and Reinert, 2017)) and Bowtie2 (v2.4.4, (Langmead and Salzberg, 2012)) and exported as a counts table.\u003c/p\u003e\n\u003cp\u003eAfter raw data processing, all further data analysis was performed using R (4.3.2, 2023-10-31 (R Foundation for Statistical Computing, 2018)) and R-Studio (2023.12.0+369 \"Ocean Storm\" (Allaire, 2012)).\u003c/p\u003e\n\u003ch3\u003eFeature and sample selection \u003c/h3\u003e\n\u003cp\u003eIn order to ensure consideration of highly stringent features only, we applied strict inclusion criteria at each omics level. For the transcriptome, from the original 61,552 features, we selected only those that had valid expression values in at least 3 samples for each stage, with the exception of the proteome, where at least 2 samples were required to accommodate the fact that of the PMN population, only two samples were successfully sequenced. Valid expression values were defined as more than 10 counts or a mean tpm value of 0.28 adjusted for library size (average 37 million) (Law et al., 2016), resulting in a total of 8,432 transcriptome features. From the total 3503 quantified proteins, we only selected features that were quantified with at least 2 peptides per protein. Furthermore, we only included proteins that were quantified in at least 2 samples per stage and excluded known contaminants (MaxQuant contaminant list and all keratins), resulting in 3,156 proteins.\u003cbr\u003e For microRNA, the raw data contained 2950 features. Filtering for features with at least 10 counts in at least 2 samples per stage resulted in 283 stringent features for downstream analysis. \u003c/p\u003e\n\u003cp\u003eAll individual data sets were analyzed for outliers using the signal intensity and principal component analysis (PCA) before unification of technical repeat samples via median expression values (sFig. 4-6). No batch effects were detected in any of the datasets and no correction methods were applied. For normalization, we applied VSN (variant stabilization normalization) to each dataset (Huber \u003cem\u003eet al.\u003c/em\u003e, 2002).\u003c/p\u003e\n\u003ch3\u003eSample analysis\u003c/h3\u003e\n\u003cp\u003eA detailed description of the analysis steps and packages used can be found in the Extended Materials and Methods as well as the provided markdown scripts. Here we provide only a short summary of the applied methods and packages. Unless specified, base-R methods were used. For clustering analysis of each omics dataset, we tested euclidean-based wardD and k-means and plotted the results on the first two principal components (Fig. 2 I-IIIa). Stage specific features were selected using PLS-DA of \u003cem\u003emicOmics\u003c/em\u003e (Rohart \u003cem\u003eet al.\u003c/em\u003e, 2017) and plotted with the \u003cem\u003eComplexHeatmap\u003c/em\u003e package (Gu, 2022). Differential expression analysis for each dataset was done with Limma (Ritchie \u003cem\u003eet al.\u003c/em\u003e, 2015). As significance thresholds, we applied a Benjamini Hochberg corrected p.value \u0026lt; 0.05 and a fold change equal or below/above +/-1. For biological-term enrichment, we employed ClusterProfiler to probe for GO terms (Gene Ontology Consortium \u003cem\u003eet al.\u003c/em\u003e, 2023) and the pathway annotations of Reactome (Milacic \u003cem\u003eet al.\u003c/em\u003e, 2024). For miRNA enrichment, we used miEAA (Aparicio-Puerta \u003cem\u003eet al.\u003c/em\u003e, 2023) and miRTargetLink 2.0 (Kern \u003cem\u003eet al.\u003c/em\u003e, 2021).\u003c/p\u003e\n\u003cp\u003eA significance threshold of adj.P \u0026lt; 0.05 and a log-fold change \u0026gt;=1 (upwards) or \u0026lt;=-1 (down) was applied for each of the three omics levels. \u003c/p\u003e\n\u003ch3\u003eTrajectory analysis \u003c/h3\u003e\n\u003cp\u003eTo derive feature trajectories over the granulopoietic developmental course, we combined weighted gene coexpression network analysis (WGCNA, (\u003cem\u003eWGCNA: an R package for weighted correlation network analysis\u003c/em\u003e, no date)) with time course analysis (timeOmics, (Bodein \u003cem\u003eet al.\u003c/em\u003e, 2022)). WGCNA enables the identification of features with correlated expression intensities in a series of groups. This approach has been successfully used to identify large clusters of features with similar expression trajectories in time series, but a post-hoc analysis of individual feature medians showed large variation of the individual patterns (data not shown). To reduce the variability of identified features, we employed timeOmics' multivariate analysis to delineate highly stringent expression clusters across time and omics layers. We found an advantageous combination of both approaches by first determining the number of correlated features showing similar trajectories using WGCNA and using the respective numbers of features per module to define the feature numbers of timeOmics components. Next, we conducted an enrichment analysis of the features of the derived components with \u003cem\u003eclusterProfiler\u003c/em\u003e and \u003cem\u003emiRTargetLink\u003c/em\u003e 2.0 (Kern \u003cem\u003eet al.\u003c/em\u003e, 2021). The resulting terms per omics layer were compared for overlaps between all 3 strata and plotted in combination to visualize the multi omics based term enrichment (Fig. 3 IV).\u003c/p\u003e\n\u003cp\u003eNext, we probed features of trajectories on transcriptome and proteome with granulopoiesis associated disease genes (severe congenital neutropenia (Donadieu and Bellanné-Chantelot, 2022) and acute myeloid leukemia (Makkar \u003cem\u003eet al.\u003c/em\u003e, 2023)) (Fig. 2 IVb).\u003c/p\u003e\n\u003ch3\u003eFunctional microRNA enrichment\u003c/h3\u003e\n\u003cp\u003eFor microRNA target enrichment, we employed the online based tool \u003cem\u003eEnrichMIR\u003c/em\u003e (Soutschek \u003cem\u003eet al.\u003c/em\u003e, 2022). It enables integration of a set of microRNA transcripts against a set of differentially expressed features to identify enrichment of targets and the motif-specific cumulative fold-change distribution of targets and non-targets. We interrogated the miRNA features of each detected trajectory (from Figure 2 III a-d) with the appropriate differentially expressed transcripts (see results text for details). First, we plotted all enriched miRNAs with their enrichment factors (x-axis), p-values (y-axis) and the fold change of the miRNA in the respective comparison group (Fig. 4 I-IV a). Next, we plotted the dose-response curves for features of the main enriched miRNA to visualize differences of expression intensity for features with (colored by binding motive) and without (black line) binding sites for the respective miRNA transcript (Fig. 4 I-IV b). Lastly, we applied GO enrichment to the miRNA targeted features to probe for their respective biological themes (Fig. 4 I-IV c) and plotted selected terms with their connected features (Fig. 4 I-IV d). \u003c/p\u003e\n\u003cp\u003eAs miRNAs generally act as mRNA transcript suppressors, we modeled the input data in a way that maintained this relationship. This required inversing the fold changes of differential expression in cases where the miRNAs decrease with maturation.\u003c/p\u003e\n\u003ch3\u003eTranscription factor enrichment \u003c/h3\u003e\n\u003cp\u003eTo derive transcription factor (TF) signatures, we used the R-package of \u003cem\u003eChEA3\u003c/em\u003e (Keenan \u003cem\u003eet al.\u003c/em\u003e, 2019) to investigate the differentially expressed features between differentiation steps. The results were first plotted as the top 10% of TF signatures per step (Fig. 5 Ia), while we secondly conducted a comparative analysis between all signatures to derive common and stage specific TF signatures using an Upset analysis (Fig. 5 Ib).\u003c/p\u003e\n\u003cp\u003eWe employed the integrated top-rank results of ChIP-X Enrichment Analysis 3 (\u003cem\u003eChEA3\u003c/em\u003e), a recently updated tool that enriches user-defined gene sets with libraries based on RNAseq, ChIPseq and published data (Keenan \u003cem\u003eet al.\u003c/em\u003e, 2019). \u003c/p\u003e\n\u003ch3\u003eOmics level shifts\u003c/h3\u003e\n\u003cp\u003eBased on our previous results of the unique proteome rank distribution in neutrophil granulocytes (Grabowski \u003cem\u003eet al.\u003c/em\u003e, 2019), we plotted the intensity based cumulative protein abundance (y-axis) per rank abundance (x-axis) of the different maturation stages (Figure 5 III a,c) on proteome and transcriptome level. Of note, we found that neutrophil granulocytes actively synthesize histones during maturation, and concluded that previously derived molecular numbers by the histone ruler (Wiśniewski \u003cem\u003eet al.\u003c/em\u003e, 2014) are not applicable, as this method assumes a stable ratio between DNA and histones. Next, we summarized the subcellular locations of all expressed features. To this end, we combined the location data of the human protein atlas (Thul \u003cem\u003eet al.\u003c/em\u003e, 2017), UniProt (UniProt Consortium, 2023), the GO consortium (Gene Ontology Consortium \u003cem\u003eet al.\u003c/em\u003e, 2023) and the neutrophil granule and membrane specific publication of Rorvig at al (Rørvig \u003cem\u003eet al.\u003c/em\u003e, 2013). For a detailed description of the process, please refer to the extended methods and the annotation preparation script. We summarized the total intensities per omics layer and subcellular location for each stage and used an alluvium plot for the location specific intensity as a percentage of the stages total (Brunson, 2020) (Figure 5 IV a,b).\u003c/p\u003e\n\u003ch3\u003eIntegration with single cell transcriptome data\u003c/h3\u003e\n\u003cp\u003eTo showcase the utility of our neutrophil atlas to enhance our understanding of scRNA-seq data, we applied Φ-Space (Mao et al., 2024) to map the neutrophil subset of the peripheral blood mononuclear cells (PMMC) scRNA-seq data from (Schulte-Schrepping \u003cem\u003eet al.\u003c/em\u003e, 2020)) to our bulk atlas. Φ-Space outputs scores measuring the similarity between each single neutrophil and each of our bulk cell stages. We first visualized the similarity scores as a heatmap. Next, we computed PC1 and PC2 of the similarity scores and visualized the PCs. Lastly, we visualized the similarity scores using uniform manifold approximation and projection (UMAP,(McInnes, Healy and Melville, 2018)). \u003c/p\u003e\n"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"Neutrophil granulocytes, Hematopoiesis, Multi Omics data integration, Transcription factors","lastPublishedDoi":"10.21203/rs.3.rs-6184761/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6184761/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eWe developed a comprehensive molecular atlas of human granulopoiesis, encompassing the transcriptome (mRNA), proteome, and non-coding RNA (miRNA, lncRNA) of seven defined stages of neutrophil maturation.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eTwo waves of mRNA transcription activity were observed. Whereas the first wave, extending from myeloblasts to metamyelocytes, showed correlated mRNA transcripts and protein translation, the second wave was characterized by uncoupling of transcription and translation. Integrated analysis revealed intricate dynamics of transcription and translation parallelization in granule subset proteins, ribosomes, and mitochondria, respectively.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eWe identified the hsa-miR-106a-5p family and hsa-miR-125a-5p as potential repressors of myeloid differentiation in myeloblasts, and hsa-miR-193a-5p as a repressor of the SWI/SNF complex. LncRNA transcripts displayed remarkable stage specificity. Whereas NEAT1 and MAILR were expressed only from the band stage onwards, ITGB2-AS1 was expressed exclusively in metamyelocytes. Novel lncRNAs with high expression were found in the terminal S-stage and PMN (ENSG00000288819, ENSG00000203644). A comprehensive analysis of transcription factor signatures, including activity on target genes, revealed previously unrecognized factors as well as a yet underappreciated transcriptional role of lactotransferrin (LTF) in early maturation stages. Characterization of single-cell-clusters via transcriptional profiles of histologically defined maturation stages enabled stage-specific cell labeling.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eTo ensure the availability of our data and the reproducibility of our analysis, our complete datasets, analysis codes and online tool for exploration are available on www.granulopoiesis.com.\u003c/p\u003e","manuscriptTitle":"A molecular atlas of human granulopoiesis","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-04-03 08:18:27","doi":"10.21203/rs.3.rs-6184761/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"nature-communications","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"NCOMMS","sideBox":"Learn more about [Nature Communications](http://www.nature.com/ncomms/)","snPcode":"","submissionUrl":"https://mts-ncomms.nature.com/","title":"Nature Communications","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Nature Communications","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"018fd020-228b-4531-a37b-b8d0976d4c77","owner":[],"postedDate":"April 3rd, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[{"id":45936772,"name":"Biological sciences/Immunology/Innate immune cells/Granulocytes/Neutrophils"},{"id":45936773,"name":"Biological sciences/Immunology/Haematopoiesis/Myelopoiesis"}],"tags":[],"updatedAt":"2025-04-03T08:18:27+00:00","versionOfRecord":[],"versionCreatedAt":"2025-04-03 08:18:27","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6184761","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6184761","identity":"rs-6184761","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}