Single-cell analyses identify dysfunctional CD16+ CD8 T cells in smokers.

Tobacco smoke exposure contributes to the global burden of communicable and chronic diseases. To identify immune cells affected by smoking, we use single-cell RNA sequencing on peripheral blood from smokers and nonsmokers. Transcriptomes reveal a subpopulation of FCGR3A (CD16)-expressing Natural Killer (NK)-like CD8 T lymphocytes that increase in smokers. Mass cytometry confirms elevated CD16+ CD8 T cells in smokers. Inferred as highly differentiated by pseudotime analysis, NK-like CD8 T cells express markers characteristic of effector memory re-expressing CD45RA T (TEMRA) cells. Indicative of immune aging, smokers' CD8 T cells are biased toward differentiated cells and smokers have fewer naïve cells than nonsmokers. DNA methylation-based models show that smoking dose is associated with accelerated aging and decreased telomere length, a biomarker of T cell senescence. Immune aging accompanies T cell senescence, which can ultimately lead to impaired immune function. This suggests a role for smoking-induced, senescence-associated immune dysregulation in smoking-mediated pathologies.

Introduction

As a risk factor for human diseases, the global disease burden attributed to tobacco smoke exposure is substantial. The World Health Organization (WHO) estimates ∼6 million deaths per year from tobacco smoke exposure, resulting from both chronic and communicable diseases., In smokers, a decline in immunity and an increased risk of inflammatory diseases, such as atherosclerosis, make the argument that smoking-associated diseases are mediated by immune dysfunction. The development and progression of atherosclerotic lesions serves as an example of a complex immune-mediated pathology because T cells, monocytes, macrophages, dendritic cells (DCs), and B cells have been reported to be involved., Refining smoking-associated changes within immune populations will enhance our understanding of how dysfunctional immune subsets arise from exposure to tobacco smoke. This will facilitate the prevention of diseases by identifying immune cells to target for clinical intervention.

Smoking alters the epigenome and transcriptome of human blood leukocytes, in addition to DNA damage.5, 6, 7 In the article by Su et al., we demonstrated that changes identified in isolated cell fractions, which correspond to major immune populations, were distinct from one another and whole blood. For example, ITGAL, expressed in T cells and involved in inflammation, had decreased methylation in smokers’ T cells but not in whole blood or isolated cell fractions. It follows that bulk data from isolated fractions, comprising multiple subtypes, would similarly mask meaningful changes, especially when differences arise in low-frequency subsets. As such, the interpretation of bulk genomic approaches is limited because changes could indicate an altered distribution of cell (sub)populations or changes in expression within (sub)populations. The recent development of single-cell methods provides the technology to resolve smoking-associated (sub)population composition changes, examine gene expression differences, and identify rare subtypes obscured by bulk fraction data. In addition, multiparameter data allow us to concordantly study multiple cell types from the same individuals.

To identify smoking-affected cell (sub)populations and connect observed immune cell changes with smoking-associated diseases, we characterized gene expression profiles and cell surface marker phenotypes from primary peripheral blood mononuclear cells (PBMCs) from four nonsmokers and four smokers by single-cell RNA sequencing (scRNA-seq) and mass cytometry. The combination of transcriptome profiling and immunophenotyping provides higher confidence in the validity of our findings than one single-cell method alone. Major population frequencies showed a strong correlation between scRNA-seq and mass cytometry. We used single-cell transcriptome profiling to further separate cell populations into multiple subsets according to differentiation, activation, or functional states. We found a population of CD16+ CD8 T cells that was increased in smokers and exhibited natural killer (NK)-like transcriptional programs. Pseudotime analysis and examination of canonical markers revealed that these NK-like CD8 T cells likely represent a terminally differentiated state. DNA methylation models demonstrated that the smoking dose was associated with accelerated immune aging and decreased telomere length (TL; a biomarker of T cell senescence) in CD8 T cells. Not limited to CD8 T cells, smokers’ other immune populations displayed senescent characteristics.

By detecting an altered abundance of a rare population, we revealed an immune target that can be isolated and explored for connections between smoking and chronic diseases. Combined with increased (pre-)senescent CD8 T cells, elevated regulatory T cells (Tregs) and induction of senescence-linked genes in multiple cell types provide evidence that smokers show signs of premature immune system aging. The potential immune function defects and inflammatory subsets demonstrated here mirror the characteristics of pathologies commonly found in smokers. Further studies of smoking-associated dysregulation of immune transcriptional programs and candidate dysfunctional T cells linked to accelerated immune system aging will lead to mechanistic insights to advance disease prevention strategies for smoking-mediated pathologies.

Results

scRNA-Seq and Mass Cytometry Profiling of Human Peripheral Blood Immune Cells in Smokers and Nonsmokers

We set out to characterize the effects of cigarette smoke on immune cells in peripheral blood using single-cell approaches to determine whether the smoking-associated gene expression changes observed within major immune cell populations resulted from an altered abundance of specific, identifiable cell subsets. We performed scRNA-seq and mass cytometry in parallel on cryopreserved peripheral blood samples from eight donors with no history of atherosclerosis, chronic obstructive pulmonary disease (COPD), or lung cancer (Figure 1A). Serum cotinine, a metabolite of nicotine and an established biomarker of recent cigarette smoke exposure, confirmed the smoking status of donors. Smokers (n = 4) used for single-cell analyses had serum cotinine levels ranging from 240 to 511 ng/mL; all nonsmokers (n = 4) had serum cotinine levels <2 ng/mL. Donors were matched based on gender and race; ages ranged from 31 to 56 years and were not significantly different between smokers and nonsmokers (p = 0.23). Demographic and smoking information is listed in Table S1. We obtained single-cell mRNA data from 45,049 cells and surface protein expression data from 990,748 cells.

Figure 1

scRNA-Seq and Mass Cytometry Profiling of Human PBMCs from Smokers (n = 4) and Nonsmokers (n = 4)

(A) Experimental design: cryopreserved PBMCs from smokers and nonsmokers were thawed for scRNA-seq and mass cytometry.

(B) Uniform Manifold Approximation and Projection (UMAP) scRNA-seq plot: ∼45,000 PBMCs colored by major cell type.

(C) Force directed layout (FDL) of mass cytometry: ∼1 million cells colored by major cell type.

(D) Canonical gene expression markers for major cell types: CD8T cells (CD3D, CD8A), CD4T cells (CD3D, CD4) natural killer T (NKT) cells (CD3D, NCR3), natural killer (NK) cells (NCR3), monocytes (CD14, FCGR3A [CD16]), dendritic cells (DCs; FCER1A), and B cells (MS4A1).

(E) Cell surface protein expression for major cell types: CD8T cells (CD3, CD8a), CD4T cells (CD3, CD4), NKT cells (CD3, CD56), NK cells (CD56), monocytes (CD14, CD16), DCs (CD123) and B cells (CD19).

(F and G) Major cell type frequency distributions by individual donor (nonsmokers, NS; smokers, SM) colored by cell type for scRNA-seq (F) and mass cytometry (G).

(H) Major cell type frequencies showed strong correlation between scRNA-seq and mass cytometry (Pearson r = 0.99, r2 = 0.98, p < 0.0001). Shapes represent matched individuals, NS (unfilled) and SM (filled), are colored by cell type.

For each single-cell approach, we assigned cells to common immune populations based on mRNA (scRNA-seq) or surface protein (mass cytometry) expression (Figures 1B–1E; Table S2). For scRNA-seq, we used Seurat, to integrate data and implement shared nearest neighbors (SNN) clustering (see STAR Methods). We then used Model-based Analysis of Single-cell Transcriptomics (MAST) to identify positive and negative marker genes for each cluster and combined cells into major immune populations (Table S2). Cells in clusters expressing CD3D as a positive marker were designated as T cells (Figure 1D). T cells were further classified into CD4 T cells, CD8 T cells, or NKT cells based on the expression of CD4, CD8A, or NCR3 (Figure 1D). NK cells were identified based on CD3D as a negative marker combined with the expression of NKG7, GNLY, GZMB, PRF1, and NCR3 as positive markers (Figure 1D; Table S2). Monocytes were positive for LYZ and either CD14 or FCGR3A (encodes CD16), characteristic of classical or nonclassical monocytes (Figure 1D). DCs were similar to monocytes but could be distinguished by the expression of FCER1A (Figure 1D). B cells were defined by MS4A1 (encodes CD20; Figure 1D).

In parallel, PBMCs from each donor were assessed by mass cytometry (see STAR Methods). Viable, single-cell events were manually gated using Cytobank (Figure S1A). We used VorteX to cluster and create a force-directed layout (FDL) graph using the X-shift algorithm (see STAR Methods). A total of 122 PBMC cell clusters were determined from the 8 donors representing 983,848 cells (Figure S1B) and shown by smoking status (Figure S1C). Cell surface protein expression profiles were used to classify the cell populations (Figures 1C and 1E). T cells displayed CD3 and were classified by CD4 and CD8 as double-negative (DNT), double-positive (DPT), CD4 T, or CD8 T cells (Figures 1C and 1E). NKT cells were identified by CD3 and CD56 with CD4 or CD8 protein expression markers. Monocytes expressed CD14 and/or CD16 and DCs had CD123 above background levels (Figure 1E). B cells were positive for CD19 (Figure 1E). NK cells were positive for CD56 but negative for CD3 (Figure 1E).

