TREM2+ and interstitial-like macrophages orchestrate airway inflammation in SARS-CoV-2 infection in rhesus macaques

Study overview

An overview of the study design is shown in Fig. 1a. We analyzed three separate cohorts of macaques: Cohort 1, Baricitinib-treated and Cohort 2. For Cohort 1 and Baricitinib-treated, a total of eight RMs (mean age 14 years old; range 11–17 years old) were inoculated intranasally and intratracheally with 1.1 × 10⁶ plaque-forming units (PFU) of SARS-CoV-2 (2019-nCoV/USA-WA1/2020). At 2 dpi, four of the eight animals started receiving baricitinib²⁰. For this study, pre-infection baseline and hyperacute time points (1–2 dpi) include n = 8 RMs, all untreated, and the remaining longitudinal time-points assessed to determine the pathogenesis of SARS-CoV-2 infection are comprised of n = 4 of the RMs that remained untreated. Inoculation with SARS-CoV-2 led to reproducibly high viral titers detectable in the upper and lower airways by genomic and sub-genomic qPCR assays (Fig. 1b). The peak of viremia in the nasal passage, throat and BAL was at 2–4 dpi (Fig. 1b). To increase the power of our scRNA-Seq and flow cytometry experiments, we analyzed an additional six macaques infected with the same dose and strain of SARS-CoV-2 (2019-nCoV/USA-WA1/2020) (mean age 10.5 years old; range 6–19.5 years old), referred to as Cohort 2. Animals from Cohort 2 served as SARS-CoV-2-infected, untreated controls in part of a larger study testing the impact of interferon blockade²⁶.

a Study design; RMs were infected intranasally and intratracheally with SARS-CoV-2 and tracked longitudinally (Cohort 1: n = 4, baricitnib cohort: n = 4, Cohort 2: n = 6). Baricitinib was administered daily to 4 RMs (baricitinib cohort) starting at 2 dpi. The 4 RMs  from Cohort 1 and 6 RMs from Cohort 2 were untreated. (Created with BioRender.com). b After SARS-CoV-2 inoculation, nasal, throat, and bronchoalveolar lavages (BAL) were collected and viral loads were quantified by qRT-PCR for total gRNA and sgRNA. c Longitudinal levels of monocytes within BAL and blood depicted as a % of CD45+ cells. p values: CD14−CD16+ Monocytes Blood: 1 dpi vs. −5 dpi = 0.03 and 2 dpi vs. −5 dpi = 0.004; CD14+CD16+ Monocytes Blood 2 dpi vs. −5 dpi = 0.004. d Longitudinal levels of plasmacytoid dendritic cells (pDCs) within BAL and blood depicted as a percentage of CD45+ cells. p value for pDCs Blood 2 dpi vs. −5 dpi = 0.02. e Longitudinal levels of NK cells expressing Granzyme B in BAL and blood. p values for Granzyme B+ NK Cells 2 dpi vs. −5 dpi BAL = 0.01 and blood = 0.008. n = 8 RM from Cohort 1 and baricitinib cohort. The red bars represent the mean. Statistical analysis was performed using one-tailed Wilcoxon signed-rank test in Graphpad Prism v7.02 comparing each timepoint to −5 dpi. *p value <0.05, **p value <0.01. Source data (b–e) are provided as a Source Data file.

SARS-CoV-2 induces a robust, but transient, expansion of pDCs during hyperacute infection

To characterize the innate immune response following SARS-CoV-2 infection, we analyzed changes in innate populations using multi-parametric flow cytometry in blood and BAL samples in the first 2 dpi, or “hyperacute” phase of infection (Fig. 1c–e and Supplementary Fig. 1), and over the full course of infection (Supplementary Fig. 2). In blood, we did not observe a significant increase in the proportion of classical monocytes (CD14+CD16−) at 2 dpi (Fig. 1c) nor at extended time-points (Supplementary Fig. 2a, d). Similar to reports in humans⁹, we observed a rapid, but transient, increase in blood CD14−CD16+ and CD14+CD16+ monocytes (Fig. 1c and Supplementary Fig. 2a, d). Using these conventional markers for blood monocyte subsets, we did not observe any significant changes in CD14−CD16+, CD14+CD16+, nor CD14−CD16+ within the BAL (Fig. 1c and Supplementary Fig. 2a).

