{'fraction_non_alphanumeric': 0.03212098487689039, 'fraction_numerical': 0.016572928383952006, 'mean_word_length': 4.90930576070901, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transcriptomic and epigenetic mechanisms underlying myeloid diversity in the lung

The lung is inhabited by resident alveolar and interstitial macrophages as well as monocytic cells that survey lung tissues. Each cell type plays distinct functional roles under homeostatic and inflammatory conditions, but mechanisms establishing their molecular identities and functional potential remain poorly understood. Here, systematic evaluation of transcriptomes and open chromatin of alveolar macrophages (AMs), interstitial macrophages (IMs) and lung monocytes from two mouse strains enabled inference of common and cell-specific transcriptional regulators. We provide evidence that these factors drive selection of regulatory landscapes that specify distinct phenotypes of AMs and IMs and entrain qualitatively different responses to Toll-like receptor 4 (TLR4) signaling in vivo. These studies reveal a striking divergence in a fundamental innate immune response pathway in AMs and establish a framework for further understanding macrophage diversity in the lung.


Introduction
Innate immune cells, including macrophages, monocytes and dendritic cells, collectively called mononuclear phagocytes (MPs), provide the first line of cellular defense against pathogens. In addition, macrophages play essential roles in diverse developmental and homeostatic processes 1 . Within the lung, MPs are situated at and between the epithelial and endothelial surfaces, and their bulk gene expression profiles change during development 2 . Alveolar macrophages (AMs) are derived from embryonic erythromyeloid progenitor cells 3,4 and from fetal liver monocytes 5 . AMs play key roles in surfactant turnover and in maintaining an anti-inflammatory state that benefits efficient gas exchange 6,7,8 . However, the lung is populated by several subsets of MPs in addition to AMs 9,10 . Interstitial macrophages (IMs) proposed to be derived from yolk sac macrophages and later to be replaced by circulating monocytes 11 are situated in the lung interstitium around airways, blood vessels and nerves 12,13 . In addition, the extensive capillary bed of the lung harbors a marginated pool of monocytes that is retained even after extensive perfusion of the lung 14 . When the lung is injured or infected, additional monocytes are recruited and can amplify the inflammatory response 15,16,17 .
The tissue environment plays a crucial role in determining the expression and function of transcription factors (TFs) that specify tissue-resident macrophage phenotypes 18 . Following entry into the developing fetal lung, embryonic macrophages destined to become AMs upregulate Pparg, Klf4, Atf5 and Cebpb 19 . Consistent with the upregulation of these TFs at the mRNA level, their corresponding DNA recognition motifs are highly enriched in AM enhancers 20 and loss of Pparg expression results in defective AM differentiation 21,22 .
Conversely, rapid changes in enhancer activity states and gene expression are observed following the transfer of resident macrophages to a tissue culture environment 23,24 . These findings imply that tissue-derived signals are constantly sensed and integrated at the level of enhancers, which function to maintain tissue-specific programs of macrophage gene expression.
Cell-specific enhancers also function as important sites of action of broadly expressed signal-dependent TFs, exemplified by NF-κB and nuclear hormone receptors 25 . In principle, the distinct enhancer landscapes of each cell type enable cell-specific transcriptional responses to a common signal, such as lipopolysaccharide (LPS). These context-dependent responses to common signals are likely to be essential to the specific functions of distinct macrophage and monocyte subsets. Here, we use genetic and genomic approaches to investigate mechanisms that regulate gene expression in distinct tissue-resident macrophages in the lung and demonstrate divergent functional responses of these cell subsets to the Gramnegative bacterial endotoxin LPS. of the lung vasculature to remove circulating blood cells, we isolated resident lung myeloid cells (CD45 + cells depleted of T and B lymphocytes) by fluorescence-activated cell sorting (FACS) from 4 mice and generated single-cell transcriptional profiles for 1,048 cells that passed rigorous quality control. Dimensionality reduction analysis (t-SNE) separated cells into 7 clusters based on their gene expression pattern ( Fig.1a and Extended Data Fig.1a). These clusters included AMs with high Ltc4, Plet1, Mt1 expression and IMs with high Il1b and S100a expression. In addition, two populations were related to the major subsets of circulating blood monocytes 26,27 ; inflammatory monocytes (iMos) with high Ly6c2 expression, and patrolling monocytes (pMos) with Nr4a1 expression.
To investigate developmental relationships between lung monocytes and macrophages we performed RNA velocity analysis, which distinguishes between spliced and un-spliced mRNA to predict the 'future' state of individual cells 28 . The directionality vectors associated with each cell suggest that a substantial proportion of iMos are predicted to become more similar to pMos with time (Fig. 1a, Extended Data Fig. 1b). In addition, a smaller subcluster of iMos is predicted to differentiate into IMs.
Based on the population structure defined by scRNA-seq, we isolated four subsets (AM, IM, iMo, pMo) using FACS (gating strategy shown in Extended Data Fig. 1c) 29 . The sorted cells were processed for bulk RNA sequencing (RNA-seq) and for chromatin accessibility by Assay for Transposase-Accessible Chromatin (ATAC-seq) 30 . Contamination by other cell types was negligible for all MP subsets and Spearman correlation values were high for RNA-seq replicates in any given MP subset (Extended Data Fig. 1d, e). To investigate the impact of genetic variation on lung MP gene signatures during homeostasis we performed these experiments using two strains of mice C57BL/6J (B6) and DBA/2J (DBA) . These strains of mice exhibit 4.4 million single nucleotide polymorphisms (SNPs) and 540,000 insertions and deletions (InDels) 31 . Phenotypically, they differ in their susceptibility to a spectrum of inflammatory lung diseases 32,33 .
Gene expression profiles of lung MP subsets were markedly different at baseline (Fig. 1b). Clustering of RNA-seq data segregated samples by cell type, with strain as a secondary determinant. Comparing lung MP subsets, AMs were most different from all other cell types in both strains. For example, in B6 mice we found 1,769 differentially expressed genes between AMs and IMs (FC > 2, FDR < 0.05) (Extended Data Fig. 1f). Diversity between transcriptional programs for the remaining three cell subsets was less pronounced but still substantial with 1,076 genes differentially expressed between IMs and iMos and 584 genes between iMos and pMos (FC > 2, FDR < 0.05) in B6 mice (Extended Data Fig. 1f). Notably, while lung iMos were very similar to circulating iMos, pMos retained in the lung following extensive perfusion exhibited 266 genes that were differentially expressed (FC > 2, FDR < 0.05) compared to pMos in circulating blood (Extended Data Fig. 1g). Ingenuity pathway analysis (IPA) showed that genes involved in migration and adhesion were more highly expressed in pMos isolated from the lung (Extended Data Fig. 1h).
To infer potential biological meaning of the differentially expressed genes between lung MP subsets, we merged the differentially expressed genes in both C57 or DBA datasets. Using the obtained gene list, we performed IPA (Fig. 1d). AMs were enriched for genes associated with cell cycle progression, supporting previous studies indicating that these cells are selfrenewing. For example, Myc, a transcription factor known to play a central role in stimulating cell cycle progression, and Id1 known to play a role in cell growth were expressed highest in AMs (Fig. 1d). Under baseline conditions, IMs were the only lung MP subset enriched in genes positively associated with inflammatory response, exemplified by high expression of Il1b (Fig. 1d). In addition, IMs were enriched in genes related to cell movement (Cxcl13, Ccl24) and angiogenesis (Mmp9) (Fig. 1d). Apoptosis related genes, including Bbc3 and Eif3f, were most enriched in iMos and pMos. Collectively, these cellspecific transcriptomes define the distinct molecular phenotypes of the corresponding lung MPs.

Strain-specific transcriptomes of lung MPs
Analysis of gene expression between B6 and DBA mice showed hundreds of mRNA transcripts expressed with at least a two-fold difference between the strains in each lung MP subset (Fig. 2a). Of note, most of the inter-strain differences in gene expression were cellsubset-specific (Fig. 2b, c). For example, Spp1 , the gene encoding osteopontin, a molecule known to recruit cells to inflammatory sites and promote cell survival by regulating apoptosis 34 was several fold enriched in AMs from DBA mice (Fig. 2c). Increased expression of Spp1 in response to cigarette smoke is linked to emphysema 35 . Thus, higher levels in AMs may contribute to the greater sensitivity of DBA mice to cigarette-smoke induced emphysema than B6 mice 33 . Contrarily, Spint1 a gene encoding a serine protease inhibitor, was much higher expressed in AMs from B6 mice. An example of a gene selectively expressed in DBA IMs is Alox15, encoding for a lipoxygenase implicated in antiinflammation 36 . Strain-and cell type-specificity could also be observed for Cd4, which was only expressed in B6 IMs. Ligation of CD4 on monocytes/macrophages modulates gene and cytokine protein expression as well as T cell responses 37 . Stab2 a transmembrane receptor involved in bacterial scavenging, cell adhesion and angiogenesis 38 was only expressed in lung MPs from DBA mice. Conversely, Trim30d, was expressed by lung MPs isolated from B6 mice but absent in DBA mice. These results thereby identify genes that may play important roles in determining the phenotypic diversity of B6 and DBA mice.

Inference of transcriptional regulators from open chromatin landscapes
To investigate mechanisms responsible for the cell-subset-specific gene programs, we identified accessible regions of chromatin in lung MPs by performing ATAC-seq (Fig. 3a) 30 . We found the biological replicates to be highly correlated (Extended Data Fig. 2a). Essential myeloid genes, such as Spi1, encoding for the master TF PU.1, showed a highly overlapping open chromatin pattern in all lung MP subsets (Fig. 3a). In contrast, chromatin regions neighboring the AM-specific Fabp1 gene exhibited accessibility only in AMs. Open chromatin near Il1b, a pro-inflammatory cytokine, showed partial overlap among lung MPs, with a clear IM-specific peak, the cell subset showing the most inflammatory transcriptional pattern at baseline (Fig. 3a). In these examples, gene expression tracked closely with chromatin accessibility.
Genome-wide comparison of ATAC-seq peak tag counts revealed several thousand genomic sites with differentially opened chromatin between lung MP subsets. For example, pairwise comparison of AMs vs IMs revealed 16,402 peaks that displayed at least a 4-fold increase in ATAC-seq peak tag counts in IMs and 7,920 peaks in AMs (Fig. 3b). Comparison of the two lung monocyte subsets identified 3,087 peaks enriched in pMos and 2,788 peaks enriched in iMos (Extended Data Fig. 2b) implying that these cells indeed use distinct gene regulatory networks and a cell-subset-specific TF repertoire.
Clustering of lung MPs based on chromatin accessibility resulted in similar segregation as observed with clustering based on gene expression. AMs were most distant from the other 3 cell subsets and pMos and iMos grouped most closely together (Fig. 3c). To identify candidate TFs that could bind to distant open chromatin regions that would include enhancers, we performed de novo motif enrichment analysis of distal ATAC-seq peaks (> 3 kb from a transcription start site (TSS)). Using GC-matched genomic background in the HOMER 39 motif analysis we identified motifs for common myeloid cell lineage factors such as PU.1, C/EBP, AP-1, KLF and RUNX (Fig. 3d). This core set of motifs was highly enriched in all lung MP in both strains (Extended Data Fig. 2c,d), consistent with general roles of TFs recognizing these motifs in establishing macrophage identity.
To increase the power to detect cell-subset-specific TFs, we performed HOMER de novo motif analysis comparing cell-specific ATAC peaks to all ATAC peaks present in any of the four lung MPs (Fig. 3e). This method normalizes out motifs for TFs that are common to the four lung MP subsets and enriches for subset-specific motifs. Using this approach, we identified both known and previously unknown cell-subset-specific TFs. For example, PPAR is one of the top motifs in AMs and NUR77 one of the top motifs in pMos. These findings are in line with previously published data showing that PPARγ is essential for AM development and that NUR77 is a master regulator of patrolling monocytes 21,40 . Notably, the RELA motif, recognized by NF-κB, showed highest enrichment in IM enhancers, consistent with the enhanced inflammatory signature of these cells. IRF and SMAD motifs were selectively enriched in iMos, which are recognized by TFs with established roles in regulating pro-and anti-inflammatory programs of gene expression, respectively.
Other enriched motifs in unique lung MP subsets are associated with TFs for which functions in MP subsets have not been described. For example, AMs were enriched in the binding motif for MITF. IM enhancer regions were selectively enriched in HINFP, EGR and ZEB motifs. The MAFF motif was selectively enriched in pMos and is recognized by MAF TFs that form heterodimers with many other TFs, notably the NFE2 family 41 .
Since TFs of the same family can bind to the same or similar motifs, motif enrichment analysis cannot differentiate between TF family members. To infer likely candidates, we associated enriched motifs with mRNA expression of highly expressed TF family members that could bind these motifs. This approach revealed further similarities and differences among lung MPs. For example, PU.1, encoded by Spi1, was the ETS-domain transcription factor family member with the highest expression in all lung MPs (Fig. 3d). Compared to Spi1, other ETS members showed a more moderate expression level and higher degree of variation between cell subsets. Elf4, a TF previously shown to be important for interferon production and innate immune response against viruses 42 , showed cell-subset-specificity with highest expression in monocytes (Fig. 3d). Members of the C/EBP, AP-1 and KLF families were highly expressed in general, with some family members showing cell-subsetspecificity. Of the five NF-κB proteins that could bind the RELA motif, Rela (also known as p65) and Nfkb2 (encoding for the p52 precursor protein also known as p100) were the most prominent ( Fig. 3e). Overall, these findings suggest the identities of TF combinations required for establishing the regulatory landscapes of each lung MP subset.

Transcriptional regulators inferred from effects of genetic variation
The functional relevance of TF motifs can be assessed by analyzing the impact of motif mutations 43,44 . To investigate the functional importance of specific motifs in establishing open chromatin, we analyzed the effects of the natural genetic variation distinguishing B6 and DBA mice 31 . Comparing open chromatin regions between B6 and DBA mice (Fig. 4a), we found 2,000-8,000 genomic regions that showed at least a 2-fold difference in ATAC-seq signal in lung MP subsets.
Effects of SNPs and InDels on open chromatin were assessed using MMARGE (Motif Mutation Analysis for Regulatory Genomic Elements) 45 . This method evaluates the significance of the relationship of mutations within TF motifs with strain-specific gain or loss of specific chromatin features, such as ATAC-seq signal. For example, if SNPs or InDels within PU.1 recognition motifs in B6 mice are statistically associated with a reduction in local ATAC-seq peaks compared to peaks in DBA mice, the PU.1 motif is scored as contributing to the formation of open chromatin. MMARGE analysis is thus qualitatively different from motif enrichment analysis in that it examines the functional consequence of a mutation within a binding motif on chromatin accessibility.
Application of this method to an annotated library of 287 distinct TF motifs yielded 51 significant motifs across the four lung MP subsets (Fig. 4b). As expected, motifs recognized by key regulators of the macrophage lineage such as PU.1, CEBP, AP-1, and KLF were highly significant in most or all lung MPs. MMARGE analysis also identified motifs in which mutations exerted subset-specific effects. For example, although a generic IRF motif was significant across all subsets, mutations in the motif for T1ISRE selectively affected ATAC peaks specific for AMs, and motifs for IRF-7,8,9 specifically affected pMos. ALX, MEF2C and ZNF423 were motifs specific for IMs. NFY and NFYB mutations were significantly associated with reduced open chromatin in iMos, but not other lung MP subsets (Fig. 4b,c). We also confirmed that mutations in the motif recognized by ZEB led to higher chromatin accessibility, consistent with its reported function as a transcriptional repressor 43 . When applying MMARGE specifically to open promoter regions in the two mouse strains, a much lower number of motifs were found to be significantly associated with chromatin accessibility. EGR was specific for AMs and MSC for IMs. Mutations in RBP-J and STAT motifs in promoter regions led to significant changes in chromatin accessibility in both monocyte subsets (Fig. 4b).
To integrate the findings of HOMER motif analyses using genomic background and lung MP common background with the functional results of MMARGE, we constructed networks of the most relevant TFs for each cell subset (Fig. 4c). In these networks the node size is proportional to the fraction of open chromatin regions containing the motif and the edge thickness represents the significance (p-value) of this motif within each analysis. Edges are colored according to the motif analysis method; HOMER using random GC background is brown, HOMER using lung MP common background is turquoise and MMARGE is purple.
These networks showed a high overlap of HOMER and MMARGE results for the critical macrophage regulators PU.1, CEBP, AP-1 and KLF. RUNX was the only macrophage lineage-determining TF (LDTF) (gray nodes) that did not show concomitant significance with motif enrichment and functional analysis in all cell subsets. In addition to the five LDTFs, this integrated approach identified additional TFs (yellow nodes) for each cell subset that reached significance in both analyses. For AMs, CTCF reached significance in HOMER and MMARGE analysis, while for IMs ZEB, FOXD, and HINFP were significant. Regarding the two monocyte subsets, IRF was significant for iMos and NUR77, HOXA, YY1 and CTCF for pMos. Taken together, these data provide functional evidence for putative TFs that shape the chromatin landscape in lung MPs and provide evidence for additional TFs that have selective roles in specific subsets of lung MPs that were not identified by conventional motif enrichment analysis.

Divergent responses of lung MPs to acute inflammation
The divergent baseline patterns of gene expression and distinct open chromatin landscapes associated with each lung MP subset suggests that they are programmed to respond to environmental perturbations in a cell-specific manner. To directly test this prediction, we performed RNA-seq analysis of AMs, IMs and iMos following exposure to LPS administered either intraperitoneally (i.p.) or intranasally (i.n.). Inflammation induced by i.p. LPS is a widely accepted model for endotoxin-mediated septic shock/acute respiratory distress syndrome and i.n. LPS for pneumonia 46 .
I.p. LPS administration resulted in a prominent inflammatory response in the lung with a significant shift in the composition of myeloid cells 2 h, 6 h and 22 h after LPS administration (Extended Data Fig. 3a). As expected, neutrophils and iMo were rapidly recruited to the lungs. The proportion of IMs gradually increased during the course of inflammation, while AMs showed a temporary relative decrease at 2 h and 6 h after i.p. LPS administration (Extended Data Fig. 3a).
IMs and iMos exhibited robust and largely overlapping changes in gene expression after i.p. LPS administration (in comparison to baseline expression), whereas AMs showed very modest changes (Fig. 5a,b and Extended Data Fig. 3b). In addition to proinflammatory cytokines (Il1b, Ccl5, Cxcl10), genes for anti-inflammatory cytokines such Il10 and reactive oxygen species scavengers such as Sod2, as well as matrix metallopeptidases (Mmp14, Mmp8), were rapidly induced in IMs and iMos after LPS administration (Fig. 5b).
IPA of differentially regulated genes in IMs and iMos compared to their own baseline indicated a pattern that was consistent with rapid activation of TLR4 and downstream signaling pathways (Fig. 5c). Induced pathways included TLR signaling, NF-κB signaling, IL-6 and interferon signaling. Enriched pathways for iMos and IMs were qualitatively similar, but some substantial quantitative differences were noted, including oxidative phosphorylation, eIF2 and ErbB2 signaling. ERK5 signaling and sphingosine-1-phosphate signaling were predicted to be activated in iMos but inhibited in IMs. I.n. delivery of LPS resulted in a robust recruitment of neutrophils to the lung, similar to the i.p. LPS model (Extended Data Fig. 3c). In contrast to i.p. delivery, AMs were the most responsive lung macrophages, with very few acute changes in gene expression observed in IMs and iMos ( Fig. 5d and Extended Data Fig. 3d). In general, the inflammatory response was more dampened in the i.n. LPS response compared to i.p. LPS with fewer number of differentially expressed genes and less pronounced changes compared to baseline values (Fig. 5e). A small set of genes was rapidly induced in AMs that overlapped with rapidly induced genes in IMs including Il1b and Tnf (Fig. 5e). However, the majority of genes responding to LPS in AMs exhibited a maximum response at 6 h, representing a delay in comparison to the maximal response of IMs and iMos, which occurred at 2 h following i.p. administration of LPS (Fig. 5d).
Comparison of genes induced at the time of maximal response in IMs (2 h, i.p. LPS) to genes induced at the time of maximal response in AMs (6 h, i.n. LPS) indicated that of the 253 up-regulated genes in AMs and 1,155 up-regulated genes in IMs, only 43 were in common (Fig. 6a). A similar pattern was observed for down-regulated genes. Remarkably, IPA analysis indicated that the predicted functional consequences of the AM and IM changes in gene expression were in opposite directions for several of the most enriched pathways (Fig. 6b). One of the most striking differences is provided by differential regulation of genes involved in oxidative phosphorylation, which were down-regulated in IMs but up-regulated in AMs (Fig. 6c).
To explore possible mechanisms, we analyzed the basal expression of genes belonging to the TLR signaling pathway (Fig. 6d). Tlr4, Ly96 (encoding MD2), Cd14, Lbp (both encoding co-receptors for LPS), and the signaling adapter protein encoded by Myd88, are expressed at similar levels in AMs and IMs. However, the majority of genes encoding additional downstream adapters, kinases and transcription factors mediating responses to TLR4 ligation exhibit lower expression in AMs. These genes include the key signaling adapters encoded by Tram1, Ticam1, Tirap and Traf6 and downstream kinases encoded by Mapk8, Ikbkb, Tab1 and Tab2. Genes encoding several major transcriptional targets of these kinases also exhibit lower levels of expression in AMs, including all five family members of the NF-κB family cJun and Irf5 and Irf7. These differences would be predicted to reduce the strength of the response of AMs to LPS, but not result in the observed divergent patterns of this response.
We next analyzed open chromatin in the vicinity (+/− 3 kb) of the transcriptional start sites (TSSs) of genes that exhibited differential responses to i.p. LPS in IMs at 2 h and responses to i.n. LPS in AMs at 6 h (Fig. 6e). In both cases, open chromatin was enriched for motifs recognized by PU.1 and SP2/KLF factors. However, the remaining significant motifs were divergent. Notably, ATAC-seq peaks associated with genes preferentially induced in IMs were enriched for NFY, NFIL3/CEBP, IRF, AP-1 and REL motifs. In contrast, ATAC-seq peaks associated with genes preferentially induced in AMs were enriched for ERG, ELF5 and HIC1 motifs.
Direct analysis of open chromatin patterns of genes exhibiting differential responses in AMs exposed to i.n. LPS and IMs exposed to i.p. LPS revealed a continuum ranging from nearly identical patterns to strongly cell-type-specific patterns (Fig. 6f). For example, while Batf and Cyb5a are induced in AMs but repressed in IMs, the open chromatin environment is nearly indistinguishable in the two cell subsets. A similar pattern of open chromatin is also observed for Ccl3, which is conversely induced in IMs compared to AMs. In contrast, open chromatin environments in the vicinity of Car4 and Ndufs8, which are selectively induced in AMs, and Ccl5, which is induced in IMs, exhibit substantial differences in their local ATACseq profiles, consistent with the results of motif enrichment analysis. Therefore, AMs and IMs exhibit qualitatively different responses to LPS in vivo associated with differences in the TLR4 signaling pathway at baseline and the open chromatin landscapes of TLR4 target genes.

Discussion
Although lung MPs reside in the same organ, they exist in distinct microenvironments and perform different functions. Consistent with this, transcriptomic analysis of major subsets of lung MPs identified by scRNA-seq indicates substantial differences in gene expression under homeostatic conditions. For each subset, natural genetic variation distinguishing B6 and DBA mice resulted in hundreds of significant changes in gene expression, the majority of which were cell-subset-specific. These differences are most likely explained primarily by effects of genetic variation on cell-specific regulatory elements 43 and contribute to the phenotypic diversity of B6 and DBA mice. The delineation of the specific patterns of open chromatin in each cell-subset and the influence of genetic variation enabled inference of major classes of TFs underlying their development and function, some of which have not previously been suggested to play roles in lung macrophages. It will therefore be of interest to perform loss of function experiments and identify signaling pathways that are necessary for the expression and activities of these factors within each macrophage subset.
The most striking findings to emerge from the present studies were the qualitatively different responses of AMs and IMs to i.p. and i.n. LPS. At one level, these divergent responses can be easily understood based on anatomic compartmentalization of the signal. Remarkably, the direct responses of IMs to i.p. LPS and AMs to i.n. LPS are mostly divergent, in many cases leading to predictions of opposite biological outcomes. Analysis of the open chromatin landscapes of genes that exhibit divergent responses to LPS in AMs and IMs suggest at least two underlying mechanisms. A subset of genes exhibit distinct patterns of open chromatin, consistent with differences in enriched motifs for differentially regulated genes as a whole. In these cases, the accessible landscapes for signal-dependent TFs, such as NF-κB, are different. Cell-specific responses to TLR4 ligation of genes with distinct patterns of open chromatin are likely to result at least in part from the corresponding genomic regions functioning as cell-specific regulatory elements.
A second subset of genes exhibiting AM-or IM-specific responses to LPS were observed to have indistinguishable patterns of local open chromatin. While it is possible that cell-specific regulation is determined by more distal regulatory elements for some of these genes, our findings are more consistent with the possibility that the mechanisms by which TLR4 ligation regulates gene expression at the signaling level are substantially different in AMs vs IMs. The reduced expression of genes encoding components of the TLR4 signaling pathway in AMs as compared to IMs is consistent with an attenuated response to LPS in comparison to IMs, but does not explain genes that are activated in an AM-specific manner. Instead, our findings suggest that TLR4 couples to alternative, as yet undefined, signaling molecules in AMs to activate a divergent set of transcriptional activators and/or repressors. This mechanism would provide an alternative explanation for cell-specific differences in responses of genes exhibiting similar patterns of open chromatin, through the recruitment of different classes of signal-dependent factors that are not involved in the selection of open chromatin itself. Moving forward, it will be of substantial interest to better understand mechanisms underlying distinct responses of AMs and IMs to TLR4 ligation and more generally to decipher their cell-type specific responses in disease models .

Animal model
Animals.-We obtained C57BL/6J and DBA/2J mice from the Jackson Laboratory. We established breeding colonies for both strains in our vivarium and housed them in adjacent cages. Male mice, 8-10 weeks old were used for experiments. All protocols were approved by the University of California, San Diego Animal Care and Use Committee. Mice were housed in cages, maximum 4 mice/cage, in a 12 h/12 h light/dark cycle with free access to food and water.
For baseline experiments we used the lungs of two mice/per sample, we performed >5 independent experiments and representative data of two biological replicates are used for anlaysis. The different lung MP subsets used for final analysis derive from multiple mice for both strains.

Lung MP isolation
Mice were euthanized with carbon dioxide. Lungs were perfused with PBS trough the right ventricle until no blood was visible in the tissue. Lungs were removed, minced and placed in a digestion buffer containing collagenase type 4 (1.6 mg/ml; Roche), DNase type 1 (50 units/ml; Roche) and Flavoperidol (1 μM, Sigma-Aldrich) in RPMI 1640 medium. Lungs were incubated at 37 °C for 15 min while on a rotator. Red blood cells were lysed using RBC lysis buffer (eBioscience) at 4 °C for a 5 min incubation. The obtained cell suspension was filtered using a 70-micron filter, washed and cells were resuspended in PBS. Cells were then incubated with anti-Fc blocking mAb (anti CD16/32, clone 93, BioLegend) and viability dye Zombie Aqua (BioLegend) for 10 min. After this, cells where stained with an antibody cocktail containing APC Cy7-conjugated CD45, PE-conjugated F4/80, BV605conjugated CD11c, v450-conjugated CD11b, PE Cy7-conjugated CD64, AF700-conjufgated Ly6c, FITC-conjugated Ly6g, AF647-conjugated SiglecF mAb for 30 min. All labeled cells were passed through a FACS Aria II flow cytometer (BD Biosciences) and after proper gating 4 subsets of lung MPs were sorted 29 . Data were analyzed with FlowJo 8.3.3 software (TreeStar Inc.). The antibodies used for FACS are listed in Supplementary Table 1.

Single-cell RNA-sequencing
Lungs of 2-week old DBA/2J mice (3 males and 1 female) were harvested and single cell suspension of lung cells were prepared as described above. Using FACS we isolated lung myeloid cells by gating on CD45 + cells and excluding T cells, B cells and γ/δT cells using antibodies against TCRβ chain, CD19 and TCR γ/δ chain. The antibodies used for FACS are listed in Supplementary Table 2.
Sorted lung myeloid cells were suspended in PBS containing 0.04% molecular biology grade non-acetylated BSA at a concentration of 1,000 cells/μL. The purity of the sorted cells as well as the sorted cell number were verified. Using the 10X GemCode technology single cells were captured with barcoded beads in a droplet emulsion for the RT reaction. After obtaining 10X barcoded cDNA, libraries were prepared following the standard protocol of 10X Single Cell 3′ v1. The libraries were sequenced on one lane of a Hi-Seq 4000 (Illumina).
Libraries were purified using 1.5 volumes of SpeedBeads in 2.5 M NaCl, 20% PEG8000, further size selected using 10% PAGE/TBE gels, eluted and quantified with Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), pooled and sequenced on a NextSeq 500 (Illumina) for 76 cycles or a Hi-Seq 4000 (Illumina) for 51 cycles.

