Predicting Lyme Disease From Patients' Peripheral Blood Mononuclear Cells Profiled With RNA-Sequencing

Although widely prevalent, Lyme disease is still under-diagnosed and misunderstood. Here we followed 73 acute Lyme disease patients and uninfected controls over a period of a year. At each visit, RNA-sequencing was applied to profile patients' peripheral blood mononuclear cells in addition to extensive clinical phenotyping. Based on the projection of the RNA-seq data into lower dimensions, we observe that the cases are separated from controls, and almost all cases never return to cluster with the controls over time. Enrichment analysis of the differentially expressed genes between clusters identifies up-regulation of immune response genes. This observation is also supported by deconvolution analysis to identify the changes in cell type composition due to Lyme disease infection. Importantly, we developed several machine learning classifiers that attempt to perform various Lyme disease classifications. We show that Lyme patients can be distinguished from the controls as well as from COVID-19 patients, but classification was not successful in distinguishing those patients with early Lyme disease cases that would advance to develop post-treatment persistent symptoms.

Introduction

Lyme disease (LD) is a tick-borne illness that has become a growing concern in the United States (US) and Canada. LD spreads exclusively by two tick species in the Northern hemisphere, Ixodes scapularis and Ixodes pacificus. The disease is caused by the spirochete bacteria Borrelia burgdorferi sensu stricto and is transmitted to humans through one of the vector's blood meals (1). Comprising 62.6% of all vector-borne diseases, and 81.19% of all tick-borne disease, Borrelia burgdorferi sensu stricto was the most common vector-borne pathogen in the United States from 2004 to 2016 (2). Approximately 30,000 diagnosed cases of LD are reported to the CDC each year, with an estimated true burden of ~300,000 cases (3, 4) and yearly healthcare cost of ~$1 billion in the US (5). Cases in the US are concentrated within the Northeast, Mid-Atlantic, Midwest, and coastal West regions (1). Due to the tick's seasonal lifecycle, transmission of the pathogen and subsequent human infection occurs at higher rates in spring, summer, and the early part of autumn when the various life-stages of the vector quest for a meal (6). Once infected, initial onset of symptoms can manifest ~7–14 days after transmission, although both earlier and later initial onset has been documented (6).

Clinical demonstration of untreated Lyme disease is divided into three stages: the early localized stage, early disseminated stage, and late disseminated stage (6, 7). The early localized stage is characterized by the erythema migrans (EM), a skin lesion that is often red and round or oval shaped that manifests at the site of the bite, and which can be accompanied by mild flu-like symptoms of fever and fatigue. However, the rash is absent or undetected in ~5–30% of cases (7). Symptoms can progress into the early disseminated stage when the patient lacks initial treatment during the early localized stage and the bacteria disseminate hematogenously to other areas of the skin or other organ systems. Multiple EMs may be present around the body, although they are typically less inflamed than the primary EM at the site of the initial tick bite. In addition, the patient may exhibit Lyme carditis or Lyme neuroborreliosis (1, 7). Approximately 6 months after initial disease onset and when left without treatment, the patient may experience signs of the late disseminated stage, specifically late Lyme arthritis (which may occur in up to 60% of patients with a history of primary EM) or neurologic disease (7). Following appropriate antibiotic treatment, a subset of patients experience a range of persistent, significant, but often non-specific symptoms which frequently include fatigue, widespread musculoskeletal pain, and cognitive difficulties, among others (8). Approximately 10–20% of patients will meet criteria for post-treatment Lyme disease (PTLD), which includes the presence of specific symptoms as well as significant impact of these symptoms on life functioning (9). Testing and diagnosis of LD has proven to be difficult or unreliable. Bacteria cultivation requires specialized medium and frequently results in low-yield cultivation (1). The universally-accepted diagnostic test for LD is a two-tier serological test: a positive enzyme-linked immunosorbent assay (ELISA) test followed by a positive Western blot test for IgM and IgG B. burgdorferi antigens (10). Unfortunately, sensitivity is low at 29% during the early localized stage, and the two-tiered serological test is not recommended for early diagnosis in the first few weeks of infection (11). In the absence of a laboratory diagnostic tool, the diagnosis of early LD is reliant on a demonstrated EM that occasionally does not present or is not observed, which can lead many patients to progress to the early disseminated stage and its debilitating symptoms before the disease is diagnosed and treated.

To further our understanding of the molecular mechanisms that lead to LD symptoms, we examined the longitudinal changes in gene expression as a tool for deep phenotyping of diagnosed LD patients and healthy controls. RNA-sequencing profiling was administered from peripheral blood mononuclear cells (PBMCs) collected during multiple patient visits over a 1-year timespan. Deep phenotyping was performed by examining the projection of RNA-seq data into lower dimensions for unsupervised clustering of the patients based on their gene expression vectors. In addition, differential gene expression analysis followed by enrichment analysis was employed to identify upstream regulatory mechanisms and disease phenotypes associated with LD. This analysis was augmented with single cell deconvolution analysis. Since we observed clear separation between controls and cases, we developed several classifiers as potential diagnostic tools and for the further identification of molecular mechanisms underlying LD phenotypes. These include predictors that distinguish healthy controls from Lyme patients, Lyme patients from COVID-19 patients, and whether Lyme patients will advance to develop persistent symptoms.

Materials and Methods

Patient Recruitment

The current study is part of a larger, ongoing prospective cohort study of patients with LD and non-LD controls. Adult patients with early LD were predominantly recruited from primary or urgent care settings in the Mid-Atlantic area of the United States. All participants were required to have a physician-documented EM of >5 cm present at the time of enrollment, thereby meeting CDC criteria for “confirmed” LD (12). Although early LD patients were eligible with up to 72 h of appropriate antibiotic exposure at the time of enrollment, the majority (63.0%) were antibiotic-naive at their first visit. Patients with a prior history of LD or those who had received the LD vaccine were excluded from the study. Patients were also excluded for a range of self-reported prior medical conditions associated with significant immunologic impact and/or subjective symptoms which may overlap with PTLD: chronic fatigue syndrome, fibromyalgia, unexplained chronic pain, sleep apnea or narcolepsy, autoimmune disease, chronic neurologic disease, liver disease, hepatitis, HIV, cancer or malignancy, major psychiatric illness, or drug or alcohol abuse. Non-Lyme infected controls were recruited from similar care settings as LD cases, as well as via community recruitment using flyers and online advertising. Controls underwent an initial screening two-tier antibody test and were required to test negative prior to enrollment. Follow-up tests were conducted at each subsequent time point, with all controls testing two-tier negative for the duration of the study. Controls were also screened for the same prior medical history conditions as cases and were also required to be free of any history of prior clinical LD. Participants with LD and non-LD controls were followed over multiple visits up to 1 year after study entry (Table 1). The Institutional Review Board of the Johns Hopkins University School of Medicine approved this study, and all participants signed written consent prior to initiation of any study activities.