To determine how well the scRNA-seq and mass cytometry corresponded with each other, we examined the individual donor proportion for each cell type. Cells colored by individual donors are shown for scRNA-seq (Figure S1D) and mass cytometry (Figure S1E). Cell-type frequencies were calculated and plotted to compare frequency distributions among individuals (Figures 1F and 1G). For both methods, all of the major populations—CD4T, CD8T, NKT, B, monocyte, and DC—were identified in all of the donors. We then compared the frequency of major populations in PBMCs by smoking status for scRNA-seq and mass cytometry using a Mann-Whitney U test. We observed no differences in the overall frequency of major cell types between smokers and nonsmokers by either scRNA-seq or mass cytometry (Figures S1F and S1G). Comparison of the major immune population percentages among individuals for scRNA-seq and mass cytometry showed significant strong correlation (Pearson r = 0.99, r2 = 0.98, p < 0.0001) between methods (Figure 1H).

scRNA-seq Reveals Increased Tregs and Altered Composition of the CD8 T Cell Population between Smokers and Nonsmokers

Single-cell transcriptome profiling can be used to separate cell populations into multiple subsets according to differentiation, activation, or functional states. Based on gene expression patterns, we clustered peripheral blood cells into 31 immune cell clusters and 1 erythroid contaminant cluster, labeled based on abundance from 0 (most abundant) through 31 (Figure 2A; Table S2). We identified 12 CD4 T cell (0, 1, 3, 4, 6, 10, 13, 14, 17, 18, 20, and 27), 7 CD8 T cell (2, 8, 11, 15, 19, 21, and 24), 3 NK cell (7, 22, and 25), 4 monocyte (5, 23, 26, and 30), and 2 B cell (9 and 12) clusters. NKTs (16), DCs (28), and MKs (31) were each contained by a single cluster. Here, clusters are referred to by major immune cell type, followed by original cluster identification (ID) (e.g., CD4T-0). To determine whether smoking altered the subtype distributions within the major cell populations, we compared the abundance of cells among clusters for each major cell type that separated into more than one cluster. Cells colored by smoking status are shown in Figure 2B. We did not observe any subset frequency shifts in B cells, monocytes, or NK cells (Figures S2A–S2F).

Figure 2

Increased CD4 Tregs and Altered CD8 T Cell Composition Observed between Smokers (SM, n = 4) and Nonsmokers (NS, n = 4)

(A) scRNA-seq UMAP colored by cluster ID.

(B) scRNA-seq UMAP colored by smoking status. Clusters CD8T-24 and CD8T-2 had more cells from nonsmokers and CD4T-17, CD8T-15, and CD8T-8 had more cells from smokers.

(C and D) Subset frequencies within CD4 T cells. Smokers (filled) had a significant increase in cluster CD4T-17 compared to nonsmokers (unfilled). Bar = median, ∗p < 0.05 by Mann-Whitney U test (C). Individual donor distribution of CD4 T cell subsets (D).

(E) FOXP3 gene expression within CD4 T subsets.

(F and G) Subset frequencies within CD8 T cells. Individual donor distribution of CD8 T cell subsets (F). Smokers (filled) had significant decreases in clusters CD8T-24 and CD8T-2 and increases in clusters CD8T-8 and CD8T-15 compared to nonsmokers (unfilled). Bar = median, ∗p < 0.05 by Mann-Whitney U test (G).

For 11 of 12 CD4 T cell subsets, frequencies were not significantly altered by smoking (Figures 2B and 2C). Although donors exhibited interindividual variation, smoking status did not appear to have a considerable impact on the distribution of CD4 T cells (Figure 2D). Only one cluster, CD4T-17, was higher in smokers than in nonsmokers (p < 0.05; Figure 2C). This cluster was relatively low in abundance among CD4 T cell subsets: a median 3.5% in nonsmokers and 5.4% in smokers. We characterized CD4T-17 cells as Tregs based on elevated FOXP3 and IL2RA (encodes CD25; Figure 2E; Table S2). No other CD4T clusters showed FOXP3 or IL2RA as strong positive markers. DNA methylation level of the AHRR gene can be used as a dose- and duration-dependent biomarker of tobacco smoke exposure; AHRR methylation is strongly negatively correlated with smoking dose and duration.,, The elevated frequency of Tregs was associated with reduced AHRR DNA methylation in whole blood (p = 0.002, r2 = 0.81; Figure S2G).

In contrast to most CD4 T cell subsets, the variation among donors in CD8 T cells depended on smoking status, as illustrated by the differences in the dominant population(s) in CD8 T cells from smokers and nonsmokers (Figure 2F). Smokers had lower proportions of 2 CD8 T cell clusters (CD8T-24 and CD8T-2) and higher proportions of 2 CD8 T cell clusters (CD8T-8 and CD8T-15) compared to nonsmokers (p < 0.05; Figure 2G). We did not observe significant differences in the proportions of the remaining 3 CD8 T cell clusters (CD8T-19, CD8T-11, and CD8T-21).

CD8 T Cell Distribution Shifts from Naive to Differentiated States in Smokers

Smoking had substantial effects on the composition of CD8 T cells. We further analyzed each of the seven CD8 T clusters to identify distinct gene expression patterns. MAST identified positive and negative markers for each CD8 T cluster relative to other CD8 T cells (see STAR Methods). We found 71, 223, 46, 144, 71, 739, and 105 positive and 86, 978, 10, 68, 59, 214, and 52 negative markers (Bonferroni adjusted p < 0.05 in both smokers’ and nonsmokers’ cells) for CD8T-24, CD8T-2, CD8T-19, CD8T-11, CD8T-21, CD8T-8, and CD8T-15 clusters, respectively (Table S3). We examined marker lists for genes associated with T cell differentiation and function to distinguish between CD8 T subsets. Genes frequently used to classify CD8 T subsets varied among CD8 T cell clusters (Figures 3A and 3B). Elevated CCR7, SELL, and IL7R combined with low CCL5 indicated that clusters CD8T-24 and CD8T-2 represented naive cells (TN). High levels of IL7R, SELL, and FOS (activated T cell proliferation,) suggested that cluster CD8T-15 exhibited characteristics of long-lived memory cells, such as central memory T cells (TCM). Reduced levels of CCR7, SELL, and IL7R and elevated CCL5 and KLRG1 in CD8T-11, CD8T-21, and CD8T-8 are indicative of later differentiation stages (e.g., TEM). The lack of CD27 expression and decrease in FOS in CD8T-21 and CD8T-8 indicated highly differentiated TEM cells (i.e., TEMRA). ZEB2, associated with terminal differentiation states,, was detected in approximately one-third of CD8T-8 cells (35.1% in smokers and 28.2% in nonsmokers).

Figure 3

CD8 T Cell Subsets Shift from Naive to Differentiated CD8 T Cell States

(A) scRNA-seq UMAP (n = 8) of the 7 CD8 T cell subsets colored by cluster ID.

(B) Expression of genes to characterize CD8 T subsets. Color indicates cluster ID (n = 8).

(C) Trajectory inference (pseudotime analysis) ordered CD8 T cell clusters into 2 lineages, which originate from cluster CD8T-24 and terminate at either cluster CD8T-8 (lineage 1) or CD8T-15 (lineage 2). The black line shows lineages and the gray line shows simultaneous principle curves (n = 8).

(D) Pseudotime for CD8 T cells by individual nonsmokers (NS, n = 4) and smokers (SM, n = 4) in lineage 1 (left) and lineage 2 (right).

(E) Gene expression data from flow-sorted human peripheral CD8 T cells (n = 6): naive (TN; CD27+CD45RA+), central memory (TCM; CD27+CD45RA−) effector memory (TEM; CD27−CD45RA−), and EMRA (TEMRA; CD27−CD45RA+). FCGR3A (false discovery rate [FDR] = 1.38 × 10−8), PRF1 (FDR = 2.69 × 10−6), and ITGAL (FDR = 5.43 × 10−3), are positive markers for CD8T-8 cells; B3GAT1 (CD57; FDR = 1.06 × 10−5) and ZEB2 (FDR = 2.16 × 10−5) are positive markers for CD8T-8 cells and indicators of senescence; PRSS23 (FDR = 8.62 × 10−4), TNFSF10 (FDR = 6.46 × 10−3), and CTSW (FDR = 2.52 × 10−3) are secreted factors associated with SASP; and CX3CR1 (FDR = 8.95 × 10−6) is a chemokine receptor.

(F and G) Frequency of FCGR3A CD8 TEMRA (F) and CD8 TN (G) cells are correlated with smoking dose as measured by AHRR smoking biomarker (n = 8).

Several clusters exhibited intermediate expression of differentiation state genes. To organize CD8 T cell clusters by their likely differentiation trajectories, we used the Slingshot algorithm to perform pseudotemporal analysis. Lineage inference ordered CD8 T cells into 2 lineages, which originate from cluster CD8T-24 and terminate at either cluster CD8T-8 or CD8T-15 (Figure 3C). Lineage 1 mostly comprised cells from CD8T-24, CD8T-2, CD8T-19, CD8T-8, CD8T-21, and CD8T-8, with minimal cells from CD8T-15. Lineage 2 mostly comprised cells from CD8T-24, CD8T-2, and CD8T-15, with minimal cells from CD8T-19. Based on the altered composition of CD8 T cell subsets—lower proportions of CD8T-2 and CD8T-24 and higher proportions of CD8T-8 and CD8T-15 cells (Figures 2B, 2F, and 2G)—we propose that tobacco smoke exposure alters CD8 T cell composition by shifting CD8 T cells toward differentiated states. Smokers’ cells were biased toward later pseudotimes in both lineages (Figure 3D), demonstrating that smokers’ CD8 T cells are skewed toward differentiated states and nonsmokers’ CD8 T cells are skewed toward naive states.

