Characterization of Peripheral Blood TCR in Patients with Type 1 Diabetes Mellitus by BD RhapsodyTM VDJ CDR3 Assay

The sequence of complementarity-determining region 3 of the T-cell receptor (TCR) varies widely due to the insertion of random bases during V-(D)-J recombination. In this study, we used single-cell VDJ sequencing using the latest technology, BD Rhapsody, to identify the TCR sequences of autoreactive T-cells characteristic of Japanese type 1 diabetes mellitus (T1DM) and to clarify the pairing of TCR of peripheral blood mononuclear cells from four patients with T1DM at the single-cell level. The expression levels of the TCR alpha variable (TRAV) 17 and TRAV21 in T1DM patients were higher than those in healthy Japanese subjects. Furthermore, the Shannon index of CD8+ T cells and FOXP3+ cells in T1DM patients was lower than that of healthy subjects. The gene expression of PRF1, GZMH, ITGB2, NKG7, CTSW, and CST7 was increased, while the expression of CD4, CD7, CD5, HLA-A, CD27, and IL-32 was decreased in the CD8+ T cells of T1DM patients. The upregulated gene expression was IL4R and TNFRSF4 in FOXP3+ cells of T1DM patients. Overall, these findings demonstrate that TCR diversity and gene expression of CD8+ and FOXP3+ cells are different in patients with T1DM and healthy subjects.


Introduction
T-cell receptors (TCRs) contain a stationary region (C region) and a variable region (V region), and the structure of the variable region acquires remarkable diversity by reconfiguration of the TCR and immunoglobulin genes. In the case of TCRβ, TCRδ, and immunoglobulin H chains (heavy chains), three complementarity-determining regions called complementarity-determining region (CDR) 1, CDR2, and CDR3 are formed by the recombination of the V, D, and J gene segments. In both TCR and B-cell receptors (BCRs), CDR1 and CDR2 are encoded in the V segment, and the same gene sequences as in the germline are used. Conversely, the CDR3 sequence can have significant diversity, even when the same V, D, and J segments are used, due to the insertion of random bases (N sequences) during V-(D)-J reconstruction.
It is widely known that abnormalities in B-and T-cell repertoires are involved in the development and progression of both autoimmune and allergic diseases, but the details have not been fully elucidated. Recently, comprehensive immunosequencing has been used to analyze the pathogenesis of type 1 diabetes mellitus (T1DM), rheumatoid arthritis, systemic lupus erythematosus, and multiple sclerosis [1]. One research group conducted a comprehensive immunosequencing analysis of TCRs and BCRs using lymphocytes isolated

Study Design and Participants
The KAMOGAWA-DM cohort study is an ongoing prospective cohort study that was approved by the Ethics Committee of the Kyoto Prefectural University of Medicine in 2013 (Kyoto, Japan, RBMR-E-466) [7]. Informed consent was obtained from all patients involved in the KAMOGAWA-DM cohort study. We randomly selected 4 patients with T1DM who visited our diabetes outpatient clinic from April to May 2021. In addition, PBMCs were collected on the day of the visit, and the experiment was conducted on the same day using fresh specimens without cryopreservation. In addition, none of those four patients had any apparent infection during the study period. T1DM was diagnosed based on the criteria of the American Diabetes Association [8]. According to the recommendation of the Committee of Experts of the American Diabetes Association, T1DM was divided into type 1A diabetes (i.e., immune-mediated), type 1 B (i.e., other forms of diabetes with severe insulin deficiency but without proof of autoimmune etiology, also known as idiopathic) [9], and slowly progressive insulin-dependent diabetes mellitus (SPIDDM) at all participating institutions in this study [10][11][12][13].
Furthermore, we used the scRNA-seq data generated from a healthy subject in a previous study from the NCBI's Gene Expression Omnibus under accession number GSE150060 [14].

Data Collection
Information regarding patients' background demographics (i.e., age, sex, disease duration, and smoking habits) was gathered from their medical records. Blood pressure was measured in an outpatient clinic. After an overnight fast, venous blood samples were collected to measure fasting plasma levels of glucose, C-peptide, triglycerides, total cholesterol, high-density lipoprotein cholesterol, low-density lipoprotein cholesterol, creatinine, and uric acid. The hemoglobin A1c level was determined by high-performance liquid chromatography and presented herein using the National Glycohemoglobin Standardization Program unit.