TABLE 1

Table 1. Breakdown of patient counts by visit.

Clinical Data Collection

A trained interviewer administered a series of detailed questionnaires regarding demographics and general medical, medication, and symptom histories to both patients with early LD and controls at each study time point. In addition, specific data were gathered from patients with early LD regarding their acute illness. At the follow-up study visits, participants self-administered a 36-item symptom list developed based on prior clinical and research experience among patients with Lyme disease. The presence or absence of self-reported fatigue, musculoskeletal pain, and/or cognitive difficulty at V3 and V4 were used to classify participants into those with persistent symptoms and those without. Given the relatively small number of participants in this analysis, we did not require participants to also meet criteria for functional impact and therefore those with persistent symptoms did not meet the full PTLD case definition. During the first study visit, a sensitive, PCR-based approach (PCR/ESI-MS), as previously described in a subset of these participants (13), was used to detect the presence of B. burgdorferi in the skin (Supplementary Figure 1A) and blood (Supplementary Figure 1B).

PBMC Isolation and Library Generation

PBMCs were isolated from fresh whole blood using Ficoll (Ficoll-Paque Plus, GE Healthcare) and total RNA was extracted from 10⁷ PBMCs using RLT Lysis Buffer (Qiagen) by following manufacturer's instructions. The NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (Cat# E7765) was used to generate RNA-seq libraries. Briefly, Poly A RNAs were isolated from total RNAs using NEBNext Poly(A) Magnetic Isolation Module (NEB #E7490) and then fragmented for cDNA synthesis. End repair is performed where 3' to 5' exonuclease activity of enzymes removes 3' overhangs and the polymerase activity fills in the 5' overhangs. An “A” base is then added to the 3' end of the blunt phosphorylated DNA fragments which prepares the DNA fragments for ligation to the sequencing adapters, which have a single “T” base overhang at their 3' end. Ligated fragments are subsequently size-selected through purification using the Sample Purification Beads included in the kit and undergo PCR amplification to prepare the “libraries.” The BioAnalyzer is used for quality control of the libraries to ensure adequate concentration and appropriate fragment size free of adapter dimers. The resulting library insert size is 200–500 bp with a median size around 300 bp. Libraries were uniquely barcoded and pooled for HiSeq2500 sequencing.

RNA-Seq Data Processing

Raw RNA-seq FASTQ files were processed by FastQC, a quality control tool for high throughput sequencing data (14). The samples were aligned to the human genome (hg38) with the STAR RNA-seq aligner (15). Picard tools were then used for manipulating the output from STAR so it can be piped into featureCounts (16) for gene, exon, and transcript quantification. This pipeline was encoded in python and it is made available on GitHub at https://github.com/lymeMIND.

Cytokine/Chemokine Assays

The levels of 38 cytokines and other immune mediators were measured using Bio-Plex cytokine arrays and the Bio-Plex 200 System (Bio-Rad Laboratories). All tests at the participant level for each specific immune mediator were run by the same system. All tests were run as recommended by the manufacturer using previously described optimized assay protocols (17). The cytokines, chemokines, and acute phase markers measured were: Eotaxin, FGF basic, G-CSF, GM-CSF, HTARC, IFN-γ, IL-1β, IL-1rα, IL-2, IL-4, IL-5, IL-6, IL-7, IL-8, IL-9, IL-10, IL-12, IL-13, IL-15, IL-17, IL-17A, IL-17F, IL-21, IL-22, IL-23, IL-25, IL-31, IL-33, IP-10, MCP-1(MCAF), MIP-1α, MIP-1β, MIP-3β, PDGF-ββ, RANTES, sCD40L, TNF-α, and VEGF. Data processing was performed using Bio-Plex manager software version 4.4.1, and serum concentrations were interpolated from standard curves for each respective cytokine. These data were transformed using the log of the “ratio to average” for ease of interpretation. This was calculated by setting all values <1–1 pg/mL, and then calculating the log base 2 of [(value)/(average value in the cohort)]. As a result, 0 represents an average value, 1 represents a value which is 2 times the average, and 2 represents a value which is 4 times the average.

Data Analysis and Visualization

The RNA-seq gene counts were log2 transformed, z-scored, and quantile normalized (18) across all patients irrespective of time-point. These features were used for the subsequent data visualizations and classifiers. For the differential expression analyses, RNA-seq gene counts were quantile normalized across all patients. Uniform Manifold Approximation and Projection (UMAP) (19) was used to perform non-linear manifold aware dimensionality reduction and is used for visualizing the features for all samples in two dimensions. Samples were clustered using the k-means clustering algorithm on the UMAP manifold. Silhouette clustering analysis (20) was applied to identify the optimal number of clusters. Comparisons by cluster were performed using Kruskall-Wallis tests for overall p-value, and Wilcoxon rank sum test for individual pairwise comparisons for continuous variables, and chi-square test for categorical variables.

The log2 transformed RNA-seq gene count features were selected in several ways and used to train and test a series of Random Forest (21) and Logistic Regression classifiers. We employed the same approach to distinguish cases from controls, as well as Lyme patients with persistent symptoms from those without at later time points Samples were randomly grouped into a stratified training and test set using 75 and 25% of the samples, respectively, while preserving the ratios of positive and negative labels in each group. The Random Forest classifier was trained using all features to distinguish patients. First, the training was done in such a way that would lead to a fully grown and unpruned set of decision trees, and then in a way that constrains the tree to utilize the most informative top 50 features across all trees. The Logistic Regression classifier was trained using all the features to distinguish patients. The top 50 features with the highest absolute slopes after training were then used to train another Logistic Regression classifier with access to only those features. Another Logistic Regression classifier was trained using the top 50 features with the highest ANOVA F-test. We computed the unit mean and standard deviations from the log2 normalized training data and applied z-score normalization to the log2 normalized training and test data. Performance of each of the models on the test sets were collected and visualized using AUROC and PR curves.