We next identified temporally associated genes for each lineage. Only 10 genes within the top 100 temporally expressed genes were shared, and most shared genes exhibited the same direction of change over pseudotime in both lineages (Figures S3A and S3B). In general, temporally expressed genes for CD8 T cell lineages were consistent with effector memory (lineage 1) and central memory (lineage 2) differentiation. For example, in lineage 1, CCL5 and NKG7 increased, while SELL and IL7R decreased over the differentiation trajectory (Figures S3A and S3C). CMC1 demonstrated a nonlinear association in lineage 1, as it peaked in CD8T-21 cells and then decreased through CD8T-8 cells (Figures S3A and S3C).

The terminal cluster in the effector memory trajectory, CD8T-8, shared many features with the penultimate cluster, CD8T-21; however, CD8T-8 increased in smokers, but CD8T-21 did not. Since the low expression of CD27 and CCR7 and elevated expression of KLRG1 in both clusters would classify these cells as highly differentiated CD8 T cells (i.e., TEMRA-like), we sought to find markers within CD8T-8 that did not occur in CD8T-21. Examination of exhaustion markers—TOX, PDCD1 (encodes programmed cell death protein-1 [PD-1]), CTLA4, and HAVCR2 (encodes TIM3)—did not distinguish the clusters. Whereas smokers’ CD8T-8 cells had elevated expression of TOX (log2FC [fold change] = 0.31, p = 1.44 × 10−9), nonsmokers’ CD8T-8 cells did not, and no other exhaustion markers were significantly elevated (Table S3). We next examined senescence-associated genes KLRG1 and B3GAT1 (encodes CD57). While KRLG1 was a positive marker for CD8T-21 (smokers [SM]: log2FC = 0.65, p = 8.11 × 10−04; nonsmokers [NS]: log2FC = 1.33, p = 2.23 × 10−22) and CD8T-8 cells (SM: log2FC = 0.94, p = 6.29 × 10−47; NS: log2FC = 1.28, p = 1.94 × 10−40), B3GAT1 was unique to CD8T-8 cells (SM: log2FC = 0.19, p = 3.83 × 10−14; NS: log2FC = 0.10, p = 0.022). We also found that 2 genes reported as having smoking-associated methylation changes, GFI1 and PRSS23, showed elevated expression in CD8T-8 cells (GFI1 SM: log2FC = 0.20, p = 5.58 × 10−7; GFI1 NS: log2FC = 0.22, p = 0.046; PRSS23 SM: log2FC = 1.07, p = 1.5 × 10−143; PRSS23 NS: log2FC = 0.91, p = 7.7 × 10−60). FCGR3A, commonly found on NK cells and nonclassical monocytes, was identified as a strong positive marker of CD8T-8 cells in smokers (log2FC = 1.63, p = 1.17 × 10−171) and nonsmokers (log2FC = 1.72, p = 5.17 × 10−119; Figure 1D; Table S3).

We reanalyzed the gene expression data from flow-sorted human peripheral naive (CD27+CD45RA+), central memory (CD27+CD45RA−), effector memory (CD27−CD45RA−), and EMRA (CD27−CD45RA+) CD8 T cells from healthy adults. We found that the TEMRA subset had an increased expression of FCGR3A (Figure 3E). Combined with our finding of high FCGR3A in our most differentiated CD8 T cluster (CD8T-8), data from Callender et al. supports elevated FCGR3A as a characteristic of CD8 TEMRA cells. In addition, CD8 TEMRA cells showed an elevated expression of senescence-associated secretory phenotype (SASP) genes (Figure 3E).

Given that human cytomegalovirus (HCMV) infection can alter T cell composition,, we tested donors for HCMV antibodies and DNA in serum (see STAR Methods). Compared to seronegative donors (immunoglobulin G [IgG]−/IgM−), seropositive donors (IgG+/IgM−) had elevated levels of CD8 TEMRA cells (CD8T-21), but not FCGR3A-expressing CD8 TEMRA cells (CD8T-8); neither naive CD8 T clusters (CD8T-24 and CD8T-2) were significantly reduced in HCMV seropositive donors (Figure S3D; Table S1).

To further support the idea that smoking elevated FCGR3A-expressing CD8 TEMRA cells (CD8T-8) and reduced naive CD8 T cells (CD8T-24 and CD8T-2), we compared the frequencies of these subpopulations against whole-blood DNA methylation of AHRR (Figures 3F and 3G; Table S1). A higher frequency of FCGR3A-expressing CD8 TEMRA cells was associated with reduced AHRR methylation (p = 0.008, r2 = 0.72), and a lower frequency of naive CD8 T cells was associated with reduced AHRR methylation (p = 0.002, r2 = 0.82). In addition, the increase in FCGR3A-expressing CD8 TEMRA cells was highly correlated with the increase in CD4 Tregs (p = 0.019, r2 = 0.63; Figure S3E).

Mass Cytometry Confirms Elevated CD16+ CD8 T Cells in Smokers

Relatively rare in nonsmokers (median: 1.8%), the FCGR3A-expressing CD8 T cell cluster (CD8T-8) comprised 7.3% of PBMCs in smokers. Reported as a low-frequency subset (∼2% of PBMCs in healthy adults), CD16+ CD8 T cells have been described.25, 26, 27 Based on these reports, we sought to ascertain whether smokers had increased levels of CD16+ CD8 T cells that expressed surface proteins for CD3, but not CD56; i.e., confirm an increase in CD16 expression within CD8 T cells that are not NKT cells. To show that smokers had increased surface protein expression of CD16 within their CD8 T cells, we ran X-shift on all CD8 T and NKT cells and visualized the FDL colored by cluster IDs (Figure 4A). We did not see any differences between smokers’ and nonsmokers’ CD8 T cells for CD56 (Figure 4B), but smokers showed an increase in the proportion of CD8 T cells expressing CD16 compared to nonsmokers (Figure 4C). We determined the frequency of CD16+ CD8 T cells in smokers and nonsmokers by manual gating (Figure S4A). CD3+CD56− negative cells were then gated to obtain single positive CD8 T cells (CD3+CD56−CD8+CD4−; Figure S4A), which were then used to determine the frequency of CD16+ CD8 T cells (CD3+CD56−CD8+CD4−CD16+; Figure S4A). Compared to nonsmokers, smokers had a significantly increased frequency of CD16+ CD8 T cells (p = 0.03; Figure 4D), confirming that smokers had elevated proportions of CD16+ CD8 T cells. In an independent group of donors, CD16+ CD8 T cells were also elevated in smokers compared to nonsmokers (p < 0.01; Figure S4B; Table S1).

Figure 4

Characterization of CD16+ CD8 T Cells

(A–D) Mass cytometry confirms elevated proportion of CD16+ CD8 T cells in smokers.

(A) FDL (n = 8) of CD8 T and NKT cells colored by cluster ID.

(B and C) Cell surface marker intensity for NK marker CD56 (B) and CD16 (C) in nonsmokers (top, n = 4) and smokers (bottom, n = 4).

(D) Frequency of CD16+ cells increases within CD8 T cells in smokers (filled, n = 4) compared to nonsmokers (unfilled, n = 4). Bar = median, ∗p < 0.05 by 1-tailed Mann-Whitney U test.

(E and F) NK-associated transcriptome characteristics (n = 8) revealed in CD8T-8 cells by gene set enrichment analysis (GSEA). CD8T-8 cells were positively enriched (FWER < 0.05) for gene signatures of T cells induced to have NK-like phenotypes; split dot plot shows 25 genes within the induced T to NK gene set (E). CD8T-8 cells were negatively enriched (FWER < 0.05) for genes higher in naive CD8 T cells relative to NK cells; split dot plot shows 25 genes within the naive T versus NK gene set (F). Color intensity represents average per-cell expression (scale shows standard deviation), and circle size indicates the percentage of cells expressing each gene.

(G and H) CD8T-8 cells show increased GZMB expression in smokers. (G) UMAP comparison of GZMB between nonsmokers’ (top, n = 4) and smokers’ (bottom, n = 4) cells. Cluster CD8T-8 is indicated by dotted circle. (H) Individual donor GZMB levels for CD8T-8 cells for nonsmokers (NS, n = 4) and smokers (SM, n = 4).

(I and J) CD8T-8 cells have increased PRF1 expression in smokers. (I) UMAP comparison of PRF1 between nonsmokers’ (top, n = 4) and smokers’ (bottom, n = 4) cells. Cluster CD8T-8 is indicated by dotted circle. (J) Individual donor PRF1 levels for CD8T-8 cells for nonsmokers (NS, n = 4) and smokers (SM, n = 4).

To phenotype the CD16+ CD8 T cell subset, we gated CD3+ T cells by a CD45RA/CD45 biaxial plot to establish an accurate CD45RA+ gate that was then applied to the CD16+ CD8 T cells (Figure S4C). Consistent with a TEMRA phenotype, the majority of CD16+ CD8 T cells were positive for CD45RA in both smokers and nonsmokers (Figure S4D).

FCGR3A-Expressing CD8 T Cells Exhibit NK-like Transcriptome Signatures

