Single-Cell RNA Sequencing Reveals Expanded Clones of Islet Antigen-Reactive CD4⁺ T Cells in Peripheral Blood of Subjects with Type 1 Diabetes

Accumulating evidence for a role of islet Ag-reactive CD4⁺ T cells in development of type 1 diabetes (T1D) has spurred efforts to use them to investigate disease mechanisms and as therapeutic targets and biomarkers for β cell destruction (1–6). Although levels of islet Ag-reactive cells may be increased in the pancreas (2, 3), biopsy of this organ is not tenable in humans. Instead, most efforts in humans have focused on peripheral blood, which is readily available for testing. Numerous studies have reported detection of islet Ag-reactive CD4⁺ T cells in blood of at-risk and T1D subjects, but these cells are often detected in healthy control (HC) subjects as well (7–9). Distinctive phenotypic properties of islet Ag-reactive CD4⁺ T cells in T1D subjects (8–11) suggest their relationship to disease. Early findings suggested that T1D was a Th1 disease (12), whereas subsequent studies suggest involvement of additional T cell subsets (13).

Another consideration in identifying CD4⁺ T cells important for disease progression is their proliferation in response to an antigenic peptide. This results in clonal expansion (14) of a population of cells with identical Ag specificity and unique, identically rearranged TCR α- and β-chains. Characterization of rearranged TCR sequence variation thus provides a measure of T cell diversity and Ag specificity, which can then be used to interrogate the role of those cells in disease.

Transcript profiling is a widely used tool for unbiased identification of phenotypic characteristics of cell populations. Increasingly, genome-wide transcriptome analysis by RNA sequencing (RNA-seq) has been extended to the single-cell level (15, 16), revealing heterogeneity that is masked in bulk profiling studies. Combining flow cytometry–based assays and single-cell RNA-seq, we have developed methods to identify TCR sequences in parallel with full transcriptome phenotypes from individual islet Ag-reactive CD4⁺ memory T cells. We have used this approach to perform an exploratory study of TCR clonotype expansion among islet T cells from HC and T1D subjects. We detected CD4⁺ memory T cells with expanded clonotypes in peripheral blood and identified their targets and transcript phenotypes.

Samples were obtained from HLA DRB1*0401 (DRB1*0401) HC and T1D subjects under informed consent (Table I). HC subjects were matched for age and gender to T1D patients, and they had no personal or family history of T1D. All protocols were approved by the Institutional Review Board at Benaroya Research Institute.

Peripheral blood (100 ml) was drawn by venous puncture using heparin as anticoagulant. PBMCs were isolated by Ficoll-Hypaque centrifugation and cultured in RPMI 1640 media supplemented with 10% commercial human serum (Gemini Bio Products, West Sacramento, CA), penicillin/streptomycin (100 U/ml, 100 μg/ml), sodium pyruvate (1 mM), and l-glutamine (2 mM). Immediately after isolation, PBMCs (10e6/ml) were stimulated with a pool of 28 islet peptides (Table II; 1.7 μg/ml each; Mimotopes, San Diego, CA) and 1 μg/ml anti-CD40 blocking mAb (Miltenyi Biotec, San Diego, CA) for 12–14 h at 37°C. As controls, PBMCs were cultured with an equal volume of DMSO (vehicle, negative control) or two influenza peptides as a positive control (MP p8 57–76 KGILGFVFTLTVPSERGLQR, MP p54 97–116 VKLYRKLKREITFHGAKEIS). Cells were harvested, labeled with PE-conjugated anti-CD154 mAb followed by anti-PE magnetic beads, and enriched using a magnetic column (Miltenyi Biotec). Enriched cells were stained with a live/dead dye (BD Via-Probe; BD Biosciences, Franklin Lakes, NJ) and Abs targeting surface markers: CD14-PerCP Cy5.5 (61D3; eBioscience, San Diego, CA), CD19-PerCP Cy5.5 (HIB19; eBioscience), CD4-Alexa Fluor 700 (OKT4; BioLegend, San Diego, CA), CD69-allophycocyanin (FN50; BioLegend), CD45RA-AmCyan (HI100; BD Biosciences, San Jose, CA), CD45RO-allophycocyanin Cy7 (UCHL1; BioLegend), CCR6-FITC (11A9; BD Biosciences), CXCR3-PE Cy7 (1C6/CXCR3; BD Biosciences), and CD38-Pacific Blue (HB7; eBioscience). Cells were gated as shown in Supplemental Fig. 2: lymphocytes, singlets, CD4⁺Via-Probe⁻CD14⁻CD19⁻, CD154⁺CD69⁺, CD45RA⁻RO⁺. CD4⁺CD154⁺CD69⁺CD45RA⁻RO⁺ T cells were sorted on a BD FACSAria II flow cytometer directly into a 96-well C1 microfluidic chip (Fluidigm, San Francisco, CA) for single-cell capture. For each sample, the CD154⁺CD69⁺ gate in islet-stimulated cultures was set based on the DMSO-treated culture. Each experiment was performed with PBMCs from a single blood draw.

