Early Detection of Peripheral Blood Cell Signature in Children Developing β-Cell Autoimmunity at a Young Age


Introduction

Family and sibling studies in type 1 diabetes (T1D) have implicated a firm genetic predisposition to a locus containing HLA class I and class II genes on chromosome 6 suggesting a role for CD4+ as well as CD8+ T cells in T1D pathogenesis (1–3). As much as 30–50% of the genetic risk is conferred by HLA class II molecules, which are crucial in antigen presentation to CD4+ T cells. Further, CD4+ cells reactive to β-cell antigen peptides are found in peripheral blood and the pancreas and typically secrete the cytokine IFNγ (4,5). CD4+ cells orchestrate adaptive immune responses, including that of antibody-secreting B cells as well as cytotoxic CD8+ T cells. Indeed, circulating autoantibodies against β-cell antigens may appear years before the clinical onset. Further, a cytolytic CD4+ subtype might directly contribute to target cell killing (6).

Although HLA class II is associated with the development of autoantibodies, HLA class I seems to be more strongly linked to disease progression (7). Histological analysis of pancreatic sections of cadaveric donors with T1D revealed that HLA class I is highly expressed in islets (8,9). Moreover, CD8+ cells are the most abundant cell type during insulitis (10), and the islets contain CD8+ cells specific for T1D autoantigens (11). Thus, the autoimmune cascade in T1D might be initiated by self-reactive CD4+ cells that activate B cells to produce autoantibodies that target the β-cells and unleash the cytotoxic activity of the autoreactive CD8+ cells. The environmental factors triggering and driving the autoimmunity in T1D are poorly defined, but the disease has been associated with viral infections (12), diet in early childhood (13), and reduced diversity of gut microbiota (14).

Currently, the appearance of T1D-associated autoantibodies is the first and only measurable parameter to predict progression toward T1D in genetically susceptible individuals. Although the disease progression rate varies considerably, children with genetic HLA risk expressing at least two T1D autoantibodies will very likely progress to clinical disease during the next 15 years (15). However, autoantibodies are poor prognostic markers for the timing of the clinical presentation of T1D. The appearance of autoantibodies indicates an active autoimmune reaction, wherein the immune tolerance is already broken. Therefore, there is a clear and urgent need for new biomarkers that predict the onset of the autoimmune reaction preceding autoantibody positivity or reflect progressive β-cell destruction. Such markers would present a window for early intervention aimed at complete disease prevention. Previously, we reported changes in whole-blood transcripts and serum proteins before the detection of diabetes-associated antibodies in children who later progressed to T1D (16,17). Therefore, we hypothesized that a comprehensive analysis of the transcriptome of longitudinal cellular samples including CD4+ and CD8+ T cells will lead to the identification of new early biomarkers.

Research Design and Methods

Study Cohort

Samples were collected as part of the DIABIMMUNE study from Finnish (n = 10) and Estonian (n = 4) participants (Supplementary Table 1). The HLA-DR-DQ genotypes were analyzed as previously described (18). A total of 836 children with HLA-DR-DQ risk allele were monitored and sampled at 3, 6, 12, 18, 24, and 36 months of age. The study protocols were approved by the ethics committees of the participating hospitals, and the parents gave written informed consent. Autoantibodies against insulin (IAA), glutamic acid decarboxylase (GADA), islet antigen-2 (IA-2A), and zinc transporter 8 (ZnT8A) were measured from serum with specific radiobinding assays (19). Islet cell antibodies (ICAs) were analyzed with immunofluorescence in autoantibody-positive subjects. The cutoff values were based on the 99th percentile in children without diabetes, which were 2.80 relative units (RU) for IAA, 5.36 RU for GADA, 0.78 RU for IA-2A, and 0.61 RU for ZnT8A. The detection limit in the ICA assay was 2.5 Juvenile Diabetes Foundation units (JDFU). A sample was considered seropositive when any of the autoantibodies exceeded the thresholds.

Sample Collections

At each study visit, 8 mL blood was drawn in sodium-heparin tubes (368480, Vacutainer; BD Biosciences). Peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll-Paque centrifugation (17-1440-03; GE Healthcare) and were suspended in RPMI-1640 medium (42401-018; Gibco) supplemented with 10% DMSO (cat. no. 0231, 500 mL, Thermo Fisher Scientific), 5% Human AB Serum (cat. no. IPLA-SERAB-OTC; Innovative Research), 2 mmol/L L-glutamine (G7513; Sigma-Aldrich), and 25 mmol/L gentamicin (G-1397; Sigma-Aldrich). After overnight incubation at −80°C, samples were stored in liquid nitrogen (−180°C). For fractionation, PBMC samples were thawed quickly in a 37°C water bath and quantitated for cell numbers and viability. On average, 90% of cells were viable. Magnetic antibody-coupled beads were used for sequential positive enrichment of CD4+ and CD8+ cells (11331D and 11333D; Invitrogen). RNA was isolated from the samples with an AllPrep kit (80224; QIAGEN), and quantity and quality were determined using a Qubit RNA assay (Q32852; Invitrogen) and Bioanalyzer 2100 (Agilent), respectively.

