Viral infection-related gene upregulation in monocytes in children with signs of β-cell autoimmunity.

OBJECTIVE: The pathogenesis of type 1 diabetes (T1D) is associated with genetic predisposition and immunological changes during presymptomatic disease. Differences in immune cell subset numbers and phenotypes between T1D patients and healthy controls have been described; however, the role and function of these changes in the pathogenesis is still unclear. Here we aimed to analyze the transcriptomic landscapes of peripheral blood mononuclear cells (PBMCs) during presymptomatic disease.
METHODS: Transcriptomic differences in PBMCs were compared between cases positive for islet autoantibodies and autoantibody negative controls (9 case-control pairs) and further in monocytes and lymphocytes separately in autoantibody positive subjects and control subjects (25 case-control pairs).
RESULTS: No significant differential expression was found in either data set. However, when gene set enrichment analysis was performed, the gene sets "defence response to virus" (FDR <0.001, ranking 2), "response to virus" (FDR <0.001, ranking 3) and "response to type I interferon" (FDR = 0.002, ranking 12) were enriched in the upregulated genes among PBMCs in cases. Upon further analysis, this was also seen in monocytes in cases (FDR = 0.01, ranking 2; FDR = 0.04, ranking 3 and FDR = 0.02, ranking 1, respectively) but not in lymphocytes.
CONCLUSION: Gene set enrichment analysis of children with T1D-associated autoimmunity revealed changes in pathways relevant for virus infection in PBMCs, particularly in monocytes. Virus infections have been repeatedly implicated in the pathogenesis of T1D. These results support the viral hypothesis by suggesting altered immune activation of viral immune pathways in monocytes during diabetes.

INTRODUCTION

In type 1 diabetes (T1D), functional pancreatic β‐cells are lost due to an autoimmune reaction. The destruction of β‐cells seems to happen in a T cell‐mediated manner after self‐antigen presentation, but several immune cell populations within both the adaptive and innate compartments are thought to take part in the process.

Activated cytotoxic CD8+ T cells and macrophages are the major contributors in active insulitis, in which they infiltrate the pancreatic islets of Langerhans. While β ‐cell specific CD8+ cells are found at similar frequencies in the peripheral circulation of healthy donors and patients with T1D, they display markers of antigen‐driven expansion in patients with newly diagnosed T1D. Additionally, the cytotoxic reactivity against islet autoantigens from human samples has been demonstrated. Specific subsets of CD4+ T helper cells have long been known to contribute to the differentiation of B cells into antibody‐secreting plasma cells. Since the most prominent genetic risk for T1D is mediated by the HLA locus, encoding for the class II MHC molecules, and as CD4+ T helper cells are also found in insulitis, CD4+ cells are an attractive candidate for facilitating the emergence of humoral immunity in T1D. A potential model for follicular and peripheral CD4+ T helper cell involvement was recently suggested. In addition, CD4+ T helper cells have been shown to play a critical role in autoreactive CD8+ T cell maintenance.

B cell derived plasma cells produce β‐cell specific autoantibodies which to date are the most important biomarkers of islet autoimmunity before clinical diagnosis of diabetes. , Monocytes are precursors to both macrophages and myeloid dendritic cells and have a role in antigen trafficking and presentation. Their subpopulation compartment sizes have been observed to be altered in T1D patients , , and the cytokine milieu of monocyte populations has been reported to favor more proinflammatory phenotypes.

Despite these discoveries, the exact mechanism underlying T1D development remains largely unknown and heterogeneity in disease pathogenesis is strongly suspected. In this study, to infer biological events during T1D development, we set out to analyze transcriptional differences in peripheral blood mononuclear cells (PBMCs) among subjects with HLA‐conferred risk for childhood T1D and signs of advanced β‐cell autoimmunity and autoantibody negative control subjects.

MATERIALS AND METHODS

Study subjects

