Diverse homeostatic and immunomodulatory roles of immune cells in the developing mouse lung revealed at single cell resolution

At birth, the lungs experience a sudden transition from a pathogen-free, hypoxic, fluid-filled environment to a pathogen-rich, rhythmically distended air-liquid interface. While many studies focus on adult tissue, the heterogeneity of immune cells in the perinatal lung remains unexplored. Here, we combine single cell transcriptomics with in situ hybridization to present an atlas of the murine lung immune compartment during a critical period of lung development. We show that the late embryonic lung is dominated by specialized proliferative macrophages with a surprising physical interaction with the developing vasculature. These macrophages disappear after birth and are replaced by a complex and dynamic mixture of macrophage subtypes, dendritic cells, granulocytes, and lymphocytes. Detailed characterization of macrophage diversity revealed a precise orchestration of five distinct subpopulations across postnatal development to fill context-specific functions in tissue remodeling, angiogenesis, and immunity. These data both broaden the putative roles for immune cells in the developing lung and provide a framework for understanding how external insults alter immune cell phenotype during a period of rapid lung growth and heightened vulnerability.


38
Prior to birth, the lung is maintained in a fluid-filled, immune-privileged, hypoxic 39 environment. Upon birth, the tissue quickly transitions to an air-filled, immune-challenged, oxygen 40 rich environment following the infant's first breath 1 . At this point, the distal, gas-exchanging 41 alveoli are suddenly subjected to the mechanical forces of spontaneous ventilation and exposed to 42 diverse pathogens present in the external environment 2 . Adaptation to this rapid environmental 43 shift is necessary for perinatal survival and is mediated by complex physiologic processes 44 including reduced pulmonary arterial pressure, an exponential increase in pulmonary blood flow, 45 establishment of air-liquid interface, and surfactant production 3 . The immune system is essential 46 for lung homeostasis, wound-healing and response to pathogens 4 . Although the development of 47 the murine immune system begins during early embryogenesis, little is known regarding how the 48 dynamic physiologic changes at birth alter the lung immune cell landscape, and whether specific 49 immune cell subpopulations influence lung growth and remodeling in addition to serving 50 established immunomodulatory functions. 51 Fifteen cell clusters were identified via Leiden community detection 21 and verified by t-113 distributed stochastic neighbor embedding (t-SNE) 22 ( Figure 1B)  2B and C). After birth, Mac I disappeared while Mac II abundance peaked to 35% of the total 134 before decreasing again and disappearing by P21. Mac III and Mac V abundance increased steadily 135 with Mac V being most abundant at all postnatal timepoints. The Mac IV population was relatively 136 stable over time at ~10% of the total. 137 The Mac I-V clusters broadly separated into two groups based upon expression of the 138 disabled 2 gene (Dab2), which regulates macrophage polarization 23 and the placenta-specific 8 139 gene (Plac8), which is related to bacterial immunity 24 . Mac I-IV cells expressed Dab2 but not 140 Plac8 while Mac V showed the opposite pattern ( Figure 2D). We confirmed by multiplexed FISH 141 that all Cd68+ cells in the lung expressed either Dab2 or Plac8 at both E18.5 and P7 ( Figure 2E). 142 Dab2 and Plac8 expression was not consistent with previously reported markers of macrophage 143 lineages derived from yolk sac and fetal liver 25, 26 (Supplemental Figure 2A). 144 145 Embryonic macrophages are proliferative and encircle developing vessels prior to birth. 146 Mac I cells are the predominant immune population at E18.5, hence we aimed to 147 understand their function and localization. Differentially expressed genes (DEGs) in Mac I 148 included the proliferation markers Mki67 and Mcm5 ( Figure 3A). Across macrophages and 149 monocytes, proliferation decreased from 60% of cells at E18.5 to only 10% by P21. Most 150 proliferating cells were Mac I prior to birth and distributed across Mac II-V after birth 151 (Supplemental Figure 2B). These data are consistent with prior reports of "bursts" of proliferation 152  after recruitment of macrophages into embryonic tissues, followed by low-level self-renewal by 153 adulthood 27 . 154 155 Localization of Mac I cells within the E18.5 lung revealed Cd68+ cells scattered 156 throughout the lung parenchyma but also, surprisingly, forming almost complete rings around 157 blood vessels of 20-30 µm in diameter found adjacent to large, conducting airways ( Figure 3B). 158 In contrast, Mac I cells were not found encircling small airways ( Figure 3C). Many of the vessel-159 surrounding macrophages expressed Mki67 ( Figure 3D) and Dab2 ( Figure 3E) but not Plac8 160 ( Figure 3F). Given that Dab2+ cells at E18.5 include Mac IV cells, we aimed to distinguish these 161 from the Mac I cells by detecting the expression of Gal and C1qa, which were determined to be 162 Tek Mmp9 specific markers by scRNA-Seq (see below). These studies demonstrated that the majority of 163 Cd68+ macrophages surrounding the vessels were Gal+ with an occasional C1qa+ cell, indicating 164 a predominance of Mac I and a small number of Mac IV cells comprise the perivascular 165 macrophage population ( Figure 3G). We then asked whether any of our Mac I-V clusters are 166 related to previously reported perivascular macrophages that promote vascular remodeling in the 167 developing hindbrain and retina 28 . Mac I cluster expressed Cxcr4 and Nrp1 but failed to express 168 many other genes characteristic of these previously reported macrophages, suggesting a distinct 169 phenotype ( Figure 3H). Moreover, no concentric perivascular macrophages were observed at P1, 170 consistent with a function specific to Mac I prior to birth ( Figure 3I). Taken together, these data 171 suggest that within the embryonic lung, Mac I macrophages are highly proliferative and localize 172 to small vessels, suggesting a potential role in pulmonary vascular growth or remodeling.  in addition to those found encircling small vessels ( Figure 4B). 194   Table  219 4), suggesting a role in the localization of CCR2 expressing monocytes 42 . The Mac IV cells also 220 expressed Cxcl16, an IFNg regulated chemokine that promotes T cell recruitment through 221 CXCR6 43 . Mac IV also highly expressed Mrc1 (CD206) (Supplemental Fig. 3A). However, within 222 the Mac IV cluster there were a group of cells with lower Mrc1 and high MHC II gene expression 223 (H2-Ab1, H2-Eb1, and Cd74). These transcriptional differences within Mac IV are similar to two 224 recently reported interstitial macrophages (IMs) in adult lung that can be differentiated by 225 expression of Mrc1 and MHC II genes 44 . However, the other genes reported to distinguish those 226 two clusters (e.g Cd68, Cx3cr1, Mertk, Cclr2) were diffusely expressed throughout the Mac IV 227 cluster (Supplemental Fig. 3B). In situ imaging to localize C1qa+Cd68+ cells in the postnatal lung 228 found Mac IV cells remained abutting small blood vessels as well as the abluminal side of large 229 airways ( Figure 4E). The characteristic location and the expression of numerous genes important 230 for leukocyte recruitment and pathogen defense suggest that Mac IV may serve as patrollers, 231 playing a key role in immune-surveillance, innate pathogen defense, and antigen presentation.  Table 5). This dual 254 immune signature was similar to that seen in Mac II and III, emphasizing the importance of a finely 255 tuned inflammatory response in the developing lung. 256 Within Mac V we found a gradient of cell states between two distinct phenotypes, with 258 some genes (e.g. Ccr1 and Ccr2) expressed at one end and other genes (e.g. Ace and Slc12a) 259   Table 10). We also identified a neutrophil cluster distinguished by expression of S100a8 and 325 S100a9, which are released during inflammation 75 , and Stfa1 and Stfa2, cysteine proteinase 326 inhibitors important for antigen presentation 76 ( Figure 6F and Supplemental Table 11). Our data 327 agree with prior work demonstrating that lung resident basophils express signaling molecules 328 important for interaction with neighboring cells 74 , and that mast cells and neutrophils are primed 329 for rapid innate immune responses upon pathogen challenge. 330 331 Naïve lymphocytes populate the lung at birth. 332