Bulk RNA Sequencing of PBMCs and Other Fractions

At least 80 ng total RNA was processed for RNA sequencing (RNA-seq) with the TruSeq Stranded mRNA Library Prep kit (RS-122-2101; Illumina). The sequencing was carried out with the Illumina HiSeq2500 instrument using TruSeq v3 (2 × 100 base pairs [bp] chemistry). The average sequencing depth was ∼51 million reads. Quality control was performed using FastQC (version 0.10.0). All the samples passed the quality criteria. The reads were aligned to the human reference transcriptome, GRCh37 assembly version 75, using TopHat (version 2.0.10) (20). Average mapping percentage was 93. The concordant pairs percentage was ∼89. The aligned reads were counted with htseq-count (HTSEq. 0.6.1; overlap mode of “intersection-strict”) (21). The read counts of genes were normalized using the trimmed means of the M values (TMM implemented in edgeR [22]). Coding, noncoding information were taken from Ensembl. Differential expression analyses were conducted separately for coding and noncoding genes, using edgeR (22). The variance of the data was estimated using the trended dispersion method. A further filtering step retained only those genes as differentially expressed (DE) that had |median log2 fold change (FC)| >0.5 and had >65% samples across all individuals regulated in the same direction (i.e., up- or downregulated). These filtering steps were added to discard false positives that may arise due to the heterogeneity of the samples resulting from normal variation, which is unrelated to T1D as well as to discard the outliers. A flowchart of the scheme of analysis has been shown in Supplementary Fig. 1.

Single-Cell RNA-seq

The concentrations of the PBMC samples varied from 0.55 to 1.80 × 106 cells/mL. From each sample, we aimed at the recovery of 5,000 single cells, loading ∼9,000 cells on the Chromium Controller using Single Cell 3′ Solution v2 reagents and following the manufacturer’s instructions (CG00052, Rev B; 10x Genomics). Single-cell RNA-seq (scRNA-seq) sample processing was carried out in three batches on consecutive days using the same lot of reagents and chips for all samples. The cDNA was further amplified using a Veriti Thermal Cycler (Applied Biosystems/Thermo Fisher Scientific), followed by cleanup (SPRIselect kit; Beckman Coulter). Finally, enzymatic fragmentation, end repair, A-tailing, adaptor ligation, and PCR were performed to produce indexed libraries, which were sequenced with Illumina HiSeq 3000 (one sample per lane) using paired end sequencing and 26 + 98 bp read-length configuration. The data were processed using the Cell Ranger pipeline, version 2.0.0, yielding on average 2,546 viable cells per sample and 114,309 reads per cell.

The reads were aligned to the human reference genome (hg19) using STAR (23). The mean raw reads per cell varied from 57,000 to 200,000. Quality control analysis and further exploration were done using Seurat (24). After filtering steps, 18,396 cells expressing 20,830 genes were retained. For details on the filtering steps please see Supplementary Data. The data were normalized using Seurat’s default. Highly variable genes were selected for principal component analysis. The top 20 PCs were used in the graph-based clustering. To identify marker genes for each cluster, cells of a single cluster were compared with the cells of all other clusters combined. A gene was considered a marker of a cluster if it was expressed in at least 25% of the cells of either of the two groups and the logFC between the cluster and all other clusters was at least 0.25.

For trajectory analysis, the pooled cells were ordered in pseudotime (i.e., placed along a trajectory corresponding with a type of biological transition, such as differentiation) using Monocle 2 (25). The analysis was performed on cells specifically from CD4+ and CD8+ T-cell clusters. For the details on the trajectory analyses, please see Supplementary Data.

RT-PCR Analysis

For PBMC samples, 50 ng total RNA was treated with DNaseI (Invitrogen), and cDNA was synthesized with the Transcriptor First Strand cDNA Synthesis Kit (Roche Life Science). For isoform-specific (IL32α, β, and γ) assay, quantitative PCR (qPCR) analysis was performed in triplicate runs using SYBR Select Master Mix (Applied Biosystems). ΔCt values was calculated relative to EF1α. For CD4+ T cells and pancreatic islets, RNA was isolated using the RNeasy Mini Kit (74106; QIAGEN) and RNeasy Plus Mini Kit (74134; QIAGEN), respectively. Purified RNA was treated with DNaseI and cDNA was synthesized with SuperScript II Reverse Transcriptase (18064014; Invitrogen). For the detection of global IL32, qPCR reactions were run using a custom TaqMan Gene Expression Assay reagent (AJ5IQA9; Thermo Fisher Scientific) in duplicate and in two separate runs. ΔCt values were calculated relative to GAPDH. The amplification was monitored with the QuantStudio 12K Flex Real-Time PCR System under the following PCR conditions: 10 min at 95°C, followed by 40 cycles of 15 s at 95°C and 60 s at 60°C, and analysis with QuantStudio software on Thermo Fisher Cloud.