The study subjects were participants in the Finnish Type 1 Diabetes Prediction and Prevention (DIPP) study and carried HLA class II genotypes associated with an increased risk for the development of T1D. Subjects with a disease‐predisposing HLA genotype were invited to a prospective follow‐up for signs of β‐cell autoimmunity and dysglycemia. At study visits, the participants were screened for signs of humoral β‐cell autoimmunity: during the early study, for islet cell antibodies (ICA), and if ICA were detected, for biochemical autoantibodies including insulin autoantibodies (IAA), antibodies against the 65 kD isoform of GAD (GADA), and antibodies against the protein tyrosine phosphatase‐related IA‐2 molecule (IA‐2A), from all available samples, including those obtained before seroconversion to ICA positivity. At later stages of the study, all participating children were screened for all four antibodies in samples collected during visits. Diabetes was diagnosed according to WHO criteria. The study protocol was approved by the local ethical committees and an informed consent was obtained from the guardians of the study participants.

The current analysis consists of two parts: a pilot cohort with nine case–control pairs and a confirmation cohort with 25 case–control pairs (Table S1). In the pilot cohort, all case subjects were positive for ICA and at least one biochemical autoantibody (IAA, GADA and/or IA‐2A) at the time of sample collection and had developed T1D during later follow‐up. The PBMC‐samples were stored frozen after sample collection. The first nine subjects for whom such a sample and a healthy control, matched for age at sampling, gender, HLA‐DR/DQ genotype and length of freezing time, were available, were selected in the cohort.

The confirmation cohort comprised 25 case–control pairs. The case subjects tested positive for at least two of the autoantibodies except for five cases having one biochemical autoantibody and ICA. The controls were matched for age at sampling, gender, HLA‐DR/DQ genotype and date of sample collection. In the confirmation cohort, fresh samples were used for cell separation and criteria meeting subjects and controls were selected among children taking part in regular follow‐up visits.

Autoantibody analysis and HLA genotyping

The analysis of the major HLA‐DR‐DQ haplotypes conferring T1D risk was performed using sequence‐specific oligonucleotide probes as described earlier. The protocol for determining ICA, IAA, GADA and IA‐2A in the DIPP study has been described previously. ,

PBMC isolation and sample handling

PBMCs were collected from lithium heparin blood using Ficoll‐Paque Plus density gradient centrifugation and resuspended in RPMI 1640 medium. In the pilot cohort, the cells were stored frozen at −150°C (cryopreserved in 10% DMSO) until analysis. Before RNA isolation, the cells were thawed and lysed in Buffer RLT Plus (Qiagen, Hilden, Germany).

PBMC fractionation

Samples from the confirmation cohort were fractionated and lysed fresh and immediately after PBMC isolation fractionated into monocytes and the remaining PBMC fraction with EasySep Human CD14 positive selection kit II (STEMCELL Technologies, Vancouver, BC, Canada) according to the manufacturer's instructions. The purity of the monocyte and remaining lymphocyte fractions was confirmed by flow cytometry (Table S2). Both cell fractions and unfractionated PBMC were immunostained with anti‐CD3 PE (SK7, BD Biosciences, San Jose, CA, USA) and anti‐CD19 APC (SJ25C1, BD Biosciences) for 30 min at +4°C. The PBMC and remaining fraction were furthermore stained with anti‐CD14 FITC (M5E2, BD Biosciences) for 30 min at +4°C to assess the initial and remaining amounts of monocytes in the sample. After immunostaining the cells were washed twice with phosphate buffered saline (PBS) for 5 min at 2500 rpm with Sorvall MC 12 V (Thermo Fischer Scientific, USA). The cells were fixed with 0.1% formaldehyde in PBS. The samples were analyzed using an Accuri C6 flow cytometer (BD Biosciences). Fractionated cells were lysed in Buffer RLT Plus (Qiagen) and stored at −80°C prior to RNA extraction.

RNA isolation

RNA was extracted from the PBMCs using the RNeasy Plus Mini Kit (Qiagen) in the pilot cohort and RNeasy Plus Micro Kit (Qiagen) in the confirmation cohort according to the manufacturer's instructions. RNA quality and quantity in these cohorts were assessed using the Agilent RNA 6000 Nano Kit (Agilent, Santa Clara, CA, USA) and Agilent RNA 6000 Pico Kit (Agilent), respectively, on a 2100 Bioanalyzer (Agilent). RNA integrity number (RIN) ≥8 was used as RNA‐quality cut‐off for inclusion.

RNA library preparation and sequencing

