High-throughput and high-dimensional single-cell analysis of antigen-specific CD8+ T cells

Samples and material

Human whole blood from patients diagnosed with T1D and T2D was obtained at Seton Family of Hospitals at Austin with informed consent. The use of whole blood from these patients was approved by the institutional review board of the Ascension Seton University Physicians Group under institutional review board number 2013-10-0140 and is compliant with all relevant ethical regulations. Human PBMCs from healthy donors were purchased from ePBMC.

Generation of DNA-barcoded fluorescent streptavidin

The conjugation of DNA linker (Supplementary Table 8) to PE- or APC-labeled streptavidin was performed as previously described, with slight modifications¹⁸. During S-HyNic modification of PE- or APC-labeled streptavidin, 2 mol equivalent of S-HyNic was used. Following the conjugation of DNA linker, peptide-encoding DNA barcodes (Supplementary Table 4) were annealed to the complementary DNA linker on the DNA-linker PE or APC streptavidin conjugate in the presence of 1× NEBbuffer2 (NEB) with the following program: 60 °C for 30 s, then −1 °C/cycle for 35 cycles. The final DNA-barcoded fluorescent streptavidin conjugate was stored at 4 °C.

IVTT

Peptide-encoding DNA oligonucleotides were purchased from Sigma-Aldrich. DNA templates (50 nM) were first amplified by PCR as described previously with modifications¹⁸. IVTT_r and IVTT_f primers (1 µM; Supplementary Table 8) were used in the following reaction conditions: 95 °C for 3 min; then 22 cycles of 95 °C for 20 s, 59 °C for 30 s and 72 °C for 30 s; then 72 °C for 5 min. The PCR product was then diluted with 50 µl nuclease-free water before proceeding to the IVTT reaction.

Generation of the pMHC tetramer library

IVTT-generated peptides were mixed with biotinylated pMHC monomers containing a UV-labile peptide. The UV-labile peptide-loaded pMHC monomers were provided by the National Institutes of Health tetramer core. The final concentration of biotinylated pMHC is 0.2 mg ml⁻¹. Individual pMHC was formed through UV exchange as described previously¹⁹. Confirmation of the quality and concentration of UV-exchanged pMHC monomer was assessed by an ELISA assay as described previously¹⁹. Individual pMHC tetramers and the tetramer library pool were generated and tested as described previously¹⁸. pMHC tetramer library should only be pooled together immediately before cell staining.

Customization of CD2 SampleTag, custom AbSeq and custom CD50 SampleTag

Anti-CD2 antibody was purchased from Biolegend (clone RPA-2.10, Biolegend). Amine-modified oligonucleotide was purchased from Sigma-Aldrich (Supplementary Table 8). The conjugation between the oligonucleotide and anti-CD2 antibody followed the CITE-seq protocol⁴.

Corresponding antibodies and used oligonucleotides are listed in Supplementary Table 9.

Twelve CD50 antibody SampleTags⁴⁰ were customized by BD Biosciences using the commercial SampleTag oligonucleotides.

Sorting and culture of antigen-specific CD8⁺ T cell polyclones

Seven types of tetramers with peptides chemically synthesized and UV-exchanged to MHC were used to raise antigen-specific polyclonal T cells (Supplementary Table 1). For each antigen specificity, 20 tetramer-positive CD8⁺ single T cells were sorted into each well of the 96-well plate and cultured for 3 weeks. Polyclonal T cell expansion and culture were performed according to a previously published protocol⁴¹.

pMHC tetramer staining and sorting of primary human CD8⁺ T cells

PBMCs from T1D whole blood were isolated using Ficoll-Paque density-gradient centrifugation (GE Healthcare). CD8⁺ T cells were then enriched from PBMCs of T1DM and healthy donors using the EasySep Human CD8⁺ T cell isolation kit (Stemcell Technologies).

