Evolving epigenomics of immune cells in type 1 diabetes at single nuclei resolution

TL;DR

This study tracks immune cell changes in children who later develop type 1 diabetes, revealing early epigenomic differences that suggest immune system divergence before diagnosis.

Contribution

The study provides new insights into the evolving epigenomic profiles of immune cells during the early stages of type 1 diabetes development.

Findings

Early monocyte lineage differences suggest heightened pro-inflammatory responses in children who develop T1D.

Autoantibody emergence varies across lymphoid and myeloid cells, indicating divergent immune responses.

Epigenomic changes in immune cells are most pronounced in early life, before antibody detection.

Abstract

The appearance of diabetes-associated autoantibodies is the first detectable sign of the disease process leading to type 1 diabetes (T1D). Evidence suggests that T1D is a heterogenous disease, where the type of antibodies first formed imply subtypes. Here, we followed 49 children, who subsequently presented with T1D and 49 matched controls, profiling single-cell epigenomics at different time points of disease development. Quantitation of cell and nuclei populations as well as transcriptome and open-chromatin states indicated robust, early, replicable monocyte lineage differences between cases and controls, suggesting heightened pro-inflammatory cytokine secretion early among cases. The order of autoantibody emergence in cases showed variation across lymphoid and myeloid cells, potentially indicating cellular immune response divergence. The strong monocytic lineage representation in peripheral blood immune cells before seroconversion and the weaker differential coordination of these gene networks close to clinical diagnosis emphasizes the importance of early life as a critical phase in T1D development.

Figures