We observed a significantly elevated level of pDCs in blood at 2 dpi and similarly, a trend of elevated pDCs in BAL samples (Fig. 1d and Supplementary Fig. 2c, f). This expansion was transient, as pDC numbers returned to baseline by 4 dpi. While the overall frequencies of natural killer cells (NK) were not changed in blood or BAL (Supplementary Fig. 2b, e), the fraction and absolute number of Granzyme B + NK cells increased significantly at 2 dpi in blood, from 4 to 25% (Fig. 1e) and remained elevated throughout the course of infection (Supplementary Fig. 2b, e). Similarly, increases in NK cell activation were also observed in the BAL, rising from 12 to 33% at 2 dpi (Fig. 1e), and persisting at this level until the study termination at 10/11 dpi (Supplementary Fig. 2b, e). Collectively, these data indicate that during the hyperacute phase of SARS-CoV-2 infection, there is a significant mobilization of innate immune cells capable of initiating and orchestrating effector responses of the Type I IFN system.

SARS-CoV-2 infection drives robust, but transient, upregulation of IFN responses in blood and lower airway

To understand the extent of immunological perturbations induced by SARS-CoV-2 infection, we performed extensive gene expression profiling of PBMC and BAL samples. During the hyperacute phase, the BAL had widespread induction of pathways associated with innate immunity and inflammation (Fig. 2a). Notably, we observed a rapid and robust induction of interferon-stimulated genes (ISGs) in the PBMC and BAL compartments starting at 1 or 2 dpi (Fig. 2b and Supplementary Fig. 3a). The ISG response, although widespread, had largely returned to baseline by 10/11 dpi (Fig. 2b and Supplementary Fig. 3a). We also detected a trend of elevated IFNα protein in 4/6 and 5/8 animals in BAL and plasma, respectively (Fig. 2c, d) and a significant increase in RNA-Seq read counts mapping to IFNA genes at 2 dpi in BAL (Fig. 2e), which coincided with the expansion of pDCs in the airway and blood (Fig. 1d). A significant enrichment of genes representing NK cell cytotoxicity (Fig. 2a) was observed at 2 dpi in BAL, consistent with our observation of elevated Granzyme B + NK cells by flow cytometry (Fig. 1e). Taken together, these data demonstrate the presence of primary cells able to produce Type I IFNs (i.e., pDCs), coincident with detectable IFNA transcripts and protein, and with downstream IFN-induced effector functions (ISGs, NK cell activation) following SAR-CoV-2 infection, and that these responses were transient, having largely subsided by 10/11 dpi.

n = 8 RM from Cohort 1 and baricitinib cohort for −5 dpi and 2 dpi except for (c) where n = 6 (3 RM Cohort 1 + 3 RM baricitnib cohort) at 2 dpi. n = 4 RM from Cohort 1 starting from 4 dpi. a Dot plots showing normalized enrichment scores and nominal p values for gene sets. Enrichment is indicated by dot color (red: positively enriched vs. −5 dpi; blue: negatively enriched), dot size indicates significance. Normalized enrichment scores and nominal p values were determined using GSEA⁶⁵. The exact nominal p values are included in Supplementary Data 1. b Heatmap of longitudinal responses for the ISG gene set. The color scale indicates log2 expression relative to the median of the −5 dpi samples. c Cytokines evaluation (Mesoscale) in BALF p values for 2 dpi vs. −5 dpi: IL6 = 0.02 and IP-10 = 0.03 and d Plasma; only significant cytokines are shown (p values for 2 dpi vs. −5 dpi: IP-10 = 0.008, MCP-1 = 0.008, MCP-2 = 0.004). e Sum of normalized expression of all IFNA genes in BAL. *p value = 0.04. f GSEA enrichment plot showing negative enrichment for AM gene signature (derived from SingleR) when comparing bulk BAL RNA-Seq samples from 4 dpi to −5 dpi (Nominal p value = 0). The red bar represents the mean. Statistical analysis for (c)–(e) was performed using one-tailed Wilcoxon signed-rank test in R v4.2.2 comparing each timepoint to −5 dpi. *p value <0.05, **p value <0.01. Source data (c–e) are provided as a Source Data file.

SARS-CoV-2 infection drives a shift in airway macrophage populations