Lymphocytes including ILC2s (expressing Areg), NK cells (Gzma), B cells (Ms4a1) and 333
T cells (Cde3) were present in the lung at low frequencies (approximately 2% of immune cells) 334 prior to birth and increased in frequency after birth. By P21 lymphocytes comprised 60% of total 335   After birth, B cell abundance increased, but the proliferating fraction decreased ( Figure 7B). The 340 cumulative distribution of V and J loci germline similarity revealed no somatic hypermutation 341 ( Figure 7C) and most B cells expressed Ighm and Ighd indicative of an IgM isotype ( Figure 7D). To test the clonality of the B cell repertoire at birth, we assembled the heavy and light chain loci 345  (E-H) At P1 (n=7 mice) and P21 (n=8 mice), lungs were processed to a single-cell suspension, and flow cytometry was used to assess frequencies of (F) CD4+CD3+ and CD8+CD3+ T cells, (G) Gata3+, Tbet+, Rorγt+ and Foxp3+ CD4+ T cells, and (H) IL-4-producing Gata3+CD4+ T cells or IFNγ-producing Tbet+CD4+ T cells. Data shown as mean ± SD, **P < 0.01, ****P < 0.00001 by Student's t test.
and performed t-SNE on a feature-selected transcriptome limited to over-dispersed genes in B 346 cells. Proliferating cells clustered together but candidate clonal families did not, suggesting 347 primarily homeostatic B cell proliferation rather than clonal expansion 81 . 348 The majority of T cells expressed Trac, suggesting ab identity, with a few Trac-cells 349 expressed Tcrg-C4, suggesting gd T cell identity (Supplemental Figure 7C). T cell receptor 350 diversity showed no sign of clonal expansion (data not shown). Outside the thymus, ab Τ cells are 351 usually either CD4+ (with subgroups Tbet+ or T helper (Th) 1, Gata3+ or Th2, and FoxP3+ or 352 Th17) or CD8+, however we characterized T cell heterogeneity in neonatal lungs by flow 353 cytometry and found that 85-90% of CD3+ cells were CD4-CD8-at both P1 and P21 ( Figure 7E), 354 confirming an earlier report suggesting this is a neonatal-specific phenotype. mRNA expression 355 analysis qualitatively confirmed this finding (Supplemental Figure 7D). The frequency of total T 356 cells was similar at P1 and P21 (Supplemental Figure 7E), however several T cell subsets (CD8 + , 357 CD4+ Th1, and Treg cells) increased by P21 while other subsets remained constant ( Figure 7F and 358 G). In response to stimulation with phorbol myristate acetate (PMA) and ionomycin, a greater 359 number of CD4 + Th2 cells at P21 produced IL-4 as compared to P1 (28.3 ± 19.7 vs. 8.6 ± 5.7, 360 molecules and spanned a gradient between two extreme phenotypes, one that expressed high levels 377 of homeostatic genes during early postnatal development, and a second with immunomodulatory 378 function that resembles previously reported nonclassical monocytes. Lymphocytes increased in 379 abundance from almost zero before birth to more than half of lung immune cells by P21, but 380 maintained a naive phenotype skewed toward type II immunity and with predominantly Cd4-Cd8-381

T cells. 382
This comprehensive study has far-reaching implications for lung biology.  Mmp9). Although the perivascular macrophages we observed in 420 the embryonic lung expressed low levels of Nrp1, they appear distinct from the macrophages that 421 influence retinal and hindbrain angiogenesis, expressing a unique set of ECM remodeling and 422 angiogenic genes, including genes that may modulate vascular tone and permeability. The 423 distinctive location of these macrophages and their gene signature imply a role in vascular 424 development. Furthermore, these encircling macrophages disappeared after birth, suggesting a 425 function temporally restricted to prenatal development. Future studies to selectively target this 426 subpopulation will be required to further establish their function and to delineate the signals 427 responsible for the cessation of the macrophage-vascular interaction after birth. 428 During embryonic development, the lung is populated by separate populations of erythroid-429 myeloid progenitors originating from the yolk-sac and fetal liver, prior to the emergence of 430 circulating monocytes and hematopoietic stem cells. Some existing data suggest that the early 431 yolk sac derived macrophages are eventually entirely replaced by fetal liver derived macrophages 432 capable of self-renewal 86 . In our study, we observed a broad division of the five macrophage 433 populations based upon expression of Dab2 and Plac8, evident both before and after birth.