a](#F1). The first time point in the cases was before seroconversion at the mean age of 1.6 years (range 0.3 to 6.1 years), the second time point soon after the seroconversion at the mean age of 3.0 years (range 0.8 to 8.2 years) and the third time point close to the time of T1D diagnosis at the mean age 6.0 years (range 1.5 to 13.0 years). The time points in the controls corresponded to similar ages in the cases. Genetically unique samples (N=98) were pooled prior to capture allowing us to minimize potential batch variation in downstream single nuclei or single cell assays. Three different ca

a](#F2)). We note that as expected correlations show 9:1 overrepresentation of high access chromatin and high gene expression as opposed to negative correlations between chromatin accessibility and gene expression. Similarly, most genes have multiple open chromatin regions correlating to its expression ([Figure 2b](#F2)).

a](#F3),[b](#F3)) predicts replication strongly in snRNA and scRNA and throughout data types, which are also correlated with the median read depths of the features and the number of replicated data layers ([Figure 3c](#F3)). The highest discovery rate and largest number of replicated genes is observed in snRNA data, which also shows the highest read depth per feature.

b](#F4)) the strongest evidence for coordinated differences were observed for monocytes and *IRF1*-motif (IRF1.MA0050.2) showing a fourfold increase in occurrence among case enriched open chromatin (*P*=6.22E-19) ([Figure 4c](#F4)), whereas among T cells numerous TFs show relative overrepresentation among cases and B cells or NK cells yielded no significant differences of directional TF binding. The pathway analyses based on differentially regulated (case vs. control) peak – gene links ([Figure 2](#F2)) for chromatin peaks per cell lineage show monocyte predominance in pathway enrichments simi

a](#F5)). Consequently, the most prominent differences in all cases vs. controls in pathway analyses ([Figure 4](#F4)) seen for early monocyte overexperessed genes were recapitulated in “GAD first” subgroup data ([Figure 5b](#F5)), whereas the “IAA first” subgroup shows enrichment for adaptive immune system networks in CD4 T-cells ([Figure 5c](#F5)) before IAA-appearance (Timepoint 1).

Full text

INTRODUCTION

Type 1 diabetes (T1D) is perceived as a chronic immune-mediated disease characterized by selective loss of insulin-producing β cells in the pancreatic islets of genetically susceptible individuals. The lifetime risk for people in the North American general population for developing T1D is around 0.3%, but it increases to 5% by the age of 15 years for children carrying risk associated HLA genotypes (DR3/DR4-DQB1*03:02)^1–3^. In addition to HLA there are over 100 other genetic loci described in the literature^4^, but collectively these genetic biomarkers do not reach high positive predictive value^5^. Symptomatic disease is preceded by an asymptomatic period of highly variable duration during which T1D-associated islet autoantibodies (IAA) appear in the peripheral circulation as markers of emerging β-cell autoimmunity^6,7^. In natural history studies, positivity for two or more autoantibodies signals a risk of approximately 70% for the development of clinical T1D over the subsequent 10 years^8^. Evolving data also suggest that T1D is not a homogenous disease but rather that there are different endotypes^9^. Such endotypes can be defined based on disease characteristics such as which autoantibody specificity appears first during the disease process^10,11^, or the age at the manifestation of clinical T1D disease^12^. These insights relied on a unique study design within The Trial to Reduce IDDM in the Genetically at Risk (TRIGR) study, including intense longitudinal follow-up starting at an early age (3 months onwards) coupled with rigorous biobanking of peripheral blood mononuclear cells (PBMCs) along with other biospecimens^13^. Therefore, TRIGR provides a unique opportunity to observe immune-cell phenotypes in an at-risk population.

We extend these efforts here, leveraging TRIGR and a time-resolved sample series of 98 participants (49 cases who progressed to clinical T1D and 49 age-matched controls), and profiled single-cell epigenomics early, before the first signs of autoimmunity (autoantibody seroconversion), soon thereafter, and close to the manifestation of clinical T1D. We performed population cell and nuclei analysis including transcript and open-chromatin quantitation across immune cell lineages in both T1D cases and controls. We integrated the order of autoantibody emergence in the case group and showed variation in immune cells, potentially indicating cellular immune response divergence among T1D endotypes.

ONLINE METHODS

The Trial to Reduce IDDM in the Genetically at Risk (TRIGR) was a randomized clinical trial designed to assess whether it is possible to prevent β-cell autoimmunity and clinical T1D by weaning high-risk infants to an extensively hydrolyzed formula. The primary outcomes of the trial turned out to be negative^14,15^. Out of the 2159 participants 173 (8.0%) progressed to clinical T1D during follow-up to the median age of 11.5 years.

Study subjects

The current study included 98 European TRIGR participants, out of whom 49 progressed to clinical T1D during the follow-up (cases). Forty-nine autoantibody-negative control subjects were matched with the cases for date of birth (± 1 year) and geographical region. All TRIGR participants had a family member affected by T1D and carried an HLA genotype conferring increased risk for T1D [the high risk genotype DQB102/DQB103:02, moderate risk genotypes DQB103:02/x* (x not DQB102, DQB103:01, or DQB106:02),* mild risk genotypes DQA105-DQB102/y (y not DQA102:01-DQB102, DQB103:01, DQB106:02, or DQB106:03)* or the rare mild risk genotype DQA103-DQB102/y (y not DQA102:01-DQB102, DQB103:01, DQB106:02, or DQB106:03)]*. The TRIGR participants were recruited in 2002–2007 before or immediately after birth followed by randomization to two weaning groups. The participants attended the study centers at the age of 3, 6, 9,12, 18 and 24 months and subsequently annually until the youngest child turned 10. A sample for the isolation of peripheral blood mononuclear cells (PBMC) was collected at each visit and the sample from European participants was shipped overnight for isolation in the TRIGR Core Laboratory in Helsinki. Seroconversion to autoantibody positivity was observed in the 49 cases at the mean age of 2.4 years (range 0.7 to 7.1 years) and T1D was diagnosed at the mean age of 6.6 years (range 1.0 to 13.7 years). The controls remained autoantibody negative throughout the follow-up. Each participant in the current study contributed three samples; the first collected before the appearance of the first autoantibody (mean age 1.6 years, range 0.3 to 6.1 years), the second after the seroconversion to autoantibody positivity (mean age 3.0 years, range 0.6 to 8.2 years) and the third close to the diagnosis of T1D (mean age 6.8 years, range1.0 to 13.6 years). The samples from the control subjects were obtained at the corresponding ages. The detailed characteristics of the participants are presented in Table 1. The TRIGR study was approved by the Ethical Committees of all participating centers.

Analysis of diabetes-associated autoantibodies

Diabetes-associated autoantibodies, i.e. insulin autoantibodies (IAA), GAD autoantibodies (GADA) and autoantibodies to islet antigen 2 (IA-2A) and zinc transporter 8 (ZnT8A), were quantified with the use of specific radiobinding assays in the Scientific Laboratory, Children’s Hospital, University of Helsinki, Helsinki, Finland in all samples obtained from the TRIGR participants during their follow-up^16^.

Endpoints

The primary endpoint was progression to overt T1D and included comparison between the cases and controls at the three time points analyzed. Secondary analyses comprised comparison of the “IAA first” and “GADA first” endotypes among the cases. The strict criterion for “IAA first” marked that IAA was the only detectable autoantibody in the first autoantibody-positive sample, while the looser criterion indicated that the first autoantibody-positive sample tested positive for multiple autoantibodies including IAA but not GADA. Similarly the strict criterion for “GADA first” marked that GADA was the only detectable autoantibody in the first autoantibody-positive sample, while the looser criterion indicated that the first autoantibody-positive sample tested positive for multiple autoantibodies including GADA but not IAA. For endophenotype analyses we utilized strict criterion “IAA first” (n = 21) and looser “GADA first” (n = 11) definitions and generated differential expression / peak analyses using the paired controls against each major cell lineage and timepoint. SnRNA signals replicated at qv = 0.1 in min of two data layers were kept for comparison to overall case – control analyses to explain the endophenotypes’ gene network contributions of to overall case – control differences.

Isolation of PBMCs from fresh blood samples

Peripheral blood mononuclear cells (PBMCs) were separated from fresh heparinized venous blood samples by Ficoll (GE Healthcare, Uppsala, Sweden) isogradient centrifugation. Prior to the isolation of PBMCs, blood plasma was separated with centrifugation and the separated plasma was substituted by an equal volume of phosphate buffered saline (PBS). The diluted blood was then placed on top of 5mL of Ficoll in 15mL conical Falcon tubes. The tubes were centrifuged for 25 minutes at 800 × g at room temperature (RT) with no brake during the deceleration step. The upper plasma layer was removed by aspirating with a sterile individually packed plastic Pasteur pipette. Then, the white layer containing the mononuclear cells was collected and directly transferred into a new 15mL Falcon tube containing 10mL PBS. The collected cells were washed three times with PBS. After the washing steps, the supernatant was carefully removed and the cell pellet was suspended in RPMI1640 cell culture medium supplemented with 5% of heat inactivated human AB serum, 2mM L-glutamine and gentamicin (Sigma-Aldrich, 25ug/mL). The cell density and viability was then analyzed using a hematocytometric chamber (Burker) and Trypan blue staining.

Freezing the isolated PBMCs

Isolated PBMCs were transferred into a new conical Falcon tube and the cells were centrifuged at room temperature (RT) for 5 minutes at 400g with centrifuge brake on. The supernatant was pipetted off and the cells were resuspended in cold (+ 4°C) RPMI1640 with 5% of heat inactivated human AB serum, 2mM L-glutamine and gentamicin. At this point the cell density was between 4–16×10^6^/mL. Then an equal volume of cold (+ 4°C) RPMI1640 cell culture medium with 5% of heat inactivated human AB serum, 2mM L-glutamine, gentamicin (25ug/mL) and 20% of dimethyl sulphoxide (DMSO, Sigma-Aldrich) was added dropwise slowly to the cell suspension while swaying the tube gently. Cell suspension was then aliquoted into the cryovials. One ml (1ml) of the cell suspension was transferred to cryovial tubes. The tubes were closed and put immediately in a MrFrosty container with isopropyl alcohol for freezing the cells in a temperature-controlled way. The MrFrosty container with the cells was then immediately transferred into a −80°C freezer overnight. The cryovials were transferred the next day into the automated gas phase liquid nitrogen freezer for storing at −180°C until thawed and analyzed.

Batching the thawed PBMCs

Cell samples were thawed in a group of 5 samples at a time in a 37°C water bath until a small ice crystal remained in the cryovials. 1 mL of pre-warmed thawing medium (RPMI 1640 culture medium (Life Technologies, Paisley, United Kingdom) containing 2 mmol L-glutamine, 25 mmol/l HEPES, 25 μg/ml gentamycin (Sigma-Aldrich) and 10% of heat-inactivated AB serum (Innovative Research, Novi, MI, USA), were added to the cryovials slowly. The cell suspensions were transferred to the 15 mL conical tubes containing 10 mL of thawing medium. The cells were centrifuged at 300 × g for 8 minutes and the cell supernatant was discarded. The cell pellet was gently resuspended in a volume of 0.5 mL of room-temperature Thawing Medium. An evaluated volume of cell suspension containing approximately 1 × 10^6^ cells per sample of 42–46 cell samples were pooled together. After filtering the pooled cell sample with a 40 μm cell strainer, the cells were counted, and the viability was assessed. On average, the viability of the cells was 89% with an SD of 9.1%. After viability analysis, the cell suspension was divided between two 15 ml conical tubes and centrifuged at 300 × g for 8 minutes at 4°C. The supernatant was discarded, and the cell pellets were resuspended in 1 mL of cold DMSO-containing cryopreservation medium and transferred to cryovials. The cryovials were placed into a slow-freeze cryopreservation chamber in a −80°C freezer. On the next day the frozen cells samples were relocated into a liquid nitrogen storage. The frozen cryovials were then transported from Helsinki to Kansas City in a cryoshipper.

DNA genotype analyses and genotype demultiplexing

The high-throughput microarray genotyping data (from Axiom when available) was merged with variant calls from WGS samples (TRIGR_1045832, TRIGR_1061181, TRIGR_1062514, TRIGR_1062558, TRIGR_1071138, TRIGR_1085671). The donor-multiplexed multiome data was deconvoluted using demuxlet v2^17,18^ that utilizes the merged imputed and WGS genotype matrix (436683 positions) of the pool to accurately assign each cell barcode (CID) to an individual donor and importantly identify cell doublets to be removed from the data.

Single-cell RNA-seq (scRNA), single nuclei RNA and ATAC (snRNA/snATAC) methods

Thawed cells were centrifuged at 300 × g for 8 minutes and resuspended in 10 mL of Thawing Medium consisting of IMDM (ATCC Cat. No. 30–2005) supplemented with 10% fetal bovine serum, 100 Units/mL Penicillin, and 100 μg/mL of Streptomycin. Cells were washed once more in Thawing Medium then twice in PBS + 0.04% BSA and filtered through a 40-μm nylon mesh cell strainer. A Countess II automated cell counter was used to assess the cell count and viability of the pool. Approximately 25,000 cells were loaded into each of four wells of a 10x Genomics Chromium Chip B, and scRNA-seq was performed using 10x Genomics Chromium Single Cell 3’ Reagent Kits v3 according to the manufacturer’s protocol. SnATAC was then performed in parallel with scRNA. Nuclei were isolated from an aliquot of the pooled cells following the 10x Genomics Demonstrated Protocol: Nuclei Isolation for Single Cell ATAC Sequencing with an optimized lysis time of 3 minutes. Approximately 15,300 transposed nuclei were loaded into each of four wells of a 10x Genomics Chromium Chip E, and snATAC-seq was performed using 10x Genomics Chromium Single Cell ATAC Reagent Kits v1 according to the manufacturer’s protocol. All scRNA and snATAC libraries were loaded on an Illumina NovaSeq 6000 sequencer using S4 Reagent Kit v1.5 (200 cycles) and sequenced to the recommended depths.

10X Genomics Chromium Single Cell Multiome ATAC + Gene Expression protocol was deployed in independent replicate cell pools for simultaneously assessing nuclear RNA expression and open chromatin by ATAC profiles on the Chromium Instrument with Next GEM technology. Pooled cells were thawed and washed as described above, and nuclei were isolated following the 10x Genomics Demonstrated Protocol: Nuclei Isolation for Single Cell Multiome ATAC + Gene Expression Sequencing with an optimized lysis time of 3 minutes. For each pool, approximately 16,100 transposed nuclei were loaded into each of eight wells of a 10x Genomics Chromium Next GEM Chip J to target at least 1,000 nuclei per donor retained after data filtering for cell quality and doublets using conservative estimates (see rationale for this strategy in more details below). Simultaneous snRNA/ATAC-seq was performed using 10x Genomics Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kits according to the manufacturer’s protocol. All Multiome libraries were loaded on an Illumina NovaSeq 6000 sequencer using S4 Reagent Kit v1.5 (200 cycles) and sequenced to the minimum depth of 50000 GEX and 25000 ATAC reads per nuclei, respectively.

Processing and differential gene / peak expression in scRNA/snATAC and snRNA-seq analyses

Primary analysis of NovaSeq6000 genome sequences, to produce sequence reads, was performed with Illumina Real Time Analysis (RTA) software and bcl2Fastq2-20. Single cell RNA, snRNA and ATAC 10X Genomics data was processed by cellranger-4.0.0 and cellranger-atac-1.2.0. The Cell Ranger Single Cell Software was used to perform sample demultiplexing, barcode processing, and single cell 3′ gene counting for RNA data, along with fragment-mapping and counting for ATAC data. The R-based Seurat 4.0 package^19^ was used for QC, analysis, and exploration of the data.

Processing and observed/expected read depths in single-nuclei Multiome augmented with snATAC data

Seurat 4.0 package was used for QC, analysis, and exploration of the scRNA-seq and Multiome data and to identify and interpret sources of heterogeneity from the measurements. Quality filtering was based on the number of unique molecular identifiers (UMIs) and/or mitochondrial content in each cell. For snATAC, we applied Signac (https://CRAN.R-project.org/package=Signac) workflow^20^, which is based on term frequency-inverse document frequency (TF-IDF) normalization and dimensionality reduction using singular value decomposition (SVD). Quality filtering based on number/fraction of reads in peaks, transcriptional start site (TSS) enrichment score, ratio of reads in blacklist (non-unique genomic segments contaminating chromatin datasets). ATAC signal was evaluated in peaks taken from lymphoid and myeloid/erythroid DNase Hypersensitivity Sites^21^ showing 34,555/99,779 snATAC peaks overlapping lymphoid, -and 22,707/99,779 peaks overlapping myeloid DHS (total overlaps 51,208/99,779).

We used unsupervised clustering and discovery of cell types/states and dimensionality reduction by the Uniform Manifold Approximation and Projection (UMAP) method as implemented in Seurat^19^. We manually curated all cell annotations provided by Seurat based on supervised literature mining. For Multiome datasets, the datasets were normalized using SCTransform and then integrated, using Seurat v4. Information from both scRNA and scATAC modalities were integrated independently for each cell using a weighted-nearest neighbor (WNN) approach. This incorporates information from both RNA and ATAC data robustly, without regard to the relative quality of each dataset.

Differential analysis was performed by merging the cell types into five main groups guided by Seurat cell label transfers where “CD4”, “CD8”, “NK” and “B-cell” included all cells in any cluster with lineage label (e.g. “NK” included nuclei from NKT and NK CD56bright/dim labelled clusters), “Monocyte” included all Monocyte (e.g. CD14, CD16) and DC labelled nuclei. Following collapse of sublineage clusters we used data of nuclei from five parent lineage clusters and applied Fisher’s exact test to compare the total number of reads in a particular peak/gene (“observed”) and the total number of reads in all other peaks/genes (“expected”) for one group and comparing to the observed/expected total reads for the second group.

In both the snRNA dataset and the integrated snATAC, we computed correlation of RNA and ATAC signals near 141 known T1D-related single nucleotide polymorphism (SNP) loci. We extracted the imputed genotype value for each sample at the SNP, represented as an alt allele count of 0, 1 or 2 for genotypes 0/0, 0/1 and 1/1, respectively. SCT normalization was used for RNA merging, macs2peaks were used for the ATAC, using the Seurat and Signac packages. The read counts were normalized per 1M total gene expression or ATAC reads for each individual and time point. For each of the five parent lineage clusters, we looked at genes (for the multiome dataset) or peaks (for the ATAC dataset) within 1Mbp of a SNP locus and used Spearman rank correlation to compute the rank correlation coefficient rho and P-value of the normalized signal in a particular peak/gene for each individual – total number of reads in a particular peak/gene divided by the total number of reads in all peaks/genes – against the matched genotype value of the SNP for each individual. As with other analyses we utilized five aggregated cell types CD4T, CD8T, B-cell, NK-cell and monocytes. Multiple testing was accounted for by Bonferroni correction using the number of tests applied per cell type and data modality.

RESULTS

Longitudinal Single-Cell and Single-Nuclei Profiling of Peripheral Blood Cell Types

To determine the immune cellular and regulatory landscape underlying T1D, we analyzed 98 individuals (Table 1, 49 progressing to T1D and 49 non-progressors) enrolled in the TRIGR study each having peripheral blood mononuclear cells collected at three time points. Overall study design is depicted in Figure 1a. The first time point in the cases was before seroconversion at the mean age of 1.6 years (range 0.3 to 6.1 years), the second time point soon after the seroconversion at the mean age of 3.0 years (range 0.8 to 8.2 years) and the third time point close to the time of T1D diagnosis at the mean age 6.0 years (range 1.5 to 13.0 years). The time points in the controls corresponded to similar ages in the cases. Genetically unique samples (N=98) were pooled prior to capture allowing us to minimize potential batch variation in downstream single nuclei or single cell assays. Three different captures were performed for each pool including 1) single cell RNA (scRNA), 2) single nuclei open chromatin (snATAC) and 3) Multiomics snATAC and snRNA. A total of 187,995, 186,449 and 263,308 singleton (cells or nuclei) passed QC in scRNA, snRNA and snATAC assays, respectively, yielding on average 2,278 independent single cells or nuclei measurements per sample and time point (Figure 1b–d; Supplementary Table 1). A total of 307,746 case and 330,006 control cells or nuclei passed the QC. The clustering (Figure 1e–g; Supplementary Figure 1a-b) of scRNA, snRNA and snATAC followed by inference of major cell types revealed no statistically significant differences (after multiple testing adjustment) in the frequency of major cell types between cases and controls at any given time point (Supplementary Figure 1c-f). However, strong calendar age-dependent trends were observed for proportional changes in B-cell and NK-cell subtypes (Supplementary Figure 2).

