The pediatric PreventCD cohort has been extensively described previously (17). Briefly, the at-risk children included in this study were newborns positive for HLA-DQ2 and/or -DQ8 from families in which at least one direct relative has been diagnosed with CeD. We excluded children with trisomy 21 or Turner syndrome and premature infants born at <36 weeks. The PreventCD samples used for this study were randomly selected based on donor disease status, sex and age. We used PBMC samples from 11 patients (confirmed diagnosis based on histopathology of small intestinal biopsy) before and after seroconversion for CeD and from 10 HLA genotype-, age- and sex-matched controls who did not develop CeD during the course of the PreventCD study. After selection of these samples, all experiments were performed in a blinded way. This study was conducted according to the guidelines laid down in the Declaration of Helsinki and it was approved by all medical ethics committees of the participating centers.

PBMCs were obtained from blood samples by gradient separation in Ficoll-Paque plus™, according to the manufacturer’s recommendations. Isolated PBMCs were cryopreserved in medium containing DMSO 10% and FCS 40% at -80°C, as described in Nazapour et al. (18), until use. Samples were thawed and prepared in batches of 10, consisting of 5 different donors sampled at two different timepoints. Upon thawing, cell concentration and viability were determined using a hemocytometer and Trypan Blue staining, respectively. Cells were then pooled into samples that contained 5 samples of different individuals (~1000 cells per donor).

Single-cell cDNA libraries of PBMC pools were generated using the 10X Genomics platform following the manufacturer’s instructions (document CG00026) and as previously described (16). Briefly, each sample pool was loaded in one lane of Single Cell A Chip (10X Genomics, 120236), and the single-cell RNA was captured by beads, retro-transcribed into single-cell cDNA and amplified using the Single Cell 3’ Library & Gel Bead kit version 2 (10X Genomics, 120237) and the i7 Multiplex kit (10X Genomics, 120262). These libraries were sequenced according to 10X Genomics guidelines using the Novaseq 6000 S2 reagent kit (100 cycles) following a paired-end configuration. The final average sequencing depth was ~50 per cell. Raw sequencing data was processed using CellRanger v1.3 software with default settings. Sequencing reads from single cells were aligned to the hg19 reference genome.

Donor DNA was genotyped using the Infinium GlobalScreeningArray-24v1.0. Next, we demultiplexed the data of the pooled samples based on their genotype using the Demuxlet algorithm, as previously described (16, 19). We applied the tool Opticall v.0.7.0 (20) to call genotypes using default settings. All samples had a call rate > 0.99. We did not identify any related individuals within the study group but did identify one individual of admixed ancestry who was kept in order to maximize our sample size. Cell doublets mapped to multiple donors were removed. We used the R package Seurat v.3.2 (21) to perform further analysis. We included cells expressing >200 and < 3000 unique genes and with < 15% mitochondrial counts for the downstream analysis.

The filtered single-cell dataset was analyzed to obtain the most variable genes across samples and cells, and the first 12 principle components were used for an unsupervised clustering with a resolution of 0.6 using the R package Seurat v.3.2 (21). Thereafter, to achieve an accurate annotation of cell cluster, we took a combined approach based on both well-annotated PBMC data and the known marker genes list of PBMC subsets. Specifically, we firstly applied a canonical correlation analysis to project the data onto a larger dataset of well-annotated PBMCs (16). The annotation of cell-types in the reference was then used to identify the cell populations in the query dataset. Second, the resulting cell types were validated by checking the expression of previously characterized marker genes of PBMC subsets (22).

The generated data are from the following four categories of samples: data obtained from patients before and after seroconversion (CeD T0 and CeD T1, respectively) and from HLA genotype-, age- and sex-matched controls (CTR T0 and CTR T1). Using the “MAST” approach (23), we performed differential gene expression analyses comparing: 1) CeD T0 vs. CTR T1, 2) CeD T1 vs. CTR T1, 3) CeD T0 vs. CeD T1 and 4) CTR T0 vs. CTR T1. We only characterized transcripts that were expressed in > 10% of the cells of each subset and filtered in DEGs with an absolute log2FC > 0.25 and an adjusted p-value < 0.05. We then checked the data for the presence of 118 genes previously genetically associated with CeD and prioritized to play a possible causal role in CeD (24).