After confirming an increase in CD16+ CD8 T cells in smokers, we further examined how the transcriptomes of FCGR3A-expressing CD8 T cells differed from other CD8 T cells. In addition to FCGR3A, CD8T-8 cells exhibited elevated expression NKG7, GNLY, FGFBP2, GZMB, and PRF1 (Tables S2 and S3). While these genes are considered cytotoxic T or NK cell-expression signatures, the presence of CD16 led us to suspect that this subset may express genes indicative of NK-like attributes. To gain insight into the functional relevance of gene expression profiles for CD8T-8 cells, we performed gene set enrichment analysis (GSEA). Consistent with an NK-like transcriptional program, “LI INDUCED T TO NATURAL KILLER UP” had a positive normalized enrichment score (NES = 1.88, FWER < 0.05) and “GSE22886 NAIVE TCELL VS NKCELL UP” had a negative NES (−2.13, FWER < 0.05). The “LI INDUCED T TO NATURAL KILLER UP” gene set encompasses expression patterns for T cells reprogrammed to have NK-like phenotypes—“induced T to NK” (iTNK) cells. Here, we found 67 and 71 genes from the iTNK gene signature as positive markers (p < 0.05) of FCGR3A-expressing CD8 T cells in smokers and nonsmokers, respectively. Figure 4E shows the 25 highest ranked iTNK genes, based on p value in CD8T-8 cells. CDT8-8 genes were negatively enriched (FWER < 0.05) for genes that are higher in naive CD8 T cells relative to NK cells (“GSE22886 NAIVE TCELL VS NKCELL UP”). Smokers’ and nonsmokers’ CD8T-8 cells had significantly reduced expression (i.e., negative markers) of 55 and 45 naive CD8 T versus NK genes. Figure 4F shows the 25 highest ranked naive CD8 T versus NK genes, based on p value for CD8T-8 cells.

To examine whether smokers’ “NK-like” CD8 T cells differed from those of nonsmokers, we compared the average expression of genes from smokers’ CD8T-8 cells to nonsmokers’ CD8T-8 cells (see STAR Methods). We found that 63 genes had an increased and 74 genes had a decreased average per cell expression in smokers compared to nonsmokers (p < 0.05; Table S4). Figures S4E and S4F show 25 genes with increased and decreased per cell expression, ordered by the difference in the percentage of cells expressing each gene between smokers and nonsmokers. Although cellular mRNA levels for effector molecules granzyme B (encoded by GZMB) and perforin (encoded by PRF1) exhibited interindividual variation, both increased in frequency of expression and average per-cell expression in smokers compared to nonsmokers (Figures 4G–4J, S4E, S4G, and S4H).

CD8 T Bulk Transcriptomes Reflect Differentiation State Shifts Observed at the Single-Cell Level

To assess the overall impact of smoking on CD8 T cells, we identified differentially expressed genes (DEGs) between smokers and nonsmokers by comparing the average per-cell expression for all of the cells in the 7 CD8 T clusters combined. Of 2,163 genes evaluated in the pseudobulk analysis, we found that 1,817 genes had higher expression and 344 genes had lower expression in smokers versus nonsmokers (q < 0.05; Table S5). To examine the interindividual variability in response to smoking, we performed hierarchical clustering using smoking scRNA-DEGs, which separated individual donors by smoking status (Figures 5A and S5A). To validate altered CD8 T gene expression profiles in smokers, we used RNA-seq and microarray on isolated CD8 T cells to examine the differences in bulk RNA expression. Isolated CD8 T cells for bulk RNA-seq included the 8 donors used in scRNA-seq (bulk RNA was from a previous visit) and 7 additional donors (Table S1). We identified 1,268 genes as differentially expressed—692 increased and 576 decreased (q < 0.05; Figure 5B). With the exception of F061, principal-component analysis (PCA) of bulk RNA-seq data separated smokers from nonsmokers (Figure 5C). We evaluated microarray data from isolated CD8 T cells from 19 donors (9 smokers and 10 nonsmokers). Isolated CD8T cells included 4 donors used in scRNA-seq, 4 donors used in bulk RNA-seq, and 11 additional donors (Table S1). We identified 51 DEGs (see STAR Methods)—46 increased and 5 decreased (Figure 5D).

Figure 5

CD8 T Bulk Transcriptomes Reflect Differentiation State Shifts at Single-Cell Level and DNA Methylomes Demonstrate Accelerated Aging and Shortened Telomeres

(A) Top 350 upregulated and downregulated differentially expressed genes (DEGs) between smokers (SM, n = 4) and nonsmokers (NS, n = 4) from CD8 T cells (scRNA-seq). Donors were separated by smoking status using smoking scRNA-seq-DEGs for hierarchical clustering. Labeled genes were also significantly upregulated in the CD8 T cell microarray results.

(B) Comparison of bulk RNA-seq expression of isolated CD8 T cells between smokers (n = 7) and nonsmokers (n = 8). DEGs with higher expression in smokers are labeled in red and DEGs with lower expression in smokers are denoted in blue. The gene names in black were altered in bulk RNA-seq and scRNA-seq. The gene names in gray were altered in bulk RNA-seq and microarray. The gene names in boldface were altered in microarray, bulk RNA-seq, and scRNA-seq.

(C) Principal-component analysis (PCA) of bulk RNA-seq. Donors used in both scRNA-seq and bulk RNA-seq are labeled with “F” and “V” patient codes. See Table S1 for donor information.

(D) Comparison of microarray expression of isolated CD8 T cells between smokers (n = 9) and nonsmokers (n = 10). The genes with higher expression in smokers are red and the genes with lower expression in smokers are blue. The gene names in black were altered in microarray and scRNA-seq. The gene names in gray were altered in microarray and bulk RNA-seq. The gene names in boldface were altered in microarray, bulk RNA-seq, and scRNA-seq.

(E) Immunological signatures and chemical and genetic perturbations that were significantly enriched (FWER < 0.05) by GSEA in scRNA-seq, bulk RNA-seq, and microarray.

(F–I) In CD8 T cells, age acceleration (F, H) and shortened telomeres (G, I) are correlated with the smoking dose (AHRR DNA methylation) in scRNA-seq donors (F–G) and a larger group (H–I) of nonsmokers (n = 71, open blue circles) and smokers (n = 60, closed red circles).

We compared results from the 3 RNA analysis platforms used to identify smoking-associated DEGs in CD8 T cells. Bulk RNA-seq confirmed 241 smoking DEGs from the scRNA-seq, 162 with increased and 79 with decreased expression in smokers compared to nonsmokers (Figures 5B; Table S5). Microarray analysis confirmed 24 smoking DEGs that were identified by scRNA-seq as increased in smokers compared to nonsmokers (Figures 5A and 5D; Table S5). Two genes, FAM129A and CD58, increased in smokers’ CD8 T cells in all 3 platforms (Figures 5A, 5B, 5D, S5A, and S5B; Table S5). Genes found to be altered by at least 2 methods include LGALS1, ADAM8, and CLDND1, which were increased in bulk RNA-seq and scRNA-seq data; GPR15, which was increased in the bulk RNA-seq and microarray data; and NDFIP1, which was decreased in the bulk RNA-seq and scRNA-seq data (Figures 5A, 5B, 5D, and S5A–S5F). ITGAL, a smoking methylation biomarker, was only found to be significantly increased by scRNA-seq (log2FC = 0.36, p = 1.70 × 10−26), and was also elevated in the NK-like subset (SM: log2FC = 0.64, p = 6.53 × 10−29; NS: log2FC = 0.78, p = 1.50 × 10−20).

Although each platform identified smoking DEGs not found by other platforms, we expect the overall changes observed to represent a similar shift in the functional states of CD8 T cells. We used GSEA to determine that CD8 T pseudobulk and bulk transcriptomes were enriched for similar functional annotations. Pseudobulk scRNA-seq, bulk RNA-seq, and bulk microarray shared 7 positively and 8 negatively enriched gene sets (FWER < 0.05; Figure 5E). CD8 T cells were positively enriched for genes with a higher expression in memory T, TCM, TEM, PD1lo (CD8 TEM), and PD1hi (CD8 TEM) cells relative to naive T cells (Figure 5E). CD8 T cells were negatively enriched for genes that have a higher expression in naive T cells relative to memory T, TCM, TEM, PD1lo (CD8 TEM), and PD1hi (CD8 TEM) cells (Figure 5E). Therefore, GSEA of CD8 T smoking DEGs identified immunological signatures indicative of an increased expression of genes associated with effector memory and central memory functions and decreased expression of genes associated with naive T cells.

CD8 T DNA Methylomes Indicate Smoking Dose-Dependent Accelerated Aging and Decreased TL

In CD8 T cells, the shift from naive to differentiated states has been associated with aging and senescence. To determine whether tobacco smoke exposure was associated with immune aging and senescence, we analyzed DNA methylation data from bulk CD8 T cells (see STAR Methods). We used DNAm PhenoAge to calculate age acceleration (see STAR Methods). In the donors used for scRNA-seq and mass cytometry, we found a strong association between reduced AHRR DNA methylation and age acceleration (p = 0.005, r2 = 0.76; Figure 5F). Since shortened telomeres (TL) are indicative of T cell senescence, we used a DNA methylation estimator of TL (DNAmTL) to estimate the TL of CD8 T cells (see STAR Methods). We found a strong association between lower AHRR DNA methylation and shorter DNAmTL (p = 0.0004, r2 = 0.89; Figure 5G). To confirm these effects in a larger cohort, we analyzed CD8 T cell DNA methylation in a group of 131 individuals (Table S1). Age acceleration was highly significantly correlated with smoking dose (p < 0.00009), and ∼11% of the variance in age acceleration in this population could be attributed to smoking dose (Figure 5H). Similarly, methylation-derived CD8 T cell TL was highly significantly associated with smoking dose (p < 0.0002, r2 = 0.11; Figure 5I).

