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.

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).

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.

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.

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.

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).

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.

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.

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.

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.

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.

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.

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.