Linking open chromatin with gene expression at single nuclei resolution across immune cell types

The full snATAC dataset including 263,308 nuclei was used to call 99,779 peaks or regions with open chromatin across immune cell types. These open chromatin regions show 51% overlap with regions mapped by DNAseI hypersensitivity (DHS) in lymphoid and myeloid lineages (see methods) published by ENCODE^22^. We evaluated correlations of gene expression and open chromatin among the 186,449 nuclei up to 1Mb distance from each transcript start site (TSS) (Figure 2a). We note that as expected correlations show 9:1 overrepresentation of high access chromatin and high gene expression as opposed to negative correlations between chromatin accessibility and gene expression. Similarly, most genes have multiple open chromatin regions correlating to its expression (Figure 2b).

In total, we identified 31,162 nominally significant (P-value <0.05, |r|>0.05) correlations between open chromatin and gene expression. As expected, the highest density of open chromatin regions was observed proximal to TSS (20kb flanking) of the associated gene with a density corresponding to 0.9 links per kb. Distal open chromatin and gene expression correlations (20–500kb) were seen at a rate of 0.02 links per kb but even very distal associations (500kb-1Mb) were detected (0.01 links per kb). However, given the width of regulatory space explored the most abundant open chromatin and gene expression associations were seen at distal sites (Table 2). We leveraged the independent gene expression and open chromatin data from same individuals (corresponding to the unpaired scRNA and snATAC datasets) to validate the gene expression and open chromatin associations. Altogether >50% of links were independently verified showing same direction of effect and strong correlation (pearson |r|>0.65) (Table 2). These data provide core immunoregulome in the developing, circulating immune system extending to even distal regulatory sequences of common immune system genes. All open chromatin and gene expression correlations are listed in Supplementary Table 2.