DRB1*0401 MHC class II tetramers (Tmrs) labeled with PE were produced at Benaroya Research Institute Tetramer Core Laboratory and loaded with exogenous islet peptides (Table II) as described previously (17). As an irrelevant control, Tmrs were loaded with an influenza hemagglutinin (HA) peptide (HA 306–318 PRYVKQNTLKLAT). CD4⁺ T cell clones or transduced 5KC murine hybridoma cells (18) expressing human CD4⁺ (provided by M. Nakayama) were incubated with Tmrs at 37°C for 1–2 h, then surface stained with anti-CD4 and analyzed by flow cytometry. Tmr staining was assessed in gated CD4⁺ T cells.

The GAD65-specific DRB4-restricted T cell clone BRI4.13 was described previously (19). Cells were used directly (unstimulated) or were stimulated before use. Polyclonal stimulation by mAbs was achieved by incubation with immobilized anti-human CD3 plus soluble anti-human CD28 mAbs (eBioscience). Stimulation by Tmr (Ag-specific stimulation) was achieved by incubation in 96-well flat-bottom plates coated with class II Tmrs loaded with GAD555–567 (20) at 20 μg/ml. Cells were stimulated at 37°C for 12 h before use.

T cell clones were established from islet Ag-reactive CD4⁺ memory T cells from visit 3 of subject T1D2 using successive rounds of nonspecific activation with PHA and irradiated PBMCs in the presence of IL-2 (10 U/ml; Roche Applied Science, Mannheim, Germany). Clones were screened for expression of TRBV6-6 using a specific mAb (JU74.3; Beckman Coulter, Brea, CA) by flow cytometry, and clones testing positive were sequenced to confirm expression of the expanded TCR pair from T1D2.

Oligonucleotides (Genscript, Piscataway, NJ) encoding codon-optimized rearranged TRAV and TRBV sequences from expanded clonotypes were cloned into the modified “TCR flex” pMP71 retroviral backbone upstream of the murine Trac and Trbc genes (21). Recombinant retroviruses were packaged using Phoenix-AMPHO (CRL-3213; ATCC, Rockville, MD) by transfection of 5 μg of retroviral vector DNA with Lipofectamine 3000 transfection reagent (Thermo Fisher Scientific, Waltham, MA). Viral supernatants were collected at 48 and 72 h posttransfection. Purified human CD4⁺ T cells from peripheral blood (10⁶) were cultured in ImmunoCult-XF T Cell Medium (STEMCELL Technologies, Cambridge, MA) and activated with CD3/CD28 T cell activator (STEMCELL Technologies) in the presence of 100 IU/ml IL-2 and 5 ng/ml recombinant human IL-15 (BD Biosciences) for 48 h. Activated CD4⁺ T cells (0.2–0.5 × 10⁶) were suspended in 1 ml of retroviral supernatant diluted 1:2 and polybrene (final concentration 10 μg/ml) and transduced by spin inoculation (2000 rpm, 90 min). To maximize transduction efficiency, we repeated spin inoculation after 24 h. After 3 d, transduction efficiency was determined flow cytometry using a murine TCRβ C region mAb (H57-597; BD Biosciences). 5KC murine T cell hybridoma cells were transduced with recombinant TCR retroviruses as described earlier and sorted by flow cytometry to yield homogenous populations of human TCR-expressing cells for Tmr binding experiments.

Peptide-induced proliferation was detected by [³H]thymidine incorporation for T cell clones, CellTrace Violet (Thermo Fisher Scientific, Grand Island, NY; plus Ki57 staining [clone B56]; BD Pharmingen) for PBMCs, and CFSE dye dilution for transduced CD4⁺ T cells (22). T cell clones were screened in triplicate for proliferation to the original islet peptide pool in the presence of irradiated DRB1*0401 PBMCs for 96 h, followed by deconvolution to smaller pools and individual peptides. Peptide vehicle (DMSO) and anti-CD3/anti-CD28 mAbs were used as negative and positive controls, respectively. T cell clones with a stimulation index >3-fold over the DMSO control were considered to have proliferated. For CFSE proliferation assays, CFSE-labeled CD4⁺ T cells (10⁴) were mixed with 2 × 10⁴ irradiated APCs (Priess lymphoblastoid cells; ATCC), previously loaded with antigenic peptides (5 μg/ml, 1 h), and were cultured for up to 5 d. Priess cells (DRB1*0401, *0401) have been used to present T1D Ags in the context of that class II molecule (23). Cells were stained with mAbs specific for CD4 and the murine TCR β-chain. CFSE intensity was measured by flow cytometry and quantified by gating on the murine TCR⁺ or TCR⁻ populations in the CD4⁺ population.

Because islet Ag-reactive CD4⁺ memory T cells are rare (a median of ∼700 cells recovered in our experiments), it was important to optimize recovery of single cells. A major advance in our procedure was to sort islet Ag-reactive T cells directly into microfluidic chips (Fluidigm C1), which decreased the input cell number required, increased the cell capture yield, and resulted in better quality libraries (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables). This direct sorting procedure resulted in a median of ∼33 high-quality libraries per 96-well chip (n = 3 samples each for T1D and HC subjects). Neither the number of cells captured nor the number of high-quality libraries recovered differed between T1D and HC subjects (p > 0.4 by two-sided t test). After capture, cells were lysed, followed by reverse transcription (SMART-Seq v1 Ultra Low Input RNA Kit; Takara, Mountain View, CA) and cDNA amplification according to the manufacturers’ protocols (Fluidigm). Sequencing libraries were prepared using Nextera XT DNA kits (Illumina, San Diego, CA).