CIBERSORT was used along with the benchmarked LM22 signature file, which contains signatures for 22 immune cell types to derive cell type ratios from our patient's RNA-seq gene counts (22). The resulting cell type ratios were investigated for statistically significant feature correlations and the features themselves were used with the Random Forest classifiers and the Logistic Regression classifiers to distinguish Lyme patients from healthy controls, and Lyme patients with persistent symptoms from those without. Our RNA-seq counts were further contrasted with the RNA-seq counts from a study of the immune responses in several RNA-seq analyzed PBMCs of patients with COVID-19 (GSE152418). The data were merged by independently log2-transformed, z-score normalized, and quantile normalized before being combined into a unified set of gene expression and corrected for batch effects which distinguished the datasets using ComBat (23).

Analyses and figures were produced using python scikit-learn and the scipy ecosystem, SAS Software (version 9.4; SAS Institute Inc., Cary, NC, USA), and GraphPad Prism (version 8.1.0; GraphPad Software, San Diego, CA, USA).

Results

RNA-Seq Processing, Clustering, and Visualization

Seventy-two cases (41 males and 31 females) diagnosed with early localized and early disseminated LD along with 44 controls (19 males and 25 females) were enrolled in the longitudinal study. After aligning the RNA-seq data, transcript counts were converted to the gene level and then counts were log transformed and normalized (see methods). We then applied Uniform Manifold Approximation and Projection (UMAP) (19) to estimate similarities and differences between the samples, each representing a patient at a specific visit. Interestingly, cases were mostly separated from the controls, even after 6 months and 1 year follow up visits (Figure 1A). Such separation cannot be attributed to batch, gender, or seasonal effects (Supplementary Figure 2). Next, we aimed to automatically identify clusters based on the samples' RNA-seq data to see if the cases and controls further cluster into subtypes. The Silhouette score identified an optimal number of 3 clusters (Figure 1B). The automatic clustering placed the controls in cluster 2. LD cases are divided into two distinct clusters, clusters 0 and 1, with a smaller group of cases clustering with the controls in cluster 2 (Figure 1C). There are two notable observations. First, most cases and controls remain in one cluster over the course of study (Figures 1D–G). Secondly, some cases that are mainly in cluster 0 move around to cluster 2 (mostly controls) (Figure 1F), but there are no patients from cluster 1 that return to the control cluster 2 even after 1 year (Figures 1G–H). This is surprising because most LD patients that receive timely treatment are expected to return to a normal state after several weeks. We observe that none of the clinical measured variables can clearly explain the overall long-lasting immune activation in the LD cases. However, it should be noted that the data collected for this study was processed in batches over several years (Figures 1I–J). If we plot the data by batch or by year, we see clear patterns that separate the controls from most of the cases (Supplementary Figure 2). While it is possible that there are some batch effects within the data, such separation cannot explain the clear long-term immune activation observed for the cases. We also attempted to remove these batch effects with the ComBat method (23) and the results regarding long-term immune activation remain.

FIGURE 1

Figure 1. (A) UMAP projections of all samples collected for the study based on RNA-seq expression where each sample, representing a patient at a specific visit, is colored by cases (blue) vs. controls (red). (B) K-means is fit to the quantile normalized gene counts matrix and silhouette analysis is performed to identify an optimal number of clusters. (C) UMAP projections of all samples collected for the study based on RNA-seq expression where each sample, representing a patient at a specific visit, is colored by automatic cluster assignment. (D) UMAP projections of all samples collected for the study based on RNA-seq expression where each sample, representing a patient at a specific visit, is colored by cases (blue) vs. controls (red) and line trace the trajectory of patients over time. (E) UMAP projections of controls that change clusters (in color), cluster membership (color) and line trace of patients over time. (F) UMAP projections of cases that change clusters from cluster 0 (in color), cluster membership (color) and line trace of patients over time. (G) UMAP projections of cases that change clusters from cluster 1 (in color), cluster membership (color) and line trace of patients over time. (H) Membership of cases (left) and controls (right) for each visit in each automatically detected cluster. (I) UMAP projections of all samples colored by the season/year when the data was collected. (J) UMAP projections of all samples colored by one of six batches.

Lyme Serology Projection and Cytokine Profiling

Control participants generally remained in cluster 2 over time, with only 9 of the 132 (6.8%) total control visits found in clusters 0 or 1. We compared those control study visits in cluster 2 with those in clusters 0 or 1 and did not find statistically significant differences in Lyme serology (ELISA values, or the number of reactive IgM or IgG bands), or the percent reporting a recent viral infection such as a cold or flu in the past 10 days (p > 0.12 for each comparison). In addition, there were no new diagnoses or tick bites in the preceding interval reported by controls at any of the non-cluster 2 study visits. Direct evidence of infection was obtained using a sensitive PCR/ESI-MS approach that detected Borrelia burgdorferi in the skin or blood of most of our patients. Interestingly, there was a subgroup of cases that were PCR/ESI-MS negative and who failed to seroconvert on standard two-tier antibody tests [Figure 2, as previously identified in (13)]. Remarkably, these patients generally cluster within the other cases, suggesting that these patients may have LD, but current assays are not capable of detecting the presence of Borrelia in their skin or blood. It should be noted that two patients were clinically diagnosed with LD and considered cases but have negative detection of Borrelia antigens and clustered with the controls. These patients may represent a Lyme look-alike group. Thirty-eight cytokines and chemokines were measured in the sera of cases and controls and 5 were statistically significant by cluster group (p < 0.05, see Figure 3), while 7 were borderline statistically significant (0.10 < p < 0.05: IL-23, IL-7, IL-13, IL-15, PDGF-ββ, IL-21, and IL-17F). The remaining n = 26 was not statistically significant by cluster group. Individual pairwise comparisons were also conducted between clusters.

FIGURE 2