We observed that SARS-CoV-2 infection induced significant enrichment of several inflammatory cytokine signaling pathways, namely IFNA, IL4, IL6, IL10, IL12, IL23 and TNF, and the chemokine pathways CXCR4 and CXCR3, in both PBMCs and BAL of RMs, with higher magnitude in the BAL (Fig. 2a and Supplementary Data 1–3). For many of these pathways, we were able to quantify significant increases in the upstream regulator at either the protein, or mRNA level, or both: IL6 protein levels were significantly increased in the BAL fluid (BALF) (Fig. 2c), as were RNA transcripts in BAL (Supplementary Fig. 3b). Similarly, the induction of CXCR3 pathways signaling was consistent with detection of increased IP10/CXCL10 protein in BALF and RNA at 2 dpi in BAL (Fig. 2c and Supplementary Fig. 3b). The appearance of inflammatory pathways in the blood and airway have been reported in a multitude of human studies (reviewed in ref. ²⁷). However, we noted that SARS-CoV-2 infection also drove early expression of several immunoregulatory/immunosuppressive pathways in the BAL, namely: PD1 and CTLA4 signaling, and negative regulators of MAP kinase and DDX58/RIG-I signaling (Fig. 2a). Previously, we reported that the myeloid fraction in BAL was primarily responsible for the production of pro-inflammatory mediators, however the specific immunophenotypes were not defined. To further investigate the presence of different macrophage subsets within the lower airway after SARS-CoV-2 infection, we performed GSEA on bulk BAL data using AM gene signature (obtained from SingleR²⁸) specific for RM pulmonary macrophages. We observed that genes specific for alveolar macrophages (AMs) were significantly enriched at baseline (−5 dpi) relative to 4 dpi, indicating a downregulation of this gene set after SARS-CoV-2 infection (Fig. 2f). Collectively, these bulk RNA-Seq data indicate a rapid and significant shift in the balance of macrophage populations in the lower airway following SARS-CoV-2 infection.

SARS-CoV-2 infection induces an influx of two subsets of infiltrating macrophages into the alveolar space

In our prior work in RMs, we demonstrated that cells of myeloid origin were the predominant subset responsible for production of inflammatory cytokines in the lower airway following SARS-CoV-2 infection²⁰. While our prior scRNA-Seq analyses determined the majority of cells in the BAL after infection to be of monocyte/macrophage origin, with relatively few neutrophils or granulocytes, the precise immunophenotypes of the myeloid cells driving inflammation in the lower airway have not been precisely delineated.

Cell classification based on cell-surface marker genes is typically problematic in scRNA-Seq data due to gene dropouts inherent to the technology. Accurate classification is further complicated in the rhesus model system, in which genomic references have incomplete annotation, and markers from other model species may not phenocopy. Several significant advances have been made recently elucidating the resident tissue macrophage subsets in the lung and their function during viral infection and inflammation^29,30,31,32. However, analysis of scRNA-Seq data from RM lung suspensions and BAL during steady state condition indicated that several key markers used to differentiate macrophages in the murine lung (e.g., Lyve1) were not expressed at levels sufficient to distinguish populations in the rhesus pulmonary myeloid populations (Supplementary Fig. 7). Therefore, we used two overlapping alternative strategies to accurately classify tissue macrophages and monocyte-derived/infiltrating macrophages in the RM airway after SARS-CoV-2 infection in our scRNA-Seq data. The first strategy was based on using existing lung scRNA-Seq data from uninfected RMs as a reference to map and annotate the BAL cells. We processed lung 10X data from three uninfected RMs (NCBI GEO: GSE149758)³³ through the Seurat pipeline³⁴ and reproduced the four reported macrophage/monocyte subsets: CD163+MRC1+, resembling alveolar macrophages; CD163+MRC1+TREM2+ macrophages, similar to infiltrating monocytes; CD163+MRC1−, similar to interstitial macrophages; and CD16+ non-classical monocytes (Supplementary Fig. 4a–d). We used Seurat to map BAL macrophages/monocytes from SARS-CoV-2 infected RMs and transfer annotations from the lung reference. The second strategy involved using bulk RNA-Seq on sorted AM and IM from the lungs of three uninfected RMs³⁵, according to the phenotype defined by Cai et al.³⁶, based on expression of CD206/MRC1 and CD163, to annotate cells using SingleR²⁸ Supplementary Fig. 4e, f). In total, 2069 genes were found to be differentially expressed between IMs and AMs (FDR < 0.05, fold-change > 2) (Supplementary Fig. 4e). Of note, CX3CR1 was upregulated in the IMs, consistent with both murine and human definitions of this subset (Supplementary Fig. 4f). APOBEC3A, an RNA-editing cytidine deaminase, was also upregulated in IMs along with PTGS2, a pro-inflammatory COX-2 cyclooxygenase enzyme, TIMP1, which enables migration of cells via the breakdown of connective tissue, VCAN, an immunosuppressive regulator, and PDE4B, which regulates expression of TNFα (Supplementary Fig. 4f). We annotated the lung macrophage/monocyte subsets using the bulk sorted AM and IM datasets and found that almost all of CD163+MRC1+ cluster and some CD163+MRC1+TREM2+ cells were annotated as AM and the remaining as IM (Supplementary Fig. 4g). Thus, benchmarking our lung scRNA-Seq based reference against rudimentary bulk transcriptomic signatures demonstrated their accuracy in resolving the AM phenotype from non-AM in steady state conditions.

