To dissect direct and indirect responses to transfected immunogenic DNA stimulation with single cell resolution, we designed the following experiment. We used single cell RNA-Seq to define the heterogeneity of cultured BMDMs (i) at steady state, (ii) in response to dsDNA alone, or (iii) in response to dsDNA + IFNAR Ab inhibition of secondary signaling, or (iv) in dsDNA-treated IFNAR KO ( Figure 2A ). BMDMs were derived by flushing hematopoietic cells from the femurs of adult C57BL6 mice and differentiating them into macrophages in culture using monocyte colony stimulating factor (m-CSF) as we previously described (20). We then subjected the cultures to 24 hours of stimulation with the conditions above, after which BMDMs were collected, stained with DAPI, purified by flow cytometry, and subjected to droplet-based microfluidic single cell barcoding (Chromium, 10X Genomics), library preparation, and next generation sequencing (NovaSeq, Illumina). After demultiplexing and mapping to a reference murine genome, we performed shared nearest neighbor (SNN) clustering on the integrated data from all three treatment groups to define a universal embedding of shared BMDM subtypes ( Figure 2B ). Unbiased clustering revealed 6 coarse BMDM clusters, 2 of which were interpreted to be proliferating based on expression of canonical cell replication marker genes and labeled them Repl(G2/G1) and Repl(S/G1) ( Supplementary Figure S1 ). Proliferating BMDMs were stratified by cell cycle phase cell based on their expression of Mcm6 and Mcm5 (S phase), Top2a and Cenpa (G1 phase), Birc5 and Nusap1 (G2 phase) ( Figures 2C, D ). Non-proliferating BMDMs were also readily stratified into 3 clusters based on their expression of Fabp4 and Fabp5 and MHCII genes such as H2-Aa and H2-Ab1 ( Figures 2C, D ). We labeled the non-proliferating BMDM subsets (Fabp4-,MHCII-), (Fabp4+, MHCII-), and (Fabp4+, MHCII+). The final cluster was defined as a monocyte-like BMDM Progenitor characterized by expression of S100a4 and S100a6.

Having defined the BMDM subtypes at steady state, we turned our attention to the response to dsDNA. Experiments based on bulk measurement techniques (e.g., qPCR) have firmly established that dsDNA complexed with a transfection reagent, when delivered to macrophages in culture, induces expression of type I IFNs and ISGs (21). However, underlying single cell heterogeneity remains incompletely defined. Consistent with results from bulk experiments, single cell RNA-Seq profiling of dsDNA-stimulated BMDMs markedly induced expression of ISGs (Irf7, Ifit2, Oasl2, Rsad2, Isg15) compared to unstimulated controls as well as compared to dsDNA in the context of IFNAR Ab, thus demonstrating successful inhibition of secondary amplification ( Figure 3A ). Violin plots confirmed dsDNA-dependent and IFNAR-dependent induction of ISGs at the single gene level while broader macrophage marker genes such as Ms4a7 remained highly expressed across all three conditions ( Figure 3B ). To evaluate dsDNA responses independent of clustering, we created an ISG Score based on the summed expression of a collection of ISGs derived from previously published studies ( Table 1 ). dsDNA-treated BMDMs had the highest ISG Scores, which were almost completely abrogated by IFNAR Ab inhibition ( Figure 3C ). Interestingly, we did not find that dsDNA induced ISGs preferentially expressed by one or more of the BMDM subtypes. Instead, they were uniformly elevated across all proliferative and non-proliferative subsets compared to control ( Figure 3D ). Although we included IFNAR KO BMDMs, they were found to have significant ISG expression in scattered cells with or without DNA stimulation, which limited their utility for identifying direct DNA-stimulated cells ( Supplementary Figure S2 ). Taken together, these data demonstrated that dsDNA stimulation induces ISG expression in both proliferating and non-proliferating BMDMs, and that secondary amplification of the response via type I IFN secretion plays a dominant role in the ISG response since it is markedly abrogated by simultaneous blockade with a neutralizing IFNAR Ab.

Since IFNAR Ab treatment markedly inhibited secondary amplification of DNA-induced IFN responses, we reasoned that any residual ISG expression must result from IFNAR-independent responses to DNA. To identify direct DNA-sensing cells, we bioinformatically combined BMDMs from control and IFNAR Ab treated cells and reclustered ( Figure 4A , Supplementary Figure S3 ). This revealed a small group of BMDMs with high ISG scores despite inhibition of IFNAR-dependent amplification ( Figure 4B ). We interpreted these ISG-expressing BMDMs to be direct responders to dsDNA stimulation. Unbiased clustering of IFNAR Ab-treated Condition 3 cells revealed 8 distinct clusters including one expressing ISGs embedded within the broader Fabp4+MHCII- cluster, as illustrated in a feature plot ( Figure 4C ) and heatmap of top marker genes expressed by each cluster ( Figure 4C ). The top marker genes for directly DNA-stimulated BMDMs were Ccl5, Irf7, Ifit1, Isg15, Rsad2, Oasl2, Bst2, and Ifi204 ( Figure 4D ). To confirm the result independent of clustering, we plotted the ISG Score for each DNA- and IFNAR-Ab-treated BMDM of Condition 3 and identified the same small population of Fabp4+MHCII-ISG+ cells ( Figure 4E ). Together, these results suggest that direct DNA-sensing BMDMs are preferentially Fabp4+MHCII- macrophages.