For EndoC-βH1 cell data, cDNA was synthesized using the Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) as per the manufacturer’s recommendations. All reactions were performed in duplicates on at least three biological replicates. Cyclophilin-A was used as an endogenous control. Primer sequences are presented in Supplementary Table 2.

Elisa

To measure secreted IL32 levels, we used IL-32 DuoSet ELISA (DY3040-05 and DY008; R&D Systems) following the manufacturer’s instructions.

Intracellular Staining and Flow Cytometry

The cells were fixed for 10 min in Fix Buffer I (557870; BD Biosciences), followed by 45 min permeabilization using ice-cold permeabilization buffer III (558050; BD Biosciences). The cells were stained using allophycocyanin-conjugated IL32⍺ antibody (IC30402A; R&D Systems) and FITC-conjugated IFNγ antibody (MHCIFG01; Invitrogen) in PBS containing 0.5% FCS. The data were acquired in BD Fortessa and analyzed using FlowJo (version 10.4.2).

EndoC-βH1 Cell Culture

The EndoC-βH1 human β-cell line was obtained from Univercell Biosolution S.A.S., Toulouse, France. The cells were cultured as previously described (26). EndoC-βH1 cells were stimulated with IL32γ either alone (100 ng/mL; R&D Systems) or in combination with a cocktail of IL1β (5 ng/mL; R&D Systems) and IFNγ (50 ng/mL; R&D Systems) for 24 h. RNA samples were collected at the end of each treatment and analyzed by RT-qPCR.

Human CD4 T-Cell Isolation and Culturing

CD4+ T cells were isolated from cord blood collected from neonates born in Turku University Hospital and were cultured in Iscove’s modified Dulbecco’s medium containing 1% AB serum in absence (Th0) or presence (Th1) of 2.5 ng/mL IL12 (R&D Systems). Cells were activated with plate-bound CD3 (0.5 μg/well of a 24 well-plate) and soluble CD28 (0.5 μg/mL), both from Immunotech, with or without 50 ng/mL rIL32γ (R&D Systems). 12 ng/mL IL2 was added at 48 h. For IFNγ neutralization, anti-IFNγ antibody (MAB285, 10 μg/mL; R&D Systems) was used. For reactivation, cells were treated with 5 ng/mL phorbol myristic acid (Calbiochem) and 0.5 pg/mL ionomycin (Sigma-Aldrich) for 5 h.

Human Pancreatic Islets and Their Infection With Coxsackie B Virus

Human islets were isolated from pancreases obtained from brain dead organ donors and purified by handpicking to a purity of >90%. Islet culturing and virus infection with Coxsackie B virus-1 (CBV-1-7-10796 [CBV-1-7]) was performed as previously described (27). Islets were collected at the day 4 time point, and RNA was extracted using the RNeasy Plus Mini Kit or the AllPrep DNA/RNA Mini Kit (QIAGEN). For RNA-seq, 100 ng total RNA from three donors was used for library preparation according to Illumina TruSeq RNA Sample Preparation v2 Guide (part no. 15026495). The high quality of the libraries was confirmed with the Agilent Bioanalyzer 2100 and Qubit Fluorometric Quantitation (Life Technologies). The libraries were pooled in two pools and run in two lanes on the Illumina HiSeq 2500 instrument using 2 × 100 bp.

Data and Resource Availability

All the raw data will be deposited to European Genome-phenome Archive for access. The study does not involve any noncommerical reagents or tools.

Results

Fractionation of PBMC Sample Into CD4+, CD8+, and CD4CD8 Cellular Subsets Reveals Distinct and Overlapping Gene Expression Signatures

We performed RNA-seq of 306 longitudinal samples including unfractionated PBMCs, as well as CD4-enriched (CD4+), CD8-enriched (CD8+), and CD4 and CD8 cell–depleted (CD4CD8) cell fractions from seven case-control pairs (Table 1). The seven case children who developed T1D-related autoantibodies (Aab+) were selected from the DIABIMMUNE birth cohort (18), where HLA-susceptible children are sampled at 3–36 months of age (Fig. 1A). All seven children developed T1D-associated autoantibodies by the age of 2 years (Table 1), and four of them developed clinical T1D between the ages of 2.4 and 3.7 years. For each case subject, an autoantibody-negative control child was matched for sex, date and place of birth, and HLA-conferred risk category.

Table 1

Summary of the case and control children sampled at the age of 3–36 months

Figure 1
Figure 1

Fractionation of PBMC sample into CD4+, CD8+, and CD4CD8 cellular subsets reveals distinct and overlapping gene expression signatures. A: Outline of the sample collection and cell fractionation. B: t-SNE (t-distributed stochastic neighbor embedding) visualization of the log2-transformed expression data (without any filtering steps) colored according to cell fraction information. C: Number of DE genes when CD4+, CD8+, and CD4CD8 fractionated samples were compared with their original PBMC aliquots. The functionally important fraction-specific upregulated genes are highlighted in red. Analysis was restricted to healthy control subjects only. For the gene lists, see Supplementary Table 3.