The R package clusterProfiler v.3.14.3 was used to investigate whether the DE transcript sets were enriched for genes involved in specific pathways (P adjusted value < 0.05) (25).

CellPhoneDB (26) was applied to the identified DEGs to infer immune ligand–receptor interactions in our dataset using the default parameters for the statistical analysis (receptor and ligand expressed in at least 10% of cells, 1000 permutations, FDR < 0.05). The R package ggplot2 was used to generate a heatmap visualizing the significance and the mean expression levels of these ligand–receptor gene pairs in different cellular subsets.

To study the differences at the onset of CeD at the transcriptomic level, we isolated PMBCs from 21 children included in the PreventCD cohort (17), who all were at high-risk of developing CeD. These samples were classified into cases or controls based on whether the donors developed CeD during the study, confirmed by histopathology of small intestinal biopsy. Samples (n=22) were collected from the same patients (n=11) before (T0) and after (T1) seroconversion for CeD as determined by the anti-TG2 antibody test (U/milliliters > 6). Control samples (n=20) from HLA genotype and sex matched donors (n=10) were taken at similar ages to T0 (n=10) and T1 (n=10) from CeD patients ( Figure 1A and Supplementary Table 1 ). The average time between T0 and T1 was 16 months. In total, 19,663 single cells were profiled. After quality control by filtering based on possible doublets, the number of genes expressed (included cells with >200 and <3,000 genes) and low quality cells (included cells with <15% mitochondrial transcript reads)( Supplementary Figure 1 ), we retained 9,559 cells for subsequent analyses.

To annotate the cell types in our dataset, we applied both a standard clustering method and a reference-based approach using an external dataset of PBMCs from healthy controls that were processed in a similar manner (16). In total, we identified 10 cell subsets in the PBMC population: CD4+ T cells, CD8+ T cells, B cells, plasma cells, natural killer (NK) cells, classical monocytes (cMonocytes), non-classical monocytes, myeloid dendritic cells (mDCs), plasmacytoid dendritic cells (pDCs) and megakaryocytes ( Figure 1B ). The annotation of these subpopulations was validated by comparing the DEGs per cluster with known marker genes for each cell type (22).

After cell cluster annotation, we compared the cell frequencies of each subpopulation between disease states (cases vs controls) and between time points (T0/before seroconversion vs T1/after seroconversion) ( Supplementary Figure 2 ). We tested whether the cell subset composition was significantly changed prior to or after CeD onset in CeD patients, excluding samples with fewer than 30 cells, and found that cell proportions were relatively stable in the blood of CeD patients compared to controls ( Supplementary Figure 2 ). It is important to note that the main tissue affected in CeD is the small intestine and that the cell composition of this compartment will change upon CeD onset, for instance due to an expansion of lymphocytes (8, 27). However, such a compositional change does not seem to be reflected in the blood of our pediatric samples.

To investigate the changes in gene expression profiles within each cell population, we performed differential expression analysis across disease status and across time of sampling ( Supplementary Figure 3 ). DEGs were determined as specified in the Methods (FDR < 0.05 and absolute log2-fold change (log2FC) > 0.25) upon pooling each cell type in all samples from the same disease condition and sampling time, ending with four different comparisons: CeD T0 vs. T1, controls T0 vs. T1, CeD T0 vs. controls T0 and CeD T1 vs. controls T1 (as exemplified for CD4+ T cells in Figure 2A ). We observed a varying number of DEGs in each cell type and category. CD4+ T cells showed the highest number of DEGs (n=467 genes), including 250 markers that were unique to only one of the comparisons, while the other DEGs were present in at least two of the comparisons ( Figure 2A ). The high number of DEGs in CD4+ T cells could be related to the increased power in this compartment as it is the most abundant cell type in our samples.