Data analysis
Single-cell RNA-sequencing.-We used the R package Seurat 49 (version 2.3.4) to perform quality control and clustering analysis. Unique molecular identifiers (UMI) that barcode each individual mRNA molecule within a cell during reverse transcription were used to remove PCR duplicates. Cells expressing <300 or >2000 genes were filtered out to exclude non-cells or cell aggregates. We also filtered out cells that contained a percentage of mitochondrial UMI reads 5%. After quality control, 1,048 cells were included in subsequent analysis. Seurat standard pipeline was used for normalization and scaling, and first 11 PCs from PCA were used for t-Distributed Stochastic Neighbor Embedding (t-SNE) dimensionality reduction. The shared nearest neighbor (SNN) modularity optimization based clustering algorithm was used to group cells into 8 clusters with resolution=0.8. Velocyto 28 a python package was used to predict the migration of cells.
RNA-sequencing analysis.-FASTQ files were processed to assess quality by determining general sequencing bias, clonality and adapter sequence contamination. RNA sequencing reads were aligned to the individual genomes for the mouse strain (based on mm10 reference genome) of origin using STAR 50 . Gene expression levels were calculated using HOMER 39 by counting all strand specific reads within exons. Only the most abundant transcripts, including multiple alternative variants, were selected for each gene, and the genes whith a length smaller than 250bp were removed. Transcripts per million (TPM) were used to evaluate the correlation among replicates. Differential gene expression was calculated using DESeq2 51 to assess both biological and technical variability between experiments. Unsupervised hierarchical clustering was used to cluster the gene expression in the heatmaps. To identify functional enrichment of genes we used Ingenuity Pathway Analysis (IPA) software (Qiagen).
ATAC-sequencing analysis.-FASTQ files were mapped to individual genomes for the mouse strain (based on mm10 reference genome) with Bowtie2 52 and shifted back to reference coordinates for DBA/2J strain. Peaks were called with HOMER (findPeaks) using parameters "-style factor -minDist 200 -size 200". After merging these peaks, correlations among replicates from the same cell subset (and the same LPS treatment time point, if applicable) were evaluated by correlation using tag counts. The two most highly correlated samples were used for identifying the most robust peaks using the irreproducible discovery rate (IDR) method. For this step, peaks were called with HOMER's findPeaks, using parameters "-L 0 -C 0 -fdr 0.9 -minDist 200 -size 200". For baseline analysis, peaks were merged with HOMER's mergePeaks and annotated with HOMER's annotatePeaks.pl using the merged IDR tag directories. For the LPS experiments, peaks were merged with HOMER's mergePeaks using parameters "-d 100". The raw tags of all samples which had reasonable correlation were quantified with HOMER (annotatePeaks.pl) using parameter "noadj". For each cell subset, these raw tags which contained at least 16 tags in at least 2 samples were used to identify differentially bounded peaks by DESeq2. De novo motif analysis was performed on peaks 3kb away from closest TSS using a modified version of HOMER's findMotifsGenome.pl from MMARGE. For LPS analysis, peaks 3kb away from the closest TSS were used and motif analysis was performed with HOMER's findMotifsGenome.pl with default parameters. 45 .-To model the impact of a motif on chromatin accessibility, MMARGE's pairwise comparison was used. In brief, peaks from both strains per cell type were merged with HOMER's mergedPeaks and annotated with normalized tag counts. These files were scanned for motifs and the distribution of chromatin accessibility is calculated for all peaks missing a motif in either C57BL/6J or DBA/2J mice. To calculate significance, a Student's t-test is performed between these two distributions.