The samples clustered according to the cell fraction (Fig. 1B), and the clustering was not affected by case-control status or sampling age, indicating that cell fraction–specific differences dominated over variation derived from other factors (Supplementary Fig. 2A and B). When CD4+, CD8+, and CD4CD8 samples from control subjects were compared with the unfractionated PBMC samples (also referred to as a fraction henceforth), 889, 399, and 1,002 genes were DE specifically in CD4+ (e.g., CD28, CTLA4), CD8+ (e.g., CD8A, CD8B, KLRK1), and CD4CD8 (e.g., IL1A, IL1B, IL6) fractions, respectively (Fig. 1C and Supplementary Table 3). CD4+ and CD8+ fractions shared 1,815 DE genes, of which 1,803 genes (99%) were concordant (either up or down in both fractions) (Supplementary Fig. 2C and Supplementary Table 3). In summary, fractionation of the PBMC population based on the T-cell phenotype allowed improved detection of DE genes and enabled identification of cell subset–specific gene expression signatures.

RNA-seq Analysis Identifies Transcriptomic Changes Associated With β-Cell Autoimmunity

Comparison of case samples with their respective controls identified 51, 69, 143, and 85 genes as DE (false discovery rate <0.05) in CD4+, CD8+, CD4CD8, and PBMC fractions, respectively (Supplementary Table 4), with a total of 278 unique DE genes in one or more fractions (Fig. 2A). Six genes, AMICA1, BTN3A2, IL32, RPSAP15, RPSAP58, and WASH7P, were upregulated in the case subjects in all four fractions (Fig. 2A). Only 16% of the DE genes have previously been reported as DE in genetically susceptible children with prediabetes, using microarrays (16,28,29) or RT-PCR (30–32), confirming dysregulation of these genes in children progressing to T1D. Besides protein-coding genes, 54 noncoding genes, including 3 antisense, 2 sense intronic, 7 enhancer, and 18 promoter-associated lncRNAs, were DE. To our knowledge, none of these lncRNAs have been linked to the etiology of T1D (16,28–32).

Figure 2
Figure 2

RNA-seq analysis identifies transcriptomic changes associated with β-cell autoimmunity. A: Number and overlap of DE genes between case and control subjects identified in cell fractions analyzed. Genes shared between all four fractions are highlighted. B: Heat map of the genes DE in CD4+ T cells between the case and control subjects. Values are presented as log2FC (truncated between [−2, 2]) between each case-control pair at each time point (3–36 months) and standardized to the mean of each gene. Genes coregulated with IL32 (<2.5 Euclidean distance) are marked with red box and text. Sample collection time with respect to (w.r.t) seroconversion, sample pairing information, and clinical status have been indicated with colors on top of the heat map. “Before/After SC” informs whether the case sample was collected before or after seroconversion. “Pair Info” provides the case-control pair information. The “SC / T1D” annotation indicates whether the case subject has progressed to clinical T1D diagnosis (T1D) or not (SC). C: Number and overlap of IL32 coclustered genes in indicated cell fractions. Genes regulated at least in two fractions are highlighted. D: Profiles of IL32, AMICA1, and BNT3A2 in CD4+ samples, presented in log2 RPKM (reads per kilo base per million mapped reads) scale. For individual profiles, see Supplementary Fig. 4. The case-control pairs are grouped according to the diagnosis of the case subjects. T1D, case subject has been diagnosed with clinical T1D; SC, case has seroconverted to autoantibody positivity.

Hierarchical Clustering Identifies Coregulated Gene Expression Clusters Associated With T1D Autoimmunity

Gene- and sample-wise hierarchical clustering for each cell fraction, including PBMCs, identified a cluster upregulated in the case samples in all four fractions (Fig. 2B and Supplementary Fig. 3AD). Interestingly, this cluster consistently contained IL32 and BTN3A2, along with other fraction-specific genes (Fig. 2C). In the CD8+ fraction, expression of a distinct cluster, including IFNG, was lower in most of the case samples than in control samples (Supplementary Fig. 3B). Surprisingly, in the PBMC fraction, we detected case-specific upregulation of a cluster, including insulin (INS), glucagon (CGC), and regulin 1α (REG1A), transcripts (Supplementary Fig. 3D), which are predominantly expressed in the pancreas.

To explicitly define coregulated genes in these clusters, we calculated Euclidean distances for IL32 (in each fraction), IFNG (in CD8+ fraction), and INS (in PBMC fraction) and considered the genes with a median Euclidean distance <2.5 across all case-control pairs to be coclustering with the gene of interest (Supplementary Table 5). In three of the four fractions, the IL32 cluster included BTN3A2, AMICA1, LARS, and RSU1 (Fig. 2C). IL32, AMICA1, and BNT3A2 show concerted gene expression profiles in CD4+ samples (Fig. 2D). In at least two of four fractions, this cluster also comprised TRBV4-1, TMEM14C, UROS, WASH7P, BTN3A3, CARD8, CCDC167, and LINC01184. The profile of these and other interesting genes is shown in Supplementary Fig. 4. Upon examination of the overrepresented transcription factor binding sites on the promoters of IL32 cluster genes, the V$IK_Q5_01 motif bound by Ikaros (IKZF1) was revealed to be among the enriched transcription factor binding sites shared in both the CD4+ and PBMC fractions (Supplementary Table 5). IKZF1 has been genetically associated with T1D (33). The T1D-associated risk allele rs10272724 (T) increases IKZF1 transcript level (34).