We next analyzed changes in the myeloid populations within the BAL of RMs after SARS-CoV-2 infection by applying these signatures to two independent scRNA-Seq datasets from rhesus macaques infected intranasally and intratracheally with 1.1 × 10⁶ PFU of the USA-WA1/2020 strain of SARS-CoV-2: Cohort 1, comprised of a dataset of n = 3 (baseline and 4 dpi)²⁰ and Cohort 2, comprised of n = 5 (baseline) and n = 6 (4 dpi)²⁶. Using the lung/scRNA-Seq reference, we found that most of the BAL macrophages/monocytes belonged to the AM-like CD163+MRC1+ macrophage subset at −5 dpi along with some cells from the CD163+MRC1+TREM2+ macrophage subset (Fig. 3a, b and Supplementary Fig. 5a). At 4 dpi, there was an influx of both CD163+MRC1+TREM2+ macrophages and the IM-like CD163+MRC1− macrophages with few cells annotated as CD16+ non-classical monocytes. The expression of gene markers such as MARCO, FABP4 and CHIT1 further supported the cell subset annotations (Fig. 3c). We also observed a similar increase in APOBEC3A and decreases in the alveolar macrophage-associated genes MARCO and CHIT1 expression in BAL samples analyzed by bulk RNA-Seq, indicating a loss of CD163+MRC1+ cells (Fig. 3d). By analyzing the sc-RNA-Seq datasets, the percentage total of CD163+MRC1+ macrophages at baseline compared to 4 d.p.i. reduced from 93.3% to 55%, and 89.3% to 63.1% of all macrophages/monocytes in BAL in Cohorts 1 and 2, respectively (Fig. 3e, f). Estimates of cellular frequencies in pooled scRNA-Seq datasets can be driven by unbalanced cell counts from individual samples. To account for this potential bias, we examined the changes in myeloid populations by individual animals. We observed that in Cohort 1, the CD163+MRC+ cells decreased from mean 90.2% (sd = 5.2%) to mean 65.4% (sd = 32%), p = ns, and in Cohort 2, they decreased significantly from mean 92% (sd = 5.1%) to mean 65.7% (sd = 16.4%) p = 0.0087 (Fig. 3g, h). Conversely, we saw an overall increase in the percentage of CD163+MRC1+TREM2+ macrophages from 6.5 to 36.8% (Cohort 1) and 10.3 to 19.8% (Cohort 2) (Fig. 3e, f). At the individual level, we observed that levels of CD163+MRC1+TREM2+ cells increased from mean 9.7% (sd = 5.3%) to mean 28.6% (sd = 25.2%) in Cohort 1 (p = ns), and mean 7.7% (sd = 4.9%) to mean 20.0% (sd = 10.4%) (p = 0.03) (Fig. 3g, h). Additionally, we observed increases in the IM-like CD163+MRC1− macrophages: in the pooled scRNA-Seq data from 0.2 to 8% in Cohort 1 and 0.4 to 17% in Cohort 2 (Fig. 3e, f). Considering individual animals, the CD163+MRC1− cells increased from mean 0.15% (sd = 0.19%) to mean 5.9% (sd = 6.9%) in Cohort 1 and significantly from mean 0.28% (sd = 0.18%) to mean 14.2% (sd = 9.1%) (p = 0.004) (Fig. 3g, h). Thus, SARS-CoV-2 infection resulted in an influx of monocyte-derived and IM-like macrophages in BAL at 4 dpi.