Smoking-Associated Gene Expression Changes in PBMC Populations

Since we did not observe substantial changes in subset distribution for the majority of PBMC populations, we compared average per-cell gene expression for all cells within each cell type (CD4 T, NKT, NK, monocyte, DC, and B) to identify smoking DEGs. For CD4 T cells, we found 1,563 DEGs—1,278 showed increased and 285 showed decreased expression (Table S6). Hierarchical clustering of donors by CD4 T smoking DEGs clustered individuals by smoking status (Figure 6A). Bulk gene expression of isolated CD4+ cells identified 2 upregulated (LRRN3 and GPR171) and 1 downregulated (APBA2) gene in common with the CD4 T pseudobulk analysis (Figures 6A and S6A). NKT cells had 89 smoking DEGs, 45 with increased and 44 with decreased expression, and NK cells had 238 smoking DEGs, 129 with increased and 109 with decreased expression (Table S6). Hierarchical clustering of DEGs separated donors by smoking status for NKT but not for NK cells (Figures 6B and 6C). CD56+ cells (bulk), which contain NKT and NK cells, shared 3 upregulated genes (PROK2, MX1, and TRAT1) and 1 downregulated gene (KLRB1) with NKT cells and 3 upregulated (MX1, DOCK5, and CSGALNACT1) and 3 downregulated genes (CD160, XCL2, and KLRB1) with NK cells (Figures 6B, 6C, and S6B). For monocytes, we found 488 DEGs between smokers and nonsmokers by scRNA pseudobulk analysis (Figure 6D; Table S6). Of the 290 DEGs with higher expression in smokers’ monocytes, 15 increased in smokers’ isolated CD14+ cells by microarray (Figures 6D and S6C). For DCs, we identified 21 smoking DEGs—5 showed increased and 16 showed decreased expression (Figure 6E; Table S6). In B cells, we found 190 DEGs, 111 with increased and 79 with decreased expression in smokers (Table S6). Using microarray in CD19+ cells, we confirmed the decreased expression of HLA-DQA1 in B cells from smokers compared to nonsmokers (Figures 6F and S6D). Table 1 lists the biological relevance of smoking DEGs found in PBMC populations by pseudobulk scRNA-seq and confirmed by microarray.

Figure 6

Major Cell Types Show Gene Expression Changes in Peripheral Blood of Smokers

(A) Hierarchical clustering of the top 350 increased and decreased DEGs between smokers (n = 4) and nonsmokers (n = 4) from CD4 T cells (scRNA-seq).

(B–F) Heatmaps of smoking DEGs (scRNA-seq) in NKT cells (B), NK cells (C), monocytes (D), DCs (E), and B cells (F).

Table 1

Smoking DEGs Altered in Both scRNA-Seq and Microarray

Gene Name	DEG		Biological Process/Function/Pathway
CD8 T
granzyme B	↑	GZMB	PD-1low44
natural killer cell granule protein 7	↑	NKG7	PD-1low, PD-1high44
perforin 1	↑	PRF1	PD-1low, PD-1high44
killer cell lectin like receptor D1	↑	KLRD1	PD-1low, PD-1high44
sterile α motif domain containing 3	↑	SAMD3	CD8 T tolerance58,59
G protein-coupled receptor 171	↑	GPR171	negative regulation of myeloid differentiation60
calpastatin	↑	CAST	PD-1low, PD-1high;44 effector memory, central memory61
chromosome 1 open reading frame 21	↑	C1orf21	cytotoxic T by scRNA-seq in Crohn's disease59
family with sequence similarity 129, member A	↑	FAM129A	PD-1low, PD-1high;44 memory, effector memory, central memory61
TNF superfamily member 10	↑	TNFSF10	marker for atherosclerosis plaque formation62
C-X3-C motif chemokine receptor 1	↑	CX3CR1	chemokine receptor associated with terminally differentiated effector CD8 cells63
pyrin and HIN domain family member 1	↑	PYHIN1	PD-1low, PD-1high,44 memory61
AT-rich interaction domain 5B	↑	ARID5B	central memory61
zinc finger E-box binding homeobox 2	↑	ZEB2a	PD-1low, PD-1high44
CD58 molecule	↑	CD58	PD-1low, PD-1high;44 memory, effector memory, central memory61
protein kinase cAMP-dependent type I regulatory subunit α	↑	PRKAR1A	enhanced cytoxicity of effector T cells64
acyloxyacyl hydrolase	↑	AOAH	effector memory phenotype, increase in chronic lymphocytic leukemia CD3 T cells65
proline rich 5-like	↑	PRR5L	PD-1low, PD-1high;44 effector memory, central memory61
MYB proto-oncogene like 1	↑	MYBL1	PD-1low, PD-1high;44 memory, effector memory, central memory61
ATPase plasma membrane Ca2+ transporting 4	↑	ATP2B4	PD-1low, PD-1high;44 memory, effector memory, central memory61
synaptotagmin 11	↑	SYT11	PD-1low, PD-1high44
killer cell Ig-like receptor, 3 Ig domains, and long cytoplasmic tail 2	↑	KIR3DL2a	T cell aging66
death domain containing 1	↑	DTHD1	unknown function
fibrinogen like 2	↑	FGL2a	replicative senescence,67 procoagulant68
CD4 T
leucine-rich repeat neuronal protein 3	↑	LRRN3a	smoking methylation,56 T cell replicative senescence69
G protein-coupled receptor 171	↑	GPR171	negative regulation of myeloid differentiation60
amyloid β precursor protein binding family A member 2	↓	APBA2	smoking methylation38
NKT
prokineticin 2	↑	PROK2	inflammatory response in smoking70 and diabetes71
MX dynamin like GTPase 1	↑	MX1a	TNF-α induced cellular senescence37
T cell receptor associated transmembrane adaptor 1	↑	TRAT1	leukocyte activation72
killer cell lectin like receptor B1	↓	KLRB1a	downregulated in human aged CD4+ memory T cells66
NK
MX dynamin like GTPase 1	↑	MX1a	TNF-α induced cellular senescence37
dedicator of cytokinesis 5	↑	DOCK5a	upregulated in CD57+ NK73
chondroitin sulfate N-acetylgalactosaminyltransferase 1	↑	CSGALNACT1	positively correlated with smoking74
CD160 molecule	↓	CD160a	increased in exhausted T cells75
X-C motif chemokine ligand 2	↓	XCL2	overexpressed in lung cancer tumors, increases with prognosis76
killer cell lectin like receptor B1	↓	KLRB1a	downregulated in human aged CD4+ memory T cells66
Monocytes
phospholipid scramblase 1	↑	PLSCR1	expression increased with cytokine treatment77
GTPase, IMAP family member 4	↑	GIMAP4	regulates INF-γ in CD4 T cell differentiation78
MX dynamin like GTPase 1	↑	MX1a	TNF-α-induced cellular senescence37
MX dynamin like GTPase 2	↑	MX2a	TNF-α-induced cellular senescence37
poly(ADP-ribose) polymerase family member 9	↑	PARP9	silences pro-inflammatory genes, found in atherosclerotic plaques41
cystatin C	↑	CST3a	atherosclerosis,38 cellular senescence39
sterile α motif domain containing 9	↑	SAMD9	anti-inflammatory factor79
TNF superfamily member 13b	↑	TNFSF13B	regulates proliferation and differentiation of atherogenic B cells40
Cysteine-rich secretory protein LCCL domain containing 2	↑	CRISPLD2	regulates anti-inflammatory effects80
erythrocyte membrane protein band 4.1-like 3	↑	EPB41L3	positively correlated smoking gene and associated with cancer74
heat shock protein family A (Hsp70) member 1B	↑	HSPA1B	overexpressed in advanced atherosclerosis42
thrombospondin 1	↑	THBS1	regulation of monocyte mobility, vascular inflammation81
IFN-induced protein with tetratricopeptide repeats 5	↑	IFIT5a	up in TNF-α-mediated cellular senescence37
formyl peptide receptor 2	↑	FPR2	increased in atherosclerotic lesions and plaque stability43
IFN induced with helicase C domain 1	↑	IFIH1a	up in TNF-α-mediated cellular senescence37
B
major histocompatibility complex, class II, DQ α 1	↓	HLA-DQA1	presents peptides from extracellular protein in T cells82

cAMP, cyclic AMP; CLL, XX; DEG, differentially expressed gene; GTPase, guanosine triphosphatase; INF, interferon; TNF, tumor necrosis factor.

Gene associated with cellular senescence.

Discussion