Bulk measurement techniques do not allow separation of direct DNA-sensing cells from the more abundant secondary IFNAR-dependent cells. This limits separation of cell-specific roles in initiation and amplification of IFN responses in vivo. To identify genes that may discriminate direct and indirect cells, we identified the most highly differentially expressed genes between direct DNA responsive cells and secondary IFNAR-dependent cells. First, we isolated Fabp4+MHCII- macrophages, which contained the majority of directly responding cells. We reclustered them, which revealed an ISG+ and ISG- population. Next, we combined the direct macrophages from condition 3 (dsDNA + IFNAR Ab) with the Fabp4+MHCII- macrophages from condition 2 (dsDNA Alone) and reclustered ( Figure 5A , Supplementary Figure S4 ). We displayed expression of 45 ISGs as a heatmap. This revealed Ccl5 to be the most differentially expressed between the direct cells from condition 3 and combined indirect and direct cells from condition 2 ( Figure 5B ). Feature plots show the cells colored by ISG Score compared to the more stringent by Direct DNA-sensing Score ( Figure 5C ). We selected the top 4 genes and created a Direct Score (Ccl5, Cdk8, Cxcl2, and Cd74) ( Figure 5D ). Based on these results, we propose that the top gene, Ccl5, may represent a useful marker gene for identifying direct DNA sensing cells.

Mouse experiments were approved and conducted under the oversight of University of California San Diego Institutional Animal Care and Use Committee (IACUC #17144) using Adult C57BL/6J (3 WT mice, stock 000664, and 2 IFNAR -/- mice, stock 028288). All experiments were performed with 10 to 14-week-old animals and were carried out using age and gender matched groups without randomization. All mice were maintained in a pathogen-free environment of UC San Diego.

Bone marrow derived macrophages (BMDMs) were isolated by flushing femurs of adult mice and culturing the resulting cells in 10% FBS 1% Pen/Strep-containing DMEM supplemented with 10ng/mL recombinant m-CSF (Peprotech) for 7 days. 5μg of immunogenic HT DNA (In vivogen) was complexed with a Lipofectamine transfection agent (ThermoFisher) at a ratio of 1:1.5 in serum free-media and added to 1 million BMDMs in a 6-well multi-well plate with serum- and mCSF-containing media for each experiment. Two independent experiments were performed, yielding n = 5 for each condition. For inhibition of secondary signaling via the IFNAR receptor, cells were treated with 20μg/ml of MAR1-5A3 IFNAR neutralizing antibody or isotype control (BioXCell).

Isolated cells were stained at 4°C in FACS buffer (PBS supplemented with 2.5% bovine serum albumin). Cell suspensions were labeled with DAPI just prior to flow cytometric analysis to allow exclusion of dead cells. Doublets and dead cells were excluded by forward scatter and Dapi. Data was acquired by Sony sorter MA900 at UCSD and analyzed with FlowJo software.

Single cell RNA-Seq was performed by microfluidic droplet-based encapsulation, barcoding, and library preparation (10X Genomics) as previously described (27). Paired end sequencing was performed on an Illumina NovaSeq instrument. Low level analysis, including demultiplexing, mapping to a reference transcriptome (Ensembl Release 85 - GRCm38.p5), and eliminating redundant UMIs, was performed with the CellRanger pipeline.

To account for variations in sequencing depth, total transcript count for each cell was scaled to 10,000 molecules, and raw counts for each gene were normalized to the total transcript count associated with that cell and then natural log transformed. Cells with between 200 and 2,000 uniquely expressed genes and < 1% mitochondrial counts were retained for further analysis. Highly variable genes across individual datasets were identified with the FindVariableFeatures method from the Seurat R package (version 3.0) by selecting 3,000 genes with the highest feature variance after variance-stabilizing transformation. Integration of multiple single-cell RNA-seq datasets was performed in Seurat to enable harmonized clustering and downstream comparative analyses across conditions (28–30). Anchoring cell pairs between datasets were identified by Canonical Correlation Analysis (CCA) and the mutual nearest neighbors (MNN) method using the Seurat FindIntegrationAnchors function.

After scaling and centering expression values for each variable gene, linear dimensionality reduction was performed on integrated data using principal component analysis (PCA). Clustering was performed using the shared nearest neighbor (SNN) clustering algorithm with the Louvain method for modularity optimization, as implemented in the Seurat FindNeighbors and FindClusters functions. To visualize data in two-dimensional space, Uniform Manifold Approximation and Projection (UMAP) dimensional reduction was performed. Differentially expressed genes (DEGs) between clusters were determined using a Wilcoxon Rank Sum test. Where specified, subsets of cells were isolated and reclustered to identify new DEGs.

ISG Scores were measured as the sum of the raw reads for the ISGS: Ifi204, Isg15, Irf7, Ifit1, Rsad2, Oasl2, Oas3, Gpb5, Bst2, Ifit2, Ifit3, Oas2, Isg20, Cxcl10, Oasl1, Gpb2, Ccl5, Ccl3. Direct Scores were measured as the sum of the raw reads for: Ccl5, Cdk8, Cxcl2, and Cd74. ISG and Direct Scores were normalized to reads per cell and scaled by 10⁴.

Statistical analysis was performed using GraphPad Prism software. All data are represented as mean values +/- standard error of mean (S.E.M.) unless indicated otherwise. A statistical method was not used to predetermine sample size. All analyses were unpaired. P values are indicated by P values less than 0.05 were considered significant and are indicated by asterisks as follows: *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001.