a Projection of single-cell macrophages/monocytes from −5 dpi (green) and 4 dpi (magenta) 10X BAL samples obtained from three SARS-CoV-2 infected rhesus macaques (Cohort 1) onto the reference UMAP of lung macrophage/monocytes from uninfected rhesus macaques (NCBI GEO: GSE149758³³). b UMAP projections showing the predicted cell type annotations based on the uninfected lung reference split by time of sample collection (Cohort 1). c DotPlots showing the expression of marker genes for the different macrophage/monocyte subsets in SARS-CoV-2 infected BAL samples (Cohort 1). d Log2 fold-changes compared to −5 dpi for APOBEC3A, CHIT1 and MARCO in bulk BAL RNA-Seq data (Cohort 1). Significance was determined using DESeq2. The p values corrected by the default Benjamini and Hochberg method were used *adj p value <0.05, ***adj p value <0.001. The exact p values are included in Supplementary Data 2. Percentage of a given subset out of all macrophage/monocyte subsets at −5 dpi and 4 dpi from all three animals pooled (Cohort 1) (e) and at −7 dpi and 4 dpi from all six animals pooled (Cohort 2) (f). Percentage of a given subset out of all macrophage/monocyte subsets for each animal at −5 dpi and 4 dpi from Cohort 1 (n = 3) (g) and at −7 dpi (n = 5) and 4 dpi (n = 6) for Cohort 2 (h). The black lines indicate the median. p values for Cohort 2 (f) for 4 dpi vs. −7 dpi: CD163+MRC1+ = 0.009, CD163+MRC1+TREM2+ = 0.03, and CD163+MRC1− = 0.004. Contribution of each macrophage/monocyte subset toward the production of the pro-inflammatory genes and ISG—Cohort 1 pooled (i), Cohort 2 pooled (j), Cohort 1 individual (k) and Cohort 2 individual (l). p values for (l): * = 0.03. The percentage contribution was calculated by dividing the sum of normalized expression of a given gene in a macrophage/monocyte subset by the sum of the normalized expression of the gene in all macrophage/monocyte subsets. a–c, e, g, i, k Cohort 1: n = 3 for both −5 dpi and 4 dpi, d Cohort 1: n = 8 for 2 dpi and n = 4 for 4 dpi, f, h, j, l Cohort 2: n = 5 for −7 dpi and n = 6 for 4 dpi. The black lines indicate the median. Statistical analysis was performed using two-tailed Wilcoxon singed rank test for (g, k, l) and two-tailed Mann–Whitney U test for (h) in R v4.2.2. *p value <0.05, **p value <0.01. Source data (d–l) are provided as a Source Data file.

To further validate our cell classification and support the observation that it is the infiltrating cells that increase in numbers and predominantly produce inflammatory mediators, we used the second strategy of using gene expression of bulk sorted AM and IM cells to classify the BAL macrophages/monocytes. Using this definition, we confirmed that there is an increase in the percentage of non-AM population with a corresponding decrease in the AM population (Supplementary Fig. 5b, d). The non-AM population was also found to show higher expression of pro-inflammatory cytokines (Supplementary Fig. 5e).

Differential expression analysis of 4 dpi and −5 dpi BAL macrophages/monocytes showed that CHIT1, MARCO and MRC1 were among the top-ranking genes exhibiting downregulation in the BAL, while genes such as ADAMDEC1³⁷ and S100A8³⁸ that are associated with monocyte-derived macrophages were among the most upregulated (Supplementary Data 4). These data demonstrate that our observation of an influx of infiltrating macrophages into the BAL at 4 dpi was consistent across multiple definitions of this phenotype.

Infiltrating macrophages produce the majority of lower airway inflammatory cytokines during acute SARS-CoV-2 infection

