CD8⁺ lymphocytes are critical for early control of tuberculosis in macaques

Prior to the 2020 COVID-19 epidemic, tuberculosis (TB) was the leading infectious cause of death in the world. With roughly 1.6 million deaths and over 10.6 million new cases of active TB in 2021, TB remains a global health crisis (WHO, 2022). It is widely agreed that there is a critical need for an effective preventive vaccine against Mycobacterium tuberculosis (Mtb) infection and/or TB disease. However, the development of such a vaccine has been hindered by our limited understanding of the complex host response to this pathogen.

Although CD4 T cells have been the focus of several TB vaccines, they are unlikely to be sufficient for protection against infection or disease (Lewinsohn et al., 2017a; McShane, 2019). Recent studies in macaque models have demonstrated robust protection against infection or disease through intravenous (IV) or mucosal delivery of the TB vaccine strain Bacille Calmette-Guérin (BCG) or by vaccination with a CMV viral vector expressing several Mtb antigens (Darrah et al., 2020; Dijkman et al., 2019; Hansen et al., 2018). IV BCG induces strong CD4 and CD8 T cell responses as well as innate CD8 cell responses, and CMV induces robust CD8 T cell responses, supporting a role for CD8⁺ lymphocytes in protection against TB (Darrah et al., 2020; Hansen et al., 2018). Chen et al. (2009) showed that CD8 depletion in intradermal BCG-vaccinated or unvaccinated rhesus macaques reduced protection against high-dose Mtb challenge. In addition, Mtb infection induces both CD4 and CD8 T cell responses in humans (Caccamo et al., 2015; Lancioni et al., 2019; Lewinsohn et al., 2017b) and macaques (Gideon et al., 2015; Lin and Flynn, 2015; Lin et al., 2012), and we demonstrated previously that primary infection provides robust protection against reinfection in macaques (Cadena et al., 2018). In humans, CD8 T cells expressing several cytotoxic markers have been linked to the control of mycobacterial infection, with granulysin implicated in the direct killing of Mtb (Balin et al., 2018; Busch et al., 2016; Stenger et al., 1998). Our recent single-cell RNA sequencing (scRNAseq) analysis of macaque granulomas revealed several innate and adaptive cytotoxic subsets that were associated with granuloma-level bacterial killing (Gideon et al., 2022). These CD8⁺ lymphocyte subsets expressed a variety of cytotoxic molecules including several different granzymes, granulysin, and perforin (Gideon et al., 2022).

Lymphocytes expressing CD8 molecules can be adaptive CD8⁺ T cells or innate-like lymphocytes, including NK cells and mucosal-associated invariant T cells (MAITs), γδ T cells, and natural killer T (NKT) cells, which recognize non-peptide antigens in the context of MHC-like molecules. Enhanced γδ T cell responses were associated with a modest improvement in TB disease in macaques (Shen et al., 2019) and were enriched in diseased human lung (Ogongo et al., 2020), while NK cells in human blood were associated with controlled (latent) Mtb infection (Roy Chowdhury et al., 2018). These innate CD8⁺ cells have been suggested for TB vaccine targeting due to their donor-unrestricted nature (Joosten et al., 2019; Ruibal et al., 2021).

Despite mounting evidence that CD8⁺ lymphocytes are associated with protection, a causal role for innate or adaptive CD8⁺ cells in early control of Mtb infection has not been established. Here, we sought to dissect the importance and functions of CD8-expressing immune subsets in early Mtb infection. Using cynomolgus macaques, we used an in vivo depletion strategy to define the roles of innate and adaptive CD8⁺ cell types in early control of Mtb infection. CD8 is expressed as a dimer and, depending on cell type, can be expressed as a CD8αα homodimer or a heterodimer of CD8α and CD8β chains. In humans and non-human primates (NHPs), many of the innate-like subsets express CD8αα, including γδ T cells, NK T cells, MAITs, and NK cells. However, adaptive, conventional CD8 T cells (TCRαβ) predominantly express the CD8αβ heterodimer, allowing for distinction and differential depletion of immune cells by CD8 chain expression.

Cynomolgus macaques were depleted with anti-CD8α or anti-CD8β monoclonal antibody (Ab) prior to and for 6 wk following infection with the intention that anti-CD8α Ab would deplete all CD8 expressing cell types (innate-like and adaptive), while anti-CD8β Ab would selectively deplete conventional CD8 T cell populations, leaving innate CD8 compartments intact (Fig. S1 A). We demonstrate that depletion of all (innate and adaptive) CD8⁺ lymphocytes prior to infection enabled Mtb better establishment of infection, increased dissemination, and worsened TB disease; depletion of adaptive CD8⁺ lymphocytes worsened bacterial control in the lymph nodes. Through flow cytometry and scRNAseq, we found that in the absence of CD8⁺ cells, CD4 and γδ T cells adopt enhanced but incomplete cytotoxic signatures, suggesting environmental drivers of cytotoxic functions which when expressed in CD4 and γδ T cells apparently do not recapitulate the full cytotoxic potential of CD8⁺ lymphocytes nor substitute for them in control of infection. Thus, this study demonstrates that CD8⁺ lymphocytes play a critical protective role in early Mtb infection.

To assess the role of innate or adaptive CD8⁺ lymphocytes during early Mtb infection, cynomolgus macaques were treated with either a control IgG, anti-CD8α, or anti-CD8β depletion antibody prior to and during 6 wk of Mtb infection (Fig. S1 A). This is a very early time point in Mtb infection and allows us to assess the role of CD8 T cells in the earliest events of infection—establishment, granuloma formation, and early dissemination; however, most macaques have a high bacterial burden at this time point (Gideon et al., 2022; Lin et al., 2014), making it difficult to detect substantial differences in bacterial burden upon depletion of immune effectors. Efficient depletion of targeted CD8 populations in blood was observed prior to and during Mtb infection as both CD8αα- (which include γδ T cells of which ∼50% express CD8α, MAITs, NK, and NKT cells) and CD8αβ-expressing T cells fell below ∼10 cells/μl (Fig. S1 B). Similarly, in bronchoalveolar lavage (BAL), anti-CD8α Ab depleted both CD8αα⁺ cells and CD8αβ T cells, while anti-CD8β Ab only depleted CD8αβ T cells with the exception of one animal (Fig. S1 C). Thus, we successfully depleted CD8 populations based on α- or β-expression in the periphery and airways.

Serial positron emission tomography (PET) with computed tomography (CT) (PET CT) scans of control and both groups of CD8-depleted macaques indicated similar total lung inflammation (measured by lung 2-deoxy-2-¹⁸F-D-deoxyglucose [FDG] activity) and numbers of granulomas at 3 wk after infection (Fig. S1 D and Table S1). However, by 6 wk after infection (pre-necropsy scan), CD8α-depleted animals had significantly higher total lung FDG activity and granuloma numbers compared with IgG controls (Fig. 1 A). Total lung FDG activity and number of granulomas varied substantially in the CD8β depletion group but were not significantly different from controls at 6 wk (Fig. 1 A).

Previously, we reported that peak bacterial burden within granulomas occurs at 4–6 wk after infection, with bacterial killing occurring by 10–12 wk after infection (Gideon et al., 2022; Lin et al., 2014). Thus, the total thoracic bacterial burden is generally quite high at early time points before cynomolgus macaques begin to control infection (Grant et al., 2022). Nonetheless, CD8α-depleted animals had significantly higher overall thoracic bacterial burden (which includes all lung CFU and thoracic lymph nodes; P = 0.0378) and higher thoracic lymph node CFU (P = 0.0224) with a trend toward higher lung CFU (P = 0.0717) compared with IgG controls (Fig. 1 B and Table S1). There were significantly more CFU⁺ lung lobes and CFU⁺ lymph nodes and significantly more overall pathology and extrapulmonary disease in the CD8α-depleted group compared with controls, indicating extensive disease and dissemination (Fig. 1, C and D). CD8β-depleted animals had a significantly higher bacterial burden in the thoracic lymph nodes compared with IgG controls (P = 0.0441; Fig. 1 B) and a trend toward more CFU⁺ lymph nodes (P = 0.0694; Fig. 1 C). These data indicate that depletion of either all (innate and adaptive) or only adaptive CD8⁺ lymphocytes leads to increased Mtb dissemination and/or increased bacterial replication in thoracic lymph nodes and support that CD8-expressing cells play a critical role in controlling early Mtb infection.

Using barcoded Mtb technology, where each bacillus in the inoculum has an individual DNA sequence tag (Martin et al., 2017), we previously showed that each granuloma is established by a single Mtb bacillus (Lin et al., 2014) and demonstrated tracking of dissemination in infected macaques (Martin et al., 2017). The macaques in the current study were infected with barcoded Mtb strain Erdman, which allows determination of the number of unique barcodes in an animal as a measure of how many individual bacteria establish infection and the extent of barcodes shared across tissues as a measure of dissemination. Mtb barcodes were sequenced from every CFU⁺ tissue sample including each individual granuloma and all other tissues categorized by anatomical site (lung, lymph nodes, and extrapulmonary). Timing of granuloma formation (3 versus 6 wk) was determined by serial PET CT scans. Of note, the sequencing of samples from one CD8β-depleted animal did not pass QC and is not included in this analysis. Circos plots depicting the complexity of bacterial barcode dissemination dynamics in each animal by depletion group showed shared barcodes among granulomas, lung lobes, lymph nodes, and extrapulmonary sites (Fig. 2, A and B) (Krzywinski et al., 2009). All animals in the CD8α depletion group showed substantial dissemination of Mtb from lung granulomas to other sites in the lung, lymph nodes, and extrapulmonary organs. Few extrapulmonary or lymph node sites have barcodes that are also not found in the lung, indicating a pulmonary origin for the bacilli. The plots qualitatively demonstrate differences in the extent of infection and dissemination among groups, with CD8α-depleted animals having notably more sharing of barcodes across tissues and CD8β-depleted animals showing an intermediate dissemination pattern.

