Single-cell analyses identify dysfunctional CD16 + CD8 T cells in smokers
 SUMMARY 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 (T EMRA ) 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 remains substantial. The World Health Organization (WHO) estimates approximately six million deaths per year from tobacco smoke exposure, resulting from both chronic and communicable diseases. In smokers, a decline in immunity and increased risk of inflammatory diseases, such as atherosclerosis, argues 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 prevention of diseases by identifying immune cells to target for clinical intervention. 
In addition to DNA damage, smoking alters the epigenome and transcriptome of human blood leukocytes. In Su et al., we demonstrated that changes identified in isolated cell fractions, which correspond to major immune populations, were distinct from each other 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, comprised of multiple subtypes, would similarly mask meaningful changes, especially when differences arise in low frequency subsets. As such, interpretation of bulk genomic approaches is limited because changes could indicate 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. Additionally, multiparameter data allows 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 (scRNAseq) 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 strong correlation between scRNAseq 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 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 smoking dose was associated with accelerated immune aging and decreased telomere length (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 provides evidence that smokers show signs of premature immune system aging. The potential immune function defects and inflammatory subsets demonstrated here mirror 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. 
scRNAseq 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 smoking-associated gene expression changes observed within major immune cell populations resulted from altered abundance of specific, identifiable cell subsets. We performed scRNAseq and mass cytometry, in parallel, on cryopreserved peripheral blood samples from eight donors with no previous history of atherosclerosis, chronic obstructive pulmonary disease (COPD), or lung cancer (Figure 1A). Serum cotinine, a metabolite of nicotine and 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 – 511 ng/ml; all nonsmokers (n=4) had serum cotinine levels below 2 ng/ml. Donors were matched based on gender and race; ages ranged from 31 – 56 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. 
For each single-cell approach, we assigned cells to common immune populations based on mRNA (scRNAseq) or surface protein (mass cytometry) expression (Figures 1B–1E and Table S2). For scRNAseq, 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 expression of CD4, CD8A, or NCR3 (Figure 1D). NK cells were identified based on CD3D as a negative marker combined with expression of NKG7, GNLY, GZMB, PRF1, and NCR3 as positive markers (Figure 1D and Table S2). Monocytes were positive for LYZ and either CD14 or FCGR3A (encodes CD16), characteristic of classical or nonclassical monocytes (Figure 1D). Dendritic cells were similar to monocytes but could be distinguished by 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). 122 PBMC cell clusters were determined from the eight 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 dendritic cells 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 scRNAseq 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 scRNAseq (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 major populations—CD4T, CD8T, NKT, B, Monocyte, and DC—were identified in all donors. We then compared the frequency of major populations in PBMCs by smoking status for scRNAseq 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 scRNAseq or mass cytometry (Figures S1F and S1G). Comparison of the major immune population percentages among individuals for scRNAseq and mass cytometry showed significant strong correlation (Pearson r=0.99, r2=0.98, p<0.0001) between methods (Figure 1H). 
scRNAseq 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 thirty-one immune cell clusters and one erythroid contaminant cluster, labeled based on abundance from 0 (most abundant) through 31 (Figure 2A and Table S2). We identified twelve CD4 T cell (0, 1, 3, 4, 6, 10, 13, 14, 17, 18, 20, and 27), seven CD8 T cell (2, 8, 11, 15, 19, 21, and 24), three NK cell (7, 22, and 25), four monocyte (5, 23, 26, and 30), and two B cell (9 and 12) clusters. NKTs (16), DCs (28), and MKs (31) were each contained by a single cluster. Herein, clusters are referred to by major immune cell type, followed by original cluster 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). 
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: median 3.5% in nonsmokers and 5.4% in smokers. We characterized CD4T-17 cells as regulatory T cells (Tregs) based on elevated FOXP3 and IL2RA (encodes CD25; Figure 2E and 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/duration. 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, variation among donors in CD8 T cells depended on smoking status as illustrated by differences of the dominant population(s) in CD8 T cells from smokers and nonsmokers (Figure 2F). Smokers had lower proportions of two CD8 T cell clusters (CD8T-24 and CD8T-2) and higher proportions of two CD8 T cell clusters (CD8T-8 and CD8T-15) compared to nonsmokers (p<0.05, Figure 2G). We did not observe significant differences in proportions of the remaining three CD8 T cell clusters (CD8T-19, CD8T-11, CD8T-21). 
CD8 T Cell Distribution Shifts from Naïve 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 (padj)<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 naïve cells (TN). High levels of IL7R, SELL, and FOS (activated T cell proliferation) suggested 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). 
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 two lineages, which originate from cluster CD8T-24 and terminate at either cluster CD8T-8 or CD8T-15 (Figure 3C). Lineage one was mostly comprised of cells from CD8T-24, CD8T-2, CD8T-19, CD8T-8, CD8T-21, and CD8T-8, with minimal cells from CD8T-15. Lineage two was mostly comprised of 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 smokers’ CD8 T cells are skewed toward differentiated and nonsmokers’ CD8 T cells are skewed toward naïve 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 one) and central memory (lineage two) differentiation. For example, in lineage one, CCL5 and NKG7 increased, while SELL and IL7R decreased over the differentiation trajectory (Figures S3A and S3C). CMC1 demonstrated a nonlinear association in lineage one, 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 PD-1), CTLA4, and HAVCR2 (encodes TIM3)—did not distinguish the clusters. Whereas smokers’ CD8T-8 cells had elevated expression of TOX (log2FC=0.31, padj=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, padj=8.11×10−04; nonsmokers (NS): log2FC=1.33, padj=2.23×10−22) and CD8T-8 cells (SM: log2FC=0.94, padj=6.29×10−47; NS: log2FC=1.28, padj=1.94×10−40), B3GAT1 was unique to CD8T-8 cells (SM: log2FC=0.19, padj=3.83×10−14; NS: log2FC=0.10, padj=0.022). We also found that two genes reported as having smoking-associated methylation changes, GFI1 and PRSS23, showed elevated expression in CD8T-8 cells (GFI1 SM: log2FC=0.20, padj=5.58×10−7; GFI1 NS: log2FC=0.22, padj=0.046; PRSS23 SM: log2FC=1.07, padj=1.5×10−143; PRSS23 NS: log2FC=0.91, padj=7.7×10−60). Surprisingly, 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, padj=1.17×10−171) and nonsmokers (log2FC=1.72, padj=5.17×10−119, Figure 1D and Table S3). 
We reanalyzed gene expression data from flow-sorted human peripheral naïve (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 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. Additionally, CD8 TEMRA cells showed 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 (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 naïve CD8 T clusters (CD8T-24 and CD8T-2) were significantly reduced in HCMV seropositive donors. (Figure S3D and Table S1). 
To further support that smoking elevated FCGR3A-expressing CD8 TEMRA cells (CD8T-8) and reduced naïve 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 and Table S1). Higher frequency of FCGR3A-expressing CD8 TEMRA cells was associated with reduced AHRR methylation (p=0.008, r2=0.72) and lower frequency of naïve 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. 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. That is, 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 and Table S1). 
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 Natural Killer-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 lead us to suspect that this subset might 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 (padj<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 naïve 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 naïve CD8 T vs NK genes. Figure 4F shows the 25 highest ranked naïve CD8 T vs 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 average expression of genes from smokers’ CD8T-8 cells to nonsmokers’ CD8T-8 cells (see STAR Methods). We found that 63 genes had increased and 74 genes had decreased average per cell expression in smokers compared to nonsmokers (padj<0.05, Table S4). Figures S4E and S4F show 25 genes with increased and decreased per cell expression, ordered by difference in 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 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 cells in the seven CD8 T clusters combined. Of 2163 genes evaluated in the pseudobulk analysis, we found that 1817 genes had higher expression and 344 genes had lower expression in smokers versus nonsmokers (q<0.05, Table S5). To examine 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 RNAseq and microarray on isolated CD8 T cells to examine differences in bulk RNA expression. Isolated CD8 T cells for bulk RNAseq included the eight donors used in scRNAseq (bulk RNA was from a previous visit) and seven additional donors (Table S1). We identified 1268 genes as differentially expressed; 692 increased and 576 decreased (q<0.05, Figure 5B). With exception of F061, principal component analysis of bulk RNAseq 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 scRNAseq, 4 donors used in bulk RNAseq, and 11 additional donors (Table S1). We identified 51 DEGs (see STAR Methods); 46 increased and 5 decreased (Figure 5D). 
We compared results from the three RNA analysis platforms used to identify smoking-associated DEGs in CD8 T cells. Bulk RNAseq confirmed 241 smoking DEGs from the scRNAseq, 162 with increased and 79 with decreased expression in smokers compared to nonsmokers (Figures 5B and Table S5). Microarray analysis confirmed 24 smoking DEGs that were identified by scRNAseq as increased in smokers compared to nonsmokers (Table S5; Figures 5A and 5D). Two genes, FAM129A and CD58, increased in smokers’ CD8 T cells in all three platforms (Figures 5A – 5B, 5D, and S5A – S5B and Table S5). Genes found to be altered by at least two methods include LGALS1, ADAM8, and CLDND1, which were increased in bulk RNAseq and scRNAseq data; GPR15, which was increased in the bulk RNAseq and microarray data; and NDFIP1, which was decreased in the bulk RNAseq and scRNAseq data (Figures 5A, 5B, 5D, and S5A – S5F). ITGAL, a smoking methylation biomarker, was only found to be significantly increased by scRNAseq (log2FC=0.36, padj=1.70×10−26), and was also elevated in NK-like subset (SM: log2FC=0.64, padj=6.53×10−29; NS: log2FC=0.78, padj=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 scRNAseq, bulk RNAseq, and bulk microarray shared seven positively and eight negatively enriched gene sets (FWER<0.05,Figure 5E). CD8 T cells were positively enriched for genes with higher expression in memory T, TCM, TEM, PD1lo (CD8 TEM), and PD1hi (CD8 TEM) cells relative to naïve T cells (Figure 5E). CD8 T cells were negatively enriched for genes that have higher expression in naïve 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 increased expression of genes associated with effector memory and central memory functions and decreased expression of genes associated with naïve T cells. 
CD8 T DNA Methylomes Indicate Smoking Dose-dependent Accelerated Aging and Decreased Telomere Length
In CD8 T cells, the shift from naïve 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 scRNAseq 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 are indicative of T cell senescence, we used DNAmTL to estimate telomere length 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 group of 131 individuals (Table S1). Age acceleration was highly significantly correlated with smoking dose (p<0.00009) and approximately 11% of the variance in age acceleration in this population could be attributed to smoking dose (Figure 5H). Similarly, methylation-derived CD8 T cell telomere length 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 1563 DEGs; 1278 showed increased and 285 showed decreased expression (Table S6). Hierarchical clustering of donors by CD4 T smoking DEGs clustered individuals by smoking status (Figures 6A). Bulk gene expression of isolated CD4+ cells identified two upregulated (LRRN3 and GPR171) and one 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 three upregulated genes (PROK2, MX1, and TRAT1) and one downregulated gene (KLRB1) with NKT cells and three upregulated (MX1, DOCK5, and CSGALNACT1) and three 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 and 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 and 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 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 scRNAseq and confirmed by microarray. 
DISCUSSION
Our study reveals CD16+ CD8 T cells and other signs of immune cell dysfunction as elevated in smokers. These cells, uncovered by scRNAseq 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 scRNAseq 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 TNFα-mediated senescence, increased in smokers in NKT cells, NK cells, and monocytes. 
Alluding to 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 naïve and gain TCM-like and TEMRA-like cells. Differentiation-state shifts in CD8 T cells were supported by bulk analysis methods. Interestingly, all three 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. (2011) demonstrated that PD-1hi CD8 T cells obtained from healthy adults had similar gene expression profiles to 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 scRNAseq, bulk RNAseq, or microarray were represented by these gene signatures (Figure 5E and Table 1). Notably, bulk RNAseq data was acquired from prior visits from donors used for scRNAseq, 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 IL-7Rα that drives terminal differentiation of CD8 T cells, was elevated in NK-like CD8 TEMRA cells. Loss of naïve 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 loss of immune function either due to replicative senescence or functional exhaustion. Gene expression changes in low-frequency subsets may not be discernable 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, 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 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 non-telomeric 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 telomere length 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 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. PRSS23 methylation level is a reproducible biomarker of tobacco smoke exposure and altered methylation persists up to at least 30 years after smoking cessation. Taken together, 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 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 illustrates 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 impact global health. 
Limitations of the 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. 
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 has been deposited at NCBI GEO (GSE138974). Bulk and single-cell RNAseq data have been deposited at NCBI GEO (GSE138851 and GSE138867). Raw idat files for DNA methylation have been deposited at NCBI GEO (GSE147430). Mass cytometry data are available from the corresponding author on request. 
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, sex and ethnicity for whole blood collection, cotinine levels were measured. See Table S1 for additional donor information. 
PBMC Isolation for scRNAseq 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 (300g, RT, 10 min). After supernatant removal, leukocytes were resuspended in the same volume of Buffer EL (5 min) before spinning (300g, 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 1×107 cells/mL. Cryopreserved cells were thawed in nonsmoker/smoker pairs following the10X 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 scRNAseq 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 (~3×106 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 60×106 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 ~5×105 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
scRNAseq, libraries were prepared with the Chromium™ Single Cell 3’ Library & Gel Bead Kit v2 (10X Genomics) according to the manufacturer’s guidelines. For bulk RNAseq, 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. 
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–C 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 (983848 cells retained). 
scRNAseq Processing and Analysis
scRNAseq 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 scRNAseq 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 RNAseq Analysis
For antibody-purified CD8 T cell fraction bulk RNAseq, 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 RNAseq, 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. (2018) 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. 

Table S1. Donor Information, Related to STAR Methods 
Table S2. Markers for 32 clusters (scRNAseq) relative to other PBMC clusters, Related to Figures 1 and 2. 
Table S3. Markers for CD8 T cell clusters (scRNAseq) relative to other CD8 T cells, Related to Figure 3. 
Table S4. Smoking DEGs (scRNAseq) for NK-like CD8 T cells (cluster CD8T-8), Related to Figure 4. 
Table S5. Smoking DEGs (bulk RNAseq) for CD8 T cells, Related to Figure 5. 
Table S6. Smoking DEGs (scRNAseq) for CD4 T cells, NKT cells, NK cells, monocytes, DCs, and B cells, Related to Figure 6. 