Given our observation of the dynamics of pulmonary macrophages within the alveolar space during early SARS-CoV-2 infection, we characterized the transcriptional changes in each macrophage/monocyte population (Fig. 3), and also in conventional myeloid dendritic cells. Myeloid DCs were present at very low frequencies (<2%), and not more than 70 total cells were detected to be harboring mRNA from TNF, IL6 or IL10 after infection. IL1B expression was slightly higher, but accounted for <2.5% of IL1B expressing cells in the BAL at 4 dpi (Supplementary Fig. 6a–e). Several chemokines (CCL4L1, CCL3, CXCL3, CCL2), multiple ISGs, NFKB1A, S100A8, and GZMB were among the most upregulated genes at 4 dpi in BAL populations (Supplementary Data 4). Elevated expression of multiple inflammatory genes, including IL6, TNF, IL10, and IL1B, were observed in the CD163+MRC1+TREM2 Mac and CD163+MRC1− subsets in both Cohort 1 and 2 (Fig. 3i–l and Supplementary Fig. 7) after infection. The infiltrating macrophages were also observed to upregulate multiple chemokines, including those specific for recruiting neutrophils (CXCL3, CXCL8), macrophages (CCL2, CCL3, CCL5, CCL4L1), and activated T cells (CXCL10) as well as multiple ISGs (Supplementary Figs. 7 and 8). When we examined CD163+MRC1+ macrophages, many of the same inflammatory cytokines and gene sets seen in the infiltrating macrophages were elevated at 4 dpi, albeit at much lower magnitude (Supplementary Figs. 7 and 8). Having observed a significantly higher average expression of inflammatory cytokines in infiltrating macrophages compared to CD163+MRC1+ macrophages, we compared the fractions of sequencing reads detected from each of the subsets to assess the overall contribution to inflammatory cytokine production (Fig. 3i). In Cohort 1, at 4 dpi, we observed that the CD163+MRC1+TREM2+ macrophages accounted for 55% of IL6, 57% of TNF, and 86% of IL10 expression while the CD163+MRC1− macrophages accounted for 20% of IL6, 21% of TNF, and 6% of IL10 expression (Fig. 3i). In Cohort 2, we also observed that the infiltrating macrophages contributed more to the expression of most inflammatory cytokines than alveolar macrophages: the expression in CD163+MRC1+TREM2+ and CD163+MRC1− cells was 39.2% and 47% for IL10, 17.4% and 74.9% for IL6, and 22.4% and 46.6% for TNF, respectively (Fig. 3j). To account for potential bias in cell counts in our pooled data, we also examined the contributions of cytokines to the BAL expression in individual animals (Fig. 3k, l). In Cohort 1, the overall trend of higher contribution to inflammatory expression in the infiltrating macrophage populations seen in the pooled data was only observed for IL10 and CCL4L1, largely due to imbalances in cell counts and low statistical power (Fig. 3k). However, in Cohort 2, we observed consistently elevated levels in the infiltrating CD163+MRC1+TREM2+ and CD163+MRC1− macrophage populations (Fig. 3l). For IL10, the mean ± sd expression in CD163+MRC1+ was 17.3 ± 12.4%, compared to 40.4 ± 12.1% in CD163+MRC1+TREM2+ cells (p = 0.03) and 42.3 ± 21.2% for CD163+MRC1− (ns) (Fig. 3l). For IL6, the mean ± sd expression in CD163+MRC1+ was 9.4 ± 12.2%, compared to 12.2 ± 24.5% in CD163+MRC1+TREM2+ cells (ns) and 78.4 ± 28.8% for CD163+MRC1− (p = 0.065). Additionally, while our observation for IL6 trended to significance, we found a wide variability in the percentages in CD163+MRC1− cells; to address this, we analyzed another set of data²⁶ in which we obtained scRNA-Seq expression of BAL macrophage after 2 dpi of SARS-CoV-2 infection—these data also trended to much higher expression of IL6 in the CD163+MRC1− cells (67.1 ± 32.4%) compared to CD163+MRC1+ macrophages (27.7 ± 29.5%)(Supplementary Fig. 9). For CCL4L1, CD163+MRC1+ contributed 9.6 ± 6.3% of expression, in comparison to 23.7 ± 13.5% in CD163+MRC1+TREM2+ cells (p = 0.03) and 66.7 ± 18.2% in CD163+MRC1− cells (p = 0.03). The contribution of different subsets towards CXCL10 expression was 19.2 ± 16.7% from CD163+MRC1+ compared to 15.3 ± 9.3% from the CD163+MRC1+TREM2+ and 65.4 ± 15.2% from the CD163+MRC1− populations. Overall, these data indicate that the infiltrating macrophage populations are responsible for the majority of lower airway inflammatory cytokine production during acute SARS-CoV-2 infection.

To validate the increase in infiltrating myeloid populations, we quantitated by flow cytometry the frequency of CCR2+ myeloid populations in peripheral blood and BAL from an additional six SARS-CoV-2-infected rhesus macaques (Supplementary Fig. 10). CCR2 has been demonstrated to regulate monocyte infiltration into the lung parenchyma of SARS-CoV-2 infected mice³⁹ and its expression is upregulated in the BAL of infiltrating macrophages in NHPs infected with influenza virus³⁵. Consistent with our observation of elevated infiltrating myeloid cells by scRNA-Seq, we found that there was a concomitant increase in the frequency of CCR2+ CD14−CD16+ and CD14+CD16+ monocytes at 2 dpi. These results further support the infiltration of inflammatory monocytes in BAL after SARS-CoV-2 infection.

Identification of pro-inflammatory subsets in human SARS-CoV-2 infection corresponding to NHP immunophenotypes