Figure 2. UMAP projection of all samples collected for the study based on RNA-seq expression where each sample, representing a patient at a specific visit, is colored by detection status of Borrelia antigens in their blood or skin. Samples colored in gray did not have data on PCR/ESI-MS status either because they represent control samples, or because they represent case follow-up time points where these data were not obtained.

FIGURE 3

Figure 3. Statistically significant immune mediator differences among transcriptome-defined clusters. Levels of the 5 immune mediators found to be statistically significant by cluster are shown.

Differential Gene Expression Analysis

Differential expression analysis was performed using the Characteristic Direction (CD) method (24) on all quantile normalized case and control samples. Interestingly, three of the top-ranked up-regulated genes in LD cases are the human leukocyte antigens HLA-A and HLA-B, which suggest inflammation via adaptive immune activation (Figure 4). The top 200 up-regulated and down-regulated genes were subjected to enrichment analysis with Enrichr (25) to identify biological processes and pathways that are altered in the LD cases. As expected, enrichment terms associated with the immune response are associated with the top 200 up-regulated genes in LD cases (Supplementary Figure 3). Specifically, NFkB/RelA is consistently detected as the most enriched transcription factor based on ENCODE (26), TRRUST (27), TRANSFAC (28) and JASPAR (29), and Transcription Factor protein-protein interactions (PPIs). The most enriched pathways are those associated with Influenza, Salmonella, and Ebola infections. Interestingly, arthritis is the top enriched term using the Jensen DISEASES (30) library. Arthritis is a known symptom for LD patients (31). LINCS (32) L1000 (33) ligand perturbations up-regulated genes are enriched for signatures for IL-1 and TNF alpha, supporting immune system activation; and dbGAP (34) enrichment analysis points to leprosy and gout as the top terms. The mycobacterium tuberculosis is the most enriched term for Microbe Perturbations from GEO up library, suggesting that Borrelia effects on gene expression might be most like the effects of this pathogen.

FIGURE 4

Figure 4. Top ranked up-regulated and down-regulated differentially expressed genes in LD cases based on the characteristic direction method.

The most striking enrichment result is the top enriched term returned from the “Rare Diseases AutoRIF ARCHS4 Prediction” gene set library which is Erythema elevatum diutinum. This library contains 3,725 gene sets created using the Geneshot (35) tool to search the names of various rare diseases on PubMed. Geneshot uses AutoRIF, which is a resource containing PubMed IDs and the genes mentioned in the title and abstract of these publications. For each rare disease term, the associated PubMed IDs are cross referenced to the AutoRIF resource and ranked by publication frequency. Geneshot then converts these generated gene lists into predicted gene lists using co-expression data from ARCHS4 (36) to identify genes that may be associated with disease terms but not yet studied or published in the literature. All 44 genes associated with Erythema elevatum diutinum were contained within the top 200 positive top ranked genes from the differential expression analysis. These results suggest that erythema disease genes are associated with an immune module that is highly relevant to the molecular mechanisms of LD.

We also examined the differentially expressed genes (DEGs) that were driving the separation of the LD cases into two clusters. Not surprisingly, each of these two clusters of cases displayed unique and defining sets of DEGs (Supplementary Figure 4). The top 200 up-regulated and down-regulated genes in each cluster was subjected to enrichment analysis to identify biological processes and pathways that are unique to each cluster. Surprisingly, HLA-E, a class Ib molecule that can serve as regulatory ligand for NK cells was prominently down-regulated in cluster 0 when compared to cluster 1. These and other enriched terms, derived from the up-regulated genes in clusters 0 vs. 1, and account for the majority of Lyme cases, are visualized as bar charts (Supplementary Figure 5). These two clusters share features that include T cell receptor signaling, the involvement of monocytes and CD4+ T cells, enrichment of Rel and NFKB transcription factors and association with arthritis and leprosy. Cluster 0 is distinct in that monocytes, dendritic cells, and CD8+ T cell involvement are implied. Distinct features of cluster 1 include neutrophils, IL-4 signaling and IgE/IgA synthesis. As above, the most remarkable enrichment came with results from the “Rare Diseases AutoRIF ARCHS4 Prediction” gene set library with each cluster associated with distinct diseases. The rare diseases associated with cluster 0 include those with autoimmune or inflammatory features while cluster 1 rare diseases are largely driven by non-immune mechanisms. This highlights the differences between patients within these two clusters. Such differences can potentially be explained by different types of cytokine CD4+ helper T cells (Th) response.

Cell Type Deconvolution With CIBERSORT

To further explore and extract information from the RNA-seq data, we next applied deconvolution algorithms to identify potential changes in cell type composition between the controls and cases. Specifically, we applied the CIBERSORT algorithm with the LM22 reference signatures (22) which identifies the proportion of 22 immune cell types in a sample (Figure 5A). We observe that the Lyme disease cases have significantly more Tregs (Figure 5B), more monocytes (Figure 5C), and less resting T memory cells (Figure 5D). These observations suggest a general activation of an immune response consistent with the conclusion made via the differential expression analysis. Interestingly, the increase in Tregs compared with the controls is more pronounced in cluster 1 compared with cluster 0 across all visits (Figure 5E). This further supports the overall potentially more robust immune response for the patients within cluster 1. Tregs are known to increase with infection and inflammation. They suppress auto-immune host tissue damage and assist establishing tolerance when the inflammation resolves. Next, we asked whether there are gene modules that return to normal level over time and approach the control subjects, and conversely others that do not. For this, we created two Venn diagrams that highlight the unique genes that are up/down early and those that are up/down late when comparing the cases to the controls. We observe that many cell cycle genes are uniquely up regulated at early visits and down regulated at late visits (Figure 5F). This is consistent with the known rapid proliferation of T cells that occurs at the initiation of innate immune response, and apoptosis and other regulatory mechanisms during adaptive immunity stages.

FIGURE 5

Figure 5. (A) Cell-type ratios were derived from gene expression with CIBERSORT. Green bar denotes Lyme disease patients. (B–D) Cell-types ratios with significant differences across all visits in healthy patients (type_subject = 0) vs. those with Lyme disease (type_subject = 1) Significance computed with Welch's unequal variance t-test for independence. (E) Predicted levels of Tregs at each visit for each cluster based on CIBERSORT LM22. (F) Comparing the identified DEGs between cases and control and early and late time-points. Enrichment analysis with Enrichr was applied to the unique genes that are only up early (left), or only down (dn) late (right) with the gene ontology biological processes gene set library.