Differential multiomics immune signatures among T1D cases vs controls

Initial differential expression analyses indicated high proportion of associations among highly expressed genes in small cell populations observed in only few samples. To mitigate this outlier bias in small cell clusters, we utilized a mapping strategy collapsing related clusters to Monocyte, B-cell, CD4 T-cell, CD8 T-cell and NK cell “parent lineages” for differential analyses between cases and controls. At a false discovery rate (FDR) of 10% across time points and cell lineages, we observed a total of 23,903 snRNA, 15,625 scRNA and 5,487 snATAC associations, respectively. We investigated the replication of cell lineage specific signals across the orthogonal data layers (time point, snRNA, scRNA, snATAC) for each individual molecular feature (Figure 3). Overall, original discovery significance (Figure 3a,b) predicts replication strongly in snRNA and scRNA and throughout data types, which are also correlated with the median read depths of the features and the number of replicated data layers (Figure 3c). The highest discovery rate and largest number of replicated genes is observed in snRNA data, which also shows the highest read depth per feature.

Early transcriptome differences in T1D cases vs controls linked to proinflammatory monocyte signatures

The primary snRNA case-control associations (at 10% FDR) were marginally (c^2^=10.2, df=1, P=0.0014) more abundant before (T1 = 8684 genes, 6% of tests) than soon after seroconversion (T2 = 8363 genes, 5.7% of tests). However, markedly fewer genes (c^2^=161, P<0.000001) were differentially expressed close to T1D diagnosis (T3 = 6856 genes, 4.7% of tests) as compared to T1. These associations were unequally distributed across cell lineages with monocytes showing greatest number of associations (N=6681 genes across 3 time points), followed by B cells (N=5312 genes), CD4 T cells (N=4557 genes), CD8 T cells (N=4162 genes) and NK cells (N=3192 genes). While read depth per lineage was correlated with discovery rate (r = 0.72) the strongest contributor to this was low read depth in NK-cells (r = 0.44 without NK-cells) and the high discovery rate in monocytes was observed at lower overall read depth (91%) compared to CD4 T-cells.