CD8⁺ T cells were resuspended in FACS buffer containing 0.05% sodium azide and 50 nM dasatinib. CD8⁺ T cells were then incubated at 37 °C for 30–60 min. Approximately 10,000 cells from an HCV peptide-binding clone used previously¹⁸ were prestained with BV510 anti-CD8a antibody (clone RPA-T8, Biolegend) and spiked into the primary CD8⁺ T cells. Following dasatinib treatment, the tetramer pool, together with BV421 anti-CD8a antibody (clone RPA-T8, Biolegend), was directly added into the cells. Cells were incubated at 4 °C for 1 h with continuous rotation. After washing, cells were further stained at 4 °C for 20 min with the presence of 5 µg ml⁻¹ mouse anti-PE (clone PE001, Biolegend) and/or mouse anti-APC (clone APC003, Biolegend). AbSeq staining mastermix was prepared by pooling 1 µl of each AbSeq together (Supplementary Table 9). Cells were washed in FACS buffer once and stained with the AbSeq mastermix. Additional dump-channel antibodies (AF488-anti-CD4, AF488-anti-CD14 and AF488-anti-CD19), 7-aminoactinomycin D and 2 µl anti-CD50 SampleTag were mixed in cells. Cells were incubated at 4 °C for 40 min prior to washing in FACS buffer twice and then sorted.

During cell sorting, approximately 50,000 tetramer-negative CD8⁺ T cells were also sorted and then later spiked into tetramer-positive T cells. FlowJo V10 was used to process FACS data.

BD Rhapsody sequencing library preparation and sequencing

Prior to BD Rhapsody processing, tetramer-negative CD8⁺ T cells were first stained with 2 µl CD2 SampleTag at 4 °C for 30 min. Cells were washed in FACS buffer three times and resuspended in 100 µl BD Sample Buffer. Sorted tetramer-positive CD8⁺ T cells and tetramer-negative CD8⁺ T cells were counted using BD Rhapsody. Tetramer-positive and tetramer-negative CD8⁺ T cells were pooled and processed on a BD Rhapsody cartridge following the user’s manual. Single-cell mRNA, AbSeq barcodes, tetramer barcodes and SampleTag barcodes were all captured by BD Rhapsody beads coated with poly(T) oligonucleotide, with a unique cell barcode and molecular barcode on each bead. Single-cell cDNA synthesis and library amplification were performed following the manufacturer’s protocol, with some modifications. Briefly, in PCR1, 1.2 µl tetramer PCR1 primer was added to the PCR reaction in addition to primers for gene expression panel, AbSeq, SampleTag and universal oligonucleotides (Supplementary Table 9). Nine and ten PCR cycles were used for 5,000–10,000 and 10,001–20,000 cells, respectively. Double-sided size selection with AMPure beads was performed to purify short amplicons (AbSeq, SampleTag and tetramer DNA-barcodes) and long amplicons (target genes and TCRα/β) separately. In PCR2, five separate PCR reactions with 15 reaction cycles were carried out to amplify gene panel, SampleTag, TCRα, TCRβ and tetramer DNA barcodes. AbSeq, tetramer and TCRα/β libraries were gel extracted for the desired band before proceeding to PCR3. Finally, eight cycles of PCR reactions were performed for all six elements following the manufacturer’s instructions. All PCR libraries were quantified using Bioanalyzer 2100 and pooled. Fifteen percent PhiX was used in all sequencing runs. Pooled libraries were sequenced on HiSeq X with PE150.

BD Rhapsody sequencing preprocessing

Sequencing reads from target gene expression, AbSeq, SampleTag, TCRα/β and tetramer DNA barcodes were processed as described below (Supplementary Note).

For target gene expression and AbSeq sequencing, reads were processed with BD Targeted Multiplex Rhapsody Analysis Pipeline Version 1.5 on the Seven Bridges platform following the manufacturer’s instructions. For tetramer and SampleTag sequencing, reads were processed with custom codes and are available in GitHub. True cell barcodes were converted to oligonucleotide sequences according to BD cell barcode indexing rules. Then, sequencing data of tetramer, TCRα and TCRβ were processed using umitools⁴² to extract a cellular barcode and unique MID for each read. Reads that are mapped to true cell barcodes were obtained.

For tetramer DNA barcodes, only reads that were an exact match for the tetramer DNA-barcode reference were retained. The number of reads of the same MID-tagged tetramer DNA barcode (unique tetramer DNA barcode) was counted for each cell. The distribution of the reads of unique tetramer DNA barcode follows a bimodal distribution as reported previously¹⁸. The first peak corresponds to PCR and sequencing errors, and thus, reads falling under the first peak were filtered out. Further, the number of MIDs aligned to each tetramer DNA barcode in each cell was determined to construct a tetramer DNA-barcode count matrix.