These qualitative analyses were supported by quantification of barcode numbers and spread among the IgG control and CD8 depletion groups. The number of unique barcodes recovered per animal was significantly higher in CD8α-depleted macaques compared with the IgG-treated control macaques (Fig. 2 C), indicating that more bacilli established infection in the context of CD8α depletion. CD8α-depleted animals had a higher frequency of disseminated barcodes (i.e., percent of unique barcodes shared at two or more sites) (Fig. 2 D) and significantly increased percentage of tissues sharing a barcode (Fig. 2 E), as well as increased sharing of lung barcodes with both lymph nodes and extrapulmonary tissues (Fig. 2, F and G). In most of these measures, CD8β depletion resulted in an intermediate phenotype between CD8α-depleted and IgG control animals; although these were not significant, there were only four monkeys in this group for these analyses. These data indicate that CD8α⁺ lymphocytes limit the establishment and early dissemination of Mtb infection.

In light of the Mtb barcode data indicating that CD8⁺ lymphocytes limit the establishment of infection, we sought to understand what cell types and states are present in the airways of immunologically naïve animals to delineate which features may contribute to an infection bottleneck. Flow cytometry on BAL from uninfected and non-depleted cynomolgus macaques revealed that greatest proportion of cells were myeloid and that both CD4 and CD8αβ T cells were present at higher frequencies than CD8α⁺ innate T or NK cells (Fig. S2, A and B), with ∼20–40% of the CD8αα and CD8αβ lymphocytes expressing some cytotoxic proteins and almost no cytokine production from either the CD8⁺ or CD4⁺ lymphocytes (Fig. S2 C). Stimulation with mycobacterial antigens did not enhance cytokine or cytotoxic effector molecule expression, confirming that these were not Mtb antigen–specific T cells (Fig. S2 D). Thus, even in uninfected macaques, airways have both CD8αα and CD8αβ lymphocytes, which could impose a modest bottleneck for infection of the lungs.

At necropsy (6 wk after infection), PET CT scan–matched granulomas were harvested individually from each animal as previously described (Cadena et al., 2018; Martin et al., 2017). Using spectral flow cytometry, we characterized lymphocytes in individual granulomas (Fig. S3, A–C; and Table S2). In total, 188 granulomas were analyzed (five monkeys per group; granuloma numbers: IgG n = 57; αCD8α n = 64; αCD8β n = 67) of which individual and median values of lymphocyte populations in granulomas are represented in Fig. 3 A. In IgG control granulomas, the primary lymphocyte populations were CD4 T cells (30.86%) and CD8αβ T cells (22.08%). We identified a variety of other lymphocytes, including CD8αα T cells (2.79%), γδ T cells (comprised of both CD8α⁺ [3.32%], and CD8α⁻ [3.72%] subsets), CD4⁺CD8α⁺ (double positive) T cells (8.17%), CD4⁻CD8α⁻ (double negative) T cells (2.09%), B cells (18.58%), and innate lymphocytes (CD3⁻CD8α⁺ 4.19%, CD3⁻CD8α⁻ 4.20%). We considered CD3⁻CD8α⁺ to be NK cells, with the understanding that some NK populations fall within the ILC definitions (CD3⁻CD8α) (Fig. S3 B). Thus, control granulomas contain diverse lymphocytes at 6 wk after infection.

CD8α depletion resulted in a robust reduction in the number of lymphocyte populations represented (Fig. 3 A). Of the 10 subsets defined in IgG granulomas, five cell subsets were reduced to <0.2% following CD8α depletion and the granulomas were dominated by five lymphocyte subsets: CD4 T cells, CD8α^- γδ T cells, B cells, CD4⁻CD8α⁻ (double negative) T cells, and CD3⁻CD8α⁻ innate lymphocytes. In the CD8α depletion treatment group, the frequencies of CD4 (P = 0.0324) and CD8α⁻ γδ T cells (P = 0.0392) were expanded, while the frequencies of CD3⁻CD8α⁺ innate lymphocytes (P = 0.0806) and CD8αβ⁺ T cells (P = 0.0392) were reduced compared with IgG controls (Fig. 3 A). In the CD8β depletion treatment group, frequencies of lymphocyte cell types remained similar to IgG controls, with the exception of CD8αβ⁺ T cells, which were significantly reduced (P = 0.0094) compared with IgG controls (Fig. 3 A).

Despite the effective depletion of CD8⁺ lymphocytes and changes in the frequencies of lymphocyte subsets, it was somewhat surprising that total lymphocyte numbers per granuloma were similar between IgG control and both depletion groups (Fig. S3 D). Total numbers of CD8αβ T cells in granulomas were significantly reduced in both anti-CD8α and anti-CD8β antibody-treated animals compared with controls, while most other cell subsets were not significantly different. This suggests that other cell types, such as modest increases in CD8α⁻ γδ T cells and double negative T cells, are taking the place of the depleted cells in granulomas (Fig. S3 D).

To assess the functions of each lymphocyte subset, we performed flow cytometry and intracellular staining on granulomas in the IgG-treated (control) animals for cytokines (IFN-γ, TNF, IL-2, IL-17, and IL-10) and cytotoxic effectors (granzyme B, granulysin, and CD107a) (Fig. 3 B and Fig. S3 E). For a simplified comparison of effector functions, we combined all cytokine or all cytotoxic markers using an AND/OR Boolean gate in Flow Jo for each cell type in control granulomas (Fig. 3 C). As in our previously published data (Gideon et al., 2015), a large percentage of T cells in granulomas did not express cytokines. A subset of CD4 T cells expressed cytokines or cytotoxic molecules. All other cell types had significantly higher cytotoxic marker expression compared with cytokine expression. Of note, CD4⁺CD8α⁺ T cells were similar in cytokine expression to CD4⁺CD8α⁻ T cells but had slightly higher frequencies of granzyme B than CD4⁺CD8α⁻ T cells (Fig. 3, B and C). Thus, at 6 wk after Mtb infection, most lymphocytes in granulomas have a predominantly cytotoxic immune signature.

To assess the multifunctionality of cytotoxic T cell subsets, limited combination gating was performed on the dominant cytokines (IFN-γ and TNF) and cytotoxic (granzyme B and granulysin) effectors observed in granulomas. We focused on CD4 T cells and γδ T cells as these are the dominant lymphocyte populations after anti-CD8α antibody treatment (Fig. 3 A). In control granulomas at 6 wk after infection, CD8αβ T cells were predominantly granzyme B single expressors with a subset co-expressing granzyme B and granulysin (Fig. 4 A), while the majority of CD4 T cells (∼80%) did not express any of the effector molecules assessed here. γδ T cells in granulomas from control animals had a range of expression of IFN-γ, TNF, granzyme B, and/or granulysin with variability across granulomas and animals, although most γδ T cells in granulomas had relatively low expression of these effector molecules. In CD8α-depleted granulomas, there was a significant increase in the frequency of CD4 and γδ T cells expressing granzyme B (Fig. 4, A and B) and a trend towards increased frequencies of granzyme B⁺ γδ T cells in the CD8β-depleted granulomas (P = 0.0954) compared with IgG control granulomas. Granulysin expression was variable but not significantly different across groups (Fig. 4, A and B). γδ T cells also had increased coexpression of granzyme B and granulysin in many CD8α- and CD8β-depleted granulomas (Fig. 4 A). There were no significant differences in frequencies of cytokine-producing T cells in granulomas from CD8-depleted animals compared with IgG controls (Fig. 4 C). These data suggest that depletion of CD8⁺ lymphocytes may reprogram CD4 T cells and γδ T cells in granulomas to adopt cytotoxic phenotypes.

The flow cytometry data suggested, surprisingly, that CD8 depletion leads to cytotoxic molecule expression from other T cells. To interrogate the lymphocyte subsets in depleted granulomas in more detail, scRNAseq was performed on a subset of granulomas from these macaques (41 granulomas from 4 anti-CD8α [n = 17 granulomas], 3 anti-CD8β [n = 12 granulomas], and 3 IgG [n = 12 granulomas] NHPs) using the Seq-Well S³ platform (Gideon et al., 2022; Gierahn et al., 2017; Hughes et al., 2020) (Table S3).

We performed iterative clustering of 42,277 high-quality single cells and identified 13 major cell-type clusters: macrophages/monocytes (example marker genes: LYZ, IFI30), neutrophils (FCGR2B, CSF2RB, LMK2), classical dendritic cell (cDCs) (CPVL, PEA15), plasmacytoid dendritic cells (pDCs) (PLAC8, SULF2, CCDC50), T/NK cells (ETS1, CD3D, IFITM1), B cells (MS4A1, BANK1, CD19), plasma cells (MXB1, JCHAIN, TXNDC5), mast cells (CPA3, CLU, MS4A2), type 1 pneumocytes (AGER, CLDN18, EMP2), secretory epithelial cells (PIGR, WFDC2, SCGB3A1), type 2 pneumocytes (SFTPB, SFTPD, NPC2), endothelial cells (PTRF, SPTN1, EPAS1), and fibroblasts (DCN, COL1A2, COL3A1) (Fig. 5, A and B; Fig. S4 A; and Tables S4 and S5). Regardless of the depletion group, the dominant cell type clusters in each granuloma were macrophages/monocytes and T/NK cells (Fig. 5, C and D; and Tables S4 and S6).

To understand the impact of depletion on the composition of granulomas at the level of these general cell types, we compared the proportional abundance of clusters expected to be directly depleted under this perturbation, i.e., those containing cells that have CD8A and/or CD8B transcript expression (Fig. S4 B). No general cell type cluster was comprised entirely of cells expressing either of these genes (Fig. S4 B), so it was expected that, at least, the CD8α depletion would reduce the proportional abundance of CD8A-expressing T/NK cells, cDCs, and pDCs, and CD8β depletion would reduce the proportional abundance of T/NK cells. Neither expression of CD8A or CD8B in CD8α-depleted granulomas nor CD8B in CD8β-depleted granulomas was observed (Fig. S4 B).