These observations including pronounced case-control associations in early time points as well as monocyte-lineage specificity were concordant in the orthogonal scRNA data sets. In fact, the primary scRNA case-control associations (at 10% FDR) show more pronounced differentiation (T1 vs. T2, c^2^=154, P<0.000001; T1 vs. T3 c^2^=428, P<0.000001) between cases and controls at the first time point (T1 = 6302 genes, 5.5% of tests) before seroconversion as compared to soon after seroconversion or close to clinical T1D (T2 = 5070 genes, 4.4% of tests; T3 = 4253 genes, 3.7% of tests). For the lineage-specific associations in the scRNA dataset 5098 genes across 3 time points were identified for monocytes, 4720 genes for CD4 T cells, 2374 genes for B cells, 1823 for CD8 T cells and 1704 genes for NK cells. We then focused the analysis to pursue GeneSet, pathway and gene ontology enrichments analyses among differential expressed genes replicated in at least 2 data layers (Figure 3). Greatest number of enriched gene networks were observed in monocytes, and a subset of these enrichments were replicated in B -and NK-cells. CD4 T cells show less coordination of genes throughout the progression of autoimmunity preceding T1D (Figure 4a). We also performed upstream regulator analyses of genes using TRRUST (Transcriptional Regulatory Relationships Unraveled by Sentence Based Text Mining)^23^ and found monocytes with NFKB1, STAT1, RELA, SP1 and IRF1 signatures particularly upregulated (early) whereas B-cells demonstrate MHC regulators (RFX, RFXANK, CIITA, RFX5) enriched for in controls (Figure 4b).