Statistical analyses
Differences in gene expression between different cell subsets and between the two mouse strains were calculated by DESeq2 and multiple-testing was corrected with the Benjamini-Hochberg procvedure. Data are presented as mean ± SD. Differences in percentage of different subsets of myeloid cells in the lung after i.p. and i.n. LPS administration were tested using the non-parametric Wilcoxon signed-rank test with Bonferroni correction. Data are presented as mean ± SD.

Data availability statement
The sequencing data reported in this article are deposited in GEO with the following accession number GSE136916. Additional data supporting the presented findings are available in the manuscript and upon request from the corresponding author.

Extended Data
Extended Data Fig. 1. Single-cell RNA-seq clusters, FACS sorting strategy and quality control of FACS and RNA-seq results a. Single-cell RNA-seq gene expression in clusters used to determine cluster identities. Heat map is representing expression values for the most significant genes in each cluster. Cells from the lungs of 3 male and 1 female DBA/2J mice were pooled. b. Zoomed-in view of velocity analysis for single-cell RNA-seq from Fig. 1a. c. Flow cytometry analysis and sorting strategy to obtain subsets of lung mononuclear phagocytes (MPs).  Bars represent the average percentage of given cell subset out of CD45+ leukocytes ± SD. Three independent experiments were pooled, 0h N=11, for all other time points N=8. Nonparametric Wilcoxon signed-rank test with Bonferroni correction was used, * p < 0.05, ** p < 0.01, *** p < 0.001.
b. Venn diagram of 4,499 differentially expressed genes in AM, IM and iMo after i.p. LPS administration. c. Flow cytometry analysis of lung MP subsets and neutrophils after i.n. LPS administration. Bars represent the average percentage of given cell subset out of CD45+ leukocytes ± SD. Two independent experiments were pooled, 0h N= 11, for all other time points N=4. Nonparametric Wilcoxon signed-rank test with Bonferroni correction was used, * p < 0.05, ** p < 0.01, *** p < 0.001. d. Venn diagram of 806 differentially expressed genes in AM, IM and iMo after i.n. LPS administration.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. c. Examples of strain specific genes included in the heat map in Fig. 2b. Bar graphs represent mean log2 transcripts per million (TPM+1) ± SD, N=2. Asterisks indicate statistical significance, * p-adj < 0.05, ** p-adj < 0.01, *** p-adj < 0.001 calculated by DESeq2 and multiple-testing corrected with Benjamini-Hochberg.  a. Inter-strain comparison of ATAC-seq peak tag counts for the four lung MP subsets. Scatter plots are showing log2 tag counts of ATAC-seq peaks, colored dots show ATAC-seq peaks with fold change (FC) 2 or higher (blue for AM, orange for IM, bordeaux for iMo and purple for pMo). b. MMARGE analysis of lung MP. Distribution of chromatin accessibility was calculated for all peaks missing a motif in either B6 or DBA mice. Student's t-test was used to test for significant differences between these two distributions. Heatmap showing p-values for all