To translate our findings in the NHP model to human SARS-CoV-2 infection, we used a similar bioinformatic approach to that employed to define rhesus myeloid. We used macrophages/monocytes from publically available scRNA-Seq dataset of lungs from six healthy human donors (GEO: GSE135893⁴⁰) and classified these based on a recent classification into FABP4hi, SPP1hi, FCN1hi and proliferating macrophages⁴¹ (Fig. 4a). When the canonical marker genes were compared between the healthy lung macrophages/monocytes of human and rhesus macaque, we found that there were comparable populations between the two (Fig. 4a–c). Namely the CD163+MRC1+ rhesus subset was highly similar to the FABP4hi human subset; the CD163+MRC1+TREM2+ rhesus subset was congruent with the SPP1hi human subset; and the CD163+MRC1− macrophages and CD16+ monocytes rhesus were transcriptionally similar to the FCN1hi human subset. Next, we combined data from all macrophage/monocyte cells from the six healthy human samples with the three healthy rhesus samples and applied reference-based integration in Seurat using the human samples as reference (Supplementary Fig. 11a, b). We looked at the distribution of different human and rhesus cell types in each cluster and found that the earlier observations regarding the similarity of subsets based on canonical markers was further supported by the global gene expression of these cells (Supplementary Fig. 11c). Finally, to test the robustness of our cellular classifications to identify macrophage subsets between species accurately, we generated gene signatures for each subset/species combination and tested for enrichment in the opposite species; for all comparisons, a signature scored highest with its corresponding opposite subset (Supplementary Fig. 11d).

a UMAP of macrophage subsets in lungs from six healthy human donors (GEO: GSE135893). b UMAP of macrophage subsets in lungs from three healthy rhesus macaques (GEO: GSE149758). c DotPlots showing the expression of marker genes. The color gradient represents the level of expression and the size of the dot represents the percentage of cells expressing a given gene. d UMAP of BAL samples from human donors that are healthy (n = 3), or suffering from moderate (n = 3) or severe (n = 6) COVID-19 (GEO: GSE145926) mapped to the healthy lung reference using the Seurat MapQuery function. e Percent of predicted cell types out of all macrophage/monocytes in each human BAL sample. The black bar indicates the median. Statistical analysis was performed using pairwise two-tailed Mann–Whitney U test in R v4.2.2. p value * = 0.02. Contribution of each predicted macrophage/monocyte subsets in human BAL toward the production of the pro-inflammatory genes and ISG—pooled (f) and individual (g). The percentage contribution was calculated by dividing the sum of normalized expression of a given gene in a macrophage/monocyte subset by the sum of the normalized expression of the gene in all macrophage/monocyte subsets. The black bars represent the median. Statistical analysis was performed using two-tailed Wilcoxon signed rank test in R v4.2.2. *p value 0.03 for all except FABP4 vs. Proliferating for IL6, IL1B, TNF, CXCL3, CXCL8 and SPP1hi vs. Proliferating for CCR2: p value = 0.04. Source data (e–g) are provided as a Source Data file.