IFNG cluster of the CD8+ cells included TBX21 (codes for TBET), BHLHE40, and ZEB2, transcription factors expressed in CD8+ T cells (35), as well as NKG7, OASL, and KLRD1 (Supplementary Table 5). ZEB2 has been reported to drive terminal effector CD8+ cell differentiation together with T-bet (36). In the PBMC fraction, GCG and REG1A were coregulated with INS (Supplementary Table 4 and Supplementary Fig. 5).

Transcriptional Changes Preceding the Appearance of T1D-Related Autoantibodies Are Enriched in the CD8+ T-Cell Fraction

To identify changes that occur immediately before the first detection of T1D-related autoantibodies (i.e., seroconversion), we performed a separate differential expression analysis for the samples drawn at most 12 months before seroconversion. Altogether, 121 coding and noncoding genes were DE in case subjects compared with matched control subjects (Supplementary Table 4 and Supplementary Fig. 6). Notably, more than half of these (58%) were detected only in the CD8+ fraction. Besides IL32, only two other genes were common to all fractions, RPSAP58 and RPSAP15—both being the pseudogenes with unknown functions with very similar expression profiles (Supplementary Fig. 4MT).

Higher IL32 expression in case subjects was validated using qRT-PCR. Interestingly, genes encoding all three major isoforms (IL32 α, β, and γ) were upregulated in PBMC samples in all the case children at each of the time points, including 3 months (Fig. 3A and Supplementary Fig. 7). Among these isoforms, the gene encoding IL32γ was expressed at the highest level, followed by IL32β and IL32α.

Figure 3
Figure 3

scRNA-seq of PBMCs identifies T and NK cells as IL32 high populations. A: Expression of IL32γ isoform in longitudinal PBMC samples of case subjects and their control subjects (n = 7 + 7), assayed by qRT-PCR (for α and β isoforms, see Supplementary Fig. 7). B: t-SNE (t-distributed stochastic neighbor embedding) clusters from the pooled data from all scRNA-seq samples (4 case and 4 control subjects [in total 18,396 cells]). Clusters are named according to the expression of classical marker genes, such as CD8A (for details and marker gene list, see Supplementary Fig. 8; for contribution of each sample per cluster, refer to Supplementary Figs. 9 and 10). C: Expression of IL32 in the 12 cell clusters (natural logarithm transformation with addition of 1). For case-control comparison, please see Supplementary Fig. 11 DF: Trajectories emerging when using the data from CD4+ cells and the precursor cells. G-I: Trajectories emerging when using the data from CD8+ and the precursor cells. Here, “precursor cells” refer to cells from the naive and RGCC+ T-cell clusters. For the trajectory analysis of all the cells from all clusters as well as the breakdown of each individual cluster, see Supplementary Fig. 12. In D and G, cells are colored based on the contributions from different t-SNE clusters. In E and H, cells are colored by case (orange) or control (gray) status. In F and I, cells are colored by the intensity of IL32 expression (log10 transformation with addition of 0.1). Act., activated; prolif., proliferating.

scRNA-seq Identifies T and NK Cells as the IL32-High Population

To specify the cell populations responsible for the IL32 and INS signatures, we performed scRNA-seq on four selected case and their nearest matched control PBMC samples where the expression of IL32 or INS was high (or low) based on the bulk RNA-seq data (Supplementary Table 6). Unsupervised clustering of 18,396 single cells from all eight PBMC scRNA-seq runs identified 13 clusters (Fig. 3B and Supplementary Fig. 8). The 2 largest clusters expressing high CCR7 were merged as one cluster of naive T cells, reducing the number of clusters to 12. Clusters named as “RGCC+ T cells,” “CD62L+ T cells,” and “activated Th cells” expressed lower levels of CCR7. Activated CD8+ T cells cluster expressed high levels of CD8A and CD8B as well as NKG7, and two separate clusters of CD8+ T cells expressing either granulysin or granzyme A were observed (“activated GNLY+ CD8+ T cells” and “activated GZMA+ CD8+ T cells,” respectively). A subcluster of activated GZMA+ CD8+ cells had higher expression of cell-cycle genes (e.g., STMN1, TUBA1B) and was named “activated proliferating GZMA+ CD8+ T cells.” An NK cell cluster was positive for expression of CD56, NKG7, and GNLY and negative for CD8A and CD3E. A B-cell cluster was identified by the expression of MS4A1, CD79A, and CD79B, whereas the monocyte/dendritic cells (DC) cluster was composed of cells expressing CD14 or FCGR3A, LYZ, and TYROBP. Interestingly, the expression of many HLA class II molecules was as high in B cells as in monocytes, suggesting high antigen-presentation potential.