However, compared with IgG control, lower proportions of several non-CD8A–expressing clusters were observed in granulomas from CD8α-depleted animals including mast cells (P = 0.004), plasma cells (P = 0.005), and B cells (P = 0.009), as well as larger proportions of secretory epithelial cells (P = 0.01) and type 1 pneumocytes (P = 0.02) (Fig. 5 D and Table S6), indicating that this depletion may have an impact on the wider granuloma microenvironment beyond the direct targets of depletion. None of these non-CD8A-expressing clusters showed decreased abundances in CD8β-depleted samples. The functional significance of these changes remains to be determined.

The expectation was that CD8α depletion would reduce the abundance of cell clusters that had CD8A gene expression in control granulomas. Of the expected target clusters of CD8α depletion, the pDC cluster (P = 0.0019) was the only one whose proportional abundance decreased in CD8α- and CD8β-depleted granulomas (Fig. 5 D and Table S6). Although subsets of cDCs cells can also express CD8A (i.e., cDC1 cells), cDC subclustering indicated no change in the frequencies of either cDC1s (CLEC9A) or cDC2s (CD1C) between treatment groups; intriguingly, despite a shared gene expression signature, cDC1s from the IgG and depletion conditions differed in their expression of CD8A, suggesting cDC functional plasticity in response to depletion (Fig. 5, J–M) (Villani et al., 2017). Meanwhile, consistent with the flow cytometry analysis, the scRNAseq data did not reveal differences in the proportion of the T/NK cluster, which had the largest number of cells expressing CD8A and/or CD8B in control granulomas, after either CD8α or CD8β depletion (Fig. 5 D and Fig. S4 B). This suggested that other T/NK cell subset types (i.e., those not expressing CD8A or CD8B) were increased in proportion in the T/NK cluster under depletion.

To gain further insight into the presumed increase of non-CD8A/CD8B expressing cells in the large T/NK cluster which led to its stable abundance after depletion, the T/NK cells identified in all experimental groups were subclustered (Fig. 6 A). In flow cytometry analyses, T/NK cell subsets were defined by their phenotypic expression of surface antigens (CD3, CD4, CD8α, CD8β, and γδ TCR). In contrast, unsupervised clustering of scRNAseq T/NK cells from all granulomas resulted in cells grouped by functional gene expression, as previously published (Gideon et al., 2022). Here, T/NK cells were first classified into seven major subclusters: cytotoxic T/NK (GZMK, GZMH, GZMB, GNLY, or PRF1), IFN-responding (CD8A, CD8B, IFI6, and MX1), non-cytotoxic CD8 (CD8A, CD8B, LAT2, TCF7, lacking cytotoxicity, and cytokine markers), proliferating T cells, stem-like (TCF7, S1PR1, IL6R), T1–T17 (IL7R, RORA, BHLHE40, NFKBIA, FURIN, IFNG, TNF), and Tregs (TIGIT, FOXP3) (Fig. 6, A and B; Fig. S4, C–F; and Tables S4 and S7). The dominant subclusters identified within 6-wk granulomas included cytotoxic T/NK (mean and standard deviation: IgG: 32.6 ± 9%, αCD8α: 33.8 ± 15%, αCD8β: 28.8 ± 9%), proliferating T cells (IgG: 17.7 ± 14%, αCD8α: 36 ± 20%, αCD8β:, 22.1 ± 20%), and T1T17 (IgG: 23.8 ± 9%, αCD8α: 20 ± 9%, αCD8β: 31.5 ± 16.6%) (Fig. 6 C and Table S4). However, neither CD8α nor CD8β depletion resulted in statistically significant changes to the proportions of cells in the large cytotoxic T/NK or proliferating T/NK subclusters as compared with control granulomas (Fig. 6 C).

Within the cells classified as cytotoxic T/NK based on their expression of GZMK, GZMH, GZMB, GNLY, and/or PRF1, expression of the CD4 gene and the genes for γδ TCR (TRGC2 and TRDC) were significantly increased in CD8α-depleted granulomas compared with control granulomas (CD4, P = 1.25E⁻⁸; TRGC2, P = 0.00027; TRDC, P = 2.2E⁻²²) (Fig. 7, A and B). This is consistent with the flow cytometry data where CD8α depletion resulted in higher frequencies of CD4 and γδ T cells expressing granzyme B and/or granulysin, indicating that these cell types acquire some cytotoxic features in the absence of CD8 cells (Fig. 4, A and B). Two CD8A/CD8B expressing subclusters were largely specific to the granulomas from control animals and not present in the depleted granulomas: the IFN-responding subcluster (IgG versus αCD8α, P = 6.38e⁻³; IgG versus αCD8β, P = 5.65e⁻²) and the non-cytotoxic CD8 subcluster (IgG versus αCD8α, P = 3.24e⁻⁸; IgG versus αCD8β, P = 9.12e⁻⁶) (Fig. 6 C).

The flow cytometry data suggested that in the absence of CD8⁺ lymphocytes, CD4 and γδ T cells could adopt cytotoxic effector expression. To further explore the observed changes in cytotoxic functions between phenotypically distinct cell types (i.e., CD4, CD8α, CD8β, γδ T cells), cytotoxic T/NK cells were subclustered from all granulomas into eight cytotoxic subclusters (Fig. 6, D and E) defined by potential functional marker genes as well as canonical phenotypic and lineage markers (CD8A, CD8B, CD4, TCRG, TCRD, and CD3) (Table S7). Unsupervised clustering of these eight cytotoxic subclusters revealed four putative functional cluster pairs (Fig. 6 E). Each of these pairs includes one member comprised primarily of cells from control granulomas that express CD8A and/or CD8B (odd number clusters) while the other is comprised primarily of cells from granulomas of CD8-depleted animals that express higher levels of CD4, TCRG, and/or TCRD (even number clusters) (Fig. 6 E). These latter types are henceforth called “moonlighting” members. The moonlighting member of each pair emerged in the context of CD8 depletion, and cells with these phenotypes were generally less prevalent in control granulomas (Fig. 6 F). Of the control-enriched clusters, two (cytotoxic 1, cytotoxic 5) are canonical CD8αβ T cells and two (cytotoxic 3 and cytotoxic 7) appear to be adaptive CD8 T cells and γδ T cells, expressing a combination of CD8A, CD8B, TRCG, and TRCD (Fig. 6 E).

The four cluster pairs with similar general phenotypes (with odd number clusters found primarily in control granulomas and even number clusters in CD8-depleted granulomas), as described by marker genes and gene ontology (GO) enrichments, include: infiltrating effector-like cytotoxic clusters 1 and 2 (marker genes: GZMK, CD28, TOX, EOMES, CCR5; GO terms: positive regulation of IL-2 production, T cell receptor signaling pathway); innate-like cytotoxic 3 and 4 (marker genes: GZMA, GZMH, GNLY, KLRD1; GO terms: positive regulation of tumor necrosis factor superfamily cytokine production, regulation of natural killer cell–mediated cytotoxicity); tissue resident effector-like cytotoxic 5 and 6 (marker genes: GZMB, GZMH, CLEC12A, CX3CR1, CD226, KLF3, ZEB2, TGFBR3); and T1/T17-like cytotoxic 7 and 8 (marker genes: GZMB, FURIN, RORA; GO terms: positive regulation of T cell activation, positive regulation of IFN-γ production) (Fig. 6 E, Fig. 7 C, Fig. S4, C–F; and Tables S4, S7, and S8).

Of these pairs, the moonlighting member of the infiltrating effector-like subclusters (cytotoxic 2) and the tissue-resident effector-like subclusters (cytotoxic 6) were dominated by CD4-expressing cells, and the moonlighting members of the T1T17-like (cytotoxic 8) and innate-like (cytotoxic 4) subclusters were enriched in γδ T cells (expressing TRGC and TRDC) compared with the control granuloma counterparts (Fig. 6, E and F). Most cytotoxic cells in control granulomas are dominated by CD8A/B T cells in the infiltrating effector-like cytotoxic 1 subcluster or T1T17-like cytotoxic 7 subcluster. Most CD8α-depleted granulomas are dominated by the tissue-resident effector-like cytotoxic 6 subclusters or one of the γδ-expressing subclusters cytotoxic 4 or T1T17-like cytotoxic 8 (Table S4).

Pair-wise pseudobulk differential expression between each of these pairs of cytotoxic subpopulations revealed that despite the shared expression of many markers, the cytotoxic NK/T cell clusters dominated by cells from control granulomas have more cytotoxic potential compared to their moonlighting counterpart (Table S9). Using cytotoxic and cytokine gene signatures described by Szabo et al. (2019), with the exception of T1T17-like cytotoxic clusters, each control cytotoxic T/NK cluster had a higher mean score on the cytotoxic gene signature than its corresponding moonlighting pair; two of the cytotoxic pairs also had higher mean scores across the cytokine gene signature (Fig. 6 G) (Szabo et al., 2019).