BD Rhapsody Single Cell Analysis System
Heparin was added to the syringe when peripheral blood was collected. Peripheral blood mononuclear cells (PBMCs) were isolated by density gradient centrifugation, counted, and resuspended in 650 µL of cold sample buffer for loading on a BD Rhapsody cartridge. Targeted scRNA-seq with TCR analysis was performed using the BD Rhapsody Express system (BD Biosciences, Piscataway, NJ, USA). Cell capture and library preparation were performed using the BD Rhapsody Targeted mRNA and AbSeq Reagent Kit (BD Biosciences), according to the manufacturer's instructions. Briefly, 1 × 10 4 cells from each fresh blood sample were captured in a microwell plate with beads. This was followed by cell lysis, bead retrieval, cDNA synthesis, template switching, Klenow extension, and library preparation (a targeted gene library using a human T-cell expression panel and a TCR gene library) following the BD Rhapsody VDJ CDR3 protocol. The final pooled libraries were spiked with 20% PhiX control DNA to increase the sequence complexity and subsequently sequenced (75 bp × 225 bp paired-end) on a HiSeq X Ten sequencer (Illumina, San Diego, CA, USA). In this study, HLA typing was not tested because the number of cells in each sample was not sufficient.

Data Analyses
The FASTQ files obtained from the sequences were processed using the BD Rhapsody Targeted Analysis Pipeline with V(D)J processing incorporated (BD Biosciences) in the Seven Bridges Platform (https://www.sevenbridges.com/d, accecessed on 20 October 2021). First, low-quality read pairs were removed based on read length, average base quality score, and highest single-base frequency. High-quality R1 reads were analyzed to identify cell labels and unique molecular identifier (UMI) sequences. The high-quality R2 reads were aligned with the reference panel sequences (Human T cell Expression panel) and TCR gene segments from the international ImMunoGeneTics information system ® (IMGT.org) using the program Bowtie2. IGBlast was utilized for CDR3 determination. Reads with identical cell labels, identical UMI sequences, and identical genes were folded into a single molecule. The obtained counts were subjected to error correction algorithms (recursive substitution error correction (RSEC) and distribution-based error correction (DBEC)) developed by BD Biosciences. The DBEC-adjusted number of molecule data obtained from the Rhapsody pipeline was imported into SeqGeq version 1.6.0, and quality control was then performed to gate out cells that were significantly smaller and with low numbers of expressed genes (dead cells). Subsequently, dimensional reduction and unbiased clustering in SeqGeq were performed using the Seurat plug-in. Briefly, Seurat was set up to include all genes used, and the QC function, log normalization, and UMAP (uniform manifold approximation and projection) were used for dimensionality reduction. These plug-ins output data, including UMAP, lists of upregulated and downregulated genes, and annotation information, using the PBMC gene model. Further clustering analysis was completed with manual curation. Integration of the cluster information and TCR CDR3 information in each cell was performed using the VDJExploler plug-in of SeqGeq. As a parameter for the structural diversity of the genes, the Shannon index H was calculated [15].  PBMC dataset generated using v3 chemistry (7865 cells passing QC, average reads per cell for mRNA libraries: 35,433) and the 5k PBMC dataset generated using NextGEM chemistry (5527 cells passing QC, average reads per cell for mRNA libraries: 30,853). See https://support.10xgenomics.com/single-cell-gene-expression/datasets/, accessed on 20 October 2021). CD8 + or FOXP3 + T cells were defined as cells expressing one or more copies of CD8 or FOXP3 (UMI).

Recombination of TRAV and TRAJ
The genetic recombination of 41 TRAVs and 50 TRAJs (excluding pseudogenes, ORFs, and low-expressed genes) resulted in a total of 2050 AV-AJ recombinants. Notably, the AV1-1-AV6 gene did not preferentially combine with the AJ50-AJ58 gene, and similarly, little recombination was observed between the AV35-AV41 gene and AJ3-AJ16. For TRB, 650 genetic recombinations occurred in 50 BV genes (excluding 11 pseudogenes and 5 ORFs) and 13 BJ genes (excluding pseudogenes). There were no restrictions on the combination of TRBV and TRBJ, as observed in the TRA. The Shannon index H' was used as a diversity index to evaluate the diversity of TRA and TRB. Shannon-index H' of TRA and TRB in S1 was 10.80 and 10.83, that in S2 was 11.62 and 11.69, that in S3 was 10.26 and 10.57, and that in S4 was 11.37 and 11.25, respectively ( Table 2).

TCR Clonotypes of CD8 + T Cells and FOXP3 + T Cells
Next, we investigated the genetic recombination of TRAVs and TRAJs and that of TRBVs and TRBJs in CD8 + T cells and FOXP3 + T cells, and we have shown the top ten TCR clonotypes in Tables 3 and 4. The most frequently observed CDR3 of TRA and TRB in CD8 + T cells in S1 were AGAISNNDMR and ASSVVGSGTDEQF; those of S2 were VVRARPPLPWSGGGADGLT and ASTPPSSPGYEQY, those of S3 were AFSGGYQKVT and ASSLAGEGSGTGELF, and those of S4 were VVSAFFSGGSYIPT and ASSSSRDRGNYEQY (Table 3).

TCR Clonotypes of CD8 + T cells and FOXP3 + T Cells
Next, we investigated the genetic recombination of TRAVs and TRAJs and that of TRBVs and TRBJs in CD8 + T cells and FOXP3 + T cells, and we have shown the top ten TCR clonotypes in Tables 3 and 4. The most frequently observed CDR3 of TRA and TRB in CD8 + T cells in S1 were AGAISNNDMR and ASSVVGSGTDEQF; those of S2 were VVRARPPLPWSGGGADGLT and ASTPPSSPGYEQY, those of S3 were AFSGGYQKVT and ASSLAGEGSGTGELF, and those of S4 were VVSAFFSGGSYIPT and ASSSSRDRGNYEQY (Table 3). In the FOXP3 + T cells of S1 were AMRFKSGYNKLI and ASSPPTSGASYEQY; those of S2 were ALSSNDYKLS and ASTLDGPGSPLH, those of S3 were GFSSGSARQLT and ASSFGRYEQY, and those of S4 were AAGRGNNRLA and ASSRTGGGYGYT (Table 4).

T-Cells in Blood of Patients with T1DM Have Phenotypic Hallmarks
Next, we performed an unbiased analysis of gene expression using Seurat and identified T-cell clusters in four T1DM patients ( Figure S1). Heatmaps of gene expression in each cluster are shown, with 10 clusters in S1, 11 clusters in S2, 19 clusters in S3, and 11 clusters in S4 (Figure 7).

Discussion
We used adapter ligation-mediated PCR, a bias-free PCR technique, for TCR repertoire analysis using NGS. This method uses a single set of primers to avoid PCR bias due

Discussion
We used adapter ligation-mediated PCR, a bias-free PCR technique, for TCR repertoire analysis using NGS. This method uses a single set of primers to avoid PCR bias due to primer competition. This method is, therefore, suitable for accurately estimating the abundance of each TCR gene in a wide variety of samples. In the present study, we comprehensively investigated the TRA and TRB repertoires from four patients with T1DM at the clonal level and evaluated a large amount of sequence data. This is the first study to reveal the TRA and TRB repertoires of patients with T1DM using BD Rapsody. Moreover, this integrated analysis makes it easy to detect the preferential use of specific TRVs and TRJs, which may be useful in studying immune responses by antigen-specific T-cells.
There are several single-cell sequencing platforms that have been widely used around the world in recent years. BD Rhapsody, which was used in this study, uses microwells and magnetic beads to isolate cells and perform bead-based 3 RNA-Seq, thus realizing highly accurate sample preparation with low doublets and cross-contamination rates. Compared to the widely used 10× Genomics, the library cost is lower, and cell viability demand does not need to be as high as 50%, cell loading can be up to 40,000 cells, and the frequency of doublets is lower [17]. Therefore, information on a large number of cells can be obtained.
The expression levels of TRAV17, a variable gene known to be enriched in a population of CD1b-restricted T-cells [18,19], and TRAV21 in patients with T1DM were higher than those in healthy Japanese subjects investigated in a previous study [20]. Conversely, the expression levels of TRAJ in patients with T1DM were not different from those of healthy subjects. The majority of TRBV in the healthy Japanese subjects was TRBV 29-1, whereas, in the patients with T1DM, there was no clear majority compared to the healthy subjects. The Shannon indices of TRA and TRB in healthy subjects in the previous report were both approximately 7, which were clearly smaller than those of the patients with T1DM observed in this study. Moreover, we surveyed the diversity index of CD8 + and FOXP3 + T cells, and the Shannon index of the cells in the patients with T1DM was lower than those of healthy subjects [21,22]. TCR variability of Tregs has been proposed to be beneficial in the maintenance of self-tolerance [6]; therefore, the findings in this study indicate that the reduced TCR diversity in Tregs of patients with T1DM in this study may indicate reduced immune tolerance in patients with T1DM.
Upregulated genes in CD8 + T cells of T1DM patients included cytotoxicity-associated genes, such as PRF1, GZMH, ITGB2, NKG7, CTSW, and CST7, whereas the expression of CD4, CD7, CD5, HLA-A, CD27, and IL-32 was downregulated. Cytotoxic CD8 + T cells are considered to be the primary mediators of β-cell injury, based on the predominance of CD8 T-cells in pancreatic islet infiltration [23,24], as well as numerous studies using animal models of T1DM caused by β-cell injury by CD8 T-cells [25,26]. In addition, several human studies have reported an expanded pool of memory T-cells in the peripheral blood of patients with type 1 diabetes [27] and resistance of effector T-cells to Treg suppression, and our results are consistent with these previous reports [28].
Upregulated genes in FOXP3 + T cells in T1DM patients included IL4R. Interleukin 4 (IL-4) has been reported to be involved in several signaling pathways in the regulation of Treg cell development and function [29][30][31][32][33]. IL-4 is a cytokine that defines the type 2 immune response, while IL-4 receptor alpha (IL-4 Rα) suppresses Treg cell function during type 2 disease [34,35]. Recent reports have shown that enhanced IL-4Rα signaling by gainof-function mutations [32,35] or chronic type 2 inflammation [36] drastically reduces the number of Foxp3 + Treg cells, impairs the suppressive function of Treg cells, and promotes their reprogramming to T helper 2 (Th2)-like or T helper 17 (Th17)-like cells. This receptor is further thought to play a role in suppressing Treg cell function. Although no difference in the frequency of Tregs in peripheral blood isolated from T1D patients has been reported, defects in the phenotype and suppressive capacity of Tregs have been reported [37][38][39][40][41]. In this study, we used single-cell sequencing for the first time in the world to reveal that IL4R expression is upregulated in the Tregs of patients with type 1 diabetes. As previously reported, this result suggests that the function of Tregs in T1DM is impaired. In addition, TNFRSF4 was upregulated in FOXP3 + T cells in patients with T1DM. TNFRSF4 is one of the most highly expressed genes in Tregs [42,43]. In addition, TNFRSF4 is one of the most highly expressed genes in tumor-invasive Tregs compared to those in healthy tissues [44,45]. While there have been no reports on the expression of TNFRSF4 in relation to Tregs in patients with type 1 diabetes, this has been reported to be significantly increased in patients with relapsed acute myeloid leukemia compared to healthy donors [46]. TNFRSF4 mediates TRAF2 and TRAF5 to activate the NF-κB pathway of TNFRSF4, and the PI3K/PKB and NFAT pathways have also been identified [47,48]. The most important function of TNFRSF4 is to promote T-cell division, proliferation, survival, and cytokine production by activating the aforementioned pathways. In this study, we found that the expression of TNFRSF4 was upregulated in T1DM regulatory T-cells, suggesting that these cells may have some immune abnormalities.
The strength of this study is that it is the first to use the BD Rhapsody system, a stateof-the-art technology for single-cell sequencing of TCR repatriation and gene expression in peripheral blood T-cells from patients with T1DM. However, this study has several limitations. First, we analyzed CD8-and FOXP3 + T cells as cells with more than one read, and Tregs as FOXP3 + T cells could also be activated T-cells and not Tregs; it would have been more accurate if we sorted each positive cell by cell sorter and performed the same analysis. Second, peripheral blood of healthy subjects was used as a control, but this data was obtained by another research group and was not analyzed simultaneously in this study. Therefore, we should prepare our own samples for future studies. Third, the subjects in this study have had diabetes for many years and may not have autoreactive T-cells in the peripheral blood collected. On the other hand, Tregs have been reported to suppress GAD-responsive T-cells in patients with type 1 diabetes who have had the disease for more than 5 years [49]. Therefore, it has been reported that Tregs suppress GAD-reactive T-cells in patients with type 1 diabetes mellitus more than 5 years after onset, and it is possible that some immune abnormalities may still be present in PBMCs over time. Finally, we randomly selected patients who visited an outpatient clinic within a limited time. Therefore, we did not match background factors such as age, gender, duration of disease, and diabetes type, which are limitations of this study.

Conclusions
In conclusion, in this study, we used the latest technology, BD Rhapsody, to analyze the pairing of α and β chains that constitute the TCR of PBMCs from patients with type 1 diabetes at the single-cell level. In this study, we identified genes that are upregulated in T-cells as well as TRB repairs. scRNA-seq has greatly improved our understanding of heterogeneity in various biological processes and has led to significant breakthroughs in the fields of immunology, oncology, and developmental biology.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cells11101623/s1, Table S1: Clinical characteristics of the study patients, Table S2: Shannonindex H', Table S3: TCR clonotype of CD8 + cells, Table S4: TCR clonotype of FOXP3 + cells, Table S5: Top five CDR3 motif of TRA in CD8 + cells and clustering, Table S6: Top five CDR3 motif of TRB in CD8 + cells and clustering, Table S7: Top five CDR3 motif of TRA in FOXP3 + cells and clustering, Table S8: Top five CDR3 motif of TRB in FOXP3+ cells and clustering, Figure S1: nFeatures_RNA and nCount_RNA of Seurat object, Figure S2: Motifs of CDR and clustering, Figure S3: The heatmaps of gene expressions in CD8 + cells and FOXP3 + cells, Figure S4: Seurat analyses of FOXP3 + cells, Figure S5: The top five CDR3 motifs of TRA and TRB.