Single-cell RNA sequencing uncovers heterogenous transcriptional signatures in macrophages during efferocytosis

Efferocytosis triggers cellular reprogramming, including the induction of mRNA transcripts which encode anti-inflammatory cytokines that promote inflammation resolution. Our current understanding of this transcriptional response is largely informed from analysis of bulk phagocyte populations; however, this precludes the resolution of heterogeneity between individual macrophages and macrophage subsets. Moreover, phagocytes may contain so called “passenger” transcripts that originate from engulfed apoptotic bodies, thus obscuring the true transcriptional reprogramming of the phagocyte. To define the transcriptional diversity during efferocytosis, we utilized single-cell mRNA sequencing after co-cultivating macrophages with apoptotic cells. Importantly, transcriptomic analyses were performed after validating the disappearance of apoptotic cell-derived RNA sequences. Our findings reveal new heterogeneity of the efferocytic response at a single-cell resolution, particularly evident between F4/80+ MHCIILO and F4/80− MHCIIHI macrophage sub-populations. After exposure to apoptotic cells, the F4/80+ MHCIILO subset significantly induced pathways associated with tissue and cellular homeostasis, while the F4/80− MHCIIHI subset downregulated these putative signaling axes. Ablation of a canonical efferocytosis receptor, MerTK, blunted efferocytic signatures and led to the escalation of cell death-associated transcriptional signatures in F4/80+ MHCIILO macrophages. Taken together, our results newly elucidate the heterogenous transcriptional response of single-cell peritoneal macrophages after exposure to apoptotic cells.

Scientific RepoRtS | (2020) 10:14333 | https://doi.org/10.1038/s41598-020-70353-y www.nature.com/scientificreports/ Cd11c + MHCII + SPMs require IRF4 for differentiation from monocytic precursors highlighting disparities among peritoneal macrophage subsets 13 . Moreover, these distinct populations rely on disparate molecular mediators to conduct and respond to efferocytosis 16 . Consequently, these MΦ populations would display heterogeneous transcriptional signatures that could not be assigned to phagocyte subpopulations after quantitative PCR or bulk sequencing of mixed macrophage cultures. Indeed the transcriptomics of bulk macrophages during efferocytosis has been characterized in detail 4,17 . Such analyses have defined conserved activation profiles including mobilization of the cytoskeletal network and transcriptional induction of target mRNAs by nuclear receptors, such as Liver X Receptor (LXR) 18 . Yet, our appreciation of how distinct phagocyte subpopulations may uniquely reprogram in response to apoptotic cells has not yet been resolved at the single cell level. Herein, we describe the transcriptomic signature of primary single cells during efferocytosis and after exposure to apoptotic cells. To control for apoptotic cell-derived transcripts 19 , primary macrophages were co-cultivated with apoptotic cells of a separate species. With this experimental design, we are able to discriminate mRNAs derived from the phagocyte versus mRNAs derived from apoptotic cells. In addition, apoptotic cells that expressed a transgenic reporter gene were utilized as a supplemental reporter to track apoptotic cell catabolism. Thus, our results newly elucidate the heterogenous transcriptional response of single-cell peritoneal macrophages, after exposure to apoptotic cells. These data and analyses will be useful for future investigations that seek to determine the unique functional roles of phagocyte subsets during efferocytosis and after exposure to apoptotic cells.
We next performed pathway analysis by directly comparing differentially expressed genes between the F4/80 + MHCII LO MΦs to other clusters individually. Compared to F4/80 + MHCII LO MΦs, F4/80 + MHCII INT MΦs were enriched for pathways associated with efferocytosis such as phagosome, actin cytoskeleton organization, and regulation of cell death (Fig. 1F). In contrast, F4/80 + IL1β + MΦs were enriched in inflammatory networks such as response to lipopolysaccharide, NF-κB signaling pathway, and TNF signaling pathway. Separately, F4/80 − MHCII HI MΦs exhibited upregulated pathways associated with cellular locomotion and chemotaxis, as well as oxidative phosphorylation and electron transport chain pathways (Fig. 1F). Thus at steady state, we newly define 4 resident macrophage populations derived from contrasting ontogenies that perform distinct immune functions within the peritoneum.
Internalized apoptotic RNA dissipates after 6 h. To ascertain the degree of transcriptional heterogeneity of the efferocytic response, we employed a strategy whereby murine peritoneal MΦs were co-cultivated with human Jurkat T cells expressing the GFP transgene under the native Survivin promoter ( Fig. 2A). The difference in species between the two cell types was chosen in order to track human-specific transcripts, such as APOL1, that would mark "passenger" transcripts from the apoptotic cell and therefore permit us to catalog specifically the macrophage transcriptional response. The GFP transgene, which is not naturally expressed in either the mouse or human genome, provided an additional reporter of the apoptotic cell transcriptome. Moreover, the GFP fluorescence also served as a reporter to confirm apoptotic cell engulfment ( Supplementary Fig. S2). As depicted in Fig. 2B, the percent total aligned reads uniquely mapping to the human genome (hg38) within peritoneal immune cells increased 2 h after co-cultivation of human apoptotic cells, with a commensurate decrease in the percentage of aligned reads mapping to the mouse genome (mm10). These reads represent the captured and analyzed portion of the transcriptome during acute exposure to apoptotic cells. Unaligned reads, which can indicate the fraction of a library with lower quality sequencing or increased RNA degradation 22 , were then separately mapped to hg38. After 2-h of co-culture, the percent of initially unaligned reads were enriched for sequences uniquely mapping to hg38 (66.1%). This unaligned portion of the library is typically excluded from standard transcriptomic analyses of gene expression. We next mapped this same bin of unaligned reads to GFP and APOL1 sequences in order to assess apoptotic transcripts in the phagocyte that would not have overlapping reads aligning to mm10. The number of raw reads aligning to these sequences was highest (113 raw counts) after 2-h of co-culture (Fig. 2B). According to our defined parameters (as detailed in "Methods"), macrophages without GFP or APOL1 reads were defined to have reached a stage of efferocytosis in which apoptotic cell mRNAs    www.nature.com/scientificreports/ would likely be degraded. It was at this time point at 6 h that we focused the majority of our transcriptional analysis.