Compared with control cytotoxic cluster 1, the depletion-dominated moonlighting cytotoxic 2 cluster expressed lower levels of effector genes including GZMH, GNLY, PRF1, CCL5, and CTSW and higher levels of CTLA4, IL2RA, SPOCK2, MAF, SLAMF1, ICOS, and the cytokine IL10, which has a potential immunosuppressive role (Gideon et al., 2015; Huynh et al., 2018; Wong et al., 2020; Fig. 6 E and Table S9). Compared with cytotoxic 5 in control granulomas, cytotoxic 6 (appearing in CD8-depleted granulomas) expressed lower levels of GZMH, GNLY, PRF1, and CCL5, suggesting reduced killing potential (Table S9). The innate-like cytotoxic 3 subcluster, primarily seen in control granulomas, appeared to contain cells with an NK or NKT phenotype due to low expression of CD3 and high expression of CD8A and NK markers (FCGR3A, KLRC2), while the CD8α-depleted dominated moonlighting innate-like cytotoxic 4 subcluster appeared to contain more γδ T cells (Fig. 6 E and Table S9). Finally, the moonlighting T1T17-like cytotoxic 8 cluster is made up of primarily TRCG and TRCD-expressing cells (i.e., γδ T cells) while the T1T17-like cytotoxic 7 cluster contains a combination of CD4-, CD8A-, and CD8B-expressing cells (Fig. 6 E and Table S9). Like the larger T1T17 cluster, both of these T1T17-like cytotoxic clusters express cytokines IFNG and TNF, which have previously been implicated in the control of Mtb infection (Cooper et al., 1993; Flynn et al., 1993, 1995; Gideon et al., 2015; Keane et al., 2001; reviewed in Shen and Chen, 2018; Zeng et al., 2018).

In short, these data support the notion that other lymphocytes, primarily CD4 and γδ T cells, gain some, but not all, cytotoxic functions in the absence of CD8⁺ lymphocytes. We speculate that the incomplete cytotoxic potential of these moonlighting populations makes them incapable of fulfilling all the activities of CD8⁺ lymphocytes.

We hypothesized that the increase in cytotoxic non-CD8 T cells in CD8-depleted granulomas could arise from the proliferating T/NK cell subcluster. We subclustered the proliferating T/NK cells and identified three cytotoxic subclusters with marker genes similar to the non-proliferating cytotoxic subclusters, as well as a proliferating subcluster with marker genes similar to the non-proliferating T1T17 cluster (Fig. S4, E and F; and Table S4). Two proliferating subclusters expressed markers specific to other non-proliferating, moonlighting cytotoxic subclusters (cytotoxic 4 and T1T17-like cytotoxic 8) and were more abundant among T/NK cells in CD8α-depleted granulomas compared with control granulomas (proliferating cytotoxic 4, P = 3.12e⁻⁶; proliferating cytotoxic 8, P = 1.99e⁻⁴) and compared with CD8β-depleted granulomas (proliferating cytotoxic 4, P = 4.82e⁻⁴; proliferating cytotoxic 8, P = 1.77e⁻⁵) (Tables S4 and S6). Similarities between proliferating clusters and moonlighting cytotoxic clusters suggest that moonlighting cytotoxic cells in depleted granulomas may be increasing due to proliferation.

The appearance of non-CD8⁺ granuloma lymphocytes with cytotoxic gene expression in the context of CD8 depletion led us to consider which immunological factors might contribute to the emergence of these unique cytotoxic populations. To address this, we sought to identify signals that these non-CD8 cytotoxic lymphocytes might be receiving from other cells in the granuloma under CD8 depletion conditions. To have sufficient cells to assess for this analysis, we focused on a comparison of IgG versus CD8α-depleted granulomas composed of >20% T cells. In each granuloma, we defined an interaction potential for each sender cell, receiver cell, and ligand combination using gene expression values of known ligand and receptors in the putative sender and receiver cells as well as the proportion of expression of these genes in each granuloma originating from each sender cell type (see Materials and methods). We reasoned that external ligands in the granuloma microenvironment may, in the absence of CD8 cells as “receivers,” drive CD4 or γδ T cells to exhibit a more cytotoxic phenotype. To identify potential differences in targets of ligand signaling, we performed differential expression analysis between all CD4 T cells in IgG granulomas and those in CD8α-depleted granulomas. We then used a NicheNet-based cell–cell interaction analysis to compare the potential ligands differentially affecting gene expression of CD4-expressing T cells between the two conditions to ask which ligands could influence the development of cytotoxicity (Fig. S5, A and B) (Browaeys et al., 2020). We further refined these lists of putative ligand drivers to include only ligands involved in cell–cell interactions, which significantly change their interaction potential between CD8α-depleted granulomas and IgG granulomas (Fig. S5, C–E). We performed the same analysis on cells expressing both TRGC and TRDC.

This analysis of both CD4-expressing T cells and TRGC and TRDC cells (γδ T cells) identified several sender–receiver pairs involving multiple ligands sent from macrophage clusters (Fig. 8, A and B; and Table S10). The ligands CCL3 and IL15 were shared across both target lineages and had increased differential interaction potential in three moonlighting cytotoxic T/NK cell populations: CD4-expressing cytotoxic 6 and TRGC- and TRDC-expressing cytotoxic 4 and 8. Among these two ligands, the differential interaction potential of IL15 was specific to the moonlighting clusters while differential interactions via CCL3 also appeared broadly across non-cytotoxic T/NK cell subclusters (Fig. 8 A). It is possible that CCL3, a chemokine involved in T cell recruitment, is promoting more recruitment of CD4 or γδ T cells (or other cells) to the granulomas in the context of depletion since overall T cell numbers remain the same across control and depletion groups.

IL-15 is a known T cell survival cytokine that promotes memory development and effector functions, including cytotoxic molecule expression in CD8 T cells and NK cells (Berard et al., 2003; Huntington et al., 2009). The signaling of IL15 is complicated, where myeloid populations express IL-15 and/or IL-15RA, with IL-15RA required to present IL15 in trans to lymphocytes. Lymphocytes in turn express IL-15RB (IL2RB), which pairs with the common γ chain (IL2RG) to form the IL-15 receptor (Ring et al., 2012).

Since our interaction analysis implicated IL-15 as a differential myeloid:T cell interaction ligand and that it has the potential to enhance cytotoxic function and myeloid cells are abundant in granulomas, we further explored the differences in IL-15 expression in our data set. As expected, the primary expressors of IL15 and IL15RA were in the macrophage/monocyte cluster with some expression by mast cells and NK/T cells (Fig. 8 C and Fig. S5 F) (Balin et al., 2018; Ring et al., 2012; Santana Carrero et al., 2019). The T/NK cluster was enriched for IL15RB (IL2RB) and IL2RG as expected, and the former was detected in the majority of the cytotoxic subsets (including those in controls and the moonlighting clusters in the CD8-depleted granulomas) (Fig. S5 G). These data indicate that both IL15 and IL15RA are expressed in the granuloma at 6 wk after infection and that they may be influencing the changes in cytotoxic phenotypes observed during depletion.

To assess IL-15/IL-5RA in granulomas at the protein level, we performed immunohistochemistry (IHC) staining for IL15RA and IL15 on granulomas from the animals in this current study to determine the location and identity of cells expressing these proteins. IL15-RA expression was seen on myeloid cells in granulomas (Fig. 8, D and E). Moreover, there was substantial colocalization between IL-15RA and IL-15 in these cells (Fig. 8 E). This pattern of IL-15RA and IL-15 expression was seen consistently in epithelioid macrophages surrounding the caseous centers of necrotic granulomas and alveolar macrophages in the granuloma-adjacent lung tissue. Macrophages in the granuloma lymphocyte cuff also expressed both proteins although there was substantial variability in the abundance of these macrophages in different granulomas. In the epithelioid macrophage areas staining for IL-15RA and IL-15, there was evidence of CD3⁺ T cells, albeit at lower levels than in the lymphocyte cuff, suggesting that these T cells could experience IL-15 signaling (Fig. 8 E). Similar expression patterns were observed in granulomas from IgG-treated control animals (Fig. 8, D and E) and animals that were depleted of CD8α⁺ or CD8β⁺ lymphocytes (Fig. 8 F). These data support our transcriptional data suggesting that macrophage IL-15 expression promotes lymphocyte function and cytotoxicity in granulomas from control and depleted animals.

To look for the putative source of IL-15, we further subclustered the “macrophages/monocytes” subset and identified four subpopulations: Mac 1 expressing alveolar macrophage marker genes (MARCO, PPARG, and FABP4), Mac 2 expressing genes associated with immunosuppressive TREM2⁺ macrophages (TREM2 and TCN2), and Mac 3 expressing inflammatory markers (CXCL11), as well as a monocyte cluster (Fig. 5, E–G). The composition of the macrophage/monocytes was quite variable between granulomas, but mean IL15 expression was comparable across all subclusters and depletions (Fig. 5, G and H). The Alveolar Mac 1 and monocyte clusters were both significantly increased in proportional abundance among macrophage/monocyte clusters in CD8α-depleted granulomas compared to IgG with substantial heterogeneity in composition between granulomas (Fig. 5 F); however, we note that depletion did not result in a statistically significant shift in the overall proportion of macrophages/monocytes, only a trend toward an increase (Fig. 5 D).

When comparing granulomas from the two depletions with those from our IgG controls, we also noted that the cDC, macrophages/monocytes, and mast cell clusters showed altered gene expression (Table S11). In the macrophage and monocyte cells, genes related to response to IL1 and response to type 1 IFN were enriched in genes up in CD8α-depleted granulomas but not CD8β-depleted granulomas compared with IgG (Table S11). These pathways were increased in CD8α-depleted granulomas in each subcluster (P < 0.05 Mann–Whitney) (Fig. 5 I). Similarly, genes specifically up in cDC cells from CD8α-depleted granulomas but not CD8β-depleted granulomas compared with IgG were enriched for cytokine production pathways including IL-6 and IL1β production (Table S11), and gene scores for these pathways were significantly increased (P < 0.05 Mann-Whitney) in CD8α-depleted samples compared with IgG in both cDC subclusters (Fig. 5 N). Increased IL1β and IL-6 signaling and response in cDC and macrophage cells point to potential changes in the pathological inflammatory response after CD8α depletion.