Single-cell libraries were sequenced on HiScanSQ or HiSeq2500 sequencers (Illumina) using single-read 100-bp dual-indexed reads (T cell clones) or single-read 58-bp dual-indexed reads (Ag-reactive T cells) to an average read depth of ∼1.8 million raw reads, a value that yields saturating numbers of genes detected in single-cell assays (24). Bulk RNA-seq libraries were sequenced to target depths of ∼10 million reads. RNA-seq pipeline analysis methods have previously been described (25). Quality metrics for aligned reads were obtained using the Picard (v.1.56) suite of tools (https://broadinstitute.github.io/picard/). For transcriptome analysis, reads were processed to remove reads with identical genome coordinates, which likely resulted from PCR amplification during library construction and comprised a large fraction of the raw reads (mean ∼42% of T cell clone reads).

To ensure the highest data quality for islet Ag-reactive CD4⁺ T cells, we sequentially examined the distribution of values for several unrelated quality control metrics (https://broadinstitute.github.io/picard/) and eliminated outlier libraries (cutoff values for elimination in parentheses): PF_ALIGNED_BASES (≤3.5e6); MEDIAN_CV_COVERAGE (<0.4 or >2.0); PCT_USABLE_BASES (<0.25); and MEDIAN_3PRIME_BIAS (log10 value+1 < 0.1 and >0.4). The combined use of these quality-control metric filters eliminated ∼24% of initial libraries (89/364 initial libraries eliminated). We also eliminated libraries that did not yield at least one in-frame rearranged TCR junction (∼10% of libraries that had passed quality-control metric filters). Altogether, the use of these conservative metrics led us to consider ∼68% of initial profiles (246/364) as having sufficient quality for subsequent analyses. Data from retained profiles were deposited in the Gene Expression Omnibus repository (Gene Expression Omnibus accession number: GSE96569; https://www.ncbi.nlm.nih.gov/geo/). We analyzed a total of 93, 35, and 25 cells from subjects T1D2, T1D4, and T1D5, and 37, 31, and 22 from subjects HC2, HC3, and HC5, respectively. The higher number of cells for subject T1D2 resulted from pooling of profiles from three different visits. Before analysis, counts were normalized using the trimmed mean of M values method (26) and transformed to reads per kilobase of transcript per million reads mapped (RPKM) values.

To determine the sequence of rearranged TCR sequences, which include nontemplated nucleotides in the CDR3 junction not present in the reference genome, we used methods for genome-independent (de novo) assembly to construct a set of overlapping DNA segments (or an assembly of overlapping sequence reads [contigs]) (27). In initial experiments, we prefiltered RNA-seq reads from the BRI-4.13 T cell clone to identify reads aligning (28) to TCR genes (29) and assembled them de novo into contigs (27). We found that each cell yielded ∼1000–3000 reads matching TCR genes, which could be assembled (27) into TCR contigs of ∼100–1000 bp in length. Submission of these contigs to IMGT/V-QUEST (30) identified productive TCR chain rearrangements. In subsequent experiments, we found that performing de novo assembly on total reads without TCR gene prefiltering gave very similar results; we consequently omitted the TCR gene prefiltering step in later experiments. Unique TCR chains for all cells were sequentially filtered for TRAV/TRBV gene usage (i.e., no TRDV or TRGV), productive rearrangements (i.e., no in-frame stop codons), and length (7–25 aa).

Differential gene expression was performed using the MAST R package (31). Linear models for gene expression contained terms for cellular detection rate (31), group (T1D or HC), or frequency of TCR sharing. X and Y chromosome genes were removed before differential gene expression comparisons. Protein–protein interactions (PPIs) were obtained from STRING (32) (http://string-db.org/) or GeneMANIA (33) (http://genemania.org/) and visualized using Cytoscape (34).

Statistical tests were performed using the R programming language and software environment. For continuous, normally distributed variables, we used t tests; for nonnormally distributed variables, Wilcoxon tests; and for categorical variables, Fisher exact test. One-sided tests were performed when testing whether a given parameter was larger than the value given by the null hypothesis. A two-sided test was used when the test was that a parameter was simply not equal to the value given by the null hypothesis (i.e., that the direction did not matter). A false discovery rate (FDR) (35) <0.1 was used to define differential gene expression. The specific test used to derive each p value is listed in the text.

For comparing fractions of cells sharing clonotypes, we devised a permutation-testing procedure. We originally used a down-sampling approach for comparing TCR frequencies (36) but found it difficult to run statistical tests on these down-sampled data, because sample sizes were insufficient for standard nonparametric tests and the distributions violated normality assumptions. Because of these problems, we devised an alternative permutation approach that estimated the probability of recovering differences as large as those actually observed if all patients had the same distribution of TCRs. We generated a single distribution of TCRs by pooling all the TCR sequences recovered from all patients. For each replicate, we drew simulated sets of TCRs from that pooled distribution that were equal in size to the actual samples we obtained from each patient. We then quantified the percentage of shared TCRs within each patient at each sharing threshold, calculated the mean percent across patients within each group (T1D or HC), and determined the between-group difference in mean percent shared. We repeated this process 1000 times to generate a distribution of expected between-group differences if all patients had the same TCR repertoire. If the observed differences between HC and T1D were due to unequal sample sizes, sampling error, or a combination of the two, the observed values should fall within this distribution. Significance values were calculated as the proportion of replicates for which the permuted difference was greater than or equal to the observed difference.

We hypothesized that single-cell RNA-seq profiling would allow parallel determination of both the rearranged TCR chains and the transcriptome phenotypes of individual T cells. We tested this hypothesis by comparing individual cell and bulk profiles from BRI-4.13 (19), an islet Ag-reactive CD4⁺ T cell clone from an individual with T1D (2Materials and Methods). Although single-cell transcript profiling has been successfully used with several cell types (37), less is known (38) about the performance of single-cell techniques with Ag-specific T cells, which contain very limited amounts of RNA. In our experiments, we detected nonlinearity between individual cell and bulk profiles in expression of low-abundance genes, indicating that genes expressed at low levels are less likely to be detected at the single-cell level than in a bulk measurement (Fig. 1A). Median expression in single cells was ∼0 for genes expressed at less than the top quartile of expression in bulk samples [log2 (RPKM + 1) ∼3.4 or ∼5 RPKM]. We calibrated our ability to detect expression of genes of different abundance (Fig. 1B): single-copy genes (median expression ∼2 RPKM) (39) were detected in ∼35% of cells; genes expressed at ∼8 RPKM (∼4 copies per cell) were detected in ∼50% of cells; and genes expressed at ∼115 RPKM (∼60 copies per cell) were detected in ∼90% of cells. Finally, we tested the consistency of expression in single-cell profiles for genes of different abundance in bulk samples (Fig. 1C). This revealed that low-abundance genes tended to show bimodal gene expression (40, 41) in single-cell profiles. For these genes, one mode was near zero, indicating a population of cells in which the gene was either not expressed or the transcript failed to be amplified during library construction (Fig. 1C).

We next tested our ability to recover rearranged TCR chains from individual cells of the BRI-4.13 T cell clone, where rearranged TRAV-CDR3-TRAJ and TRBV-CDR3-TRBJ sequences that bind antigenic peptide were known (42). Using our RNA-seq and TCR clonotype identification pipeline (2Materials and Methods), we identified the expected rearranged TCR chains, as well as a previously undescribed productively rearranged TRAV chain in 82–98% of cells (Fig. 1D). Most cells (∼74%) yielded all three sequences (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables).

To extend our approach to analyze primary Ag-specific T cells in peripheral blood, we isolated influenza Tmr-reactive CD8 T cells from a healthy subject and subjected them to single-cell profiling (2Materials and Methods; Supplemental Fig. 1A). We recovered rearranged TCRs from ∼76% of cells (34/45), of which TRAV and TRBV chains were found in 33% (15/45) and 42% (19/45) of cells, respectively, and both TRAV and TRBV chains were found in 18% (8/45) of cells (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables). The recovery of both TRAV and TRBV chains from the same cells was lower with influenza-reactive cells than with the BRI-4.13 T cell clone or islet Ag-reactive cells (see later), likely because we used frozen PBMCs for identification of influenza-reactive cells. Most of the recovered TCR sequences (∼88% or 30/34) were expanded (Supplemental Fig. 1B). These expanded rearranged TCRs shared sequence identity with previously identified immunodominant TCR chains for influenza (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables; Supplemental Fig. 1C) (43). Using flow cytometry of influenza Tmr-reactive cells, we confirmed the expansion of the TRBV19 gene segment predicted by single-cell RNA-seq (Supplemental Fig. 1D). Together, these results validate the sensitivity and specificity of our procedures for determining transcript profiles and TCR sequences from RNA-seq profiles of individual Ag-specific T cells.

To investigate the diversity of islet-specific CD4⁺ T cells in disease and health, we extended our methods to include comparisons of islet Ag-specific T cells in blood from HC and T1D individuals (Fig. 2). We relied on CD154 upregulation (44) to identify CD4⁺ T cells that became activated when pooled islet Ag peptides were added to PBMCs. We then isolated and sorted these activated cells into microfluidic chips using flow cytometry and subjected them to single-cell RNA-seq. We next processed RNA-seq reads along two parallel paths to identify rearranged TCR chains and elucidate transcript phenotypes. From these results, we identified paired TCR chains that were found in multiple individual cells (expanded), expressed them in recombinant form, and deconvoluted the islet Ag peptide pool to identify the specific antigenic peptides recognized (2Materials and Methods and Fig. 2).

We recruited a group of well-characterized high-risk DRB1*0401 T1D subjects and matched HC subjects (Table I). All T1D subjects were adult, within 2 y of diagnosis at their initial visit, and had detectable levels of C-peptide and multiple islet autoantibodies (Table I). HC subjects had no personal or family history of T1D. To identify islet Ag-reactive CD4⁺ memory T cells, we stimulated fresh PBMCs by exposure to a pool of DRB1*0401-restricted islet Ag peptides (Table II) and enriched for cells with upregulated CD154 using magnetic bead separation (Supplemental Fig. 2A). We then sorted enriched cells for those that had upregulated expression of the CD154 (44) and CD69 activation markers (2Materials and Methods; Supplemental Fig. 2B). We included two activation markers to increase the specificity of the CD154 assay. We chose stimulation times of 12–14 h because we found that this condition yielded improved downstream library quality. CFSE labeling and Ki67 staining experiments showed that relative to DMSO-treated (negative control) cells, <1% of islet Ag- or influenza-reactive cells proliferated under these conditions (Supplemental Fig. 2C). The islet peptides we used (Table II) represent consensus immunodominant epitopes recognized by CD4⁺ T cells in DRB1*0401 T1D subjects over many functional and epitope mapping studies (45–47). We also tested additional unpublished peptides identified using similar procedures, including several derived from ZNT8, a major islet autoantigen (48, 49). From these previous studies, we expected that most DRB1*0401 T1D subjects would have cells reactive with some, but not all, of the peptides used. After flow cytometry, the isolated cells (Supplemental Fig. 2B) were predominantly effector T cells, because regulatory T cells do not strongly upregulate CD154 under these experimental conditions (50) and we did not detect genes clearly related to regulatory T cells in subsequent transcript analyses.

If islet Ag-reactive T cells expand after encountering islet Ag(s), we expected to detect expanded TCR clonotypes of shared TRAV and TRBV chains. To test this prediction, we isolated CD4⁺ memory T cells from T1D and HC subjects after stimulation with islet Ag peptides and performed single-cell RNA-seq. We applied quality metrics to restrict RNA-seq analysis to 246 high-quality libraries (2Materials and Methods).

The frequencies of islet Ag-reactive CD4⁺ T cells detected in total or CD4⁺ memory T cell populations were similar between T1D and HC subjects (p ≥ 0.35, one-sided Wilcoxon test) (Fig. 3A). Likewise, similar numbers of cells were captured and passed quality filters between T1D and HC subjects at each visit (2Materials and Methods) (p > 0.48, two-sided t test). We then examined the distribution of rearranged CDR3 junctions, TRAV and/or TRBV, detected in individual islet Ag-reactive CD4⁺ memory T cells from T1D and HC subjects (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables). From N = 246 total cells, n = 165 (67%) contained both TRAV and TRBV chains, n = 30 (12%) contained TRAV chains only, and n = 51 (21%) contained TRBV chains only (https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables). Some cells (n = 21, 8.5%) contained two TRAV chains; fewer contained two TRBV chains (n = 9, 3.7%). Many rearranged CDR3 junctions, particularly from T1D subjects, were shared between cells from the same subject (Fig. 3B). The fraction of shared junctions within cells from the same participant was higher than sharing between subjects, where we detected no shared CDR3 junctions (p = 2.5e−9, Fisher exact test). This indicates that the islet CD4⁺ T cell responses were subject-specific, or “private.” Despite similarities in frequencies of islet Ag-reactive CD4⁺ T cells between T1D and HC subjects (Fig. 3A), fractions of cells sharing CDR3 junctions were higher in T1D than HC subjects (Fig. 3C). These differences were significant (by permutation testing) over a range of thresholds (2–8 cells) (Fig. 3C). To calculate these p values, we used permutation testing (2Materials and Methods). In accord with the observation of increased CDR3 junction sharing, we also found that values for Shannon entropy, a measure of clonal diversity (36), were higher for HC than T1D subjects (means of 4.3 versus 3.4, respectively; p = 0.04, one-sided Wilcoxon test). Shannon entropy values were calculated using down-sampling (36).

Within T1D subjects, we found identical rearranged TRAV and TRBV protein sequences for the most highly shared junctions (sequences 1–8, Table III) in 4–18 individual cells per subject (frequencies ∼0.10–0.47). From each subject, nearly all cells with one of these rearranged chains also contained the other rearranged chain, thereby permitting unambiguous determination of TRAV/TRBV pairing (Table III). Expanded TCR sequences in individual cells were identical at the amino acid (Table III) and nucleotide levels (data not shown) through the V, J, and D genes (for TRBV chains) and CDR3 regions. This indicates clonal rather than convergent origin of the expanded clonotype sequences. Subjects T1D2 and T1D4 yielded a single expanded clone, whereas T1D5 yielded two. In one T1D subject (T1D2), who was sampled three times, we observed extensive sharing of expanded rearranged TCR chains (sequences 1 and 2, Table III) over a period spanning >15 mo (Fig. 3D). Approximately 13% of junctions (12/93) were shared between visits 1 and 3, and one expanded junction was shared between all three visits by subject T1D2. The reduced TCR sharing at visit 2 may have been caused by, in part, lower cell yield at this visit. In comparison, there was no sharing of junctions (0/355) between any two T1D or HC subjects (p = 1.8e−6, Fisher exact test).

This demonstrates stability of a shared clonotype over time, a feature expected in cells relevant for disease progression. Taken together, these findings illustrate more extensive clonotype sharing among islet Ag-reactive CD4⁺ memory T cells in T1D than HC subjects. This suggests the possibility of in vivo clonal expansion of T cells with certain clonotypes, as the result of repeated encounters with Ag. The higher clonotype expansion in T1D subjects may indicate that such encounters are more frequent in T1D than HC subjects.

We then identified specific antigenic peptide(s) from the pool used to trigger T cell activation to clarify whether expanded TCRs from different individuals recognize the same or different islet Ags and/or epitopes. Our procedures involved isolation and characterization of islet Ag-reactive T cell clones and retroviral expression of recombinant TCR sequences.

We first identified the Ag recognized by the TCR clonotype expanded in subject T1D2 (sequences 1 and 2, Table III) (Fig. 4). We generated T cell clones from islet Ag-reactive CD4⁺ memory T cells from this subject and screened them by flow cytometry for expression of different TRBV genes. We found that ∼23% (11/47) of clones expressed TRBV6-6, the TRBV segment expanded in T1D2 (sequence 2). We selected five of these TRBV6-6⁺ clones for RNA-seq analysis and confirmed that they yielded rearranged TCR chains identical to sequences 1 and 2. We then tested these TRBV6-6⁺ clones for proliferation in response to pooled and individual peptides. All five TRBV6-6⁺ T cell clones proliferated in response to a pool of peptides from IGRP, but not to other pooled peptides (Fig. 4A, Table II). Testing against individual peptides showed that only IGRP 305–324 (QLYHFLQIPTHEEHLFYVLS, Table II) triggered specific proliferation (Fig. 4A). Consistent with this finding, the T cell clones expressing the expanded clonotype bound class II Tmr loaded with IGRP 305–324_, but not Tmrs loaded with an influenza HA peptide (51) or other IGRP peptides (Fig. 4B). These experiments demonstrate that the expanded TCR pair comprising sequences 1 and 2 (subject T1D2) recognizes IGRP 305–324 in the context of DRB1*0401 MHC class II molecules.

Although we successfully used T cell clones to elucidate the specificity of the expanded TCR clonotype from subject T1D2, there are clear drawbacks to this approach, including difficulties in access to patients and isolating T cell clones of some specificities. We therefore developed recombinant and ectopic retroviral expression methods to demonstrate the specificity of expanded TCR clonotypes (Fig. 5). We transduced primary human CD4⁺ T cells with retroviruses expressing recombinant rearranged TRAV and TRBV sequences expanded in subjects T1D2 and T1D4 (sequences 1 and 2 and sequences 3 and 4, respectively). We then tested transduced T cells for proliferation in response to islet peptides by flow cytometry using CFSE dye dilution. Because autologous APCs were not available for these experiments, we instead used Priess lymphoblastoid cells, which have been used to present T1D Ags in the context of DRB1*0401 class II molecules (23). Although these cells are DRB1*0401, *0401, they may have other HLA mismatches with patient cells (MHC class-1, etc.). To control for alloreactivity and other potential background issues, we compared proliferation in T cells expressing the recombinant TCR clonotype (17–29% of total cells), with untransduced T cells in the same culture, as a negative control. As shown in Fig. 5A, we observed minimal proliferation in either transduced or untransduced T cells in the absence of peptide, indicating that alloreactivity was not a major concern under these conditions. As expected from Fig. 4, we found that transduced T cells expressing the expanded TCR clonotype from subject T1D2 proliferated in response to the IGRP pool and IGRP 305–324 peptide, but not with other peptides (Fig. 5A). These results confirm that the expanded TCR clonotype from subject T1D2 specifically recognizes the IGRP 305–324 peptide. Surprisingly, however, transduced T cells expressing the expanded clonotype in subject T1D4 (sequences 3 and 4) also proliferated in response to the IGRP peptide pool, but not other peptide pools (Fig. 5B and data not shown). In contrast with subject T1D2, however, proliferation with the expanded clone from subject T1D4 was specific for peptide IGRP 241–260 (KWCANPDWIHIDTTPFAGLV; Table II), a different IGRP peptide than recognized by the clonotype from subject T1D2. 5KC cells expressing these same recombinant TCRs from subjects T1D2 and T1D4 specifically bound class II Tmrs loaded with IGRP 305–324 and IGRP 241–260, respectively, but not an influenza HA peptide (Fig. 5C, 5D). Thus, the expanded clones from two subjects with T1D, comprising sequences 1 and 2 and sequences 3 and 4, respectively, recognize distinct epitopes of the IGRP protein in the context of DRB1*0401.

To examine whether transcript profiles of individual cells differed by disease status, we performed global comparisons of transcript profiles from T1D and HC subjects using an unsupervised approach (principal component analysis [PCA]). Comparison across the top three principal components showed small group differences between T1D and HC cells (Fig. 6A, 6B). Moreover, there were no robustly differentially expressed genes between the groups (FDR < 0.1, with an estimable log fold change [i.e., flagged as not applicable (NA)]; https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables).

We reasoned that group differences between T1D and HC subjects might be obscured by heterogeneity at the cellular level within subjects and groups. To examine whether clonal expansion was associated with transcriptional heterogeneity within the T1D group, we focused on comparing T1D cells having expanded TCR sequences (T1D-E cells) with T1D cells having nonexpanded TCR sequences (T1D-NE). We selected T1D cells with TCRs shared in four or more cells (T1D-E4 cells), that is, cells expressing expanded clonotypes (sequences 1–8, Table III). We compared the distribution of T1D-E4 cells, T1D cells with TCRs shared in less than four cells (T1D-NE), and HC cells (Fig. 6C, 6D) by PCA. In principal component space, T1D-E4 profiles were shifted toward lower PC3 values than T1D-NE or HC profiles (Fig. 6C, 6D). We found even more pronounced differences (Fig. 6E, 6F) when repeating this analysis using T1D cells with TCRs shared in eight or more cells (T1D-E8 cells), indicating that the greatest differences in gene expression were among cells with the most expanded TCRs. These shifts in T1D-E4 and T1D-E8 versus HC profiles were highly significant (Fig. 6G).

We found more genes showing significant differences in expression (n = 62, FDR <0.1; https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables) when including frequency of TCR sharing as a term in a linear model for gene expression (31), again demonstrating that degree of expansion is correlated with variation in gene expression. Genes positively related to degree of clonotype expansion (up in T1D-E cells) were enriched in T cell activation and leukocyte differentiation genes (http://genemania.org/; FDR = 6.5e−3 for both terms). Genes negatively associated with clonotype expansion (down in T1D-E cells) were enriched in type I IFN signaling genes (FDR = 4.0e−10). Volcano plots showed that expression of selected Th2-related genes (e.g., GATA3, CCR4, IRF4) and genes involved in IFN responses (e.g., IFNG, CD69, GBP5) was higher and lower, respectively, in T1D-E cells (Fig. 7A). Projection of differentially expressed genes onto PPI networks showed significant interconnectedness of genes down and up in T1D-E cells (Fig. 7B). Other differentially expressed genes that were not represented in the PPI networks were also noted (Fig. 7B; https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables), but these have not been investigated further.

To explore the consistency of gene expression within and between subjects, we used PCA plots to examine transcript profiles from expanded and nonexpanded cells from each of the three T1D subjects individually, compared with all other cells in the data set (Fig. 7C–E). To these PCA plots we added vectors (biplots) for expression of individual differentially expressed genes (from https://github.com/linsleyp/Cerosaletti_Linsley/tree/master/Supplemental_Tables), including GATA3, a Th2 cell marker, and IFNG, a Th1 cell marker (Fig. 7C–E). These projections showed that T1D-E4 cells were shifted relative to T1D-NE cells in principal component space for each subject, especially for subject T1D4, and had different relationships with the GATA3 and IFNG gene expression biplots. T1D-E cells from subject T1D4 were more shifted in the dimension of the Th2 marker, suggesting that these cells were the primary source of the Th2-like genes differentially expressed at the group level. Taken together, these findings demonstrate that cells expressing expanded TCR clonotypes differ from cells with nonexpanded clonotypes, but that there is heterogeneity between donors for the extent of differences.

The presence of Ag-specific T cells in peripheral blood of T1D and HC subjects, even in cells with a memory phenotype, has been a surprising finding (7–9, 52). The widespread detection of these islet Ag-reactive T cells may result from their expression of TCRs that cross-react with pathogen-derived Ags, as has been reported for the islet Ags insulin and IGRP (53, 54). Attempts to identify disease-related differences in islet T cells have yielded inconsistent results, confounding efforts to use them to elucidate mechanisms of pathophysiology or as biomarkers of disease progression and therapeutic targets.

In this study, we show that novel features of islet Ag-reactive T cells from the peripheral blood of T1D subjects can be uncovered by utilizing the power of single-cell RNA-seq profiling to identify their TCR clonotypes in parallel with full transcript profiles. One of our key findings was the demonstration of expanded clones of islet Ag-reactive T cells, particularly in T1D subjects. Because we used 12- to 14-h activation times for our CD154 assays to increase RNA yield, we were concerned that some of the TCR sharing, or expansion, we observed resulted from proliferation in culture rather than in vivo. However, direct measurements showed negligible cell proliferation occurred under our culture conditions (Supplemental Fig. 2). Even if some cells did proliferate during culture, resulting in daughter cells sharing the same clonotypes, the length of the cell cycle in human lymphocytes (10–15 h) (55) would have precluded more than a single cell division during the culture period. At most, in vitro proliferation could account for no more than two daughter cells sharing the same clonotype in an experiment. To further mitigate this concern, we focused on more highly expanded clones (more than four cells), which were most likely to have resulted from in vivo expansion. The presence of more expanded clonotypes in cells from T1D compared with HC subjects links expansion of these cells to disease progression, possibly occurring during immune destruction of the pancreas. Although we have been able to analyze longitudinal visits from only a single subject thus far, our finding of the same rearranged TCR sequences in multiple visits from subject T1D2 suggests that expanded clonotypes can be stable over time. It is important to note that we have examined only adult T1D subjects in this study, primarily because of the blood volumes needed for our technology in its current form. It will be important to determine in future studies how the present results compare with results in children with T1D.

Another of our key findings was the demonstration of the specificity of expanded clones of islet Ag-reactive CD4⁺ T cells from different individuals for distinct peptides from the immunodominant islet protein IGRP. Although numerous T cell islet Ags and epitopes have been described in T1D (Table II), it remains unclear which of these are most important in disease progression. The CD154 assay (44) provides an opportunity to compare the relative frequencies of CD4⁺ T cells recognizing different published and unpublished Ags and epitopes in side-by-side testing in the same assays. In our studies, this “competitive” approach has highlighted the importance of IGRP as a target for expanded clones of islet Ag-reactive CD4⁺ T cells, including IGRP 305–324, a previously unpublished epitope.

IGRP is a metabolic enzyme (glucose-6-phosphatase 2) that is recognized as a major CD8⁺ T cell autoantigen for T1D in the NOD mouse model (56–59). IGRP also was recognized by T1D-related CD4⁺ T cells in mouse (60) and human studies (46). Our present studies showed the presence of expanded clonotypes of IGRP-specific CD4⁺ T cells in two T1D subjects where TCR specificity was established, suggesting a pathogenic role for IGRP-specific CD4⁺ T cells. One reason for the immunodominance of this islet protein may stem from activation of IGRP-specific T cells by molecular mimicry with microbial Ags (53, 54).

Intriguingly, we observed heterogeneity of transcriptional responses in islet Ag-reactive CD4⁺ T cells. Whereas group differences between T1D and HC cells were small, differences related to clonotype frequency were larger. We do not currently know how the in vitro stimulation in our assay affects transcriptome differences between T1D and HC cells. In fact, our activation conditions may have obscured subtle differences in cells from T1D and HC subjects (9). Although, to date, our attempts to use less stimulated cells for single-cell analysis have yielded poor results, future iterations of our technology, perhaps using Tmr staining under nonactivating conditions, may be useful in addressing this limitation. However, between individual subjects, transcript phenotypes of individual cells differed both qualitatively and quantitatively, demonstrating phenotypic heterogeneity. IGRP-reactive T1D-E4 cells from one subject (T1D4) that recognized peptide IGRP 241–260 had a more Th2-like phenotype. In contrast, T1D-E4 cells from other subjects, including cells from T1D2 that recognized IGRP 305–324, did not show this phenotype as clearly. It is also worth considering whether there is intraclonal heterogeneity in gene expression profiles. Consistent with this possibility, the PCA plots in Fig. 7C–E show cells that seem to cluster separately from other cells having the same clonotypes. Although the present studies were not powered sufficiently to conclusively demonstrate whether these outlier cells represent intraclonal heterogeneity, the possibility should be considered in future studies. Importantly, neither interclonal nor intraclonal heterogeneity would be apparent in bulk analyses of islet Ag-reactive T cells.

Although early studies on autoreactive T cell responses suggested Th1-type proinflammatory polarization in CD4⁺ T cells in T1D (12), other studies have indicated a more complex scenario (13). More recent studies showed different ratios of Th1, Th2, and T regulatory type 1 cells in IGRP-reactive CD4⁺ T cells from adult- and juvenile-onset T1D subjects (61). A role for Th2 cells in T1D has been suggested by increased levels of type 2 cytokines in the serum of T1D subjects (62) and by genetic and epigenetic fine mapping studies of causal autoimmune disease genetic variants (63). Taken together with previous studies, our findings challenge a simple classification of T1D as a Th1-mediated pathology. Our results suggest either the existence of different disease subtypes (61) or changes in disease over time that were not resolved in our study, which mostly were taken from single visits.

Our results demonstrate the power of single-cell RNA-seq profiling for simultaneously determining T cell clonotypes and linking these with expression profiles. When used in conjunction with technologies for isolating Ag-specific T cells, our methods allow an unprecedented view of specific T cells likely to be involved in pathogenic responses. In contrast with other methods (64, 65), our procedures for clonotype determination do not require the use of multiple sets of PCR primers for determining TCR sequences and couple the power of unbiased, full-transcriptome analysis with TCR clonotype determination. Similar procedures for linking TCR clonotypes and single-cell RNA-seq transcriptomes were recently published (38, 66), as was a computational method to infer the CDR3 sequences of tumor-infiltrating T cells in RNA-seq profiles from tumor profiles (67). Our procedures extend these previous studies by using short single-end reads, rather than longer paired-end reads, which reduces the cost of sequencing, and by eliminating initial TCR gene filtering, which reduces the number of steps in data processing. We also confirmed the accuracy of our approach by defining the specificity of two of the TCR clonotypes for individual islet peptides. One limitation of single-cell RNA-seq is that transcripts present at low-to-moderate abundance (i.e., ≤∼4–8 copies per cell in our studies) were not uniformly detected in transcriptomes of individual cells. Bulk RNA-seq is likely to be more useful for detection of genes with low expression, albeit at the price of averaging expression over all cells in the sample.

A key feature of single-cell RNA-seq data sets is that they are typically more powered for detecting gene expression differences between individual cells than between individual subjects. For example, our exploratory studies used 246 cells, but only six subjects, three T1D and HC subjects each. The limited number of subjects in our studies means that our findings, although significant (i.e., unlikely to have happened by chance), should be confirmed in more highly powered and/or differently designed studies. These studies should include expanded cross-sectional studies to verify the extent and specificity of expanded islet Ag-reactive CD4⁺ T cell clones in T1D and HC subjects, and longitudinal studies to verify the extent and stability of expanded clonotypes and to determine their relationship to disease progression. Further studies also will be needed to compare the extent of clonotype expansion in T1D versus other autoimmune conditions.

Our findings could have implications for the treatment of T1D. CD4⁺ T cells with expanded clonotypes may provide new biomarkers for disease progression and potential targets for Ag-specific therapies. Biomarkers and therapeutics involving islet Ag-reactive T cells will likely need to be individualized, because we observed that they had mostly unique or private sequences with distinctive specificities. T cells with expanded clonotypes may also provide new and better targets for immunotherapy than islet Ag-reactive T cells that have not expanded and are therefore less likely to be involved in the disease process. Extending our approaches to include more subjects and/or longitudinal studies may reveal how levels of T cells with expanded clonotypes change during disease progression, how their levels are modified during therapeutic intervention, and which specificities will provide the best therapeutic targets.

We acknowledge the National Institutes of Health Tetramer Core Facility for provision of the influenza Tmr. The National Institutes of Health Tetramer Core Facility was supported by contract HHSN272201300006C. We also acknowledge Jane Buckner for comments on the manuscript; Chester Ni, Mark Robinson, Masanao Yajima, and Greg Finak for bioinformatics assistance; Janice Chen, John P. McNevin, Kimberly O’Brien, and Quynh-Anh Nguyen for technical assistance; and Anne Hocking for assistance preparing the manuscript.

The authors have no financial conflicts of interest.