RNA libraries for the pilot cohort were made using a modified version of the single‐cell tagged reverse transcription (STRT) method, described in detail in Reference 20 to prepare a 48‐plex Illumina‐compatible sequencing library from 10 ng of each RNA sample. Briefly, RNA samples were placed in a 48‐well plate in which a universal primer, template‐switching oligos, and a well‐specific 6‐bp barcode sequence (for sample identification) were added to each well. , The synthesized cDNAs from the samples were then pooled into one library and amplified by single‐primer PCR with the universal primer sequence. The library was sequenced on three Illumina HiSeq2000 (Illumina, San Diego, CA, USA) lanes, using the Illumina TruSeq v3 60‐bp single‐read protocol. Sequencing was performed at the Bioinformatics and Expression Analysis (BEA) core facility at Karolinska Institutet (Huddinge, Sweden).

RNA libraries for the confirmation cohort were made using 20 ng RNA as starting input and the libraries were sequenced on an Illumina NextSeq 500, High Output (75 cycles). Sequencing was done at Biomedicum Functional Genomics Unit (FuGU), University of Helsinki, Finland.

Sequencing data analysis

For the pilot cohort, sequence data was converted to fastq files using Casava 1.8.2 (Illumina), and processed using the STRTprep pipeline available at https://github.com/shka/STRTprep (also described in Reference 20).