Our study reveals CD16+ CD8 T cells and other signs of immune cell dysfunction as elevated in smokers. These cells, uncovered by scRNA-seq and confirmed by mass cytometry in human PBMCs from multiple individuals, shared transcriptomic features with iTNK cells, which acquire NK surface receptors and have increased cytotoxic potency. Combined with CD16 and CD45RA protein expression, the transcriptome of the NK-like CD8 T cells implies an innate-like, terminally differentiated CD8 T subset. CD16, commonly associated with NK cells, binds IgG antibodies to mediate antibody-dependent cellular cytotoxicity (ADCC); exogenous or endogenous CD16 expression enables T cells to execute ADCC., Consistent with a heightened cytolytic potential, the NK-like CD8 T cells had elevated mRNA expression of cytolytic effector molecules GZMB and PRF1; these transcripts were also higher in smokers than nonsmokers within this subset. Granzyme B and perforin-expressing CD8 T cells contribute to the development of atherosclerotic plaques in mice., As such, our results highlight a potential link between smoking-induced functional changes in human CD8 T cells and atherosclerosis.

Applying scRNA-seq and mass cytometry to PBMCs from tobacco smoke-exposed individuals, we show that major immune populations can be discerned and disparate subsets can be identified among CD4 T cells, CD8 T cells, NK cells, monocytes, and B cells. Pseudobulk analysis of immune cell populations revealed smoking DEGs, several of which were confirmed in bulk cell-type fractions. The increase in smokers’ Tregs is likely masked in bulk data from isolated CD4 T cells because Tregs are a low-frequency subset. In a larger cohort of healthy women (N = 75), smoking was found to be a positive predictor of Treg levels in peripheral blood. Notably, Tregs have been shown to induce T cell senescence, highlighting a potential role for the increase in Tregs observed here. Senescence-related genes were altered in smokers in multiple cell types (Table 1). MX1, induced in tumor necrosis factor-α (TNF-α)-mediated senescence, increased in smokers in NKT cells, NK cells, and monocytes.

Alluding to the shared regulation of pro-senescent and pro-atherosclerotic signaling, TNF-α-induced senescence genes are enriched for atherosclerosis signaling genes. CST3, associated with subclinical atherosclerosis and cellular senescence, increased in smokers’ monocytes. We also identified TNFSF13B, a critical regulator of atherogenic B cell proliferation and differentiation. Other notable genes connected to atherosclerosis that increased in smokers’ monocytes include PARP9, HSPA1B, and FPR2.

Two approaches, established markers and trajectory inference, demonstrate that smokers lose naive and gain TCM-like and TEMRA-like cells. Differentiation state shifts in CD8 T cells were supported by bulk analysis methods. All 3 transcriptomic platforms identified changes associated with PD-1hi CD8 T cells. With persistent antigen stimulation, the inhibitory effect of the PD-1 pathway contributes to pathologies associated with T cell dysfunction during chronic viral infection and tumor evasion of host immune response. Duraiswamy et al. demonstrated that PD-1hi CD8 T cells obtained from healthy adults had gene expression profiles similar to those of PD-1lo CD8 T cells and did not show either exhausted gene signatures or phenotypes characteristic of PD-1hi cells obtained from humans or mice with chronic infections. Smoking DEGs found in CD8 T cells in pseudobulk scRNA-seq, bulk RNA-seq, or microarray were represented by these gene signatures (Figure 5E; Table 1). Notably, bulk RNA-seq data were acquired from prior visits from donors used for scRNA-seq, indicating a chronic or recurring state of activation in smokers’ CD8 T cells. Chronic activation suggests that CD8 T cells obtained from healthy smokers here may represent a dysfunctional phenotype.

Consistent with an end-stage TEMRA phenotype, NK-like CD8 T cells had the latest pseudotimes. GFI1, a transcriptional repressor of interleukin-7 receptor α (IL-7Rα) that drives the terminal differentiation of CD8 T cells, was elevated in NK-like CD8 TEMRA cells. The loss of naive and accumulation of terminally differentiated T cells, observed here in smokers, mimics the altered distribution of T cell subsets reported in aging and chronic infection that is proposed to result from repeated or persistent stimulation of immune cells, ultimately leading to the loss of immune function either due to replicative senescence or functional exhaustion.46, 47, 48 Gene expression changes in low-frequency subsets may not be discernible in pseudobulk and bulk analyses. Therefore, we looked for indicators of T cell dysfunction within the smoking-associated NK-like CD8 TEMRA subset. While TOX, a transcription factor that controls fate commitment in exhausted T cells, was elevated in smokers’ NK-like CD8 TEMRA cells compared to other CD8 T cells, it was only detected in 8.9% of cells within this cluster. Other exhaustion markers PDCD1, CTLA4, and HAVCR2 were not increased. Whereas exhausted CD8 T cells lack cytotoxic activity, the high expression of genes encoding proteins responsible for cytolytic activity in NK-like CD8 TEMRA cells suggests that these cells more likely represent a senescent or pre-senescent state. In support, compared to other CD8 T cells, the NK-like CD8 TEMRA cells from both smokers and nonsmokers had an elevated expression of KLRG1, an inhibitory receptor correlated with extensive proliferative history, and B3GAT1 (CD57), a marker of limited proliferative potential and shortened telomeres. Senescence is induced as the result of telomere shortening or nontelomeric DNA damage, both of which have been reported to occur in smokers.,

To further examine the effects of smoking on immune aging and senescence in CD8 T cells, we used DNA methylation-based models to demonstrate that age acceleration and reduced TL correlated with smoking dose. Accelerated immune system aging accompanies T cell senescence and can manifest as impaired immunological memory, which could contribute to attenuated immune responses in smokers.

In addition to impaired immune function, prolonged SASP, driven by the accumulation of senescent cells, can lead to chronic inflammation. Callender et al. showed that CD8 TEMRA cells exhibit both inflammatory and senescent phenotypes. Consistent with SASP, TNFSF10 (secreted factor) and CX3CR1 (chemokine receptor) were increased in smokers’ CD8 T cells, and NK-like CD8 TEMRA cells expressed secreted factors CTSW and PRSS23. The PRSS23 methylation level is a reproducible biomarker of tobacco smoke exposure,55, 56, 57 and altered methylation persists up to at least 30 years after smoking cessation., This indicates that epigenetic modifications likely contribute to the senescent attributes of CD8 T cells in smokers. The acquisition of CD16 and NK-like characteristics implies an underappreciated role for CD16 receptor in the maintenance of cytotoxic activity in TEMRA cells in smokers.

In conclusion, we show an association between smoking and an immune-cell subtype that can be isolated to investigate how NK-like CD8 TEMRA cells contribute to proinflammatory states in smoking-mediated chronic inflammatory conditions. Our data illustrate links between smoking-induced gene expression changes and a T cell senescent phenotype, immune system aging, and, potentially, atherosclerosis. Consequently, our use of recently developed single-cell technologies to address tobacco smoke exposure has great potential to affect global health.

Limitations of Study

The single-cell approaches have low sample numbers. Reproducing our current findings in a larger group and additional characterization of CMV status will refine the roles of smoking and CMV in CD8 T cell aging. In addition, functional studies are needed to determine the importance of the immune cell changes reported here.

STAR★Methods

Key Resources Table

Resource Availability

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Dr. Douglas A. Bell (bell1@niehs.nih.gov), Senior Investigator, Environmental Epigenomics and Disease Group, Immunity Inflammation and Disease Laboratory, National Institute of Environmental Health Sciences, Research Triangle Park, NC, USA.

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

The datasets generated during this study are available at NCBI GEO, https://www.ncbi.nlm.nih.gov/geo/. Microarray data have been deposited to NCBI GEO: GSE138974. Bulk and single-cell RNA-seq data have been deposited to NCBI GEO: GSE138851 and NCBI GEO: GSE138867. Raw idat files for DNA methylation have been deposited to NCBI GEO: GSE147430. Mass cytometry data are available from the corresponding author on request.

Experimental Model and Subject Details

Human Subjects

All donors were recruited with written informed consent under approved human IRB protocol NIEHS 10-E-0063 by the NIEHS Clinical Research Unit between March 2013 to January 2018 from the Raleigh, Durham and Chapel Hill, NC area., Whole blood was obtained from healthy (without acute disease according to self-reported medical histories) from nonsmokers, not having smoked > 100 cigarettes in their lifetime, and smokers who reported their average daily cigarette consumption for the past 3 months. Serum nicotine/cotinine levels were measured by HPLC-MS (Quest, Inc.) as an indication of their smoking exposure status. Human cytomegalovirus (HCVM) status was determined by two methods at the NIH Clinical Center in Bethesda, MD. Anti-cytomegalovirus IgG and IgM antibodies were measured in serum by a chemiluminescence immunoassay, and HCMV DNA viral load was determined using CMV-specific probes by quantitative real-time PCR. Donors were recalled matching nonsmokers/smokers on age, gender, and ethnicity for whole blood collection, cotinine levels were measured. See Table S1 for additional donor information.

Method Details

PBMC Isolation for scRNA-seq and Mass Cytometry

Whole blood was diluted 1:5 v/v with QIAGEN Buffer EL and incubated at room temperature (RT) until clarified (∼10 min) before centrifugation (300 g, RT, 10 min). After supernatant removal, leukocytes were resuspended in the same volume of Buffer EL (5 min) before spinning (300 g, 8 min). Leukocytes were then washed twice in autoMACS Running Buffer (Miltenyi Biotec), counted, and cryopreserved [20% Iscove’s Modified Dulbecco’s Medium (IMDM), 70% Fetal Bovine Serum (FBS), 10% Dimethyl sulfoxide (DMSO)] at a concentration of 1x107 cells/mL. Cryopreserved cells were thawed in nonsmoker/smoker pairs following the 10X Genomic’s protocol for “Fresh Frozen Human Peripheral Blood Mononuclear Cells for Single Cell RNA Sequencing.” Briefly, cells were serially diluted dropwise in complete media (IMDM,10% FBS) adding 50U/mL Benzonase (Millipore Sigma) for the first dilution. After centrifugation (1100 rpm, 8 min, RT), cells were resuspended in complete media and incubated with CD15 Dynabeads (Thermo Fisher Scientific) according to the manufacturer’s instructions to deplete the neutrophils from the PBMCs. PBMCs were then counted for viability and aliquoted for scRNA-seq or mass cytometry in parallel.