The contribution of different case or control samples to the cells in a given cellular population (cluster) varied from cluster to cluster (Supplementary Figs. 9 and 10A and B). The naive T-cell cluster was dominated by the cells from the control samples (P < 0.05) whereas the monocyte/DC cluster had more cells from case subjects (P < 0.005) (Supplementary Fig. 10B). Case 9, with the highest IL32 expression levels in the bulk RNA-seq data, dominated the CD62L+ T-cell cluster, activated NK cell cluster, and, most clearly, activated and proliferating GZMA+ CD8+ T cell clusters (Supplementary Fig. 10B). Conversely, control children 5 and 9 seemed to dominate the cluster of developing T cells expressing pre–T-cell receptor PTCRA, suggesting the presence of immature T cells in those samples.

Insulin, glucagon, and REG1A expression was not detected even in the INS-high samples of cases 5 and 9, leaving the origin of these transcripts in bulk RNA-seq as an open question. In contrast, IL32 expression was clear, and as expected, it was explicitly overexpressed in the case samples (Supplementary Fig. 11). IL32 was expressed at a very low level in monocyte/DC, B cell, and developing T cell clusters; however, it was expressed at higher levels by both the T cells and the NK cells (Fig. 3C).

To further define the relationship of IL32 expression and T-cell activation status, we performed separate trajectory analyses for the CD4+ and CD8+ T cells. The less activated precursor populations (naive and RGCC+ T cells), which detect CD4 and CD8 encoding transcripts in low abundances, were used as a starting point for the trajectory analyses. The results revealed three major cellular branches (I–III) in the data in both CD4+ and CD8+ T cells (Fig. 3D–I). Branch I consisted mainly of naive T cells, among which cells from the control samples were enriched (Fig. 3E and H and Supplementary Fig. 12). In contrast, the highest levels of IL32 were expressed by cells close to the end points of branches II and III, corresponding to more advanced stages of differentiation (Fig. 3F and I and Supplementary Fig. 12).

IL32 and IFNγ Are Coexpressed by Th1 Cells

To further study IL32 expression, we measured intracellular IL32 expression at the protein level in CD4+ T cells isolated from human umbilical cord blood. Cells were either activated through CD3/CD28 in the absence of cytokine (Th0) or were differentiated toward a Th1 cell lineage for 72 h. IL32 was induced upon activation and, unlike IFNγ, was expressed in both Th0 and Th1 cells (Fig. 4A). Interestingly, in Th1 cells, most IFNγ-producing cells were also positive for IL32 (Fig. 4A and Supplementary Fig. 13A) and the proportions of IL32-positive cells and the per cell IL32 levels were higher in IFNγ-producing Th1 cells than in Th0 cells (Fig. 4B and C). Furthermore, neutralization of IFNγ significantly reduced IL32 secretion by Th1 cells (Fig. 4D), confirming that IFNγ positively regulates IL32 expression. IL32 expression was also induced by IL32 itself in Th1 cells, both at the RNA level (Fig. 4E) and in the culture supernatant upon 48-h restimulation after 7 days of polarization in Th1 condition (Fig. 4F).

Figure 4
Figure 4

Virus- and cytokine-induced IL32 expression by pancreatic β-cells. A: Representative FACS dot plots showing IFNγ and IL32 double staining in Th0 and Th1 polarized CD4+ cells. Staining controls and two other replicates are shown in Supplementary Fig. 13A. Percent IL32-positive cells as well as median fluorescence intensity (MFI) data (mean ± SD) from all the three replicates are shown in B and C, respectively. Statistical significance was determined by paired two-tailed t test. D: IL32 secretion in culture supernatant as measured by ELISA. Cells were cultured in Th0/1 condition for 72 h in the presence (+) or absence (−) of anti-IFNγ. The expression plotted is relative to Th0 (−). Statistical significance was determined by paired two-tailed t test. E: IL32 expression in nonpolarized Th0 cells and cells differentiated to Th1 for 72 h in the presence (+) or absence (−) of IL32γ as measured by the Taqman assay. The expression is calculated relative to EEF1A. Statistical significance was determined by unpaired two-tailed t test. F: IL32 secretion in culture supernatant as measured by ELISA. Cells were cultured in Th0/1 condition for 7 days in the presence (+) or absence (−) of IL32γ, followed by washing and restimulation by phorbol myristic acid and ionomycin for 48 h. The expression plotted is relative to Th0 (−). Statistical significance was determined by paired two-tailed t test. G: Expression of the TNFA and IL6 or IL8 and IL32 genes when the EndoC-βH1 cells were stimulated with IL32γ alone or in combination with other inflammatory cytokines for 24 h. The fold change is calculated compared with nontreated cells. The results shown here are from four independent biological replicates (mean ± SD). Statistical significance for the effects on IL32 expression was determined by paired two-tailed t test. ns, not significant. H: IL32 expression as measured in an RNA-seq experiment where pancreatic islets were infected with CBV1-7. Statistical significance was determined by edgeR. I: IL32 expression in virus-infected pancreatic islets as measured by RT-qPCR Taqman assay. The expression is calculated as 2^−(dCt). The statistical significance is determined by paired two-tailed t test. RPKM, reads per kilo base per million mapped reads. *P < 0.05; **P < 0.01; ***false discovery rate <0.001.