498
Mouse lung cell isolation. C57BL/6 mice were obtained from Charles River Laboratories. For 499 studies using E18.5, P1, and P7 murine lungs, pregnant dams were purchased, and pups aged prior 500 to lung isolation. At E18.5, dam was asphyxiated with CO2 and pups extracted. At P1, P7, and 501 Sytox Blue (Invitrogen), was added to cells and incubated for 3 min prior to sorting into 384-well 518 plates (Bio-Rad Laboratories, Inc) using the Sony LE-SH800 cell sorter (Sony Biotechnology Inc). 519 FACS sorts were performed with a 100 μm sorting chip (Catalog number: LE-C3110). Prior to 520 cell sorting, the cell sorter and chip were calibrated with SH800 setup beads. Droplet targeting into 521 the middle of four corner and center wells of 384-well plates was manually calibrated. Single color 522 controls were used to perform fluorescence compensation and generate sorting gates. The 384-523 well plates pre-loaded with lysis buffer (Triton X-100 solution, dNTP, poly dT, RNase inhibitor, 524 ERCC, and Triton X-100) were loaded onto the Sony SH800 for gated single cell capture using 525 the ultra-purity mode. Following completion of sorting, 384-well plates containing single cells 526 were spun down, immediately placed on dry ice, and stored at -80C. 527 528 cDNA library generation using Smart-Seq2. Complementary DNA from sorted cells was 529 reverse transcribed and amplified using the Smart-Seq2 protocol on 384-well plates as previously 530 described 20, 95 . Concentration of cDNA was quantified using picogreen (Life technology corp.) to 531 ensure adequate cDNA amplification. In preparation for library generation, cDNA was normalized 532 to 0.4 ng/uL. Tagmentation and barcoding of cDNA was prepared using in-house Tn5 transposase 533 and custom, double barcoded indices 96 . Library fragment concentration and purity were quantified 534 by Agilent bioanalyzer. Libraries were pooled and sequenced on Illumina NovaSeq 6000 with 535 2x100 base kits and at a depth of around 1 million read pairs per cell. 536 537 Data analysis and availability. Sequencing reads were mapped against the mouse genome 538 (GRCm38) using STAR aligner 97 and gene were counted using HTSeq 98 . FZ has been the main 539 maintainer of HTSeq for several years. To coordinate mapping and counting on Stanford's high-540 performance computing cluster, snakemake was used 99 . Gene expression count tables were 541 converted into loom objects (https://linnarssonlab.org/loompy/) and cells with less than 100,000 542 uniquely mapped counts were discarded. Counts for the remaining cells were normalized to counts 543 per million reads. For t-distributed stochastic embedding (t-SNE) 22 , 500 features were selected that 544 had a high Fano factor in most mice, and the restricted count matrix was log-transformed with a 545 pseudocount of 0.1 and projected onto the top 25 principal components using scikit-learn 100 . 546 Unsupervised clustering was performed using Leiden (C++/Python implementation) 21 . Custom 547 Python 3 scripts were used for specific analyses and are available at 548 https://github.com/iosonofabio/lungsc. T Cell receptors were assembled using TraCeR 101 using the 549 default parameters of the Singularity image. B cell receptors were assembled using BraCeR 102 with 550 the parameter -IGH_networks, which agreed with our in-house pipeline consisting of Basic 103 and 551 Change-O 104 . Raw fastq files, count tables, and metadata are available on NCBI's Gene Expression 552 Omnibus (GEO) website (GSEXXXXX). 553 554 In-situ validation using RNAscope and immunofluorescence (IF). Embryonic and post-natal 555 mice were euthanized as described above. E18.5 lungs were immediately placed in 10% neutral 556 buffered formalin following dissection. P1, P7, and P21 murine lungs were perfused as described 557 above, and P7 and P21 lungs inflated with 2% low melting agarose (LMT) in 1xPBS, and placed 558 in 10% neutral buffered formalin. Following 20 hours incubation at 4C, fixed lungs were washed 559 twice in 1xPBS and placed in 70% ethanol for paraffin-embedding. In situ validation of genes 560 identified by single cell RNA-seq was performed using the RNAscope Multiplex Fluorescent v2 561 Assay kit (Advanced Cell Diagnostics) and according to the manufacturer's protocol. Formalin-562 fixed paraffin-embedded (FFPE) lung sections (5 μm) were used within a day of sectioning for 563 optimal results. Nuclei were counterstained with DAPI (Life Technology Corp.) and extracellular 564 matrix proteins stained with hydrazide 105 . Opal dyes (Akoya Biosciences) were used for signal 565 amplification as directed by the manufacturer. Images were captured with Zeiss LSM 780 and 566 Zeiss LSM 880 confocal microscopes, using 405nm, 488nm, 560nm and 633nm excitation lasers. 567 For scanning tissue, each image frame was set as 1024x1024 and pinhole 1AiryUnit (AU). For 568 providing Z-stack confocal images, the Z-stack panel was used to set z-boundary and optimal 569 intervals, and images with maximum intensity were processed by merging Z-stacks images. For 570 all both merged signal and split channels were collected. 571 572 Intracellular Flow Cytometry. P1 and P21 male and female murine lungs were isolated as 573 .5, P1, P7 and P21 murine lungs were processed for flow cytometry and the frequency of (B) total and (C ) CD45+ (immune) dead cells was assessed. (D-E) Frequency of dead cells of P7 murine lungs was determined by flow cytometry following 30 minute (D) enzymatic digestion with 0.38 mg/mL liberase with manual tituration (Liberase) and/or mechanical disruption using the lung program on the GentleMACS dissociator (Lung GM) or (E ) enzymatic digestion with collagenase and manual tituration (Collagenase), Liberase, or 0.38 mg/mL liberase with the Multi_D program on the GentleMACS dissociator (Liberase + Multi_D GM). (F) Number of viable P7 murine lung cells following Collagenase, Liberase, or Liberase + Multi_D GM was quantified. (G) Frequency of murine lung cells at gestational age E18.5, P1, P21 was quantified following incubation with 0.38 mg/mL liberase for 15, 25, and 30 min and manual tituration. To quantify whether the different mice contributed spurious variation to the data, a distribution level approach was chosen. For each cell type and time point, 100 pairs of cells from either the same mouse or between different mice were chosen and the distance in tSNE space calculated. The cumulative distributions for those pairs were subsequently plotted to check whether pairs from different animals had a significantly longer distance than cells from the same mouse.