Open chromatin analysis confirms early monocyte-specific regulation in T1D cases vs controls

Similar to primary transcriptome analyses in case-vs-control the open chromatin signal robustness as measured by replication (either via linked gene expression or independent time point) were strongly dependent on read depth (Figure 3c, right panel) but overall discovery significance was narrower reflecting lower power for primary chromatin associations (Figure 3a, right panel). In total, 22.7% of the 5487 unique peak / cell type/ time point associations demonstrated evidence for primary case-control difference in chromatin openness with replication in at least one data layer.

We looked for evidence of regulators of chromatin openness among associated peaks within each cell type. Prevalence of transcription factor (TF) binding sites were extracted separately for peaks showing increased access in cases vs. controls and vice-a-versa using TFmotifView^24^. These analyses can reveal directional impact of certain TF networks. Consistent with TRUSST analyses of gene regulators from scRNA/snRNA gene expression (Figure 4b) the strongest evidence for coordinated differences were observed for monocytes and IRF1-motif (IRF1.MA0050.2) showing a fourfold increase in occurrence among case enriched open chromatin (P=6.22E-19) (Figure 4c), whereas among T cells numerous TFs show relative overrepresentation among cases and B cells or NK cells yielded no significant differences of directional TF binding. The pathway analyses based on differentially regulated (case vs. control) peak – gene links (Figure 2) for chromatin peaks per cell lineage show monocyte predominance in pathway enrichments similar to scRNA/snRNA analyses with the early time point in CD8 T-cells and B-cells clustering together with monocytes (Figure 4d). We also looked for regionally clustered open chromatin signal differences between cases and controls to explore potential chromatin domain-wide differences in regulation. We observed 1,727 adjacent chromatin associations across time points (210 expected by chance) and noted many regions showing dense clustering (Supplementary Figure 4). These may reflect higher order chromatin domains acting differently among cases and controls.