Machine Learning Applications to Classify Patients

Finally, we applied machine learning methods to attempt to diagnose patients based on their RNA-seq gene expression profiles. Random Forest and Logistic Regression classifiers were trained using stratified training and testing data accounting for 75 and 25% of the RNA-seq from Lyme disease and control samples, respectively. These features are RNA-seq gene expression counts that were log2-normalized, variance-filtered, z-score and quantile normalized. Overall, five different classifiers were constructed: (1) a classifier for predicting whether a patient has Lyme disease, or is an healthy control (Figure 6A), (2) a classifier for predicting whether a patient that has Lyme disease will progress to having persistent symptoms (Figure 6B), (3) and (4) classifiers that attempt to perform the same classifications as 1 and 2, but with CIBERSORT transformed features instead of using the gene expression features directly (Figures 6C,D); and finally, (5) a classifier for predicting whether a patient has Lyme disease or COVID-19 (Figure 6E). The 17 COVID-19 RNA-seq samples were taken from a recent study that profiled patients' PBMCs (GSE152418). We evaluate the performance of each of the trained classifiers instance by computing AUROC and Precision-Recall curves (Figure 6). After five passes of this procedure, we aggregate these metrics to provide a sense of performance and generalizability of the classifiers tested. Overall, we observe that the Logistic Regression classifiers outperform the Random Forest classifiers when using all genes or when restricting the classifiers to use the top 50 most informative genes. The performance of these classifiers suggests that they could potentially aid in diagnosis but may not be good enough to replace existing methods. Additional information about the population of Lyme disease patients would be necessary to adequately assess the performance of this approach in a clinical setting. In addition, more samples and uniform processing of the data will be required. We should also note that the comparison to the COVID-19 patients has several caveats. For one, we are comparing acute viral infection to persistent bacterial infection. The data from the COVID-19 study and our LD study were applied to different populations. These studies were conducted by different groups that employed different protocols to collect and process the data.

FIGURE 6

Figure 6. (A) AUROC and PR curves for classifiers predicting whether a gene expression sample is from a Lyme disease diagnosed patient or a healthy control across all visits. (B) AUROC And PR curves for classifiers predicting whether a gene expression sample is from a Lyme disease diagnosed patient with persistent symptoms. (C) AUROC and PR curves for classifiers predicting whether a gene expression sample is from a Lyme disease diagnosed patient or a healthy control across all visits using CIBERSORT features instead of RNA-seq gene expression. (D) AUROC and PR curves for classifiers predicting whether a gene expression sample is from a Lyme disease diagnosed patient with persistent symptoms using CIBERSORT features instead of RNA-seq gene expression. (E) AUROC and PR curves for classifiers predicting whether a patient's observed gene expression, normalized and batch effect corrected across experiments, suggests they have Lyme disease or COVID-19.

Web supplement

To make the data, analysis, and results from this study accessible and reusable, we developed two web supplement components, one is a Jupyter Notebook with the code, markdown text, and figures generated for this study. The other is an interactive plot that provides enrichment analysis for the different clusters. These two web-based components can be accessed from: https://commons.lymemind.org/#/Notebook and https://commons.lymemind.org/#/Viewer.

Discussion

The molecular mechanisms underlying the disease course and outcomes of LD are still poorly understood. In this study, we applied in-depth examination of RNA-seq profiling of 73 LD patients over a period of 1 year and compared these to 44 non-LD infected controls over a period of 2 years. The gene expression analysis clustered controls and LD cases into three distinct clusters, each with distinct clinical and immunological features. The majority of cases remained within a single cluster even up to 1 year after diagnosis. Enrichment analysis of the differentially expressed genes between the cases and controls identified up-regulation of immune response genes as well as genes specific for diseases displaying erythema. The clear separation between controls and cases enabled us to develop expression-based machine learning classifiers that may inform improved diagnosis. In fact, we developed three different types of classifiers. One to discriminate between LD cases and healthy controls, one to distinguish between LD cases and those cases that will progress to having persistent symptoms, and one that to classify LD from COVID-19 patients. Overall, the results from these classifiers are encouraging. The cost of RNA-seq is continually dropping so it could be used as a practical approach. While it is difficult to tell whether such a classifier will be able to discriminate from other similar bacterial infections, the up-regulation of immune response genes observed in the cases can be readily detected via expression profiling of their PBMCs.

Two other prior studies attempted to use genome-wide gene expression data to profile LD patients (37, 38). Both studies profiled gene expression from PBMCs extracted from LD patients and controls. Both studies had less patients compared with the number of patients profiled in our study. The size of the cohorts is critical to gain statistical insights and identify patterns in the data collected from these patients. The first study (37), analyzed acute LD patients using RNA-seq profiling and the conclusions of this study are similar to what we observed, which is an increase in immune response genes in the cases and continued dysregulation of gene expression months post infection and treatment. However, it was unknown whether patients return to normal after several months from their initial diagnosis.

The second more recent study (38), was a longitudinal study using cDNA microarrays to profile gene expression and, unlike ours and other findings (37), found that the DRGs identified during the acute phase returned to normal levels at 6 months. There are several potential reasons for this difference. First, this study was conducted using a smaller number of LD cases, with 10 analyzed at the 6-month time point. Second, the study exclusively used cases that displayed disseminated lesions, while our study included all LD patients that presented with a 5 cm or larger EM rash. Interestingly, our study includes 24 cases displaying disseminated lesions and we did not observe a relationship between lesion features and transcriptional clustering (Table 2). We suggest that the differences in outcomes between studies may be due to variation in clinical case definition, sample size and potentially differences in Borrelia burgdorferi sensu stricto strains that initiate disease in differing geographic regions.

TABLE 2

Table 2. Baseline demographic and clinical characteristics of 72 participants with early Lyme disease by v1 cluster group^a.

