A comprehensive genetic map of cytokine responses in Lyme borreliosis

Abstract

The incidence of Lyme borreliosis has risen, accompanied by persistent symptoms. The innate immune system and related cytokines are crucial in the host response and symptom development. We characterized cytokine production capacity before and after antibiotic treatment in 1,060 Lyme borreliosis patients. We observed a negative correlation between antibody production and IL-10 responses, as well as increased IL-1Ra responses in patients with disseminated disease. Genome-wide mapping the cytokine production allowed us to identify 34 cytokine quantitative trait loci (cQTLs), with 31 novel ones. We pinpointed the causal variant at the TLR1-6-10 locus and validated the regulation of IL-1Ra responses at transcritpome level using an independent cohort. We found that cQTLs contribute to Lyme borreliosis susceptibility and are relevant to other immune-mediated diseases. Our findings improve the understanding of cytokine responses in Lyme borreliosis and provide a genetic map of immune function as an expanded resource.

Fig. 1. Cytokine responses are organized around the different cytokines measured.

a Overview of the LymeProspect project. Schematic created with BioRender.com. b Heatmap of cytokine-stimulation correlation. Samples were taken at the start of antibiotic treatment. Spearman’s correlation of the cytokine concentrations. c–e Heatmap of the significance of associations. *FDR <0.05, **FDR <0.01, ***FDR <0.001. Created with BioRender.com. c Paired two-sided Wilcoxon signed-rank test of cytokine values before and after antibiotic treatment. d Two-sided Wilcoxon rank-sum test of cytokine responses in patients with disseminated Lyme borreliosis and patients with erythema migrans. e Spearman’s correlation of C6-ELISA, i.e., B. burgdorferi s.l antibody index., and cytokine responses.

Fig. 2. Genome-wide and study-wide significant cytokine QTL.

a Manhattan plot of cQTL. Colored loci indicate genome-wide significant cQTL. The colors indicate the cytokine for which the genome-wide association is found. b Table of study-wide significant loci (P < 1.5 × 10⁻⁹). Each locus is described by its top SNP and cytokine-stimulation for which it is found. The eQTL column indicates expression QTL effects, extracted from the eQTLgen consortium. Genes are labeled as “Bb-responding” if they are also responding to B. burgdorferi. All p values were calculated using a linear model, associating cytokine measurements with genetic variants. Multiple-testing correction limits were set to genomewide significant (P < 5e-8), and study-wide significant (P < 1.5 × 10⁻⁹), based on the number of independent tests.

Fig. 3. Across condition consistency of genome-wide cQTL.

a Heatmap of cQTL signed −log₁₀ p values for the different conditions studied. b, c Between conditions comparison. The top plot shows the comparison between t-statistics for both QTL mappings. The plot at the bottom shows the relationship between minor allele count and replicability, as shown by the minimum absolute t-statistic between the two conditions compared: both time points in (b) and IL-6 compared to IL-10 production in (c). d, e Two loci colocalized with blood eQTL signals. The Y-axis represents significance. The first plot is the locuszoom plot of cQTL. The second plot, locuszoom plot of blood eQTL. The third plot is the posterior probability of having a colocalized signal. P values in (a, d, e), and t-statistics in (b, c) were calculated using a linear model, associating cytokine measurements to genetic variants. P values are shown without multiple-testing correction.

Fig. 4. Exploring the TLR1-6-10 locus: Fine mapping and differences across conditions.

a Locuszoom plot of the top association found for the TLR locus. b Fine-mapping of the TLR locus per condition. The color indicates the probability of the variant being causal. PIP posterior inclusion probability. c Comparison of two TLR2 ligands. Same as Fig. 3b–d, d Highlighting the specific regulation of IL-1Ra responses to P3C. The Y-axis represents the log2 transformed cytokine concentrations; the X-axis, the different genotypes present for rs5743618 (n = 1060). e Replication of the specific regulation in the RESIST senior individual (SI) cohort (n = 100). The Y-axis represents library-size corrected gene expression, X-axis, different genotypes present. f Schematic of the specific regulation of TLR1-2 responses. Created with BioRender.com. Boxplots in d and e represent the median value (middle line), 25th and 75th quantiles (box limits), and whiskers extend to 1.5 times the interquartile range. P values in (a), and t-statistics in (c) were calculated using a linear model, associating cytokine measurements to genetic variants. P values are shown without multiple-testing correction. *P < 0.05 **P < 0.01 ***P < 0.001.

Fig. 5. The identified cQTL are relevant for different immune-mediated diseases.

a Heatmap indicating the significance of the cQTL (rows) in public GWAS results from immune-mediated diseases (IMD) (columns). b–e Colocalization of the significant associations (FDR <0.05) between cQTL and IMD GWAS. f Causal link between cQTL and IMD GWAS. All p values in (a–e) were calculated using a linear model, associating cytokine measurements to genetic variants.

Fig. 6. cQTL shed light on Lyme borreliosis susceptibility.

a Lyme borreliosis susceptibility locus previously reported. b Locuszoom plot of the same locus for the ratio between IL-6 and IL-10 production upon bbmix stimulation in whole blood at baseline. Colored by inflammation induction (higher IL-6 to IL-10 ratio). c Same as b, at the 6-week time point. d Schematic of the proposed effect exerted by the Lyme borreliosis susceptibility locus. Created with BioRender.com.