Pancreatic β-Cells Can Express IL32 in Response to Cytokine Stimulation and Viral Infection

To study how the elevated IL32 levels may influence β-cell function, we treated human EndoC-βH1 β-cell line for 24 h with recombinant IL32γ either alone or in combination with the proinflammatory cytokines IL1β and IFNγ. In agreement with earlier published data on pancreatic ductal cancer cell lines (37), IL1β and IFNγ significantly induced IL32 expression in human EndoC-βH1 cells (Fig. 4G). However, addition of IL32γ did not further enhance 1) the IL1β- and IFNγinduced IL32 expression; 2) the expression of inflammatory cytokines TNFA, IL6, and IL8 (Fig. 4G); or 3) the expression of ER stress marker genes (ATF3, ATF4, ATF6, HSPA5, CHOP, and sXBP1) (Supplementary Fig. 13B) in EndoC-βH1 cells. Furthermore, the IL32γ treatment did not affect the expression of β-cell–specific genes, such as INS, MAFA, and PDX1 (Supplementary Fig. 13C). These results suggest that, while IL32 does not appear to directly affect the survival or the differentiation status of the β-cells, β-cells actively contribute to inflammation in the islets by secreting IL32 upon stimulation by cytokines.

Coxsackie B viruses are β-cell trophic viruses that have been linked to the development of T1D (38–43). To study the possible trigger of IL32 expression in β-cells, we infected purified human pancreatic islets of three cadaveric donors with CBV1-7 strain. Infection by the virus led to the induction of IL32 expression in the islets (Fig. 4H). We further validated this finding in the three islet samples used for RNA-seq as well as one additional islet sample using qRT-PCR assays and found a consistent increase in the IL32 expression upon CBV1-7 infection (Fig. 4I). Taken together, these results suggest that upon a viral infection (Fig. 4H and I) or a cytokine rush (Fig. 4H), β-cells may upregulate IL32 secretion, contributing to inflammation.

Discussion

We identified a panel of novel molecular players detected early in children who developed T1D-associated autoantibodies or even the clinical disease at a young age. Since the immunological changes related to T1D are known to be strongest among the T1D cases diagnosed at an early age (44), focusing on this age-group should enhance the possibility of detecting aberrations in the immune system predisposing to the disease. In this study, unbiased RNA-seq of CD4+ and CD8+ cells revealed many T1D-associated DE transcripts not previously reported. Analysis of the PBMC population offers an excellent overview of stable gene expression patterns but, at the same time, appears to mask some of the subtle fraction-specific changes. Such changes included upregulation of CD52 detected only in the CD4+ cell fraction and downregulation of the IFNG and associated transcription factors ZEB2, TBX21, and ZNF683 detected specifically in the CD8+ cells. Further studies are needed to understand whether at-risk children have defects in formulating effector CD8+ response or their effector CD8+ cells have homed to the sites of inflammation in the pancreas.

We selected IL32 as our candidate for functional studies because it has not been linked to seroconversion before, it is easy to measure with available assays from clinical samples, and as a secreted molecule it can potentially affect the function of several cell types in paracrine and systemic fashion. Increased expression of IL32 in case subjects across many cell types before seroconversion suggest that IL32 is a critical member of the immunological signature characteristic for children developing β-cell autoimmunity.

IL32 is expressed by many immune and epithelial cells and has been described to be proinflammatory (45). However, to our knowledge, it has not been associated with human β-cell autoimmunity. In contrast, IL32 is downregulated in CD4+ T cells from recently diagnosed adult T1D patients (46), which, along with our findings, suggests dynamic changes in immune cell signaling during the pathogenesis of the disease. On the other hand, IL32 overexpression was observed in synovial biopsies of patients with rheumatoid arthritis (47), in inflamed mucosa of inflammatory bowel disease patients (48), and in the serum of myasthenia gravis patients (49), indicating a connection between IL32 and autoimmunity in general. In T cells, IL32 is induced by T-cell activation, and it modulates human CD4+ T-cell effector function by promoting Th1 and Th17 responses (50). Both Th1 and Th17 cells have been linked to T1D pathogenesis in both human and mouse (51). The IL32 gene has been identified only in higher mammals, excluding rodents. Nonetheless, human IL32γ transgenic mice exhibit impaired glucose tolerance and increased levels of IFNγ and other proinflammatory cytokines in the pancreas, as well as accelerated streptozotocin-induced experimental T1D (52). No specific cell-surface receptor for IL32 has been identified, but it may act through cell-surface integrins or proteinase-3 (53).