For the confirmation cohort sequence data was processed as described previously. Briefly, raw base call (BCL) files were demultiplexed and converted to FASTQ files using Picard tools (v2.10.10; http://broadinstitute.github.io/picard/), and aligned to the human reference genome hg19, human ribosomal DNA unit (GenBank: U13369), and ERCC spike‐ins (SRM 2374) with the GENCODE (v28) transcript annotation by HISAT2 (v2.1.0). The uniquely mapped reads within the 5′‐UTR or 500 bp upstream of the protein‐coding genes were counted using Subread featureCounts (v1.6.2).

After quality check, three controls and two cases were excluded from the PBMC dataset, and one control and one case were excluded from the lymphocyte dataset.

In all three datasets, differential expression analysis between the controls and cases was performed using the R (v3.6.2) package DESeq2 (v1.24.0), where gender was considered as a covariate. Gene set enrichment analysis (GSEA) was performed using GSEA (v4.0.3) using the GSEAPreranked tool, where genes were preranked based on their p‐values and fold changes.

RESULTS

To characterize the profiles of RNA expression in immune cell subsets in children with advanced autoimmunity and compare those to that of matched controls, RNA sequencing was performed from the whole PBMC compartment (pilot cohort) and later from monocytes and lymphocytes separately (confirmation cohort).

In the pilot cohort, case subjects with advanced β‐cell autoimmunity that developed into T1D during later follow‐up were compared to matched control subjects. No significant differences in gene expression between the two groups were observed (data not shown). However, in a subsequent gene set enrichment analysis (GSEA), implemented on RNA sequencing data pre‐ranked based on fold changes and significance of differential expression, differences linked to virus immunity were detected (Table 1). The upregulated genes included gene sets corresponding to the terms “defence response to virus” (FDR <0.001, ranking 2), “response to virus” (FDR <0.001, ranking 3) and “response to type I interferon” (FDR = 0.002, ranking 12; Table S3).

TABLE 1

Gene set enrichment analysis (GSEA) results of upregulated gene genes in pilot and main cohorts

	Pilot						Monocyte						Lymphocyte
Rank	Name	Size	ES	NES	NOM p‐val	FDR q‐val	Name	Size	ES	NES	NOM p‐val	FDR q‐val	Name	Size	ES	NES	NOM p‐val	FDR q‐val
1	Regulation of response to biotic stimulus	92	0.29	3.22	0.00	0.00E + 00	Response to type I interferon	71	0.28	2.88	0.00	0.02	Electron transport chain	145	0.29	4.10	0.00	0.00E + 00
2	Defense response to virus	162	0.21	3.10	0.00	0.00E + 00	Defense response to virus	178	0.18	2.81	0.00	0.01	Cellular respiration	162	0.28	4.04	0.00	0.00E + 00
3	Response to virus	215	0.19	3.10	0.00	0.00E + 00	Response to virus	240	0.15	2.60	0.00	0.04	Small molecule catabolic process	293	0.21	4.01	0.00	0.00E + 00
4	Response to molecule of bacterial origin	197	0.19	3.09	0.00	0.00E + 00	Vesicle organization	250	0.14	2.57	0.00	0.04	Cofactor metabolic process	351	0.18	3.98	0.00	0.00E + 00
5	Positive regulation of defense response to virus by host	20	0.57	3.04	0.00	2.44E − 04	Interferon gamma mediated signaling pathway	71	0.26	2.55	0.00	0.04	Generation of precursor metabolites and energy	389	0.17	3.82	0.00	0.00E + 00
6	Defense response to other organism	273	0.16	2.99	0.00	2.04E − 04	Response to interferon gamma	141	0.18	2.50	0.00	0.05	Mitochondrial translational termination	88	0.33	3.70	0.00	0.00E + 00
7	Response to bacterium	309	0.15	2.97	0.00	5.30E − 04	Ribonucleoprotein complex biogenesis	402	0.11	2.46	0.00	0.06	ATP synthesis coupled electron transport	79	0.36	3.69	0.00	0.00E + 00
8	Lipopolysaccharide mediated signaling pathway	39	0.38	2.83	0.00	1.67E − 03	Regulation of gene silencing	111	0.19	2.38	0.00	0.09	Mitochondrial translation	132	0.27	3.68	0.00	0.00E + 00
9	Positive regulation of cytokine production	282	0.15	2.81	0.00	1.76E − 03	Ribosome biogenesis	261	0.12	2.33	0.00	0.11	Respiratory electron transport chain	95	0.32	3.65	0.00	0.00E + 00
10	Positive regulation of DNA binding transcription factor activity	152	0.19	2.78	0.00	2.08E − 03	De novo protein folding	35	0.32	2.33	0.00	0.10	Anaphase promoting complex dependent catabolic process	78	0.35	3.63	0.00	0.00E + 00
11	Negative regulation of viral genome replication	38	0.39	2.78	0.00	2.00E − 03	Organic cyclic compound catabolic process	465	0.09	2.32	0.00	0.10	Cellular amino acid metabolic process	218	0.21	3.59	0.00	0.00E + 00
12	Response to type i interferon	64	0.29	2.76	0.00	1.94E − 03	ncRNA metabolic process	403	0.10	2.31	0.00	0.10	Mitochondrial respiratory chain complex assembly	89	0.33	3.59	0.00	0.00E + 00
13	Cytokine production	450	0.12	2.71	0.00	3.02E − 03	Multi organism localization	67	0.24	2.29	0.00	0.10	Oxidative phosphorylation	114	0.28	3.50	0.00	0.00E + 00
14	Regulation of defense response to virus by host	29	0.42	2.68	0.00	4.12E − 03	Response to interferon alpha	17	0.45	2.26	0.00	0.11	Aerobic respiration	75	0.34	3.48	0.00	0.00E + 00
15	Response to interferon beta	23	0.45	2.68	0.00	3.84E − 03	Organelle localization	477	0.09	2.24	0.00	0.13	NADH dehydrogenase complex assembly	57	0.38	3.40	0.00	0.00E + 00
16	Leukocyte cell adhesion	218	0.16	2.68	0.00	3.60E − 03	Vesicle targeting	82	0.21	2.22	0.00	0.13	Translational termination	100	0.29	3.40	0.00	0.00E + 00
17	Cellular response to biotic stimulus	148	0.18	2.66	0.00	4.04E − 03	ncRNA processing	343	0.11	2.21	0.00	0.13	DNA dependent DNA replication	134	0.25	3.28	0.00	0.00E + 00
18	Cytokine mediated signaling pathway	476	0.11	2.62	0.00	5.59E − 03	tRNA transport	35	0.31	2.20	0.00	0.14	Mitochondrial gene expression	156	0.23	3.24	0.00	0.00E + 00
19	Cell adhesion	344	0.12	2.58	0.00	7.55E − 03	COPII coated vesicle budding	64	0.24	2.20	0.00	0.13	Energy derivation by oxidation of organic compounds	219	0.19	3.19	0.00	0.00E + 00
20	Regulation of defense response to virus	55	0.29	2.57	0.00	7.78E − 03	Nuclear transcribed MRNA catabolic process nonsense mediated decay	114	0.18	2.19	0.00	0.13	Nucleobase containing small molecule metabolic process	278	0.16	3.15	0.00	0.00E + 00
21	Response to interferon gamma	130	0.19	2.54	0.00	8.98E − 03	ncRNA export from nucleus	39	0.30	2.19	0.00	0.12	Antigen processing and presentation of exogenous peptide antigen via MHC class i	75	0.30	3.11	0.00	3.33E − 05
22	Inflammatory response	370	0.11	2.52	0.00	0.01	rRNA metabolic process	193	0.13	2.14	0.01	0.16	Organic acid catabolic process	187	0.19	3.05	0.00	9.43E − 05
23	Negative regulation of viral process	65	0.27	2.50	0.00	0.01	Membrane fusion	111	0.17	2.13	0.00	0.16	Regulation of cellular amino acid metabolic process	53	0.35	3.02	0.00	1.20E − 04
24	Regulation of multi organism process	253	0.13	2.50	0.00	0.01	Regulation of nuclease activity	20	0.39	2.11	0.00	0.18	Detoxification	90	0.27	2.99	0.00	1.44E − 04
25	Regulation of body fluid levels	229	0.14	2.47	0.00	0.01	Endoplasmic reticulum to golgi vesicle mediated transport	165	0.14	2.11	0.00	0.17	Mitochondrial electron transport NADH to ubiquinone	44	0.38	2.98	0.00	1.38E − 04
26	Regulation of cell adhesion	227	0.14	2.46	0.00	0.01	Cotranslational protein targeting to membrane	94	0.18	2.09	0.00	0.20	Nucleoside phosphate biosynthetic process	216	0.17	2.94	0.00	1.33E − 04
27	Negative regulation of multi organism process	109	0.20	2.45	0.00	0.02	DNA recombination	208	0.12	2.09	0.00	0.19	Regulation of cellular amine metabolic process	66	0.30	2.93	0.00	1.53E − 04
28	Adaptive immune response	261	0.13	2.43	0.00	0.02	Apoptotic DNA fragmentation	19	0.39	2.06	0.01	0.22	Drug metabolic process	219	0.17	2.91	0.00	1.48E − 04
29	Cellular response to interferon beta	15	0.51	2.40	0.00	0.02	Glycosylation	164	0.14	2.06	0.00	0.21	Cellular ketone metabolic process	142	0.21	2.88	0.00	1.67E − 04
30	Regulation of immune effector process	245	0.13	2.36	0.00	0.03	Response to topologically incorrect protein	162	0.14	2.05	0.01	0.22	Regulation of cell cycle g2 m phase transition	197	0.17	2.87	0.00	2.08E − 04
31	Positive regulation of NF‐KappaB transcription factor activity	102	0.20	2.31	0.00	0.04	Positive regulation of defense response	360	0.09	2.04	0.00	0.21	Antigen processing and presentation of peptide antigen via MHC class i	91	0.26	2.82	0.00	3.60E − 04
32	Negative regulation of immune system process	252	0.13	2.30	0.00	0.04	Vesicle budding from membrane	90	0.18	2.03	0.00	0.23	Cellular detoxification	84	0.26	2.81	0.00	4.36E − 04
33	Regulation of response to external stimulus	376	0.10	2.29	0.00	0.04	Transport of virus	54	0.23	2.02	0.01	0.24	Amine metabolic process	102	0.23	2.75	0.00	8.24E − 04
34	Production of molecular mediator involved in inflammatory response	41	0.30	2.29	0.00	0.04	Synapse organization	162	0.14	2.02	0.00	0.23	Proteasomal ubiquitin independent protein catabolic process	21	0.49	2.73	0.00	9.44E − 04
35	Interleukin 6 production	86	0.21	2.28	0.00	0.04	Vesicle localization	192	0.13	2.01	0.01	0.23	Purine containing compound biosynthetic process	154	0.19	2.72	0.00	9.77E − 04
36	Positive regulation of protein kinase b signaling	68	0.24	2.28	0.00	0.04	Spliceosomal SNRNP assembly	35	0.29	2.00	0.01	0.23	Nuclear DNA replication	50	0.32	2.70	0.00	1.11E − 03
37	Interleukin 6 secretion	22	0.39	2.28	0.00	0.04	DNA catabolic process endonucleolytic	23	0.35	2.00	0.00	0.24	Meiotic cell cycle process	122	0.21	2.69	0.00	1.13E − 03
38	T cell mediated immunity	64	0.24	2.27	0.00	0.04	Lymphocyte chemotaxis	32	0.29	1,99	0.01	0.24	tRNA metabolic process	150	0.18	2.67	0.00	1.34E − 03
39	Immune response regulating signaling pathway	390	0.10	2.24	0.00	0.05	RNA catabolic process	336	0.10	1,99	0.00	0.24	Antigen processing and presentation of peptide antigen	170	0.18	2.66	0.00	1.43E − 03
40	Interferon gamma production	68	0.23	2.22	0.00	0.05	Golgi vesicle transport	285	0.10	1,98	0.01	0.24	Cellular protein complex disassembly	182	0.17	2.61	0.00	2.15E − 03
41	Cytokine metabolic process	71	0.23	2.22	0.00	0.05	Recombinational repair	95	0.17	1,96	0.00	0.26	Negative regulation of cell cycle g2 m phase transition	94	0.23	2.60	0.00	2.29E − 03
42	Response to lipid	457	0.09	2.19	0.00	0.06	Nuclear transport	279	0.10	1,96	0.00	0.26	Branched chain amino acid catabolic process	19	0.48	2.59	0.00	2.45E − 03
43	Positive regulation of cytokine secretion	68	0.22	2.19	0.00	0.06	Telomere maintenance via semi conservative replication	23	0.34	1,96	0.00	0.25	Cofactor biosynthetic process	179	0.17	2.58	0.00	2.44E − 03
44	Cytolysis	17	0.44	2.18	0.01	0.06	Golgi vesicle budding	72	0.19	1,96	0.00	0.25	Nucleobase containing small molecule catabolic process	38	0.36	2.57	0.00	2.51E − 03
45	Cytokine production involved in inflammatory response	24	0.37	2.17	0.00	0.07	Regulation of posttranscriptional gene silencing	88	0.18	1,96	0.01	0.25	Nucleoside monophosphate biosynthetic process	38	0.35	2.57	0.00	2.45E − 03
46	Regulation of lymphocyte migration	37	0.29	2.14	0.00	0.08	DNA repair	440	0.08	1,94	0.01	0.27	Tricarboxylic acid cycle	32	0.38	2.57	0.00	2.45E − 03
47	positive regulation of myeloid leukocyte mediated immunity	21	0.38	2.13	0.00	0.09	Vesicle targeting to from or within golgi	67	0.20	1,93	0.01	0.28	Ribonucleoside catabolic process	17	0.52	2.56	0.00	2.54E − 03
48	Cytokine production involved in immune response	64	0.23	2.11	0.01	0.09	RNA export from nucleus	129	0.14	1,93	0.01	0.27	DNA conformation change	199	0.16	2.56	0.00	2.53E − 03
49	Positive regulation of ERK1 and ERK2 cascade	85	0.20	2.11	0.00	0.09	Establishment of protein localization to endoplasmic reticulum	106	0.16	1,93	0.01	0.27	Monosaccharide catabolic process	32	0.38	2.56	0.00	2.54E − 03
50	Cytokine secretion	120	0.17	2.11	0.00	0.09	Defense response to other organism	310	0.09	1,93	0.00	0.26	Fatty acid beta oxidation	62	0.28	2.55	0.00	2.64E − 03

Note: GSEA was performed with RNA sequencing data pre‐ranked based on fold changes and significance of differential expression. The analysis revealed a monocyte specific upregulation of gene sets relating to viral response and response to type I interferon in autoantibody positive cases. The table details the enrichment score (ES), normalized enrichment score (NES), nominal p‐value (NOM p‐val) and false discovery rate corrected q‐value (FDR q‐val) for each term in pilot and main cohorts.

A confirmation cohort, comparing case subjects with advanced β‐cell autoimmunity and matched control subjects, was then analyzed to further investigate these findings. In this cohort, fresh PBMC samples were separated into monocyte and lymphocyte compartments to study the role of monocytes in viral and type I interferon responses observed in the pilot. Both fractions were analyzed separately.

As in the pilot cohort transcription profiles, there were no significantly differentially expressed genes when comparing cases and controls (data not shown). However, as in the pilot cohort, the GSEA analysis suggested differences in virus‐associated immune activation between case and control subjects in the monocyte compartment (Table 1). The GSEA confirmed the terms “defence response to virus” (FDR = 0.02, ranking 2), “response to virus” (FDR = 0.04, ranking 3) and “response to type I interferon” (FDR = 0.02, ranking 1) among upregulated genes (Table S3). In contrast, enrichment of these gene sets between cases and controls could not be observed in the lymphocyte compartment in the GSEA analysis.

DISCUSSION

Various immune cell populations are implicated to play a role in the β‐cell destruction leading to T1D. However, factors affecting altered immune activation are not fully described. Understanding the differences in the distinct immune cell compartment function might provide essential information about the pathogenesis of T1D. Here we explored transcriptional profiles in PBMC of children with advanced β‐cell autoimmunity and compared them with those of autoantibody negative children matched for sex, age and HLA. The study was conducted in two parts, first a pilot cohort performed with frozen PBMC and second, a confirmation cohort with fresh PBMC that were fractionated into monocytes and remaining lymphocytes. While statistically significant gene expression differences could not be observed, three gene sets associated with the terms “defence response to virus,” “response to virus” and “response to type I interferon” were consistently upregulated in PBMCs and further in monocytes of case subjects.

Viral infections have long been linked with T1D pathogenesis. Especially enteroviral infections have been found to associate with increased risk for disease onset , and this has also been seen in the DIPP cohort. Many strains are known to be able to cause chronic systematic infections as well as infect the pancreas. According to the current understanding, these conditions may drive strong inflammatory responses and autoimmunity. In our present study, the observed upregulation of genes essential in response to virus infections was detected in PBMCs but in the further analysis the finding was restricted to peripheral blood monocytes. Innate immunity is classically responsible for the acute response to viral threats, but the combination of a lack of detectable response from lymphocytes and our specific set of three significant GSEA terms also suggested that the monocytes themselves could be infected with a virus. Coxsackie virus B4, which belongs to the group of enteroviruses, has been shown to infect monocytes and monocyte‐derived macrophages, , with the potential to establish a persistent infection. Monocyte derived macrophages also produce a strong cytokine response, including IL‐6 and TNFα, to Coxsackie virus B4. Another study by Alidjinou et al reported that enteroviral RNA could be detected in monocytes of some T1D patients, although viral loads in many cases seemed low and difficult to detect with RT‐PCR. Furthermore, the presence of enteroviral RNA coincided with the presence of IFNα mRNA in most subjects. It is therefore possible that some of the cases in our study may have an ongoing enteroviral infection, reflected both by the upregulation of virus response genes and type I interferon response genes.

Innate immune function accompanied by a type I interferon signature, that is, detectable transiently starting shortly before seroconversion, has been reported in longitudinal studies investigating T1D pathogenesis. Kallionpää et al detected this signature in whole blood transcriptomics of autoantibody positive DIPP children, starting before seroconversion and persisting until diagnosis of clinical disease. Enterovirus‐associated transcriptomic profiles were also observed in a subset of these children. Similar findings to Kallionpää et al were evident in the Environmental Determinants of Diabetes in the Young study (TEDDY) among the children whose first islet autoantibody was against insulin. Interestingly, the association of Coxsackie B1 enterovirus infections and islet autoimmunity was found specifically in children with insulin autoantibodies as the first sign of autoimmunity in the DIPP study. Enterovirus‐associated transcriptomic profiles were observed in a subset of these children. Ferreira et al reported a transient type I interferon signature in genetically predisposed children before the autoantibodies were developed, but not in children with existing disease. Our observation, that a gene set corresponding to the term “response to type I interferon” is upregulated in peripheral blood immune cells, and particularly in monocytes, is in line with these previous observations.

Several studies have explored peripheral blood transcriptomic signatures in the context of T1D from various angles. Stechova et al compared the transcriptional profiles of pediatric T1D patients, their clinically healthy first‐degree relatives and healthy, unrelated controls and found that the most significant difference was between first‐degree relatives and unrelated controls. Accordingly, Elo et al did not observe differences in gene expression profiles of children positive for β‐cell autoantibodies and children who had progressed to T1D. Similarly, in a study investigating monocytes of twins discordant for T1D, healthy twin pairs and healthy singleton controls, Beyan et al saw that most of the abnormally expressed genes observed in T1D twins were also abnormal in their non‐diabetic twins. It would therefore seem like gene expression differences already exist in genetically predisposed but healthy individuals. Additionally, many previous findings of differential gene expression in PBMCs in the context of T1D have been made with T1D patients , , or a combination of presymptomatic cases and those diagnosed with the disease , compared to healthy controls. Observations concerning peripheral blood monocytes have similarly been made predominantly in patients with existing T1D , , , and could be attributed to the metabolic crisis and ongoing stress triggered by disease onset, as the loss of glucose tolerance appears only shortly before it. As a consequence, it is possible that the changes our cases have, especially in monocytes, are difficult to distinguish due to some of the strengths of this study: close genetic matching of cases and controls and using samples predating the metabolic state caused by T1D itself. Therefore, the controls in our study may also have changes in their PBMCs because of the genetic T1D‐risk they carry and immunological changes in the early phase of disease progression are likely to be relatively minute compared to those during disease onset.

Limitations of this study include the use of peripheral blood cells, limiting statistical power and in parts of the study, heterogenous populations. Additionally, there is a lack of a control group without HLA‐conferred genetic risk to T1D. It is likely that all these factors contribute to the lack of statistically significant gene expression differences in this study. This could be addressed in future studies by more detailed cell fractionation and possible additional controls. A time series could help to pinpoint the timing of monocyte activation in T1D.

CONCLUSION

Transcriptional profiles of children with advanced β‐cell autoimmunity and those of their autoantibody negative controls matched for age, sex and genetic T1D‐risk did not differ significantly in monocytes or monocyte‐depleted PBMCs. However, gene sets essential in responses to virus were consistently upregulated in PBMCs and specifically in monocytes of subjects with advanced β‐cell autoimmunity. This result supports earlier findings implicating the role of viral infections in T1D pathogenesis and the emergence of β‐cell autoimmunity.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Juha Kere, Jorma Ilonen, and Johanna Lempainen designed the study. Milla Valta, Masahito Yoshihara, Elisabet Einarsdottir, Sirpa Pahkuri, and Sini Ezer, conducted the laboratory analyses. Mikael Knip, Riitta Veijola, and Jorma Toppari, provided study material. Masahito Yoshihara and Shintaro Katayama, analysed the data and Milla Valta, Masahito Yoshihara, Juha Kere, Jorma Ilonen and Johanna Lempainen interpreted the results. Milla Valta drafted the manuscript, Masahito Yoshihara, Elisabet Einarsdottir, Mikael Knip, Jorma Toppari, Juha Kere, Jorma Ilonen and Johanna Lempainen reviewed the manuscript and contributed to the discussion.

ETHICS STATEMENT

The study protocol has been approved by the local ethical committees. The guardians of the study subjects have given informed consent to study participation.

TABLE S1Attributes of study subjects in the confirmation cohort.

TABLE S2:Monocyte fraction purity in the confirmation cohort. The monocyte fraction was separated using magnetic separation and the purity was analyzed with a simple panel where CD3+ T cells, CD19+ B cells and CD14+ monocytes were detected with flow cytometry. The monocyte fraction contained on average <5% T and B cell contamination.

Supplementary Table 1: Lists of genes from the statistically significantly upregulated gene sets in the gene set enrichment analysis. A) “response to virus” in pilot cohort, B) “response to virus” in monocytes of main cohort, C) “defense response to virus” in pilot cohort, D) “defense response to virus” in monocytes of main cohort, E) “response to type I interferon” in pilot cohort, F) “response to type I interferon” in monocytes of main cohort. The tables contain the name of the gene, rank in original gene list, test values of rank metric score and running enrichment score (running ES) and information whether the gene contributes to the enrichment result of the gene set (core enrichment).

Click here for additional data file.