In our study, dimensionality reduction followed by clustering analysis of the longitudinal RNA-seq data identified three distinct clusters. In one cluster, containing mostly cases, the patients can be observed mixing with another cluster, containing mostly controls, over multiple visits. However, patients within another cluster containing mostly cases, stay in the same region of expression space even at later visits after treatment. This pattern of remaining within a cluster over time suggests that there is a long-term alteration of genes targeted for transcription in PBMCs following acute LD diagnosis. However, it is still unclear what is driving this persistent transcriptional activity. Possibilities include the continued presence of spirochetes or foreign bacterial antigen as well as epigenetic changes, all of which could drive long term transcriptional alterations. It is interesting to point out that for a subset of 21 cases, Borrelia DNA was not detected in the skin or blood using extremely sensitive PCE/ESI-MS based approach. Initially, we categorized these cases as potential cases of STARI, or sampling or diagnostic errors. However, in this current study, most of these patients have transcriptional profiles that cluster them with other cases, suggesting that these may represent true LD cases for which the applied detection methodologies failed. This may also suggest that profiling host gene expression may be a useful and potentially a sensitive approach for the diagnosis of acute LD.

The three clusters defined by the DEGs also have distinct clinical and immunological features. Cluster 0, the largest cluster, was found to have a higher number of symptoms of acute disease, higher rates of two-tier antibody positivity, as well as a non-statistically significant trend toward higher rates of abnormal liver function tests. All together, these imply that this cluster has more severe acute disease and a robust immune response. This is consistent with the enrichment analysis, which associates cluster 0 with immune mediated common and rare diseases. Also revealed by the enrichment analysis, was an association with CD8+ T cells, a subset usually linked to the control of intracellular pathogens via cytotoxic effector function. Interestingly, CD8+ T cells have been implicated in LD (40, 41). Cluster 0 also uniquely displays elevated levels of IL-31, member of the gp130/IL-6 family of cytokines (42). This cytokine is thought to be produced by CD4+ Th2 cells but Th1 cells can be induced to express IL-31 (43). IL-31 is also known to promote skin inflammatory disorders in mouse models and in the human setting (42, 44). Interestingly, in a retrospective study, the development of systemic autoimmune joint diseases, including psoriatic arthritis was associated with prior LD (45). Also, in our study, differential gene expression analysis comparing the controls to the cases identify immune response genes that may be specific for LD. Enriched terms include LD common symptoms such as erythema and arthritis. The overlapping genes that lead to the significant overlap between the differentially expressed genes and these terms may directly suggest novel molecular mechanisms of disease. Cluster 1 is likewise distinct as these cases have milder disease. In addition, enrichment analysis associates this cluster with distinct cells, pathway and processes including neutrophils, IL-4 signaling, IgA/IgE production as well a novel and rare diseases which are thought not to involve immune mediated processes. Neutrophils are thought not to be a major effector cell in LD, although Osp A can activate human neutrophils (46). The notion that infection with B. burgdorferi may trigger such non-immune pathways is an intriguing idea that may help explain the range of disease outcomes linked to LD. Previous work found that the human immune response in LD is largely driven by monocytes and Th1 driven mechanisms (47, 48). However, these studies largely focused on T cells recovered from inflamed joints in Lyme arthritis, a late feature of untreated LD. Our studies indicate that a wide range of immune effector cells may be engaged in LD including NK cells, CD8+ T cells and neutrophils, which may help define disease subgroups and inform treatment strategies. Clearly, future studies involving detailed immune profiling of multiple LD cohorts are required.

Analysis of our gene expression data by CIBERSORT revealed an enriched signal for Treg cells in Lyme cases. T regulatory cells have well established roles in the regulation of self-immunity as well as influencing the outcomes of infection (49). In human Lyme disease, low levels of Tregs in the synovial fluid were associated with a longer duration of illness in antibiotic-refractory Lyme arthritis, a late manifestation of human borrelliosis (50, 51). This is also consistent with the finding that in vivo depletion of Tregs accelerated inflammatory arthritis in a mouse model (52), suggesting that Tregs in the context of Lyme arthritis exhibit anti-inflammatory properties. The significance of the elevated Treg signal in our study of acute Lyme disease is not clear. Several studies have shown that in human tuberculosis increased Tregs levels are associated with active disease, suggesting that in some cases Treg responses may impair pathogen clearance and disease resolution (53, 54). Whether this is the case in some stages of human Lyme disease, or if these cells contribute to the observed dysregulated or maladaptive immune response in Lyme disease (55, 56), will require additional studies on the identification and characterization of T regs in various stages and outcomes of Lyme disease. This can be potentially achieved via single cell RNA-sequencing.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/LymeMIND/LM2-study-supporting-materials/tree/main/data.

Ethics Statement

The studies involving human participants were reviewed and approved by JHU & ISMMS IRB. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

AR, JF, MWE, MJM, WR, NR, MS, and JA collected the data. DC, AR, AB, MW, MD, and AM performed the analyses. DC, AR, AB, SJ, JB, MS, JA, and AM wrote the paper. All authors contributed to the article and approved the submitted version.

Funding

The project was partially supported by funds from the Cohen Lyme & Tickborne Disease Initiative, NIH grants U54HL127624, U24CA224260, and P30AR070254.

Conflict of Interest

MWE and MJM were employees of Ibis Biosciences, an Abbott Company, which developed the PCR/ESI-MS assays and instrumentation used in these studies; assays described are for research use only.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would also like to thank Marcia Daniela Villega de Flores, Thelio Sewell, and Erica Mihm for their efforts in patient recruitment, consent, sample processing and archiving. We would like to thank Michal Caspi Tal for help with the CIBERSORT analysis. We would like to thank members of the Dudley Lab for performing the initial data analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.636289/full#supplementary-material

References