Using the healthy human dataset as reference, we classified macrophages/monocytes from the myeloid cluster (Supplementary Fig. 11e, f) in a publically available scRNA-Seq dataset of human BAL samples from healthy donors or subjects with moderate or severe COVID-19 infection⁸ (Fig. 4d and Supplementary Fig. 11g). As reported in the original study⁸, we found that there was a significant increase in the SPP1hi and FCN1hi subsets with COVID-19 infection while the FABP4hi population that is largely representative of the resident alveolar macrophages was found to be significantly reduced in both moderate and severe COVID-19 infection (Fig. 4d, e). In addition, we also looked at the contribution of these different populations toward inflammation and found that the non-resident SPP1hi and FCN1hi subsets were largely responsible for the expression of pro-inflammatory mediators in patients with severe COVID-19 (Fig. 4f and Supplementary Fig. 12). Respectively, (FABP4hi, SPP1hi, and FCN1hi) the contribution of each population to overall expression was: IL10 (1.3%, 56.7%, 42%), IL1B (6.6%, 37%, 56%), IL6 (2.4%, 40.8%, 56.7%), TNF (1.6%, 44.1%, 54.1%), CXCL10 (1.8%, 36.5%, 61.7%) and CXCL8 (2.6%, 48.1%, 49.1%) (Fig. 4f). Of note, ISG expression (IFI27, ISG15, ISG20, MX2) was significantly upregulated in the SPP1hi, and FCN1hi populations relative to the FABP4hi in severe disease but not in moderate cases. Lastly, when we examined the contribution of cytokine expression in by individual patients in these dataset, we noted that for severe COVID-19, we observed that the infiltrating populations had significantly higher expression of IL10 (mean ± sd SPP1hi 50.9 ± 8.7%, p = 0.03; FCN1hi 46.9 ± 8.3%, p = 0.03), IL1B (mean ± sd SPP1hi 31.7 ± 14.9%, p = 0.03; FCN1hi 63.3 ± 20.1%, p = 0.03), TNF (mean ± sd SPP1hi 43.4 ± 9.3%, p = 0.03; FCN1hi 53.7 ± 13.7%, p = 0.03), CXCL10 (mean ± sd SPP1hi 27.3 ± 9.5%, p = 0.03; FCN1hi 69.8 ± 11.1%, p = 0.03), CXCL3 (mean ± sd SPP1hi 55.9 ± 9.6%, p = 0.03; FCN1hi 37.6 ± 13%, p = 0.06) and CXCL8 (mean ± sd SPP1hi 42.7 ± 14%, p = 0.03; FCN1hi 53.4 ± 17.4%, p = 0.03) compared to the FABP4hi population (mean ± sd for IL10 2.2 ± 2.8%, IL1B 4.7 ± 8.2%, TNF 2.7 ± 4.5%, CXCL10 2.6 ± 2.9%, CXCL3 6.4 ± 9.5% and CXCL8 3.8 ± 5.9%) (Fig. 4g). Collectively, these analyses demonstrate that the myeloid subsets defined transcriptionally in RMs have analogous populations in the human lung, and an overall concordance in their expansion and contribution to SARS-CoV-2 induced inflammation in the lower airway.

Baricitinib treatment prevents the influx of inflammatory IM into the lower airway

Baricitinib is a JAK1/2 inhibitor approved for the treatment of active rheumatoid arthritis that was recently approved by FDA for COVID-19 treatment in certain hospitalized adults⁴², and reported to reduce mortality when administered as monotherapy⁴³ or in combination with remdesivir²⁵. In our earlier study, using data from the Cohort 1 and baricitinib-treated animals (Fig. 1a), we found that baricitinib was able to suppress the expression of pro-inflammatory cytokines in BAL of RMs infected with SARS-CoV-2²⁰. Here, we extended this study to further characterize the impact of baricitinib on the myeloid populations in the airway from five RMs before infection (−5 dpi) and at 4 dpi, with three RMs that remained untreated and two that received baricitinib. We found that 2 days of baricitinib administration virtually abrogated the influx of infiltrating macrophages into the alveolar space at 4 dpi, as we did not detect any increase in the CD163+MRC1− or CD163+MRC1+TREM2+ populations in baricitinib-treated animals (Fig. 5a–d). This observation was consistent using classifications of macrophages either based on mapping to 10X lung reference or using bulk sorted AM/IM cells (Fig. 5a–d and Supplementary Fig. 5a–d). In addition to preventing the influx of infiltrating macrophages, baricitinib treatment also resulted in significantly lower expression of inflammatory cytokines and chemokines, but the ISG expression remained comparable to untreated animals (Fig. 5e–g and Supplementary Fig. 5e). In summary, these data further elucidate the mechanism of action by which baricitinib treatment abrogates airway inflammation in SARS-CoV-2 infection²⁰, by demonstrating its ability to block infiltration of discrete pro-inflammatory macrophage populations into the alveolar compartment. Similar to our observations using flow cytometry, there was an increase in the abundance of pDC at 4 dpi in the BAL detected by scRNA-Seq data, however this increase was abrogated in baricitinib-treated animals (Fig. 5h, i).

a Projection of macrophages/monocytes from −5 dpi and 4 dpi 10X BAL samples from three untreated and two baricitinib-treated rhesus macaques on the reference UMAP of uninfected lung macrophages/monocytes (NCBI GEO: GSE149758). b UMAP split by treatment and timepoint showing predicted cell annotations based on mapping to the reference lung macrophages/monocytes. Percentage of a given macrophage/monocyte subset of all the macrophages/monocytes in the BAL samples—pooled (c) and individual (d). Violin plots showing expression of pro-inflammatory cytokines (e), chemokines (f) and ISG (g) in the different macrophage/monocyte subsets in BAL 10X samples from baricitinib-treated and untreated samples. h Absolute number of pDC in BAL samples from scRNA-Seq and i Percentage of pDC out of all cells in the BAL samples from scRNA-Seq. Source data (c, d, h, i) are provided as a Source Data file.