Single-cell trajectories reveal distinct transcriptional responses between F4/80 + and MHCII HI macrophages after efferocytosis. With our inference that each cell represented in our analysis was
at a unique stage of efferocytosis, we performed a pseudotemporal reconstruction using Slingshot, a tool that infers global lineage structures and pseudotime variables for each cell along a lineage 23 . Two distinct trajectories were observed which corresponded to unique transcriptional changes in our F4/80 + LPMs and MHCII HI SPMs (Fig. 2C,D). Gene expression was subsequently plotted as a function of pseudotime to track changes across separate MΦ efferocytic states (Fig. 2C,D). While expression of genes such as Cd36 and H2-Ab1 increased as a function of pseudotime in both MΦ subsets, expression of Il1b and Tgfbi was unique between LPMs and SPMs. In LPMs, Il1b expression is greatly decreased as a function of pseudotime, but in SPMs, Il1b expression increases marginally as a function of pseudotime. Furthermore, expression of Il1rn, a potent antagonist of the IL1 receptor, was greatly increased in LPMs but decreased in SPMs during efferocytosis (Fig. 2D). These results reveal divergent transcriptional programming of MΦ subpopulations after exposure to apoptotic cells modeled using pseudotime trajectory analyses. While visualizing gene expression as a function of pseudotime provides insights into how single cell gene expression matures temporally, determining differentially expressed genes remains a statistically difficult challenge 24 . Therefore, we leveraged our physically separated co-culture time points (Fig. 3A,B) to determine differentially expressed genes before and after apoptotic cell exposure. After 6 h of apoptotic cell co-cultivation, F4/80 − MHCII HI SPMs are enriched while F4/80 + MHCII LO MΦs and F4/80 + IL1β + MΦs represent a smaller percentage of the captured MΦ population (Fig. 3B). Genes which had distinct expression among MΦ subsets as a function of pseudotime were confirmed to have similar expressional changes after comparison with data from separate co-culture timepoints (Fig. 3C). Pathway analysis of differentially expressed genes in each MΦ subset before and 6 h after co-culture with apoptotic cells revealed distinct global transcriptional changes after efferocytosis (Fig. 3D,E). Overall, LPMs were enriched for homeostatic and developmental pathways such as tissue development, regulation of development, and cellular homeostasis (Fig. 3D). SPMs, however, uniquely downregulated inflammatory pathways such as response to LPS and positive regulation of NF-κB signaling (Fig. 3E), relative to non-treated MΦs. Interestingly, SPMs also enriched transcription of genes associated with cell proliferation (Fig. 3E), suggesting the increase in SPMs in Fig. 3B may be due to induced cell proliferation. Lastly, to investigate acute transcriptional changes of efferocytosis, gene set enrichment analyses were also performed after 2 h of co-cultivation with ACs. While F4/80 + LPMs significantly increased expression of genes for oxidative phosphorylation (Supplemental Fig. S2), MHCII HI SPMs, conversely, significantly reduced expression of genes for oxidative phosphorylation (Supplemental Fig. S2).