While CD4 T cells have long been associated with protective immune responses against TB, the roles of innate and adaptive CD8⁺ lymphocytes have been less clear. Here, we show that CD8⁺ lymphocytes are important in limiting the establishment and dissemination of infection, demonstrating for the first time the critical importance of CD8⁺ lymphocytes for controlling early Mtb infection in an NHP model that recapitulates human Mtb infection. Our results support that polyfunctional cytotoxic CD8⁺ lymphocytes are important in limiting the establishment and dissemination of Mtb infection. A surprising finding was that the depletion of CD8⁺ lymphocytes led to the incomplete adoption of cytotoxic molecule expression by CD4 and γδ T cells in granulomas, which could not completely substitute for CD8⁺ cytotoxic lymphocytes. Our data suggest that expression of IL-15 in granulomas may reprogram CD4 and γδ T cells to express some cytotoxic molecules in the absence of CD8⁺ lymphocytes. These data support that enhancing cytotoxic cell subsets in a vaccine strategy may improve efficacy against Mtb infection or disease.

Analysis of Mtb barcodes supported that CD8α⁺ cells play a critical role in the limiting establishment of lung granulomas. The increased number of unique bacteria establishing granulomas in the lungs in the absence of all CD8 lymphocytes indicates that CD8⁺ cells contribute to a modest bottleneck of infection. Even prior to infection, there are both CD4 and CD8αβ T cells as well as innate lymphocytes in the airways where the Mtb bacillus initially encounters host cells, and depletion of the CD8⁺ cells increased the number of bacteria that transit to the lungs to establish granulomas. The CD8 cells in the airways may modify the macrophage environment or directly kill intracellular or extracellular Mtb, limiting the number of Mtb that can establish a granuloma.

Our previous data indicate that each granuloma is formed from a single bacillus which then replicates over the next 4–6 wk to reach a peak of ∼10⁵ CFU (Grant et al., 2022; Lin et al., 2014). This appears to be the bacterial carrying capacity of any single granuloma, at which point bacterial replication is either controlled or the bacteria disseminate to form new granulomas. Dissemination during the early phases of infection is associated with the development of active TB in cynomolgus macaques (Coleman et al., 2014). There were more granulomas in the CD8-depleted animals, and barcode analysis demonstrated increased dissemination, not only within lung tissue but also to lymph nodes and extrapulmonary sites. Increased dissemination suggests that granulomas in the absence of CD8 cells either reach carrying capacity sooner, suggesting that CD8 cells may impose some restraint of initial Mtb replication or are organized within the lymphocyte cuff of granulomas as a barrier to bacterial dissemination. Whether this early control by CD8 cells is either due to cytotoxic effectors that result in the lysing of infected cells or direct killing of Mtb or by other effector functions that maintain the integrity of the granuloma is yet to be determined. At this time, reagents do not exist to directly address the functional cytotoxic contribution of specific canonical cell types (such as γδ T cells or NK cells) through the depletion of these cell types in macaques. However, increasing the number of γδ T cells was beneficial in Mtb-infected macaques (Shen et al., 2019). It is possible that the loss of half the γδ T cells in our study may be partially responsible for the increased disease in CD8-depleted animals.

Granulomas at 6 wk after Mtb infection have a diverse landscape of lymphocytes, with conventional CD4 and CD8 T cells and B cells being the primary lymphocytes represented, with several other populations present at lower levels including γδ T cells, and NK and NK T cells. Cytotoxic signatures were dominant in CD8⁺ cells and overall there were relatively low frequencies of cytokine-producing T cells in granulomas, as previously published (Gideon et al., 2015; Grant et al., 2022). In the context of depletion, we saw no significant differences in overall numbers of lymphocytes or T cells in granulomas, despite the efficiency of CD8α and CD8β depletion, suggesting either increased recruitment to or proliferation of non-CD8 lymphocytes in CD8-depleted granulomas. scRNAseq data support both possibilities, with increased T/NK cell proliferation and increased expression of the chemokine CCL3 in CD8-depleted granulomas. CCL3 recruits CCR5-bearing T cells, shown previously to be present in granulomas (Lin et al., 2006).

After CD8 depletion, both flow cytometry and scRNAseq analyses demonstrated that CD4 and γδ T cells in granulomas had gained expression of cytotoxic molecules. The scRNAseq data uncovered eight unique cytotoxic subclusters that upon unsupervised analysis sorted into four “pairs” related to control or depletion status. These paired subclusters bore similar cytotoxic transcriptional signatures (granzymes although with varying levels of expression upon depletion), yet the subcluster pairs present in depleted granulomas differed in lineage markers (CD4 and γδ TCR in contrast to CD8α or CD8β) and notably had absent or reduced expression of some key effector molecules: perforin, granulysin, and CCL5. We speculate that these moonlighting CD4 and γδ T cells in CD8-depleted granulomas were not fully capable of cytotoxic activity compared with CD8 T cells in non-depleted granulomas and could not compensate for the lack of cytotoxic CD8 lymphocytes.

Our interaction analysis of scRNAseq data revealed IL15, a cytokine that induces and promotes cytotoxicity, as a ligand whose downstream targets were enriched in CD4 and γδ T cells in CD8α depletion compared with control. IL-15 is presented in trans to cells expressing IL15Rβ (IL-2Rβ) and the common γ chain by IL-15Rα. IL-15Rα and/or IL-15 transcripts were expressed by myeloid subsets, including macrophages and mast cells, in granulomas. While we did not observe a significant increase in IL15 or IL15RA transcripts upon depletion, CD8 depletion effectively removed many of the cell types that would canonically interact with and respond to IL-15, potentially allowing cell types such as CD4 T cells and γδ T cells to utilize this pathway to greater efficiency and skew their functions towards cytotoxic phenotypes. IHC analysis of granulomas confirmed colocalization of IL-15RA and IL-15 on macrophages, particularly those in the epithelioid macrophage layer of the granuloma. Although most lymphocytes are in the outer layer of granulomas, we detected T cells in proximity to epithelioid macrophages in which IL-15Rα/IL15 were colocalized, suggesting that these cells could signal to T cells within this region of the granuloma. Even in the CD8β-depleted animals, where there is only a partial loss of all CD8⁺ cells, there was increased expression of granzyme B and/or granulysin by CD4s and γδ T cells, supporting our hypothesis that the loss of even a subset of CD8⁺ cells could result in other cells responding more to IL15. It is of course possible that CD8⁺ lymphocytes regulate cytotoxic molecule expression by T cells in other ways, either directly or indirectly, and that other molecules, including IL-2, could enhance cytotoxic molecule expression by CD4 and γδ T cells in the absence of CD8 cells. Future studies are required to investigate these hypotheses more fully, but the importance of cytotoxic function in granulomas suggests that IL-15 may be an important contributor to the control of Mtb infection and could be harnessed as an adjuvant in vaccines, particularly those that are poor at induced CD8 T cell cytotoxicity.

Of note, the endpoint of the current study was 6 wk after infection, an early time point in the context of Mtb infection and adaptive immune responses to Mtb. In our model, bacterial killing is not consistently observed until 10–11 wk after infection (Lin et al., 2014), a time point that coincides with the induction of the adaptive immune response and conventional CD8αβ T cell activation in macaque granulomas (Grant et al., 2022). Although we did not see a strong effect of CD8β depletion on overall bacterial burden or dissemination at the early 6 wk time point, the bacterial burden in the thoracic lymph nodes in CD8β-depleted animals was significantly higher than in control animals, and the lymph node bacterial burden is an early indication of dissemination. It is also possible that the loss of CD8αβ T cells led to increased Mtb growth in the lymph nodes. Given that the adaptive T cell response is slow to become functional in granulomas (Grant et al., 2022), we suspect that extended CD8β depletion beyond 6 wk would result in worse disease outcomes. While innate CD8 cells or the combination of innate and adaptive CD8 cells are required for early control of establishment and dissemination, we postulate that conventional CD8αβ T cells will be required to maintain this control, presumably through cytotoxic effector functions. However, we also initiated anti-CD8β depletion 6 wk prior to infection, while anti-CD8α antibodies were initiated 2 wk prior to infection. Thus it is formally possible that the longer depletion time with anti-CD8β antibodies could have resulted in compensatory effects on the immune system that were not detected in our analyses. Another limitation of this study is that anti-CD8α antibodies also target CD4⁺CD8α⁺ T cells, which are a smaller population compared with CD4⁺CD8α⁻ T cells in granulomas (8% versus 30%). It is thought that these double-positive CD4 T cells are of a more “activated” phenotype and thus depletion of these cells could contribute to the increased dissemination phenotype observed in CD8-depleted animals. Our analysis of the function of double-positive CD4 T cells indicates that they are similar to single-positive CD4 T cells in terms of cytokine production, although slightly higher in terms of cytotoxic function. It remains a possibility that the depletion of these cells influences the outcome of infection.

In summary, our data demonstrate an important role for CD8⁺ lymphocytes that primarily express cytotoxic functions in controlling early Mtb infection both in terms of limiting the establishment and dissemination of infection. Further studies will be required to address which specific cell subsets contribute to the protection, although we hypothesize that it is not a single cell subset but rather a combination of cell types and immune functions. TB is a complicated disease and the immune response has proven to be equally complex. Our study sheds light on the importance of the early immune interactions during Mtb infection, including innate CD8 responses. While targeting unconventional T cells (γδ T cells, NK, NKT, and MAITs) for vaccine therapeutics has only recently gained traction, our study highlights the potential importance of these cell types in early disease control. To combat an unconventional disease like TB, we must think outside the conventional standards of vaccine and therapeutic strategies.

All experimental manipulations, procedures, protocols, and care of the animals were approved by the University of Pittsburgh School of Medicine Institutional Animal Care and Use Committee (IACUC). The protocol assurance number for our IACUC is A3187-01. The specific protocol approval number is 18073085. The University of Pittsburgh’s IACUC adheres to national guidelines established in the Animal Welfare Act (7 U.S.C. Sections 2131–2159) and the Guide for the Care and Use of Laboratory Animals (eighth edition), as mandated by the U.S. Public Health Service Policy.