To better understand the role of DEGs in the context of CeD, we classified the DEGs into three categories: time-related genes, genes characteristic of active CeD and pre-diagnostic CeD markers. In the time-related gene category, we included all DEGs that overlapped or were unique to the T0 to T1 comparison within the controls. These genes are unlikely to be relevant for CeD and may be differentially expressed simply as a consequence of aging (these DEGs are indicated in blue in Figure 2 ). Active CeD genes were defined as DEGs from the comparison of cases vs. controls after seroconversion and from the comparison of T0 and T1 timepoints in cases, i.e. before seroconversion vs. after seroconversion (these DEGs are indicated in red in Figures 2A, B and Supplementary Table 2 ). Finally, we categorized the DEGs from the comparison of cases vs controls before seroconversion as pre-diagnostic CeD markers (indicated in yellow in Figures 2A, B and Supplementary Table 3 ).

Several cell populations showed distinct DEGs that were present before clinical onset of CeD ( Figure 2B and Supplementary Figure 3 ). The top 4 DEGs for each comparison in the five most abundant cell types are depicted in Figure 3 . These results demonstrate that CeD patients display a PBMC transcriptome profile that differs from that of healthy individuals. The existence of PBMC subsets with an altered transcriptome in pre-seroconversion conditions indicates that several biological processes may be altered in pediatric CeD patients even before they develop villous atrophy.

The CD4+ T cell population exhibited the largest number of DEGs, including 358 DEGs classified as active CeD markers and 71 classified as pre-diagnostic markers ( Figure 2A and Supplementary Tables 2 , 3 ). In Supplementary Figures 4A, B , we show that the changes were robust for each category. Some of the upregulated genes with the highest log2FC in the active CeD category included genes associated with T cell functions, such as SELL (also known as CD62L or l-selectin), S100A4 and CD52 ( Figure 3 and Supplementary Figure 4A ). SELL and S100A4 are involved in the migration of T cells (28, 29), whereas CD52, which is expressed in antigen-activated T cells, can induce suppression of other T cells (30, 31). Other top DEGs encode well-known transcription factors, including FOS (associated with a broad range of processes in the T cell activation) (32) and KLF2 (involved in cell motility) (33). The increased expression of these genes is therefore in line with a deregulation of CD4+ T cells in active CeD. In the pre-diagnostic marker category, the top genes were associated to diverse functions, such as thioredoxin interacting protein TXNIP (a regulator of cellular redox state whose function in CD4+ T cells seems to affect the proliferation of effector T cells) (34), and IL32 (previously associated with the progression of various inflammatory disorders including inflammatory bowel disease and gastric inflammation) (35, 36).

We next explored which pathways were likely affected in each cell type by performing pathway enrichment analysis using the REACTOME database ( Supplementary Figures 5 – 9 ). For CD4+ T cells, the genes upregulated after the seroconversion for CeD ( Figures 4A, B ) were significantly enriched for pathways related to immune-related functions (P adjusted value < 0.05) such as “TCR signaling” and “Generation of second messenger molecules” and mainly included genes such as CD3G, CD3E and CD3D ( Figure 4A and Supplementary Table 4 ). In contrast, the downregulated DEGs were significantly enriched for transcripts involved in “Signaling by Interleukins” and “mRNA splicing”, amongst other pathways ( Figure 4A and Supplementary Table 4 ). Some of the downregulated genes involved in interleukin signaling included important transcription factor–encoding genes such as NFKB2 and NFKBIA. Both these genes code for proteins that have a regulatory function in NFκB signaling: NFKB2 codes for p100, which can be processed into p52, a protein that act as a transcriptional repressor when it is in the homodimer form, and NFKBIA codes for IkBa, which inhibits NFκB signaling by retaining NFκB complexes in the cytosol (37). The downregulation of both inhibitory proteins may indicate an enhanced pro-inflammatory NFκB response. Additionally, some other genes found in the “Signaling by Interleukins” pathway, such as SOCS1 and SOCS3, also code for suppressor proteins in the cytokine signaling (38), indicating a loss of inhibition of this pathway. NFκB pathway is known to be upregulated in CeD and plays an important role in the regulation of the inflammatory and immune response by controlling the expression of pro-inflammatory cytokines, adhesion proteins and enzymes in the small intestine mucosa (37, 39, 40).