PBMC Preparation for Purified Cell Fractions

The mononuclear layer was isolated directly from whole blood using density gradient centrifugation with Histopaque-1077 Ficoll and ACCUSPIN Tubes (Sigma Millipore). Purified CD4+, CD8+, CD14+, CD19+, and CD56+ cell fractions were collected using antibody-coated magnetic beads (Dynabeads, Thermo Fisher Scientific; CD56, Miltenyi Biotec). Antibody-purified fractions were then extracted for DNA and RNA using the AllPrep DNA/RNA/miRNA Universal Kit according to the manufacturer’s instructions (QIAGEN).

Mass Cytometry

Thawed PBMCs (∼3x106 cells) were spun (300 g, 5 min) and resuspended in calcium magnesium-free phosphate buffered saline (PBS). 1μM Cisplatin (Fluidigm) was added for viability staining for 5 minutes before quenching the reaction with MaxPar Cell Staining Buffer (CSB, Fluidigm). After centrifugation (300 g, 5 min), cells were resuspended in CSB at a concentration of 60x106 cells/mL and incubated (RT,10 min) with Fc receptor binding inhibitor before adding 26 MaxPar metal-conjugated antibodies (Fluidigm) against immunophenotypic markers for an additional 30-minute incubation at RT. Stained cells were then washed two times before resuspension in MaxPar fix and perm buffer with 125μM 191/193Ir intercalator for either an hour at RT or 4°C overnight. Cells were then washed twice with CSB and two times with Nuclease-Free water (Thermo Fisher Scientific) followed by filtering through 40μM strainers to remove aggregates. Cells were then counted and resuspended in nuclease-free water at ∼5x105 cells/mL with 1:10 volume of four-element calibration beads (Fluidigm) and analyzed on a Helios instrument (Fluidigm) for 250,000 events for each donor at the NIEHS Flow Cytometry Center. Following the manufacturer’s instructions, downstream processing involved normalization by the calibration beads and fcs files were uploaded to Cytobank.

Preparation of Sequencing Libraries

scRNA-seq, libraries were prepared with the Chromium Single Cell 3′ Library & Gel Bead Kit v2 (10X Genomics) according to the manufacturer’s guidelines. For bulk RNA-seq, RNA from isolated CD8+ purified cell fractions were prepared using the TruSeq Stranded Total RNA Library Prep Gold (Illumina).

Microarray

Isolated RNA from antibody-purified cell fractions (CD4, CD8, CD14, CD19, CD56) from 19 individuals (9 smokers, 10 nonsmokers (See Table S1)), was amplified using NuGEN WT-Ovation Pico RNA Amplification System followed by labeling with NuGEN Encore Biotin Module according to the manufacturer’s protocol (NuGEN). Amplified biotin-cDNA was fragmented and hybridized to streptavidin/phycoerythrin-stained arrays using the GeneChip™ Hybridization, Wash and Stain Kit according to the manual protocol FS450-0004 (Thermo Fisher Scientific). The Life Technologies Human Transcriptome Arrays v2.0 were then scanned by an Affymetrix Scanner 3000 and using Transcriptome Analysis Console (TAC) Software (Thermo Fisher Scientific).

DNA Methylation arrays

Isolated DNA was bisulfite converted using the EZ-96 DNA Methylation MagPrep Kit (Zymo Research). Converted DNA was applied to either an Illumina Infinium HumanMethylation450 BeadChip (450K) or an Infinium MethylationEPIC BeadChip Kit (EPIC) according to the manufacturers’ protocol to measure methylation at ∼450,000 (450K) or 850,000 (EPIC) CpG sites genome wide.

Quantification and Statistical Analysis

Mass Cytometry Gating Strategy and Analysis

Events were gated in Cytobank to identify single viable cells (Figure S1A). Cells gated from spiked-in normalization beads were subsequently gated by Iridium (191Ir) and Cisplatin (198Pt) to obtain DNA positive cells. Single cells were identified by event length and Iridium (193Ir) and viable cells by Cisplatin-198Pt and leukocyte marker CD45. Viable cells were exported as fcs files and imported into VorteX using all events for each donor totaling 990,748 cells. Using the default parameter recommendations, all data were transformed using hyperbolic arcsin (f = 5). Applying a noise threshold of 1.0, clustering analysis was performed using a Euclidean length profile of 1.0 in X-shift and the weighted K-nearest neighbor density estimation (K). An elbow point validation was performed to determine the optimal cluster K value (K = 25) which was then used to create a Force-Directed Layout (FDL) for visualization colored by major immune cell type, expression of marker genes, cluster ID, smoking status and subject ID (Figures 1C, 1E, S1B, S1C, and S1E). 136 clusters were identified from the 990,748 events, one cluster was determined to be red blood cells (RBCs; positive expression for CD235a/b) and 13 clusters had multiple lineage markers and were determined to be doublets (e.g., positive expression profiles for CD19 and CD3) which was a total of 6900 cells that were removed prior to downstream analysis (983,848 cells retained).

scRNA-seq Processing and Analysis

scRNA-seq data was aligned to the hg19 genome and processed with CellRanger version 2.0.2. Uniquely aligned reads sharing equal barcode × unique molecular identifier (UMI) tags but annotated to multiple protein-coding transcripts (i.e., ambiguous UMIs), within each replicate were discarded from the analysis. Cells with less than 200 or greater than 3000 genes and cells with greater than 10% of UMIs from mitochondria were removed. Dataset integration, SNN clustering, and UMAP visualization of scRNA-seq data were performed with Seurat version 3.0.0.9000. To integrate data across eight donor samples, we used 2,500 genes and 50 dimensions. Clusters were identified by SNN clustering with Seurat FindNeighbors and FindClusters functions using 50 dimensions and a resolution parameter of 2.5. UMAP was used to visualize cells colored by major immune cell type (Figure 1B), cluster ID (Figure 2A), and smoking status (Figure 2B). Cell cluster marker genes (padj < 0.05) and smoking DEGs (padj < 0.05) were identified using Seurat version 2.3.4 implementing the MAST algorithm with UMI included as a latent variable. Expression of marker genes were visualized via UMAP (Figure 1D).

Slingshot version 1.2.0 was used to perform pseudotime analysis. Principal component analysis and UMAP were run on CD8 T cells using 45 dimensions. Slingshot was then used to infer cluster lineages and assign pseudotime to CD8T cells (Figure 3C). Temporally expressed genes for each CD8 T lineage were then identified by fitting a generalized additive model to each gene using loess-smoothed pseudotime as the predictor variable (Figure S3).

Bulk RNA-seq Analysis

For antibody-purified CD8 T cell fraction bulk RNA-seq, reads were aligned to the hg19 genome with STAR. Gene read counts were obtained with featureCounts from the Subread package using release 27 of the GENCODE annotation. DEGs were determined using DESeq2 with FDR-adjusted p value < 0.05 as the cutoff for significance (Figure 5B). PCA was performed with the prcomp() function in R (Figure 5C).

Microarray Analysis

For microarray, differentially expressed genes were detected using log2-transformed expression fold-change estimates with respect to the composite average of RMA-corrected fluorescence log-intensity levels (log2FC) across matched fractions (CD14, CD19, CD4, CD56, and CD8) from multiple individual female donors, both smoking and nonsmoking (N = 53 overall, with N ≥ 5 per cell fraction × smoking status group). Probe-wise log2FC values were tested across statistical groups through a resolution-weighted ANOVA; resolution weights represented relative metrics of fluorescence discrimination in the dynamic range of detection, i.e., cumulative hazard of multivariate ANOVA significance scores (probe × cell fraction × smoking status) from probe-wise generalized linear modeling of RMA-corrected fluorescence log-intensities using an exponential distribution and inverse link function. DEGs were detected from the annotation of probes with significance level p < 0.05 adjusted for multiple comparisons, then filtered against a minimum probe-wise effect size δlog2FC > 0.3 × σlog2FC and post hoc pairwise significance (Student’s t test p < 0.05) between log2FC values of at least one matched comparison between smokers and nonsmokers on same cell fraction levels. For probe-level effect size filtering, δlog2FC = 0.3 × σSSR is 5% of the 6σ-spread log2FC regression error with respect to a probe’s grand mean [where (σSSR)2 = (SSRlog2FC)/(N-1)] compared to 5% of the 6σ-spread in measurement error about the mean log2FC of each statistical group in the probe [where (σlog2FC)2 = (SSElog2FC)/(N-1)].

GSEA

We used GSEA, via GenePattern to perform gene set enrichment analysis for Chemical and Genetic Perturbations and Immunological Signatures gene sets for scRNaseq, bulk RNA-seq, and microarray using FWER < 0.05 as the cutoff for significant enrichment (Figure 5E).

DNA Methylation Processing and Analysis