All housing, care, and experimental procedures were approved by the University of Pittsburgh School of Medicine IACUC. Examination of animals was performed in quarantine to assess physical health and confirmation of no previous Mtb infection was performed. Cynomolgus macaques (Macaca fascicularis) (n = 15) were purchased for this study from Valley Biosystems. Animals were randomly assigned to the three depletion groups (IgG, n = 5; anti-CD8α, n = 5; anti-CD8β, n = 5). Rhesus recombinant depletion antibodies (control IgG1 DSPR1, anti-CD8α antibody M-T807R1, and anti-CD8β antibody CD8b255R1) were obtained from the National Institutes of Health NHP Reagent Resource and administered at 50 mg/kg/dose (IV) every 2 wk. Based on discussions with other NHP researchers who indicated that CD8β depletion was not as efficient as CD8α depletion, we chose to begin depletion with anti- CD8β antibody 6 wk prior to infection while anti-CD8α antibody was only initiated 2 wk prior to infection. The first dose of depletion was administered subcutaneously while all subsequent infusions were delivered intravenously. Each infection cohort included control, CD8α, and CD8β depletion animals to reduce any experimental bias between infection cohorts. All animals were infected via bronchoscopic instillation with 19–23 CFU of a bar-coded Erdman strain of Mtb (Martin et al., 2017). Administration of depletion antibodies continued after infection every 2 wk until necropsy (6 wk after infection) (Fig. S1 and Table S1).

PET CT imaging was performed with FDG probe as previously described. A subset of animals (IgG, n = 1; anti-CD8α, n = 2; anti-CD8β, n = 2) was scanned after depletion but prior to infection to ensure depletion had no effect on baseline inflammation. Serial scans were performed at 3 wk after infection and prior to necropsy at 6 wk after infection. Total lung FDG activity and number of granulomas at each time point were measured as previously described (White et al., 2017).

Cynomolgus macaques (n = 15) utilized for this study were necropsied as previously described (Gideon et al., 2022; Lin et al., 2014; Martin et al., 2017). In summary, a variety of tissues (granulomas, lung lobes, lymph nodes, liver, and spleen) were excised and homogenized into single-cell suspensions, with subsets of tissues sent for histopathological analysis. Granulomas were identified and excised by PET CT image-guided procedures. Individual lesions (of at least 2 mm) were cut in half with one-half processed for homogenization and the other for histological analysis. Bacterial burdens were assessed from homogenized tissues, as previously described, using serial dilutions on 7H11 agar plates and incubated for 21 d at 37°C in 5% CO₂ for CFU enumeration. A subset of tissue was also processed using enzymatic dissociation for scRNAseq.

Single-cell suspensions of granulomas were acquired as previously described. Cells were incubated with 0.1% brefeldin A (Golgiplug; BD Biosciences) in RPMI 1640 containing 1% HEPES, 1% L-glutamine, 10% human AB serum for 3 h at 37°C in 5% CO_2. After blocking, samples were stained with a fixable viability dye (Invitrogen) prior to staining surface and intracellular antigens using standard protocols (BD Cytofix/Cytoperm; BD Bioscience). Samples were acquired on a four-laser Cytek Aurora spectral cytometer (16V-14B-10 YG-8R; Cytek) using SpectroFlo software (Cytek). Final analyses were completed in FlowJo (Treestar, Inc). Table S2 provides a list of antibodies and clones for applicable samples.

To track the establishment and dissemination of Mtb, bar-coded Mtb Erdman was used for infection and barcodes were determined after necropsy as previously described (Martin et al., 2017). Serial PET CT scans were used to classify granulomas as those seen at 3 or 6 wk. Scan-matched TB lesions were harvested at necropsy, and other gross pathologies, lung lobes, peripheral and thoracic lymph nodes, spleen, and liver were homogenized and plated onto 7H11 agar plates. After 3 wk, the plates were counted for CFU. Quantitative bacterial burden measurements were performed as previously described (Maiello et al., 2018). For barcode analysis, the Mtb colonies on plates were DNA extracted using phenol–chloroform methods. After DNA purification, samples were subjected to amplicon-based sequencing to identify all the barcode tags present and shared across tissues.

A library of digitally barcoded plasmids containing 7-mer barcodes and adjacent 74-mer library identifier sequences was stably transformed into the bacterial chromosome of Mtb Erdman, as previously described (Martin et al., 2017). Each library contained ∼16,000 unique barcodes and three independently generated libraries were combined into a master library to increase barcode diversity, thereby ensuring a <2% chance that a barcode would be represented twice if 20 bacteria were randomly selected (Martin et al., 2017).