1. Bush LM, Vazquez-Pertejo MT. Tick borne illness—Lyme disease. Dis A Month. (2018) 64:195–212. doi: 10.1016/j.disamonth.2018.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rosenberg R, Lindsey NP, Fischer M, Gregory CJ, Hinckley AF, Mead PS, et al. Vital signs: trends in reported vectorborne disease cases—United States and Territories, 2004–2016. Morb Mortal Wkly Rep. (2018) 67:496. doi: 10.15585/mmwr.mm6717e1

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kuehn BM. CDC estimates 300 000 US cases of Lyme disease annually. JAMA. (2013) 310:1110. doi: 10.1001/jama.2013.278331

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Nelson CA, Saha S, Kugeler KJ, Delorey MJ, Shankar MB, Hinckley AF, et al. Incidence of clinician-diagnosed Lyme disease, United States, 2005–2010. Emerg Infect Dis. (2015) 21:1625. doi: 10.3201/eid2109.150417

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Adrion ER, Aucott J, Lemke KW, Weiner JP. Health care costs, utilization and patterns of care following Lyme disease. PLoS ONE. (2015) 10:e0116767. doi: 10.1371/journal.pone.0116767

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Sanchez JL. Clinical manifestations and treatment of Lyme disease. Clin Lab Med. (2015) 35:765–78. doi: 10.1016/j.cll.2015.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cardenas-de la Garza JA, De la Cruz-Valadez E, Ocampo-Candiani J, Welsh O. Clinical spectrum of Lyme disease. Eur J Clin Microbiol Infect Dis. (2019) 38:201–8. doi: 10.1007/s10096-018-3417-1

CrossRef Full Text | Google Scholar

8. Aucott JN, Rebman AW, Crowder LA, Kortte KB. Post-treatment Lyme disease syndrome symptomatology and the impact on life functioning: is there something here? Qual Life Res. (2013) 22:75–84. doi: 10.1007/s11136-012-0126-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Aucott JN. Posttreatment Lyme disease syndrome. Infect Dis Clin. (2015) 29:309–23. doi: 10.1016/j.idc.2015.02.012

CrossRef Full Text | Google Scholar

10. Schoen RT. Challenges in the diagnosis and treatment of Lyme disease. Curr Rheumatol Rep. (2020) 22:1–11. doi: 10.1007/s11926-019-0857-2

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Horn EJ, Dempsey G, Schotthoefer AM, Prisco UL, McArdle M, Gervasi SS, et al. The Lyme disease Biobank–characterization of 550 patient and control samples from the east coast and upper midwest of the United States. J Clin Microbiol. (2020) 58:e00032–20. doi: 10.1128/JCM.00032-20

PubMed Abstract | CrossRef Full Text | Google Scholar

12. CDC. Lyme Disease (Borrelia burgdorferi) 2017 Case Definition. National Notifiable Diseases Surveillance System (NNDSS) (2017). Available online at: https://wwwn.cdc.gov/nndss/conditions/lyme-disease/case-definition/2017/ (accessed July 1, 2020).

Google Scholar

13. Mosel MR, Rebman AW, Carolan HE, Montenegro T, Lovari R, Schutzer SE, et al. Molecular microbiological and immune characterization of a cohort of patients diagnosed with early Lyme disease. J Clin Microbiol. (2020) 59:e00615–20. doi: 10.1128/JCM.00615-20

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Andrews S. FastQC: A Quality Control Tool for High Throughput Sequence Data. Cambridge, United Kingdom: Babraham Bioinformatics, Babraham Institute (2010).

Google Scholar

15. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. (2014) 30:923–30. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Soloski MJ, Crowder LA, Lahey LJ, Wagner CA, Robinson WH, Aucott JN. Serum inflammatory mediators as markers of human Lyme disease activity. PLoS ONE. (2014) 9:e93243. doi: 10.1371/journal.pone.0093243

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bolstad BM, Irizarry RA, Åstrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. (2003) 19:185–93. doi: 10.1093/bioinformatics/19.2.185

PubMed Abstract | CrossRef Full Text | Google Scholar

19. McInnes L, Healy J, Saul N, Großberger L. UMAP: Uniform Manifold Approximation and Projection. J Open Source Softw. (2018) 3:861. doi: 10.21105/joss.00861

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Rousseeuw PJ. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J Comput Appl Math. (1987) 20:53–65. doi: 10.1016/0377-0427(87)90125-7

CrossRef Full Text | Google Scholar

21. Liaw A, Wiener M. Classification and regression by randomForest. R News. (2002) 2:18–22.

Google Scholar

22. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. (2007) 8:118–27. doi: 10.1093/biostatistics/kxj037

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Clark NR, Hu KS, Feldmann AS, Kou Y, Chen EY, Duan Q, et al. The characteristic direction: a geometrical approach to identify differentially expressed genes. BMC Bioinform. (2014) 15:79. doi: 10.1186/1471-2105-15-79

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. (2016) 44:W90–7. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

26. ENCODE Project Consortium. The ENCODE (ENCyclopedia of DNA elements) project. Science. (2004) 306:636–40. doi: 10.1126/science.1105136

CrossRef Full Text | Google Scholar

27. Han H, Shim H, Shin D, Shim JE, Ko Y, Shin J, et al. TRRUST: a reference database of human transcriptional regulatory interactions. Sci Rep. (2015) 5:11432. doi: 10.1038/srep11432

CrossRef Full Text | Google Scholar

28. Matys V, Fricke E, Geffers R, Gößling E, Haubrock M, Hehl R, et al. TRANSFAC®: transcriptional regulation, from patterns to profiles. Nucleic Acids Res. (2003) 31:374–8. doi: 10.1093/nar/gkg108

CrossRef Full Text | Google Scholar

29. Sandelin A, Alkema W, Engström P, Wasserman WW, Lenhard B. JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Nucleic Acids Res. (2004) 32:D91–4. doi: 10.1093/nar/gkh012

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pletscher-Frankild S, Pallejà A, Tsafou K, Binder JX, Jensen LJ. DISEASES: text mining and data integration of disease–gene associations. Methods. (2015) 74:83–9. doi: 10.1016/j.ymeth.2014.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Steere AC, Schoen RT, Taylor E. The clinical evolution of Lyme arthritis. Ann Int Med. (1987) 107:725–31. doi: 10.7326/0003-4819-107-5-725

CrossRef Full Text | Google Scholar