Discussion
Using single-cell RNA sequencing to unbiasedly investigate how exposure to apoptotic cells reprograms phagocytes, our findings reveal a heterogeneity of efferocytic transcriptional responses in resident macrophages not previously appreciated by traditional bulk mRNA sequencing approaches. Importantly, our approach was careful to exclude to the best of our ability the contribution of apoptotic cell-derived transcripts. This is particularly important as low abundant "passenger" RNAs from apoptotic cells could misrepresent the phagocyte transcriptional response, as previously reported 19 .
Our analysis uncovered four transcriptionally distinct phagocyte populations and defined their unique gene expression signatures after co-cultivation of apoptotic cells. A subpopulation of LPMs (F4/80 + MHCII LO ) was enriched for Klf2 21 , Timd4 26 , and Mertk 26 , implicating a predisposition for efferocytosis, as these genes encode known regulators of apoptotic cell clearance. Furthermore, the zinc finger transcription factor GATA-binding protein 6 (Gata6), a necessary regulator of peritoneum-specific gene expression program in LPMs 15 27 . Alternatively SPMs have been shown to promote inflammation particularly after thioglycolate-elicitation. However, our data newly identifies a highly inflammatory LPM population (F4/80 + IL1β + ), a population that may have been unidentifiable in previous studies of bulk MΦs. Whether this inflammatory LPM population is the key driver of inflammation within the physiological peritoneal cavity warrants further investigation. Overall, defining differences between ontogenetically divergent MΦs may elucidate mechanisms by which resident peritoneal MΦs uniquely contribute to tissue homeostasis and physiological inflammation. Upon challenge with apoptotic cells, macrophage subpopulations induced unique transcriptional responses, consistent with the premise that macrophage ontogeny predisposes their response to apoptotic cells. Previously, Cain et al. determined that LPMs are more efficient at efferocytosis than SPMs, but importantly both subsets exhibited apoptotic cell engulfment one hour after exposure 28 . Furthermore, the vast majority of apoptotic cell-derived transcripts are degraded after 6 h of co-cultivation, and thus no phagocytes were removed from the downstream transcriptional analyses at this time point. Consequently, this experimental design does not exclude "late-eater" phagocytes that may have initiated apoptotic engulfment at a later time. Of note, acute transcriptional changes in macrophages after 2 h of co-cultivation are subtle relative to transcriptional responses 6 h after co-cultivation. This observation is consistent with the paradigm that acute phagocyte responses activate post-transcriptional cytoskeletal remodeling, while phagocyte transcriptional reprogramming is a subsequent event. Therefore to better model the transcriptional landscape throughout differential phases of efferocytosis, we employed a pseudotime trajectory analysis using Slingshot. Since these algorithms rely on transcriptional differences independent of time to produce a pseudotime trajectory, the modeled response and reprogramming trajectories during efferocytosis are better synchronized to allow for comparison between macrophage clusters. Unsurprisingly, both LPMs and SPMs revealed distinct transcriptional programs upon exposure to apoptotic cells. LPMs robustly reduced their inflammatory gene expression as a function of pseudotime, suggesting that exposure to apoptotic cells polarizes LPMs towards an anti-inflammatory state. Furthermore, many of the differentially expressed genes in LPMs after apoptotic cell exposure for 6 h are enriched in homeostatic and developmental pathways. These data are consistent with apoptotic cell clearance by LPMs as an essential mechanism for tissue homeostasis. Conversely, SPMs reduce enrichment of inflammatory pathways such as NF-kB signaling and response to LPS, while increase enrichment of cell proliferation pathways after apoptotic cell exposure (Fig. 3E). These results may lead to the observation of an increase in MHCII HI SPMs after 6 h of AC co-cultivation. Our careful approach using pseudotime algorithms combined with discrete time points allow us to compare the transcriptional changes during efferocytosis between the distinct macrophage populations.
We also leveraged macrophages deficient in Mertk, a key efferocytosis receptor 25 to validate and investigate the effects of impaired apoptotic cell clearance on peritoneal macrophage subsets. Mertk is primarily expressed in LPMs and Mertk-deficient macrophages exhibited acute reductions in F4/80 + macrophages (Fig. 4D). At steadystate, Mertk-deficient LPMs exhibited significantly increased MHCII expression and pro-inflammatory cytokine expression. Upon challenge with apoptotic cells, Mertk −/− LPMs exhibited a signature of programmed cell death (Fig. 4F). These data are consistent with previous studies implicating Mertk as a mediator of cell survival after the stress of efferocytosis [29][30][31] . Intriguingly, expression of Mertk is associated with promoting phagocyte survival after exposure to oxidative stress through the upregulation of anti-apoptotic signals including pAkt and pErk 29 . Previously it has been shown that Abca1 reverse cholesterol efflux pathway dampens oxidative stress preserving phagocyte viability after exposure to apoptotic cells or oxidized phospholipids 32 . Without the beneficial role of Abca1 reverse cholesterol efflux, engulfment of apoptotic cells induced an excessive oxidative burst leading to activation of the apoptotic cell death program. Therefore, our findings suggest that Mertk-expression specifically induces cellular homeostatic responses to oxidative stress in LPMs that preserve phagocyte viability after efferocytosis.
Future studies are warranted to examine the heterogeneity of efferocytic transcriptional responses in situ and in vivo. This could be achieved by coupling single cell-sequencing approaches with tissue phagocytes that are positive for the engulfment of reporter apoptotic cells. However, under these conditions, it would be difficult to exclude the contribution of apoptotic cell-derived transcripts. One potential circumvention would be the employment of Ribotag technology 33 to capture newly formed transcripts in the phagocyte. While scRNAseq has emerged as a powerful and insightful technique, additional studies are also needed to validate the transcriptional signatures of our findings at the protein level, possibly through technologies such as CyTOF or high-dimensional flow cytometry. With these future directions in mind, our studies newly elucidate the heterogenous transcriptional response of single-cell peritoneal macrophages, after exposure to apoptotic cells. These findings will benefit future studies that seek to determine the unique functional roles of phagocyte subsets during efferocytosis.