In the case of the pre-diagnostic markers for CD4+ T cells, the upregulated genes included ribosomal proteins (e.g. RPL39, RPL37, RPS28 and RPS21) that showed enrichments in pathways related to translation and “Infectious disease” ( Figure 4B and Supplementary Table 5 ). This upregulation in the ribosomal proteins could be a response of CD4+ T cells that are already being stimulated in this early stage and could start preparing for the production of proteins. Downregulated genes were more abundant and included cytoskeleton genes (e.g. ACTG1, TUBB2A and TUBA4A) and HLA genes (e.g. HLA-C and HLA-A). Some of the enriched pathways of the downregulated genes were pathways related to antigen presentation and interferon signaling ( Figure 4B and Supplementary Table 5 ).

NK cells showed the highest number of DEGs classified as pre-diagnostic markers (a total of 125) ( Figure 2A and Supplementary Tables 2 , 3 ). The top genes with the highest log2FC were TXNIP, PRF1 and GZMA ( Figure 3 , shown in red in Figure 5A and Supplementary Table 3 ), which are involved in NK cell activation, delivery of granzymes to target cells and cytotoxicity, respectively (41, 42). Pathway enrichment analysis showed that the upregulated pre-diagnostic markers are involved in interactions between lymphoid and non-lymphoid cells (e.g. SELL, CD247 and LAIR2) and phagocytosis (e.g. RAC1, ACTB and ARPC4). Conversely, downregulated DEGs are enriched in “translation processes”, likely caused by the high number of ribosomal proteins that are downregulated ( Figure 5B and Supplementary Table 6 ). One of the top downregulated genes was ARID5A, a ribosome-binding protein that has been previously associated with inflammatory autoimmune diseases (43). Overall, our results indicate that NK cells are already dysregulated before the onset of CeD.

We investigated whether cell–cell interactions are altered in CeD based on responding genes by applying CellPhoneDB (26) to all DEGs in the cells from the CeD cases after seroconversion. CellPhoneDB is a tool that can be used to identify signaling crosstalk between immune cell subpopulations. This analysis identified seven significant ligand–receptor pairs (receptor and ligand expressed in at least 10% of cells, 1000 permutations, FDR < 0.05) expressed among the different cell subpopulations ( Figure 6 ). These include potential CD74–MIF signaling interactions from possible antigen-presenting cells (such as Monocytes, DCs and B cells expressing CD74) towards adaptive immune cells (including CD4+ T cells, CD8+ T cells and plasma cells expressing MIF). Moreover, the ligand–receptor pairs ANXA1–FPR1, HLA-DPB1–TNFSF13B, TNF–TNFRSF1B and TNFRSF1B–GRN were found mainly between monocytes, while ICOSLG–ICOS and TNF–ICOS were mainly found between monocytes and CD4+ T cells. Thus, the cell–cell interactions we identified may be disrupted as a consequence of the differential expression of these genes during the onset of CeD.

Like inflammation and environmental factors, genetic risk factors may also assert effects on gene expression in the context of CeD. To date, 43 genetic risk factors have been identified for CeD, excluding the HLA locus, which accounts for approximately 40% of disease risk (44, 45). In order to detect the cell types that are most likely affected by CeD genetic factors, we intersected our DEGs with the CeD risk/causal genes prioritized in another study (24). Among 118 genes potentially causal in CeD, nine (NCF2, RAC2, TAGAP, REL, GPR183, STAT1, IRF1, ICOS and RGS1) were significantly differentially expressed in at least one comparison in our data ( Figure 7A ). Eight of these were observed in both T cells and NK cells ( Figures 7A, B ), suggesting the importance of the genetic background in these cell types in CeD context. One example here is TAGAP, which is a marker for active CeD in CD4+ T cells and a pre-diagnostic marker in NK cells. A potential role for genetic risk factors associated with CeD in the regulation of TAGAP expression is partly supported by public data, which shows that CeD risk allele T (rs1738074) [Z-OR=0.15, P value= 1.04e-15] (46) leads to a downregulation of TAGAP in blood [Z= -9.01, P value= 2.0419e-19] (47). Indeed, eight of the nine genes (NCF2, TAGAP, REL, GPR183, STAT1, IRF1, ICOS and RGS1) are downregulated in cases compared to controls and/or after developing CeD.