The raw idat files from the 450K and EPIC methylation arrays were read into R with the Minfi package separately. Data were combined for common CpG sites on the two arrays with Minfi combineArrays function and then preprocessed with background and dye bias correction using the Noob method. Methylation β values for CpG sites were calculated and values were input into a DNA Methylation Age Calculator to predict epigenetic biomarkers for aging (DNAm PhenoAge) and aging acceleration, determined by the correlation between DNAm PhenoAge and chronological age (DNAm PhenoAgeAccel, Figures 5F and 5H). Telomere length was estimated with DNAmTL (Figures 5G and 5I). Univariable linear regression was used to determine p values for the association of AHRR methylation versus DNAm PhenoAgeAccel and DNAmTL.

Reanalysis of Microarray Data

We reanalyzed Callender et al. gene expression data (GSE98640) from isolated CD8 T cells subsets: TN, TCM, TEM;, and TEMRA using GEO2R with default parameters to identify genes with differential expression among CD8 T subsets using limma (Linear Models for Microarray Analysis) with FDR-adjusted p values.,

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
CD45 (clone HI30) 89Y	Fluidigm	Cat# 3089003; RRID: AB_2661851
CD235ab (clone HIR2) 141Pr	Fluidigm	Cat# 3141001B; RRID: AB_2651154
CD19 (clone HIB19) 142Nd	Fluidigm	Cat# 3142001; RRID: AB_2651155
CD4 (clone RPA-T4) 145Nd	Fluidigm	Cat# 3145001; RRID: AB_2661789
CD8a (clone RPA-T8) 146Nd	Fluidigm	Cat# 3146001; RRID: AB_2687641
CD7 (clone CD7-6B7) 147Sm	Fluidigm	Cat# 3147006; RRID: AB_2802104
CD66 (clone CD66a-B1.1) 149Sm	Fluidigm	Cat# 3149008; RRID: AB_2802105
CD61 (clone VI-PL2) 150Nd	Fluidigm	Cat# 3150001; RRID: AB_2661793
CD123 (clone 6H6) 151Eu	Fluidigm	Cat# 3151001; RRID: AB_2661794
CD36 (clone 5-271) 152Sm	Fluidigm	Cat# 3152007; RRID: AB_2802106
CD45RA (clone H100) 153Eu	Fluidigm	Cat# 3153001; RRID: AB_2802108
CD163 (clone GHI/61) 154Sm	Fluidigm	Cat# 3154007; RRID: AB_2661797
CD10 (clone HI10a) 156Gd	Fluidigm	Cat# 3156001; RRID: AB_2802107
CD11c (clone Bu15) 159Tb	Fluidigm	Cat# 3159001; RRID: AB_2661800
CD14 (clone M5E2) 160Gd	Fluidigm	Cat# 3160001; RRID: AB_2687634
CD16 (clone 3G8) 165Ho	Fluidigm	Cat# 3165001; RRID: AB_2802109
CD34 (clone 581) 166Er	Fluidigm	Cat# 3166012; RRID: AB_2756424
CD38 (clone HIT2) 167Er	Fluidigm	Cat# 3167001; RRID: AB_2802110
CD206 (clone 15-2) 168Er	Fluidigm	Cat# 3168008; RRID: AB_2661805
CD33 (clone WM53) 169Tm	Fluidigm	Cat# 3169010; RRID: AB_2802111
CD3 (clone UCHT1) 170Er	Fluidigm	Cat# 3170001; RRID: AB_2661807
CD20 (clone 2H7) 171Yb	Fluidigm	Cat# 3171012; RRID: AB_2802112
CD15 (clone W6D3) 172Yb	Fluidigm	Cat# 3172021; RRID: AB_2802113
HLA-DR (clone L243) 174Yb	Fluidigm	Cat# 3174001; RRID: AB_2665397
CD56 (clone NCAM 16.2) 176Yb	Fluidigm	Cat# 3176008; RRID: AB_2661813
CD11b (clone ICRF44) 209Bi	Fluidigm	Cat# 3209003; RRID: AB_2687654
Human Fc Receptor Binding Inhibitor	Thermo Fisher Scientific	Cat# 14-9161-73; RRID: AB_468582
Biological Samples
Whole blood from healthy donors	NIEHS Clinical Research Unit (CRU)	See Table S1 for details
Chemicals, Peptides, and Recombinant Proteins
Histopaque®-1077	Sigma Millipore	Cat# 10771
ACCUSPIN Tubes	Sigma Millipore	Cat# A2055
Iscove’s Modified Dulbecco’s Medium (IMDM)	Thermo Fisher Scientific	Cat# 12200036
Fetal Bovine Serum (FBS)	Gemini Bio-Products	Cat# 100-106
autoMACS Running Buffer	Miltenyi Biotec	Cat# 130-091-221
Dimethyl sulfoxide (DMSO)	Millipore Sigma	Cat# D84-18
Ultrapure BSA	Ambion by Life Technologies	Cat# AM2618
Benzonase Nuclease	Millipore Sigma	Cat# E8263-25KU
Phosphate Buffered Saline (PBS)	Thermo Fisher Scientific	Cat# 20012-027
Dynabeads CD15	Thermo Fisher Scientific	Cat# 11137D
Cell-ID Cisplatin-198Pt	Fluidigm	Cat# 201198
Cell-ID Intercalator-Ir—500 μM	Fluidigm	Cat# 201192B
Maxpar® Cell Staining Buffer	Fluidigm	Cat# 201068
Maxpar® Fix and Perm Buffer	Fluidigm	Cat# 201067
EQ Four Element Calibration Beads	Fluidigm	Cat# 201078
Macron 0754-06 Sodium Citrate, Dihydrate	Thomas Scientific	Cat# 0562N59
SPRIselect Reagent	Beckman Coulter	Cat# B23318
Dynabeads MyOne Silane	Thermo Fisher Scientific	Cat# 37002D
Dynabeads CD19 Pan B	Thermo Fisher Scientific	Cat# 11143D
Dynabeads CD4	Thermo Fisher Scientific	Cat# 11145D
Dynabeads CD8	Thermo Fisher Scientific	Cat# 11147D
Dynabeads CD14	Thermo Fisher Scientific	Cat# 11149D
CD56 MicroBeads, human	Miltenyi Biotec	Cat# 130-050-401
AllPrep DNA/RNA/miRNA Universal Kit	QIAGEN	Cat# 80224
Flowmi™ Tip Strainers	Bel-Art	Cat# H1680-0040
Nuclease-Free Water	Thermo Fisher Scientific	Cat# AM9930
Critical Commercial Assays
Chromium Single Cell 3′ Library & Gel Bead Kit v2	10X Genomics, Inc.	Cat# 120237
Chromium i7 Multiplex Kit	10X Genomics, Inc.	Cat# 120262
Bioanalyzer High Sensitivity DNA Analysis	Agilent	Cat# 5067-4626
Ovation Pico WTA System V2	NuGEN	Cat# 3302-60
Encore Biotin Module	NuGEN	Cat# 4200-60
GeneChip Hybridization, Wash, and Stain Kit	Thermo Fisher Scientific	Cat# 900720
GeneChip Human Transcriptome Array 2.0	Thermo Fisher Scientific	Cat# 902233
EZ-96 DNA Methylation MagPrep	Zymo Research	Cat# D5041
Infinium HumanMethylation450 BeadChip Kit	Illumina	Cat# WG-314-1002
Infinium MethylationEPIC BeadChip Kit	Illumina	Cat# WG-317-1003
TruSeq Stranded Total RNA Gold	Illumina	Cat# 20020598
Deposited Data
Single-cell RNA sequencing	This paper	GEO: GSE138867
Bulk RNA sequencing	This paper	GEO: GSE138851
Microarray	This paper	GEO: GSE13897
CD8 T DNA Methylation Array Profiling	This paper	GEO: GSE147430
Flow-sorted CD8 T Microarray	Callender et al.17	GEO: GSE98640
Software and Algorithms
GraphPad Prism 7	GraphPad Software, Inc.	GraphPad Prism,RRID: SCR_002798
Cytobank	Cytobank	https://cytobank.org, RRID: SCR_014043
Transcriptome Analysis Console (TAC)	Thermo Fisher Scientific	Transcriptome Analysis Console, RRID: SCR_016519
VorteX v26	Samusik et al.14	https://github.com/nolanlab/vortex, RRID: SCR_017047
CellRanger v2.0.2	10x Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation, RRID: SCR_017344
Seurat v3.0.0.9000	Stuart et al.11	https://satijalab.org/seurat
Seurat v2.3.4	Butler et al.10	https://satijalab.org/seurat
R	R Core	https://www.r-project.org/
Slingshot v1.2.0	Street et al.22	https://bioconductor.org/packages/release/bioc/html/slingshot.html, RRID: SCR_017012
GSEA	Subramanian et al.83	https://www.gsea-msigdb.org/gsea/downloads.jsp, RRID: SCR_003199
	Mootha et al.84
GenePattern	Reich et al.85	https://www.genepattern.org, RRID: SCR_003201
GEO2R	Barrett et al.86	https://www.ncbi.nlm.nih.gov/geo/geo2r/
Minfi	Aryee et al.87	bioconductor.org/packages/release/bioc/html/minfi.html
PreprocessNoob	Fortin et al.88	https://rdrr.io/bioc/minfi/man/preprocessNoob.html
DNAm PhenoAge	Levine et al.30	http://dnamage.genetics.ucla.edu/home
DNAm TL	Lu et al.31	http://dnamage.genetics.ucla.edu/home