To identify and quantify the q-tags and barcodes present in each Mtb in each tissue sample, purified bacterial genomic DNA plated from tissue samples (see above) was diluted to 10 ng/µl in dH₂O and amplified with Q5 polymerase (catalog #M0491; New England Biolabs) in two rounds of PCR. The first round of PCR (annealing for 30 s at 62°C, extension at 72°C for 30 s, 4 cycles) tagged gDNA with primers to introduce a unique mcount motif to count the original template and correct for uneven skewing due to amplification jackpot events. 5 µl of PCR product from the first round of PCR was used as the template for a second PCR reaction (annealing for 30 s at 62°C, extension at 72°C for 30 s, 25–30 cycles) to add Illumina adapter sequences. Amplified dual-indexed product was purified using AMPure XP beads (catalog #A63881; Beckman Coulter), the amount of DNA was quantified using a Qubit dsDNA Broad-Range assay kit (catalog #Q32850; Thermo Fisher Scientific), and samples were pooled to a final concentration of 8 pM (concentration checked using Kapa Library Quantification Kit (catalog #KK4824; Roche) before spiking with 20% PhiX control library (Illumina) to increase the sequencing diversity. The pool was sequenced on the Illumina MiSeq system using V2 chemistry for 300 cycles. Typical total reads from a library are around 40 million with 35 million PF reads (passed illumina filter) and each individual tissue sample should make up 1–3% of the total PF reads.

Sequencing data was demultiplexed on Basespace (Illumina), downloaded as Fastq files, and further analyzed using a custom Perl script (available upon request) to quantify individual barcodes in each sample. Briefly, the script required that a sample had more than 10,000 total reads and searched Read 1 sequences for three motifs indicative of a valid barcoded strain:

• (1)

  An mcount motif that labels each individual PCR amplification event: ([ATCG] C[ACTG]C[ACTG]C[ACTG]GCG​CAA​CGC​G);

• (2)

  A barcode region motif: (CGA[ACTG]C[ACTG]AAT​TCG​ATG​G); and

• (3)

  A strain-specific qTAG region combining the two motifs above: (MCOUNT_MOTIF) ([ATCG]+) (BARCODE_MOTIF) ([ATCG]+TGG​TGT​TCA​AGC​TT) ([ATCG]).

A wobble was inserted at the start of the mcount motif. In between the mcount motif and the barcode motif is a fixed sequence that is not explicitly defined in the regular expression. In between the barcode motif and the last 12 bases is another fixed sequence coded by variable and fixed bases in the regular expression. The last 12 bases of the match are in the variable region of the qTAG and are used to identify the library. The combination of fixed and variable sequences allows specificity in the match without being overly stringent. The searches must yield over 30 bases of exact matches to pass the threshold to be defined as a valid qTAG read.

Previous sequencing of barcoded Mtb strains yielded a database of up to 16,384 distinct barcodes present in our library of isogenically tagged Mtb strains (Martin et al., 2017). The custom script mapped the 12-base qTAG sequences detected by sequencing against this library of known barcodes, allowing up to two mismatches. This yields a list of counts for each barcode found in a given sample. To eliminate reads with putative sequencing errors and chimera sequences (in which a read contains a qTag that is associated with a different barcode than it should), reads with lower counts that differed by one base pair from another higher counted barcode were deleted. This was iteratively done for all barcodes in the list. The remaining barcodes and counts in the list were log₂ transformed, and the second derivate of rank-ordered data was calculated for each pair of normalized adjacent counts in the list (e.g., (1,2), (2,3), (3,4), etc.). The lowest value is used to separate true counts from the background. Finally, every sample’s analysis was manually reviewed to ensure consistency and check for aberrant sequences such as missed chimeras. The resulting data consisted of tables that quantified all unique barcodes in a batch of samples on the same sequencing run and the fraction of reads assigned to each barcode in a sample.

For each animal, the resulting table of identified barcodes in each tissue and associated metadata including tissue CFU burdens and time of lesion detection into text formats that were compatible with plotting by Circos software (Krzywinski et al., 2009). Each plot was then manually adjusted such that links between tissues sharing barcodes were drawn, when possible, to originate from lung sites detected early by PET-CT with the highest CFU burden; however, other lung sites with the same barcode may also be the source of dissemination between tissues. Finally, separate configuration files were created for each animal, specifying relevant tracks to plot, and plotted using Circos. All custom scripts and relevant configuration files used to generate the Circos plots are available upon request.

Data were tested for normality using the Shapiro–Wilk test. One-way ANOVA with Dunnett’s multiple comparison adjustment (for normal data) or Kruskal–Wallis test with Dunn’s multiple comparison adjustment (for non-parametric data) was used to compare the treatment groups. There were two comparisons (each treatment group against the control group: IgG). Cytotoxic and cytokine effector functions within granulomas were compared using the Wilcoxon matched-pairs signed rank test. For immunological outcomes of granulomas, super-plots show the variability within each animal. The mean (per animal) was calculated for each outcome before running a statistical test to prevent within-animal bias. Summary immunologic outcomes for granulomas and multiple Kruskal–Wallis tests (in Fig. 3 A) were calculated in JMP Pro (version 14.3.0). All other statistical tests were conducted in GraphPad Prism 9 (version 9.4.1 for macOS).

Seq-Well S³ was performed as described previously (Gierahn et al., 2017; Hughes et al., 2020). Briefly, about 15,000 cells were loaded onto each Seq-Well array preloaded with uniquely barcoded mRNA capture beads (ChemGenes). Arrays were washed with protein-free RPMI media and then sealed with polycarbonate membranes. Arrays were incubated at 37°C for 30 min to allow membranes to seal, then transferred through a series of buffer exchanges to allow for cell lysis, transcript hybridization, bead washing, and bead recovery from arrays post membrane removal. Reverse transcription was performed with Maxima H Minus Reverse Transcriptase (Thermo Fisher Scientific) and excess primers were removed using an Exonuclease I digestion (New England Biolabs). Samples were stored at 4°C to await further processing. Second-strand synthesis was performed and whole transcriptome amplification (WTA) by PCR was performed using KAPA Hifi PCR Mastermix (Kapa Biosystems). WTA product was purified using Agencourt Ampure beads (Beckman Coulter) and dual-indexed 3′ digital gene expression (DGE) sequencing libraries were prepared using Nextera XT (Illumina). Libraries were sequenced on a NovaseqS4 with a paired-end read structure (R1: 20 bases; I1: 8 bases; I2: 8 bases; R2: 50 bases) and custom-sequencing primers.

Libraries were sequenced to depth as determined by saturation curves and aligned to macFas5 using the Dropseq-tools pipeline on the Terra platform (app.terra.bio). Predicted doublets were identified using scrublet (Wolock et al., 2019) on each sample using default parameters. Count matrices were combined and analyzed using Scanpy (version 1.7.2; Wolf et al., 2018). Cells were filtered with min_umi >750 and <120,000, min_genes = 500, percent _mito <0.15, and total_genes <10,000. Genes expressed in fewer than five cells were removed. Expression was normalized to 1e4 total counts per cell and the log base 2 + 1 was taken. Number of cells filtered from each sample and quality metrics of these cells are shown in Fig. S4 A and Table S4.

After filtering, iterative subclustering was used to identify cell clusters and remove low-quality groups of cells missed by the initial filtering steps. More specifically, on the full filtered count matrix, the top 2,300 variable genes using scanpy highly_variable_genes function with batch key = monkey, to prevent inclusion of genes only variable in two or fewer monkeys, was used to identify variable genes. As is common practice, the expression was scaled to 10 and the neighborhood graph was built with the scanpy.pp.neighbors function with parameters neighbors = 15 and pcs = 28. The cells were clustered using Leiden clustering with a resolution of 0.7 and 25 clusters were identified. Marker genes for each cluster were identified using the “rank_genes_groups” function in scanpy with default parameters. Clusters were classified into five groups for further subclustering based on markers genes: NK/T cells, macrophage/neutrophil/DC clusters, B cells, non-immune cells, and mast cells. Clusters coexpressing lineage marker genes from multiple groups were classified with the related group whose subclustering was executed first. For example, a cluster containing cells that express NK/T cell markers as well as cells that express macrophage cell markers would be assigned to the NK/T cluster for further filtering as this subclustering was executed first. At this stage, it can be identified if individual cells in that cluster coexpress multiple lineage markers, indicating that they may be doublets, or if a subset of these cells should be classified as a macrophage and included in macrophage subclustering. We found that while different clustering parameters identified different numbers of clusters, the final classification of cells into groups for subclustering did not change substantially according to the choice of clustering parameters. During each stage of iterative subclustering, variable gene selection and dimensionality reduction were repeated on only the subset of cells included in that iteration to identify features most relevant to that subset of cells. Parameters for each subclustering step were selected such that the maximum number of clusters with unique marker genes were identified. Marker genes of general cell type clusters after iterative filtering were determined using the Scanpy “rank_genes_groups” function, using a Wilcoxon test, filtered for log₂FoldChange >1.0, fraction of cells in group expressing >0.6, fraction of cells outside of group expressing <0.7, and P < 0.05 (Tables S6 and S7).

Feature selection and dimensionality reduction were repeated on only cells classified as NK/T cells during the clustering of all cells, and clusters expressing genes from other cell lineages were removed. The remaining NK/T cells were then subclustered into a total of 17 clusters (Leiden resolution 0.7; Fig. S4 C). Clusters 4, 0, and 9 were labeled “T1T17” clusters due to expression of IL7R, RORA, BHLHE40, NFKBIA, FURIN, as described in Gideon et al. (2022) (T1T17 mean frequency of T cells per granuloma: 23.8 ± 9% IgG, 20 ± 9 αCD8α, 31.5 ± 16.6 αCD8β). Cluster 13 was labeled “interferon responding” due to the unique expression of genes involved in IFN response such as IFI6 and MX1 (IFN-responding mean frequency of T cells per granuloma: 2.3 ± 4% IgG, 0.5 ± 0.7 αCD8α, 0.7 ± 0.7 αCD8β). Cluster 10 was labeled “Treg” due to expression of Treg markers (TIGIT, FOXP3) (Treg mean frequency of T cells per granuloma: 2.9 ± 1.2% IgG, 3.2 ± 2.8 αCD8α, 4.9 ± 3.5 αCD8β). Clusters 2 and 16 were identified as stem-like clusters as described in Gideon et al. (2022) (stem-like mean frequency of T cells per granuloma: 14.3 ± 7.2% IgG, 5.8 ± 4 αCD8α, 11.1 ± 5.8 αCD8β). Cluster 12 was labeled “non-cytotoxic CD8” due to the expression of CD8A and CD8B without the expression of cytotoxicity markers. Marker genes and mean abundances per group are shown in Table S4.

Feature selection and dimensionality reduction were then repeated on only cells classified as cytotoxic NK/T cells during NK/T cell subclustering and they were subclustered into a total of 16 subclusters (Leiden resolution 1.0), which were then grouped into nine consensus clusters based on shared marker genes (Fig. S4 D). All cells classified as proliferating in either the NK/T cell subclustering (Fig. S4 C) or cytotoxic NK/T subclustering (Fig. S4 D) were further subclustered to identify heterogeneity in the proliferating populations (Leiden resolution 0.5) (Fig. S4 E). Proliferating NK/T subclusters were named based on the expression of non-proliferating cytotoxic subcluster marker genes. Marker genes and mean abundances per group are provided in Table S4, and full marker genes are listed in Table S5.

cDC cells were subclustered as described for NK/T cells with a resolution of 0.07 (Fig. 5 J). Three clusters were identified and the two were merged based on shared marker gene expression. Clustering of the macrophages/monocytes cluster required batch correction due to donor-level batch effects. This was performed using pyharmony (Korsunsky et al., 2019). Subclusters were then determined as described for NK/T cells using a Leiden resolution of 0.7 resulting in 12 clusters, which were merged based on shared marker genes into four clusters.

We aimed to define marker genes for each subcluster as genes that identify cells in that cluster regardless of their sample of origin. Given monkey-to-monkey variability and sample-to-sample differences, we generated pseudobulk counts for the cells in each subcluster in each granuloma by summing the raw counts for each group and used these to identify differentially expressed genes between subclusters and between conditions in a sample and monkey aware manner (Amezquita et al., 2020; Crowell et al., 2020; Tung et al., 2017). We removed subcluster/sample pools with fewer than 10 cells from these comparisons. Using the DESeq2 package, we performed differential expression analyses using the Wald statistical test with the design formula “∼Monkey + cell type” to identify cell type marker genes and “∼Monkey + treatment” to identify condition-specific marker genes (Love et al., 2014). Afterward, genes were then filtered for adjusted P value <0.05, log₂FoldChange greater than 0.2 or less than −0.2 (Tables S5 and S7).

GO analysis was run by first calculating the marker genes for each moonlighting pair of cytotoxic NK/T cell subclusters compared to other cytotoxic NK/T cell subclusters using pseudobulk differential expression analysis as described above. Functional enrichment analysis on these marker genes was performed using Enrichr using the gseapy package with the gene set GO_Biological_Processes_2021 (Ashburner et al., 2000; Fang et al., 2023). Full results of these enrichment analyses are included in Table S8 and visualization in Fig. 7 C includes manually pruned GO terms. For each GO term, the mean gene set score was calculated for each cytotoxic NK/T subcluster using the full list of genes assigned to that GO term with the scanpy function “score_genes.” Heatmap visualizations display per-subcluster mean gene set score for all genes in the GO term z-scored across subclusters using the seaborn “clustermap” function.

Analogous analysis was run on differential genes between each moonlighting pair of cytotoxic NK/T cell subclusters, and the results passing an adjusted P value threshold of 0.05 and a minimum of four genes overlapping each gene set are included in Table S9. GO analysis of genes downregulated in CD8α-depleted versus IgG control granulomas and CD8α-depleted versus IgG control granulomas for each non-T/NK cell general cell type cluster was also run and the results yielding any GO terms are included in Table S11. Results on lists of genes upregulated in CD8α-depleted but not CD8β-depleted graunlomas compared with IgG are also included in Table S11.

The differential abundance of groups of scRNAseq clusters between treatment groups was evaluated using a generalized linear model (GLM) QL test implemented in the edgeR package as described previously (Amezquita et al., 2020; McCarthy et al., 2012). Briefly, counts of the number of cells in each cell type for each sample were loaded into a matrix. The design matrix was specified as “∼0 + factor(Monkey) + factor(treatment)” and the glm was fit using the function “glmQLFit” with parameters “robust = TRUE” and “abundance.trend = FALSE” because sample numbers were too low to support a trend. Pairwise significant differences were identified from the results matrix output from the function “qlmQLFTest” with the coefficient corresponding to the pairwise comparison of interest. Results with FDR < 0.05 were reported as significant. Separate models were run for composition over general cell types (Fig. 5 D), T/NK cell clusters (Fig. 6 C), and cytotoxic T/NK cell subclusters (Fig. 6 F). Full differential abundance analysis results are provided in Table S6.

NicheNet (Browaeys et al., 2020) was used to identify putative differential ligands between CD4-expressing cells from granulomas with sufficient NK/T cell proportions (Table S4). First, qualifying single cells were identified as cells in the T/NK cell cluster (Fig. 5 A) not expressing CD8A nor CD8B and expressing CD4 where expression is defined as a value >0 in the count matrix for CD4-expressing cells analysis. For γδ-expressing cell analysis, qualifying cells identified as cells in the T/NK cell cluster (Fig. 5 A) expressing TRDC and TRGC but not CD8A were considered γδ T cells. Pseudobulk count matrices were generated for each granuloma in clusters 2 or 4 by summing the raw count values of qualifying cells. DESeq2 was run to identify differentially expressed genes between CD8α-depleted and IgG control granulomas as before with the design formula “design = ∼ monkey+ gran_cluster.” A gene list consisting of significantly (P_adjusted < 0.05) differentially expressed genes with log₂FoldChange of CD8α-depleted over IgG control granulomas >0.5 and expressed in >10% of qualifying cells was used as the geneset of interest for NicheNet analysis. NicheNet analysis is used to identify ligands whose targets are enriched in this geneset of interest which was run with ligands of interest consisting of ligands expressed in >10% of some other cell type cluster. The top 25 Pearson-correlated ligands found via the “predict_ligand_activities” function were designated as the top ligands (Fig. S5 and Table S10).

To identify potential cell–cell interactions within each granuloma using expression of ligand and receptor, we devised an interaction potential score for each potential receptor “R,” ligand “L,” sender cell type “Cs,” and receiver cell type “Cr” combination in each granuloma G using the formula: (mean R expression in Cr) * (mean L expression in Cs) * (number of cells in Cr expressing R)/(number of cells in G expressing R) * (number of cells in Cs expressing L)/(number of cells in G). This score increases with increasing mean expression level of the ligand or receptor. It takes into account the decreased theoretical potential for a receiver cell to come into contact with a ligand if more cells in a granuloma are expressing the receptor for that ligand. Ligand interaction potentials were calculated for each sender–receiver pair in each granuloma by summing the ligand–receptor potentials involving that ligand for that sender–receiver pair. To identify differential sender–receiver ligand interactions between groups of granulomas, a log₂FoldChange and P value from a t test was calculated. Interaction potentials with P < 0.05 were retained (Table S10 and Fig. S5 D).

Ligands identified by NicheNet as differentially increased in CD8α-depleted or IgG control granulomas were used to filter the statistically significant changing sender–receiver ligand scores derived above. Any significantly changing interaction potential whose ligand was identified by NicheNet as one of the top 25 ligands increased in either CD8α-depleted granulomas or IgG granulomas and whose receiver cell type was a CD4-expressing cluster are visualized (Fig. 8 A and Fig. S5, E and F).

IHC was performed on formalin-fixed paraffin-embedded lung granulomas that were cut into 5-µm thick sections and applied to SuperFrost Plus Slides (Thermo Fisher Scientific). Sections that contained individual granulomas were selected for staining and slides from each group were processed and stained in parallel in a blinded fashion. Slides were deparaffinized and antigen retrieval was performed in a Retriever (Electron Microscopy Sciences) in buffer containing 20 mM Tris/820 μM EDTA/0.0001% Tween 20 (pH 9.0) as previously indicated (Talukdar et al., 2022) followed by 30 min of blocking with Universal Fc Receptor Blocker (Innovex Biosciences) and 30 min in 1% BSA-PBS. The staining was performed in two rounds to accommodate an expanded antibody panel and multiple antibodies from the same host with the same isotype. In the first round, sections were stained for CD3 (rabbit polyclonal, 1:50 dilution; Agilent) and IL15RA (mouse IgG2a, 1:100 dilution; Novus Biologicals). A section on the same slide stained with equivalent concentrations of isotype control antibodies was processed in the same manner as the positively stained slides. The anti-IL-15RA antibody was detected with the Excel anti-mouse IgG DyLight 594 amplification kit (Vector Laboratories) according to the manufacturer’s instructions, and the anti-CD3 antibody was detected afterward with an AF488-conjugated donkey anti-rabbit secondary antibody (1:1,000 dilution, Thermo Fisher Scientific) for 1 h at room temperature. All antibodies were diluted in 1% BSA-PBS, and between each step, the slides were washed three to five times with PBS. The slides were washed and CD11c was stained for CD11c (mouse IgG2a, 1:30 dilution; Leica Microsystems) with AF647-conjugated anti-CD11c antibodies using the Zenon mouse IgG2a-labeling kit (Thermo Fisher Scientific). The Zenon-labeled CD11c conjugate was incubated overnight at 4°C, and at the completion of the first round of staining, a coverslip was applied with ProLong Gold Mounting Medium containing DAPI (Thermo Fisher Scientific). The sections were imaged with a Nikon e1000 epifluorescence microscope (Nikon) with an X-Y motorized stage (Prior Scientific) and Nikon DS-Qi2 camera operated by NIS Elements AR software version 5.002 (Nikon).

After imaging, the slides were incubated in deionized water until the coverslips fell off and then were then washed with dH2O and PBS to remove the residual mounting medium. Antibodies from the first round were stripped by repeating the antigen retrieval process in a Retriever and loss of staining and maintenance of tissue integrity was visually confirmed. The AF488 and AF647 antibodies were effectively removed by this process whereas the DyLight 594 signal (IL-15Rα) was reduced but not eliminated, so the red (TRITC) channel was not used in the second round. Slides were blocked before the application of anti-CD20 (mouse IgG2a; 1:100 dilution; Dako) and anti-IL-15 (goat polyclonal; 1:30 dilution; Novus Biologicals) antibodies. The anti-IL-15 antibody was directly labeled with AF488 with the Zenon goat IgG-labeling kit (Thermo Fisher Scientific) and slides were incubated with this conjugate and anti-CD20 overnight at 4°C. After washing, the sections were stained with AF647-conjugated goat anti-mouse IgG2a isotype-specific secondary antibodies (1:1,000 dilution; Jackson ImmunoResearch) for 1 h at room temperature. After incubation with the secondary antibodies, coverslips were mounted as before. Positively stained and isotype control sections were imaged with the same settings and the same look-up tables (LUTs) were applied to both sections. The microscope is a Nikon e1000 epifluorescence microscope and images were taken at 20× magnification with a Nikon Plan APO lens with a 0.75 numerical aperture. Stitched image fields for four colors/images were acquired with a Nikon Ds-Qi2 camera as 14-bit ND2 files with the look-up tables for the individual channels based on the isotype control image before being exported as TIFF files. Composite five-color images were created in Adobe Photoshop (Adobe Microsystems) by aligning the DAPI images of two independent images of the same tissue so that the nuclei were aligned at the individual cell level and registration was confirmed at the nucleus level at multiple points for each image.

The images are presented with antigens represented in separate colors.

Supplemental material includes additional data to supplement the figures in the main text, a table of data on the macaques used in this study, a table of flow cytometry reagents, and several tables of primary data from the scRNAseq analysis of granulomas.

All original data from this study are available at https://fairdomhub.org/studies/1134. Code for the scRNAseq data is available at https://github.com/ShalekLab/mTB_CD8_Depletion.

We acknowledge the hard work and dedication put forth by all technical and veterinary staff of the Flynn, Lin, Shalek, and Fortune Laboratories. In addition, we acknowledge that the members of these labs have contributed to this study through scientific discussions. We also thank the Unified Flow Core at the University of Pittsburgh. Antibody reagents for in vivo administration were provided by the NIH Nonhuman Primate Reagent Resource (P40 OD028116, U24 AI126683), and we are grateful to the staff there for guidance in using these antibodies.

This study was funded by NIH National Institute of Allergy and Infectious Diseases (NIAID) grant R56 AI139053 (P.L. Lin and J.L. Flynn); NIH National Heart, Lung, and Blood Institute grant K12 HL143886 (C.G. Winchell); NIH NIAID grant 75N93019C00071 IMPAcTB (J.L. Flynn, S.M. Fortune, and A.K. Shalek); the Bill and Melinda Gates Foundation (A.K. Shalek and J.L. Flynn); Harvard University Center for AIDS Research grant P30 AI060354 (J.M. Rosenberg); and The Harvard Clinical and Translational Science Center grant KL2 TR002542 (J.M. Rosenberg).

Author contributions: C.G. Winchell: conceptualization, data curation, formal analysis, funding acquistion, investigation, methodology, resources, visualization, writing—original draft. S.K. Nyquist: data curation, formal analysis, investigation, methodology, software, visualization, writing—original draft, writing—review and editing. M.C. Chao: formal analysis, software, visualization, writing—review and editing. P. Maiello: data curation, formal analysis, visualization. A. J. Myers: data curation, investigation. F. Hopkins: data curation, formal analysis, investigation, methodology, resources, software, validation, visualization. M. Chase: software. H.P. Gideon: data curation, investigation, methodology. K.V. Patel: investigation. J.D. Bromley: formal analysis. A.W. Simonson: formal analysis, investigation, visualization. R. Floyd-O’Sullivan: investigation. M. Wadsworth: data curation, methodology, resources. J.M. Rosenberg: conceptualization, formal analysis, investigation. R. Uddin: data curation. T. Hughes: data curation, formal analysis, investigation, methodology. R.J. Kelly: investigation. J. Griffo: investigation. J. Tomko: resources. E. Klein: investigation. B. Berger: formal analysis, funding acquistion, methodology, supervision, writing—review and editing. C.A. Scanga: project administration, supervision. J. Mattila: conceptualization, investigation, methodology, visualization, writing—original draft, writing—review and editing, S.M. Fortune: conceptualization, funding acquistion, supervision, writing—original draft, writing—review and editing. A.K. Shalek: conceptualization, funding acquisition, project administration, resources, supervision, visualization, writing—original draft, writing—review and editing. P.L. Lin: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, writing—original draft, writing—review and editing. J.L. Flynn: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, visualization, writing—original draft, writing—review and editing.