For the SampleTag DNA barcodes, reads were mapped to SampleTag DNA-barcode reference using bowtie2 with – norc and – local mode⁴³. Aligned reads were then processed using umitools to count the number of MIDs for each SampleTag DNA barcode in each cell. The distribution of MID counts for each SampleTag was fitted by a bimodal distribution, and the cutoff between two distributions was set as the negative threshold for the corresponding SampleTag. In addition, to recover false-negative SampleTag signals, SampleTags whose MID counts accounted for more than 50% of total SampleTag MID counts were also classified as a positive event. Cells containing CD2 SampleTag were tetramer-negative cells, whereas cells with more than two regular SampleTags were multiplets and were removed from further analysis.

For the TCR sequencing reads, we adapted a subclustering algorithm as previously described⁴⁴ to remove PCR and/or sequencing errors and identify VDJ and CDR3, with some changes. Reads were first aligned to TCR J and C region reference. Only reads that are more than 62.5% identical were retained. Reads with the same cellular barcodes and MID were grouped together. Under each group, reads within a Levenshtein distance of 15% were further clustered into a subgroup. For each subgroup, a consensus sequence was built based on the average nucleotide at each position, weighted by quality score. After ranking the consensus sequences by their abundance, the most abundant consensus sequence was selected, and other sequences with an edited distance of less than three were removed. In cases where the most abundant consensus sequence was nonproductive, the next most abundant productive sequence (if it existed) was selected as the unique consensus sequence for that cell. The second TCR chain was retained when its MID count accounted for more than 20% of total TCRα or TCRβ MID counts.

Dimensionality reduction, clustering and differential expression of single cells

All single cells were first filtered to exclude low-quality cells whose total gene and AbSeq expression MID counts were in the last 1% quantile. Then, cells identified as multiplets with SampleTag and cells with two productive TCRβ chains were also removed. Additionally, genes or AbSeqs whose expression was detected in fewer than 50 cells were filtered out. Gene expression and AbSeq data from different Rhapsody chips were pooled together and used to perform joint probabilistic modeling of RNA expression and surface-protein measurement with totalVI²³. Each donor was treated as an independent batch factor and 200 epochs were used to train the model. Other parameters were set as default in totalVI. The posterior dataset was then used for dimensionality reduction (UMAP algorithm) and clustering (Leiden algorithm), both with Scanpy⁴⁵.

Calling tetramer specificity for each cell

First, for each tetramer fluorescent color, the distribution of total tetramer DNA-barcode counts per cell was fitted to a bimodal distribution. The cutoff counts were set as the negative threshold to capture positive tetramer-binding events. Tetramer DNA-barcode counts were then ranked for each cell, and the knee point on the count-rank plot was selected. Antigens that ranked higher than the inflection point were included as putative binding antigens, and antigens that ranked below the inflection point but showed a difference of three amino acids or fewer compared with higher-ranking antigens were also included as putative cross-reactive binding antigens. For each cell, the tetramer MID signal fraction was defined as the fraction of the cumulative MID count from putative binding antigens divided by the cumulative MID count from all bound antigens:

Cells with a tetramer MID signal fraction below 0.4 were prefiltered in the preprocessing step to identify antigen specificities. Further, cells with the same TCRα/β were pooled together. The correlation coefficient of antigen binding for each cell in the pool was calculated between detected tetramer DNA-barcode counts and the corresponding median tetramer DNA-barcode counts within the pool. This correlation coefficient for each cell was used as the tetramer-binding noise. The knee point of the distribution of correlation coefficients was set as the threshold below which cells were removed due to high tetramer-binding noise.

For analysis of viral antigens, we selected antigens detected in more than five cells to ensure the capture of low-frequency antigen-specific CD8⁺ T cells while limiting nonspecific binding.

For sensitivity analysis to demonstrate the robustness of TetTCR-SeqHD, we set the negative threshold of tetramer MID to 15 to capture positive binding events. This threshold was then used for all experiments.

Precision and recall rate calculation for TetTCR-SeqHD

In the TetTCR-SeqHD clone experiment, true positive is defined as antigen-matched TCRs between MIDCIRS and TetTCR-SeqHD. Predicted condition positive is defined as antigen-specific TCRs identified by pMHC DNA barcodes. The condition positive is defined as antigen-specific TCRs identified by MIDCIRS. Precision and recall are then calculated as follows:

Prediction of pMHC class I binding

HLA-A02:01-bound T1D autoantigens were curated from the IEDB (www.iedb.org) database, while HLA-A01:01- and HLA-B08:01-bound T1D autoantigens were predicted using NetMHCpan 4.0 (ref. ⁴⁶). The half-maximum inhibitory concentration cutoff for HLA-A01:01 and HLA-B08:01 was 950 nM and 500 nM, respectively.

TCR clonality calculation

TCRs that have productive paired α and β chains were used to calculate TCR clonality, which is a score to characterize T cell expansion. Higher TCR clonality indicates that corresponding TCRs are more clonally expanded. If there is a singleton TCR, we define the TCR clonality as 0, while single TCR species with multiple copies have a TCR clonality of 1. For all other situations, the TCR clonality is defined using the following formula:

Calculation of antigen-specific T cell frequency

The absolute frequency of antigen-specific T cells for antigen a_i in each donor was calculated as follows:

.

For cross-reactive cells, especially when cells are cross-reactive with more than two antigens in the antigen panel, one cell can be identified to bind a combination of antigen specificities by TetTCR-SeqHD. Each combination is a binding pattern. The frequency for cross-reactive antigen-specific T cells was calculated for each binding pattern.

TCR transduction

We generated TCR constructs as previously described¹⁸ and cloned them into an empty pCDH (System Biosciences) vector driven by the MSCV promoter. Lentivirus was generated using the Virapower (ThermoFisher Scientific) system and concentrated 10 times using an Amicon Ultra column. Freshly thawed CD8⁺ T cells from an HLA-A2-, HLA-B8 and HLA-A1-negative donor were stimulated with Immunocult (Stemcell Technologies) and incubated with the concentrated virus for 2–3 days. The cells were expanded for a minimum of 10 days and then assessed for murine TCRβ chain expression.

Flow cytometry on transduced cells

Tetramer staining was performed as previously described¹⁸ with tetrameric MHC loaded with chemically synthesized peptides (Genscript). Briefly, the transduced cells and negative controls were stained with an anti-CD8a antibody (clone RPA-T8, Biolegend) before the addition of tetramer for 1 h on ice. Negative controls were established using nonspecific tetramer (HLA-A*02:01:HCVns3:1406-1415 – KLVALGINAV) and untransduced T cells from the same donor. Cross-TCR and cross-HLA negative controls were also included to assess the degree of nonspecific activity. After washing, the cells were stained with an anti-murine TCRβ antibody (Biolegend) and 7-aminoactinomycin D before analysis on a BD Accuri.

T2 cells (generously provided by the Mark Davis lab) were pulsed with a chemically synthesized peptide (10 µM) for 2 h at 37 °C. The cells were then washed and incubated 1:1 with the transduced cells for 4 h at 37 °C. Negative controls were performed using nonspecific peptide (HCVns3:1406-1415) and cross-TCR nonspecific peptides (for example, EBV-BLMF1 peptide was used as a negative control for T1D antigen cross-reactive TCR, TCR51), while positive control was performed using PMA/ionomycin (Cell Stimulation Cocktail, Biolegend). During incubation, anti-CD107α (Biolegend) antibody and monensin were added to detect and stabilize degranulation events. The assay was stopped via the addition of cold PBS and subsequent staining for CD107α, CD8α and murine TCRβ (Biolegend). Cells were analyzed via a BD Accuri.

Detection of autoantibodies

The presence of anti-GAD, anti-IA2 and anti-Znt8 antibodies was determined via ELISA obtained from Kronus and performed according to the manufacturer’s instructions. Whole, undiluted plasma was used in this assay. Absorbance was measured using a SpectraMax M3 plate reader, and analysis of the standard curve was performed in R using a cubic-spline fit. The antibody concentration for each sample was then interpolated, with all positive controls falling within the reported concentrations. Patients were reported as positive if the detectable antibody levels were in excess of 5 IU ml⁻¹, 7.5 U ml⁻¹ and 15 U ml⁻¹ for the anti-GAD, anti-IA2 and anti-Znt8 antibodies, respectively, according to the manufacturer’s instructions.

Statistics and reproducibility

The relevant statistical test, sample size, replicate type and P values for each figure are found in the figure and/or corresponding figure legend.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.