T1D genetic loci and sn-epigenomes

We next sought to analyze differences in snRNA expression and open chromatin among cases and controls through the lens of genetic risk using loci from genome-wide mapping of T1D risk^4^ including each of the 144 independent loci associated at FDR<0.01. We associated each index single nucleotide variant (SNV) for peak genotype with open chromatin or gene expression in the five cell lineages (CD4T, CD8T, NK, MONOCYTE, B CELL) at each time point (T1/T2/T3) using Spearman rank test and adjusted for multiple testing. We identified 244 associations (51 gene expression-SNV and 193 open chromatin-SNV associations within 1Mb of SNV, respectively) at 57 T1D loci (Supplementary Table 3). We examined these genetically linked loci in our case-control differential gene expression/open chromatin data and observed that 11/57 loci (19%) also showed significant differences between cases and controls. However, the rate of case vs. control differences was not overall enriched among T1D SNV linked loci as compared to snRNA traits overall. Separately, we queried snRNA expression differences from cases and controls for each of the named candidate genes at 144 loci^4^, i.e. irrespective if a significant snRNA QTL was detected in our data. Among these, 94 (66%) were observed at least in one cell lineage/time point to be differentially expressed (10% FDR) in snRNA expression (Supplementary Table 4), where most genes appeared on multiple occasions (total n of case-control differences = 412). Among the top 50 most significant differences in disease genes between cases and controls (Table 3) we note relative upregulation of NFKB1 and TNFAIP3 across several cell lineages and time points among cases. In contrast, cases show predominant downregulation of HLA-DQB1 and BACH2. Similarly, regions with genetic loci queried for differential chromatin activity between cases and controls revealed 155 independent peaks within 100kb of index variant including 82 loci (Supplementary Table 5). The chromatin signals are significantly more common (42% excess, c^2^=11, P<0.001) at time point 1 in genomic intervals with T1D SNV as compared to genome wide differential peak accessibility between cases and controls, but as noted above the vast majority of these peaks do not reach significant genetic association with the T1D lead SNV.

T1D endotypes and replicated snRNA traits