32. Keenan AB, Jenkins SL, Jagodnik KM, Koplev S, He E, Torre D, et al. The library of integrated network-based cellular signatures NIH program: system-level cataloging of human cells response to perturbations. Cell Syst. (2018) 6:13–24. doi: 10.1016/j.cels.2017.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, et al. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell. (2017) 171:1437–52.e17. doi: 10.1016/j.cell.2017.10.049

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Mailman MD, Feolo M, Jin Y, Kimura M, Tryka K, Bagoutdinov R, et al. The NCBI dbGaP database of genotypes and phenotypes. Nat Genet. (2007) 39:1181–6. doi: 10.1038/ng1007-1181

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Lachmann A, Schilder BM, Wojciechowicz ML, Torre D, Kuleshov MV, Keenan AB, et al. Geneshot: search engine for ranking genes from arbitrary text queries. Nucleic Acids Res. (2019) 47:W571–7. doi: 10.1093/nar/gkz393

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Lachmann A, Torre D, Keenan AB, Jagodnik KM, Lee HJ, Wang L, et al. Massive mining of publicly available RNA-seq data from human and mouse. Nat Commun. (2018) 9:1366. doi: 10.1038/s41467-018-03751-6

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bouquet J, Soloski MJ, Swei A, Cheadle C, Federman S, Billaud J-N, et al. Longitudinal transcriptome analysis reveals a sustained differential gene expression signature in patients treated for acute Lyme disease. mBio. (2016) 7:e00100–16. doi: 10.1128/mBio.00100-16

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Petzke MM, Volyanskyy K, Mao Y, Arevalo B, Zohn R, Quituisaca J, et al. global transcriptome analysis identifies a diagnostic signature for early disseminated Lyme disease and its resolution. mBio. (2020) 11:e00047–20. doi: 10.1128/mBio.00047-20

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Eshoo MW, Crowder CC, Rebman AW, Rounds MA, Matthews HE, Picuri JM, et al. Direct molecular detection and genotyping of Borrelia burgdorferi from whole blood of patients with early Lyme disease. PLoS ONE. (2012) 7:e36825. doi: 10.1371/journal.pone.0036825

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Busch DH, Jassoy C, Brinckmann U, Girschick H, Huppertz HI. Detection of Borrelia burgdorferi-specific CD8+ cytotoxic T cells in patients with Lyme arthritis. J Immunol. (1996) 157:3534–41.

PubMed Abstract | Google Scholar

41. Whiteside SK, Snook JP, Ma Y, Sonderegger FL, Fisher C, Petersen C, et al. IL-10 deficiency reveals a role for TLR2-dependent bystander activation of T cells in Lyme arthritis. J Immunol. (2018) 200:1457–70. doi: 10.4049/jimmunol.1701248

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Gibbs BF, Patsinakidis N, Raap U. Role of the pruritic cytokine IL-31 in autoimmune skin diseases. Front Immunol. (2019) 10:1383. doi: 10.3389/fimmu.2019.01383

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Stott B, Lavender P, Lehmann S, Pennino D, Durham S, Schmidt-Weber CB. Human IL-31 is induced by IL-4 and promotes TH₂-driven inflammation. J Allergy Clin Immunol. (2013) 132:446–54.e5. doi: 10.1016/j.jaci.2013.03.050

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Dillon SR, Sprecher C, Hammond A, Bilsborough J, Presnell SR, Haugen HS, et al. Interleukin 31, a cytokine produced by activated T cells, induces dermatitis in mice. Nat Immunol. (2004) 5:752–60. doi: 10.1038/ni1084

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Arvikar SL, Crowley JT, Sulka KB, Steere AC. Autoimmune Arthritides, Rheumatoid Arthritis, Psoriatic Arthritis, or Peripheral Spondyloarthritis Following Lyme Disease. Arthritis Rheumatol. (2017) 69:194–202. doi: 10.1002/art.39866

CrossRef Full Text | Google Scholar

46. Morrison TB, Weis JH, Weis JJ. Borrelia burgdorferi outer surface protein A (OspA) activates and primes human neutrophils. J Immunol. (1997) 158:4838–45.

PubMed Abstract | Google Scholar

47. Strle K, Shin JJ, Glickstein LJ, Steere AC. Association of a toll-like receptor 1 polymorphism with heightened Th1 inflammatory responses and antibiotic-refractory Lyme arthritis. Arthritis Rheum. (2012) 64:1497–507. doi: 10.1002/art.34383

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gross DM, Steere AC, Huber BT. T helper 1 response is dominant and localized to the synovial fluid in patients with Lyme arthritis. J Immunol. (1998) 160:1022–8.

PubMed Abstract | Google Scholar

49. Richert-Spuhler LE, Lund JM. The immune fulcrum: regulatory T cells tip the balance between pro- and anti-inflammatory outcomes upon infection. Prog Mol Biol Transl Sci. (2015) 136:217–43. doi: 10.1016/bs.pmbts.2015.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Vudattu NK, Strle K, Steere AC, Drouin EE. Dysregulation of CD4+CD25(high) T cells in the synovial fluid of patients with antibiotic-refractory Lyme arthritis. Arthritis Rheum. (2013) 65:1643–53. doi: 10.1002/art.37910

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Shen S, Shin JJ, Strle K, McHugh G, Li X, Glickstein LJ, et al. Treg cell numbers and function in patients with antibiotic-refractory or antibiotic-responsive Lyme arthritis. Arthritis Rheum. (2010) 62:2127–37. doi: 10.1002/art.27468

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Siebers EM, Liedhegner ES, Lawlor MW, Schell RF, Nardelli DT. Regulatory T cells contribute to resistance against Lyme arthritis. Infect Immun. (2020) 88:e00160–20. doi: 10.1128/IAI.00160-20

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Chen X, Zhou B, Li M, Deng Q, Wu X, Le X, et al. CD4(+)CD25(+)FoxP3(+) regulatory T cells suppress Mycobacterium tuberculosis immunity in patients with active disease. Clin Immunol. (2007) 123:50–9. doi: 10.1016/j.clim.2006.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Beiranvand E, Abediankenari S, Rezaei MS, Khani S, Sardari S, Beiranvand B. Increased expression of forkhead box protein 3 gene of regulatory T cells in patients with active tuberculosis. Inflamm Allergy Drug Targets. (2014) 13:330–4. doi: 10.2174/1871528113666140908112939

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Steere AC. Posttreatment Lyme disease syndromes: distinct pathogenesis caused by maladaptive host responses. J Clin Invest. (2020) 130:2148–51. doi: 10.1172/JCI138062

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Sellati TJ, Barberio DM. Mechanisms of dysregulated antibody response in Lyme disease. Front Cell Infect Microbiol. (2020) 10:567252. doi: 10.3389/fcimb.2020.567252

PubMed Abstract | CrossRef Full Text | Google Scholar