Methods
Animal studies. C57BL/6J mice were used as WT controls and bred in the Northwestern Center for Comparative Medicine facility. MerKD (referred to herein as MerTK −/− ) have been previously described 25 and were backcrossed to C57BL/6J mice for ten generations. Eight-to twelve-week old mice were utilized for experiments. Mice were bred and housed in a pathogen-free, temperature-and humidity-controlled environment with access to standard mouse chow and water ad libitum. Mice were kept on a 12:12 light/dark cycle. All studies Scientific RepoRtS | (2020) 10:14333 | https://doi.org/10.1038/s41598-020-70353-y www.nature.com/scientificreports/ were approved and reviewed by the Institutional Animal Care and Use Committee at Northwestern University (Chicago, Illinois), protocol #IS00000378. All experiments were performed in accordance with the guidelines and regulations named in the protocol.
Ex vivo efferocytosis assay, microscopy, and flow cytometry. Resident peritoneal cells were harvested after lavage with cold saline. Peritoneal cells from 3 mice were pooled together for each experimental group. Initial cell selection was achieved through adherence to non-treated, low-adherence cell culture plates for 1 h and rinsed to remove non-adherent cells. Apoptotic cells (ACs) were generated using GFP-Jurkat Human T cells (GenTarget) exposed to UV radiation for 7 min followed by a 2-h incubation at 37 °C. Apoptosis was confirmed by annexin V positive, propidium iodide negative identification affirming greater than 80% apoptosis. Adherent resident peritoneal cells were co-cultivated with ACs at a ratio of 5 ACs to 1 peritoneal cell for 2 h or 6 h as indicated. Control cells were given a media change corresponding to the 6-h timepoint. Non-engulfed ACs were removed from the co-culture through rigorous rinsing with warm saline. Adherent macrophages were removed from the plate with Accutase (StemCell Technologies) and resuspended into a single cell suspension. Apoptotic cell engulfment was confirmed with confocal fluorescence microscopy. For in vivo efferocytosis microscopy, fluorescently labelled (green calcein-AM) apoptotic Jurkat cells were injected into the peritoneum of Mertk +/+ versus Mertk −/− mice and subsequently peritoneal lavages were collected and stained for F4/80 + resident macrophages. Percent double positive F4/80 + calcein-AM + cells as function of total F4/80 + cells were enumerated. *p < 0.03.
Single cell library preparation and RNA sequencing. The single cell RNA-Seq libraries were prepared using the 10X Genomics Single Cell 5′ Gel Bead and Library Kit pipeline following manufacturer's protocols. Cell suspensions were diluted to target a recovery of 4,000 cells per sample. A total of 7,674 resident peritoneal cells were sequenced to a read depth ~ 100,000 reads per cell. The Illumina libraries were run on an Agilent Bioanalyzer High Sensitivity Chip and Kapa Library Quantification Kits for Illumina platform (KAPA Biosystems) for quality control before sequencing. In collaboration with the Northwestern University Sequencing Core (NuSeq), the libraries were sequenced on the Illumina HiSeq 4,000 with the following parameters: Read 1-26 cycles; i7 Index-8 cycles; Read 2-98 cycles.
Cell ranger, read alignment, and quantifying cells with 'unaligned reads'. The sequenced data were processed with the Cell Ranger Single Cell software suite 1.3.1 by 10 × Genomics (GEO accession number pending). Briefly, raw base-call files from a HiSeq4000 sequencer were demultiplexed into FASTQ files. FASTQ files from each of the samples were mapped and the genes were counted using cellranger count. Given that the sequencing parameters included paired-end reads and a full length (150 bp) second read, the "SC5P-PE" option was used for the chemistry parameter. The sequencing reads were aligned to both hg38 (Homo sapiens) and mm10 (Mus musculus) reference genomes using STAR aligner. Reads mapped to mm10 were used for downstream analysis with Seurat.
The unaligned reads were extracted from the 'possorted_bam.bam' file created by CellRanger using Samtools with the '-f 4' flag. Awk was then used to remove any duplicate reads with the '!seen[$0] +' distinction. The unmapped read ids were then used as a pattern to grep the FASTQ reads from the raw FASTQ file with 'seqkit' for both read one and two. These reads were then mapped to a custom genome containing only the APOL1 genes using the STAR aligner. After alignment, the mapped reads were extracted using Samtools view with the '-F 4' flag and deduplicated using 'awk' with the '!seen[$0]++' distinction. For the purpose of determining the number of filtered cells containing the human DNA; the mapped read ids were grepped by 'seqkit' to extract the reads from the raw fastq files that had been mapped to APOL1 followed by the 10 × cell barcodes being identified using 'sed' with the '2 ~ 4p' distinction for the '-n' parameter. Lastly, the 10 × barcodes from the Seurat object were extracted and compared using 'comm' to the read extracted 10X barcodes to count the number of cells containing the human APOL1 DNA.