Our results showed that IL32 was often coregulated with genes previously linked to autoimmunity. For example, the BTN3 gene cluster resides in the extended MHC class I locus. Further, BTN3 genes have been associated with T1D in a genetic screen, especially in the case of BTN3A2 (54). AMICA1 is a plasma membrane protein involved in lymphocyte migration through its interaction with Coxsackie-adenovirus receptor (CAR) expressed in epithelial cells and has been associated with multiple sclerosis (55). An analogous scenario could be envisaged for T1D: CAR is expressed by the pancreatic islet cells, including β-cells (42), and its expression is elevated in autoantibody-positive individuals and patients with T1D (56), suggesting that it might help recruit T cells to the islets. Interestingly, the findings point to human-specific phenomena not detectable in mouse models, as IL32 and the BTN3 protein family are not encoded by the mouse genome.

The strength of our study is that the children studied here comprise a homogeneous population with the early appearance of T1D-associated autoantibodies. Increasing evidence suggests that T1D can be subdivided into different phenotypes, e.g., characterized by age-dependent B-cell infiltration in the pancreas (57), defect in Coxsackievirus-induced antibody response in children with early insulin autoimmunity (58), or rapid versus slow progression to clinical disease (59). Thus, our results may not apply to “late progressors,” adolescents, and adults. Although the analysis of the global transcriptome of T-cell subsets of children with prediabetes over the period of seroconversion is unique, a limitation of the current study is the analysis of only seven Aab+ children. The results of this study need to be validated and expanded in a larger cohort of children with prediabetes but serve as a starting point for better understanding of immunological changes preceding the clinical onset of the disease. In the future, we are interested in addressing whether our findings on a cellular level are reflected also in IL32 levels in plasma as well as studying whether IL32 alone or in combination with other identified molecules would have sufficient sensitivity and specificity as an early indicator for T1D.

Article Information

Acknowledgments. The authors are grateful to the families for their participation in the DIABIMMUNE study. The DIABIMMUNE study group is acknowledged for excellent collaborations, work with the families, and collection of the samples for the study. Marjo Hakkarainen, Sarita Heinonen, Päivi Junni, and Elina Louramo (Turku Bioscience Centre, University of Turku and Åbo Akademi University, Turku, Finland) are acknowledged for skillful assistance in the laboratory. Next-generation sequencing was performed at the Finnish Functional Genomics Centre (FFGC), Turku, Finland, part of the Biocenter Finland network. The authors thank Satu Mustjoki and her team at University of Helsinki for advice in designing scRNA-seq experiments and Riina Kaukonen at the FFGC for the sample preparation.

Funding. This work was financially supported by the JDRF; the Academy of Finland (AoF) Centre of Excellence in Molecular Systems Immunology and Physiology Research (SyMMyS) 2012–2017 (grant no. 250114); the AoF Personalized Medicine Program (grant no. 292482); AoF grants 294337, 292335, 319280, and 314444; the Sigrid Jusélius Foundation; the Diabetes Research Foundation (Diabetestutkimussäätiö); the Novo Nordisk Foundation Innovative Medicines Initiative 2 Joint Undertaking under grant agreement no. 115797 (INNODIA). This Joint Undertaking receives support from the Union’s Horizon 2020 research and innovation programme and EFPIA, JDRF, and The Leona M. and Harry B. Helmsley Charitable Trust. The DIABIMMUNE study was supported by the European Union Seventh Framework Programme (grant no. 202063). T.L. was supported by the AoF (311081).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. H.K., S.T., and U.U. were responsible for the interpretation of the results. J.S. conducted bioinformatic analyses. H.K., J.S., S.T., and U.U. drafted the manuscript. H.K., J.S., and U.U. prepared the figures. H.K. was responsible for supervising E.K. R.d.A., E.K., and O.R. were responsible for the isoform-specific IL32 RT-PCR assay and the intracellular IL32 staining in T cells and interpretation of the results. T.L. provided expertise in scRNA-seq study design, sample and data analysis, and interpretation of the results. H.S., J.H., T.H., A.P., and V.T. were responsible for sample collection, sample storage, and further clinical information of the children. V.C. and T.O. carried out the experiments and interpreted the results of the studies in pancreatic β-cells. M.K.A. and G.F. were responsible for experiments on virus-infected pancreatic islets. R.Lu. was responsible for study design, cell fractionation, sample analysis, and data production. H.L. was responsible for computational data analysis, interpretation of the results, editing the manuscript, and supervising J.S. M.K. was responsible for the DIABIMMUNE study design, sample collection, sample storage, clinical information for the children, directing of the clinical study, interpreting the results, and editing the manuscript. R.La. was responsible for study design, sample and data analysis, interpretation of the results, writing the manuscript, and supervision of the study. All authors contributed to the final version of the manuscript. H.L. and R.La. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.



Source link

Leave a Reply

Your email address will not be published.


*


Shares