To gain insight into potential heterogeneity of single nuclei signals among cases we queried whether the order of appearance (IAA first strict criterion n=21 / GAD first looser criterion n=11) of autoantibodies as an endophenotype could explain any of the overall case – control patterns observed. We matched each subphenotype to its paired control for analyses to ensure age-related factors would not interfere with the analyses and compared the subgroup case against its control at each cell. We then queried if a significant (qv<0.1), replicated (replication criteria across data layers as above) snRNA expression difference overlapped and matched the pattern seen in replicated snRNA/scRNA tests (Figure 3) between cases and controls. On average cell lineage specific up -or downregulation in cases vs. controls showed overlap of 25% with either “IAA first” or “GAD first” signal. However, in early time points (T1/T2) for monocytes the patterns in “GAD first” subgroup of patients were overlapping larger fraction of genes in same direction as seen in all cases and CD4 T-cell first timepoint showed larger overlap with “IAA first” subgroup (Figure 5a). Consequently, the most prominent differences in all cases vs. controls in pathway analyses (Figure 4) seen for early monocyte overexperessed genes were recapitulated in “GAD first” subgroup data (Figure 5b), whereas the “IAA first” subgroup shows enrichment for adaptive immune system networks in CD4 T-cells (Figure 5c) before IAA-appearance (Timepoint 1).

DISCUSSION

Exploration of molecular features of human complex disease by genomics is rapidly expanding beyond DNA sequences (genotyping) alone to layers of gene regulation and gene expression^25^. Single-cell (or nuclei) experiments are most recent additions to “multiomics” tools, offering to reduce confounding impact of primary tissue heterogeneity. However, all multiomics approaches are limited by access to disease relevant tissues and generally by measuring tissue response to chronic disease state. The latter scenario, as demonstrated in a population single-cell study for systemic lupus^26^, can identify known molecular differences in cases vs. controls, but they are unable to provide prognostic or etiologic biomarkers.

The precise drivers of T1D remain unknown, but islet autoimmunity often occurs months to years before clinical disease onset. Previous studies after the appearance of islet autoimmunity or T1D have identified differences in immune cells, including expansions of islet-reactive T cells, between individuals with islet autoimmunity and healthy controls^27^. Moreover, the targets of islet autoantibodies as well as the number of epitopes targeted are associated with differential increased T1D disease risk. One challenge of studying T1D in children is separating the impact that age has on the developing immune system from true changes in children with islet autoimmunity^28^. Here we benefit from longitudinal sampling of a cohort of children at risk for T1D to pursue the first single cell multiomics case – control analyses months to years before the first clinical biomarkers appear or overt disease is present. The observation of differential interferon and cytokine signaling early before autoantibody formation and during seroconversion yields support from earlier gene expression studies in smaller, cross-sectional samples^29^. In contrast the genetically associated loci for T1D risk showed only modest enrichment for the signal in early time points.

Most strikingly, the strong monocytic lineage representation in this response in peripheral blood and the relatively weaker differential coordination of these gene networks close to clinical diagnosis are a novel observation with important implications for follow-up studies. The endophenotype subgroup analyses suggests that this phenomenon is driven more prominently in patients with GAD antibodies appearing first. In a recent pancreatic single-cell study of a rat model of autoimmune diabetes at pre-diabetic stage highlighted increased monocyte/macrophage chemokine and interferon response at onset of autoimmunity^30^.

Despite measuring millions of cells and nuclei in the experiment and recovering 670,000 high quality, singleton nuclei and cells for analyses the data for individual data layers suggests moderate discovery power correlated with reads per measured element and most prominently affecting single nuclei ATAC-seq data. However, the multiomic approach and primary analyses based on signals replicated in independent data types and/or second time point we present here was deployed to ensure robustness of data. The gene network enrichments point towards strong activation of innate immunity, interferon signaling, chemo/cytokine signaling, immune cell activation and inflammation biased to early (T1 and T2) time points and stronger in genes relatively upregulated in cases. A previous study performing longitudinal whole-blood transcriptome analysis of selected TEDDY participants before the onset of T1D also identified an innate immune cell gene signature that was associated with the individual T1D risk^31^. In that study an age-dependent natural killer (NK) cell gene signature correlated with progression to T1D. We had lower power for NK discovery and in TEDDY the majority of at-risk patients did not have familial T1D used as inclusion criteria for cases in TRIGR.

Some previously highlighted genetic risk factors (e.g. BACH2, TNFAIP3) that are also quantitative traits were differentially expressed among cases and controls. However, most case/control trait differences in proximity of disease loci did not reach significance for T1D SNV association, i.e. despite being associated with case status the differences were not driven by genetic disease association detectable in our sample size. However, the excess of early (T1) open chromatin traits in proximity of T1D loci suggests that studying genetic influences on immune regulation predisposing to T1D could require monitoring regulatory states even earlier in infancy.

While the study in TEDDY highlighted the need for longitudinal measures of the dynamic infant immune system and demonstrated that distinct immune signatures are present prior to T1D onset and could be predictive, the study lacked single-cell resolution that is added by our work. In both cases, identifying validated immune signatures that could stratify T1D disease risk will allow for early monitoring of infants at risk. Improved predictive value of biomarkers based on single nuclei in patients at risk could allow more targeted testing of interventional therapies to prevent T1D.