Single-cell RNA sequencing analysis and visualization. The barcode-gene matrices from the Cell
Ranger pipeline were further analyzed using the Seurat R package (v.3.1) 34 . Following standard practices to exclude low-quality cells, cells that expressed less than 500 genes, had a ratio of unique molecular identifiers (UMIs) per gene below 0.8, or had greater than 10% mitochondrial genes were filtered from the datasets. Furthermore, each cell was assigned a Cell Cycle Score to address potential differences in cell cycle phase.
Normalization and variance stabilization of molecular count data was carried out using the R package sctransform (v0.2.1) 35 . After transformation, each sample was integrated using Seurat's modified integration workflow to harmonize the output Pearson residuals across datasets. To reduce the dimensionality of the expression matrices, we performed a principal component analysis (PCA). The number of principle components for further downstream applications was decided to be 40, as UMAP was employed for dimensionality reduction and visualization of the data.
Clusters were identified using graph-based clustering approach implemented by the FindCluster function in Seurat at resolutions of 0.2, 0.4, 0.6, 0.8, and 1.0 to determine transcriptionally distinct populations of resident peritoneal cells. Clusters with high cell cycle scores were excluded from further downstream analyses. To unbiasedly identify resident peritoneal cells present in the dataset, SingleR (v1.0.5) was employed. Briefly, SingleR infers the origin of each individual cell by referencing transcriptomic datasets of pure cell types. We utilized the ImmGen database, which contains normalized expression values for immune cells from 830 murine microarrays to ID our peritoneal cell types. These classifications were confirmed with canonical immune cell markers.
Scientific RepoRtS | (2020) 10:14333 | https://doi.org/10.1038/s41598-020-70353-y www.nature.com/scientificreports/ For all differential expression tests, we utilized the model-based analysis of single-cell transcriptomics (MAST) test in the Seurat package. To identify unique expression profiles for each cluster, differential expression was tested between each macrophage cluster and all other macrophage clusters combined. The top 50 differentially expressed genes (based on average log fold change) unique to each cluster were visualized in a heatmap. Differential expression testing was also performed between two individual clusters, between genotypes, and between timepoints.
Pathway enrichment analysis. To identify enriched molecular pathways based on differentially expressed genes, gProfiler 36 was applied. Gene sets from Gene Ontology biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used. GOplot (v1.0.2) was employed for visualization of the GO enrichment analyses.
Single-cell trajectory inference. Slingshot (v1.4.0) 23 was applied to infer efferocytosis reprograming trajectories of macrophages which experienced apoptotic cell engulfment. Slingshot combines highly stable techniques to infer novel lineages of the global structure. Each macrophage subpopulation (F4/80 + and MHCII +) were isolated into separate single-cell experiments and reduced dimensionality via both PCA and UMAP. The UMAP trajectories were utilized for visualization purposes, while the PCA trajectories were utilized for gene expression plots as UMAP dimensionality reduction has been shown to distort the pseudotime data 37 . Changes in gene expression as a function of pseudotime were visualized using the PlotExpression command.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. The datasets are available in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) repository (GSE156234).