Antibody binding of a total of 536 human serum samples representing seropositive (LD+) and seronegative (LD−) Lyme disease and 13 diseases with symptoms similar to LD (LADs; see Table 1 for cohort breakdown) were profiled on the peptide arrays. The “seropositive” designation was assigned to the LD patients who presented with a rash of >5-cm diameter and tested positive with the Centers for Disease Control and Prevention (CDC)-recommended serologic test. Lyme disease patients who presented with a rash of >5-cm diameter but tested negative with standard two-tier test (STTT) were categorized as “seronegative” LD.

The sequences of the peptides on the array were designed randomly with a goal of equally representing 19 (cysteine was excluded due to synthesis constraints) canonical amino acids on the array. As a result, the peptide arrays can be understood as an unbiased and agnostic way to evenly interrogate the binding of the patient’s antibody repertoire to the entire combinatorial space of 10-mer peptides. Due to the random and agnostic peptide array design, any differences observed in antibody binding between the study cohorts can be attributed to a specific humoral immune system response and should not be affected by the peptide array design per se.

The distribution of binding intensities of total serum IgG to peptide array sequences shows an overall stronger binding for samples from the patient cohort with LADs compared to the LD+ and LD− disease cohorts ( Figure 1 ). All three cohorts show bimodal distributions with two distinct peaks in their binding intensity. The first peak is centered at the low end of the binding intensity range close to the background binding signal on the array. Thus, the first peak likely represents the fraction of peptides that bind antibodies non-specifically and only weakly or not at all. The second peak in all three cohorts involves stronger binding peptide sequences. Despite the similar shape of the distributions between the three cohorts, the distributions show differences. The shift of the second peak with respect to the first varies markedly with cohort. The look-alike cohort shows the largest separation between the peaks with approximately 4× higher intensity of the second peak (0.6 on the log10 scale). The seronegative LD (LD−) group of patients exhibits the smallest ratio of ~1.6× in the shift between the two peaks, with the seropositive (LD+) cohort showing an intermediate shift of ~2×. The position of the second peak presumably captures the more specific binding and thus likely contains the disease-specific response. This suggests that antibodies in the look-alike disease group show an overall stronger immune response compared to both LD cohorts. Consistent with this, the LD+ cohort where clear antibody reactivity has been measured results in overall stronger binding than the LD− cohort, presumably reflecting the larger number of reactive antibodies present ( Figure 1A ).

The binding intensities between the three cohorts were qualitatively compared using the Uniform Manifold Approximation and Projection (UMAP) method for data dimensionality reduction and visualization ( Figure 1B ). The three disease cohorts show substantial overlap with one another in this representation (also see Supplementary Figure 1 for pairwise comparison using UMAP; the look-alike and LD− cohorts show partial separation in a binary comparison; Supplementary Figure 1B ). This suggests that despite the observed differences in the binding intensity distributions between the cohorts, differentiation based on the binding intensities of the individual peptide is not as clear. Interestingly, the UMAP representation suggests the presence of at least three subclusters in the data ( Figure 1B ), indicating underlying additional complexity with differences in fractional distributions among the cohorts. For example, cluster 3 contained only 27 out of 235 (11%) of the LAD patients, with the rest of the cohort distributed approximately evenly between clusters 1 and 2. In contrast, 39 out of 93 (42%) LD− individuals were included in cluster 3, while 17 (18%) and 37 (40%) were contained in clusters 1 and 2, respectively. The LD+ patients are distributed approximately equally among the three clusters.

Because of the differences observed between the seropositive and seronegative LD groups compared with the look-alike diseases, the combined LD (LD+/LD−) cohorts were compared with the look-alike disease cohort in determining p-value distributions and building classifier models. As observed earlier when analyzing the binding intensity distributions ( Figure 1 ), the majority of the peptides showed lower binding intensities in the LD cohorts compared with the look-alike diseases regardless of whether the LD cohorts were combined or not ( Figures 2A–C ). A comparison of classification performance using the Extreme Gradient Boosting (XGBoost) algorithm ( Figures 2D–F ) showed higher Area under the curve (AUC) values for both LD+ vs. LADs and LD− vs. LADs (AUC = 0.83, 95% CI: 0.77–0.98, and AUC = 0.85, 95% CI: 0.81–0.89, respectively) compared to the combined (LD+/LD−) vs. LADs [AUC = 0.77 (0.72–0.82)]. This result indicates a somewhat better classification in terms of AUC value when the two LD cohorts are considered separately. The classification results imply that the LD+ and LD− cohorts differ not only in terms of known biomarkers being positive in LD+ but also in terms of antibody reactivity that is unique to LD−.

Due to the random nature of the peptide sequences in the array library, one cannot directly use the information about the sequence-binding relationship to determine what proteins from the B. burgdorferi proteome may be immunogenic and serve as potential candidate biomarkers. To select proteins from the B. burgdorferi proteome for further validation, ML approaches were used to model the underlying sequence-binding patterns and then project the patterns onto all B. burgdorferi proteins. The sequences of each protein from the B. burgdorferi B31 strain proteome (n = 1,219) were split into 10 AA-long tiles with nine amino acids (AA) overlapping between adjacent tiles. The proteome tiles were then one-hot encoded and used as input for the neural network (NN) models trained on the peptide array binding data of each sample, resulting in predicted binding intensities of every tile. These were assembled into a predicted binding map of each protein. Note that one separate model was trained on the binding data from each patient, resulting in B. burgdorferi binding predictions generated for each patient. A comparison of the predicted binding distributions ( Figure 3A ) shows similar characteristics to the measured binding on the peptide array. All three cohorts show distinct bimodal shapes with the first peak representing weak binders and the second peak capturing mainly the stronger interactions. With respect to the LD+ and LD− cohorts, the LAD group shows the largest separation between the two peaks. It is also shifted most toward the higher intensities (stronger interactions) compared with the other two cohorts, followed by the LD+ and LD− groups. However, there is one substantial difference between the measured and predicted distributions. The second peak (strong interactions) in the distribution of predicted values to the B. burgdorferi proteome is higher than the first peak ( Figure 3A ), in contrast to that observed in the measured array peptide binding values ( Figure 1A ). This suggests that the sequences from the B. burgdorferi proteome overall contain more peptides (tiles) that resemble antigenic targets of antibodies in each patient. Interestingly, a UMAP representation of the predictions ( Figure 3B ) revealed a distribution with fewer distinct subclusters as compared to the measured intensities ( Figure 1B ) even though the shape and overall overlap between the distributions of the three cohorts are similar. In addition, the LAD and LD− cohorts (green and blue dots in Figure 1B ) appear to be separated more than observed in the measured data ( Figure 1B ), at least visually. Nevertheless, the overall similarity between the measured and predicted distributions indicates reliable model performance with regard to projecting patterns in the learned data onto a biologically relevant context.

The model predictions of antibody binding to the B. burgdorferi proteome were analyzed in terms of p-value distributions ( Figures 4A–C ). Similar distribution characteristics were observed to those of the data measured on the peptide arrays ( Figure 2 ), with the majority of the proteome tiles showing lower predicted binding intensities to antibodies in sera of the LD+ and LD− cohorts compared to the LAD group of patients. Two arbitrary threshold ratio values of 0.8 and 1.2 were used to better highlight trends in the differential binding profiles. There are relatively few tiles that show stronger binding in the LD cohort than in the LAD cohort, consistent with the binding distributions for the cohorts ( Figure 3A ). Classifier models were trained on the predicted B. burgdorferi binding values ( Figures 4D–F ). Lower classification performance was observed compared to models trained on the measured peptide array data ( Figures 2D–F ). All three classifiers showed similar AUC values, suggesting, in contrast to the results obtained with the peptide array data, comparable differentiation between the two LD and the LAD cohorts. Unlike the models trained on the measured peptide array data, the classification of the combined LD+/LD− vs. LAD cohorts was similar to the classification using the separate cohorts.

Note that the analysis of predicted total IgG binding to the B. burgdorferi proteome, as described above, was performed at the level of separate sequence tiles and not complete proteins. This resulted in generally low effect sizes and mediocre classification performance, suggesting that there are no strong candidate biomarkers at the individual sequence tile level.

It is not very surprising that short sequence pieces of potential antigens do not provide strong differentiation. Past work has shown that the immune response to Lyme disease is very heterogeneous (36, 44), and thus, even if patients have antibodies against common proteins, they may well not bind to the same epitopes in those proteins. An alternative is to consider the aggregate predicted binding to each protein in the proteome and in this way select candidate protein biomarkers from the B. burgdorferi proteome with high predicted differentiating power between the combined LD+ and LD− cohorts and the LAD group of patients. Two different complementary methods were used for candidate protein selection. The first method was based on ranking the proteins by the lowest p-value determined from the predicted binding values of the tiles that make up its sequence. The p-values were calculated using the outlier sum statistics (45) (see Materials and Methods) selecting the proteins with a false discovery rate (FDR) of <0.05. The outlier sum statistics method was chosen to account for the long tails in the binding intensity distributions observed on the peptide array assays. It has been demonstrated that outlier sum statistics outperforms the t-test method for calculating the statistical significance of distributions with outliers (45). This method of selection resulted in a total of 19 protein candidates being selected from a total of 1,281 proteins (Supplementary Data Sheet 1) contained in the B. burgdorferi proteome ( Table 2 ) (Equations 1–4). The last two proteins in Table 2 were included because the FDR values were above the 0.05 cutoff by only a small margin.

The second method for candidate protein selection utilizes the XGBoost classifier trained on the predicted binding intensities comparing the combined LD cohorts vs. the LAD cohort ( Figure 4A ). The XGBoost algorithm is based on a decision tree structure and intrinsically performs feature selection during training. As a result, a trained classifier is based on a subset of the features (protein tiles, in this case) that contribute to the classification of the two cohorts. To identify protein candidates, first, only the tiles that were selected by the algorithm in all the cross-validation rounds (n = 10) were kept. Second, the B. burgdorferi proteins were then ranked by the number of different tiles they contained that met this criterion (Supplementary Data Sheet 2). Proteins that contained at least three tiles from the list above were selected as candidate biomarkers, resulting in a total of 34 proteins ( Table 3 ). The arbitrary threshold for the number of tiles per protein was set to 3 to both limit the number of candidate proteins and increase the likelihood of a protein showing true differentiating power.

Interestingly, neither of the lists contains any of the serologic biomarkers used in the current LD testing standard. This indicates that the differential humoral immune response between the LD and LAD cohorts may have a different set of target antigens than when comparing LD with healthy controls (36). The two lists show no overlap and represent two unique sets of proteins. The list produced by method 2 does contain several proteins [including the two top-ranked proteins in Table 3 : the outer membrane protein (O51735) and an uncharacterized protein (O51465)] that are either known to be located on the membrane of the bacterium or are predicted extracellular (secreted) proteins. The cellular location makes these proteins accessible to antibody binding and provides further support for biological inference of the protein selection method. The finding that neither of the two lists contains any of the known LD biomarkers, e.g., the VlsE protein, is likely due to the difference in cohorts being compared. The standard serologic biomarker panel was developed to differentiate between LD and healthy controls, whereas here, LD was compared with LADs. Furthermore, the finding that the LD+ and LD− groups of patients were characterized by distinct antibody reactivity profiles, as demonstrated previously (36) and in this study ( Figure 2 ), means that combining the two cohorts as conducted has the potential to result in a set of biomarkers that better differentiates the combined LD cohort from LADs.

A total of 53 candidate proteins from the B. burgdorferi proteome were selected for further consideration. However, only 44 of these were successfully synthesized; the nine remaining proteins were excluded from synthesis due to substantial transmembrane regions. In addition, the VlsE protein, a biomarker that is currently being used for standard serology testing for LD, was included as a positive control for the LD+ samples. For validation assays, the proteins were attached to carboxylated paramagnetic beads using N-hydroxysulfosuccinimide sodium salt (NHS) chemistry (Materials and Methods). The beads were incubated with serum samples, and the antibody binding was measured as fluorescence intensity using a secondary polyclonal anti-IgG antibody labeled with phycoerythrin. A total of 185 LD+, 102 LD−, and 236 LAD samples were used for validation with the bead-based assays ( Supplementary Table 1 ). All samples were assayed in duplicate, and the mean values of the duplicates were used for further analysis. For data analysis, the binding values were converted to a log10 scale ( Supplementary Figure 3A ). The background binding signal was determined as the fluorescence intensity of the protein with the lowest coefficient of variation (CV) across all three cohorts. Low variation in the binding intensity of such a protein indicates that its binding is not disease-specific and can be used as a reference. While “blank” beads, i.e., beads not conjugated to any of the proteins, were also included in the assay as negative controls, it was observed that they showed some differential binding between the cohorts. It is possible that some of the antibodies in the three groups of patients exhibit preferential binding to the carboxyl moiety on the negative control beads. Further analyses, including classifier training, were performed with log-transformed intensity values that were used directly without further normalization. For classifier training, the LD+ and LD− cohorts were separated and contrasted against the LAD cohort individually. This was conducted based on the earlier findings in this study that suggested a different humoral immune response profile in the two groups of patients ( Figures 1 - 4 ). One of the goals of the biomarker validation study was to determine if it would be possible to reduce the number of candidate proteins to a smaller subset to reduce the complexity of a potential diagnostic assay. To this end, classifier training was performed using all biomarkers and selected subsets. Feature selection was performed using two different methods: a) first, an XGBoost classifier was trained on data from all proteins. Because the algorithm performs feature selection during training, the proteins with the highest differentiation power (see Supplementary Figure 3B for a cohort-level comparison of two highest ranking by importance proteins) were used for training another XGBoost classifier on the reduced set of features. b) Proteins were selected based on the p-values of a t-test whereby a number of proteins ranked by increasing p-value (decreasing statistical significance) were chosen for classifier training. Figures 5A–C show the receiver operating characteristic (ROC) curves of classifiers trained using the XGBoost algorithm and a subset of biomarkers that resulted in the best classification accuracy (highest AUC value) in the entire range of the number of proteins selected for training. The classifiers were trained to distinguish between three different contrasts: a combined LD+/LD− cohort and LAD, LD+ vs. LAD, and LD− vs. LAD patients. Classification accuracy in terms of AUC value was determined as a function of the number of proteins included in the training dataset. A comparison of the best classification performance achieved for each contrast when varying number of proteins in the panel was slightly lower for the combined LD vs. LADs with AUC = 0.84 (95 CI: 0.81–0.88) ( Figure 5A ) than the two other contrasts with AUC values of 0.88 (95 CI: 0.84–0.94) ( Figure 5B ) and 0.87 (95 CI: 0.82–0.93) ( Figure 5C ) for LD+ vs. LAD and LD− vs. LAD contrast, respectively. Classification performance as a function of the number of features (proteins) used in training ( Figures 5C–E ) varies only slightly for the combined LD vs. LADs and LD+ vs. LADs. In comparison, the AUC values for the LD− vs. LAD differentiation showed a slight but notable upward trend with the increasing number of features ( Figure 5D ). This suggests that the differentiating information is distributed more broadly in the LD−/LAD than in the LD+/LAD contrast. A comparison of the AUC values obtained with the full set of proteins (n = 45) indicates that all three contrasts show a reduction in classification accuracy when the full set of proteins is used. Interestingly, the LD+ vs. LAD contrast shows improved performance with AUC = 0.88 (95 CI: 0.84–0.94) when using 11 proteins that rank highest by importance/differentiating power as determined by the classifier trained on the full set of proteins ( Figure 5A ). While the LD− vs. LAD classification gave a comparable performance to the LD+ vs. LAD classification when all features were used, there is one substantial difference between them. Interestingly, when comparing the biomarkers ranked by predictive power to distinguish between LD+/LADs and LD−/LADs, it is clear that panel composition differs substantially between the two comparisons ( Table 4 ). Here, the gain parameter represents the fraction of overall classification performance that a particular protein biomarker contributes to the total classifier models shown in Figures 5B, C . The higher values represent stronger differentiating power with the gain values of all biomarkers used adding up to 1. As expected, VlsE (UniProt ID G5IXI6) exhibits by far the most differentiating power (0.37) in separating LD+ from LAD patients. Consequently, the removal of VlsE from the training dataset resulted in a marked drop in the AUC value to 0.71 (data not shown). However, it is also clear that VlsE alone was not sufficient to achieve the demonstrated classification performance ( Figure 5B ) and that additional biomarkers made smaller but substantial contributions to the classification performance. In contrast, VlsE had a distinguishing power of only 0.05 in the LD−/LAD contrast with three other proteins scoring higher. This is consistent with the LD− cohort being negative on standard panels containing VlsE. Overall, the data in Table 4 suggest differences in biomarker panels between the two contrasts, with VlsE being the only protein shared among the 10 highest-ranking biomarkers.

Both feature selection methods (t-test and XGBoost) resulted in comparable classification outcomes, with the t-test-based selection method resulting in lower performance compared to the classifier-based method ( Figures 5C–E ). The classification performance as a function of the number of selected proteins demonstrates that one can substantially reduce the number of antigens in the panel without markedly affecting classification performance. For example, the data in Figure 5D show that one can achieve comparable performance for differentiating between the combined LD and LAD patients with as few as five proteins, whereas only slightly reduced performance can be reached with five and six proteins for LD+ vs. LADs ( Figure 5E ) and LD− vs. LADs ( Figure 5F ), respectively.

In summary, the validation assay data analysis suggests that one can accurately differentiate between the two LD cohorts and the look-alike diseases using a subset of protein biomarkers predicted to have differentiating power using the peptide array data. The most important finding is the fact that the LD− patients who tested previously negative in the standard serology test can be reliably differentiated from the LAD patients.

Seropositive LD (LD+) patient serum samples were obtained from the Lyme Disease Biobank Foundation, Portland, OR (3), CDC, and several commercial biobanks (Boca Biolistics, Pompano Beach, FL; Discovery Life Sciences, Huntsville, AL; and SeraCare, Milford, MA). Samples were collected from patients with signs and symptoms of LD. Samples were tested using the STTT and categorized as seropositive Lyme disease having either an EM rash greater than 5 cm in diameter or PCR confirmation combined with positive STTT serology. The seronegative Lyme disease samples were obtained from patients having an EM rash greater than 5 cm in diameter, but without positive STTT serology, and were obtained from the Lyme Disease Biobank Foundation. These patients were diagnosed with Lyme disease by a physician based on the patients’ clinical symptoms. Participants were enrolled in East Hampton, NY, Central Wisconsin, and Martha’s Vineyard, MA. Each of the three cohorts contained equivalent numbers of patients from each collection site. Cohorts and collection sites were balanced across each assay batch of microarrays. The patient samples for the diseases with similar etiology to LD (look-alike diseases) were obtained from several commercial sources (Boca Biolistics, Pompano Beach, FL; Discovery Life Sciences, Huntsville, AL; Creative Testing Solutions, Tempe, AZ; and SeraCare, Milford, MA). Note that these samples were obtained from commercial biobanks, with no data about their previous exposure to LD provided. Given that the samples were collected outside of the LD endemic areas, it is unlikely that these patients have been exposed to B. burgdorferi infection.

Peptide microarrays containing diverse peptides were synthesized in a commercial production facility (Cowper Sciences, Chandler, AZ), following a previously described library design and photo-lithography-based manufacturing process (34, 36, 37). Briefly, the microarrays used contained 126,051 diverse peptides plus a set of 6,203 control peptides of varying lengths ranging from 6 to 13– amino acids. The standard serum Ab profiling assay protocol described by Arvey et al. (34) was used and modified for a modular research use assay system as described by Kelbauskas et al. (36). Samples were thawed from single-use aliquots and diluted to 1:625 in assay buffer (Phosphate-buffered saline/Tween (PBST) with 0.05% Tween 20, 0.1% ProClin 950, and 1% mannitol). Diluted samples (90 μL) were applied to arrays and incubated for 1 h at 37°C with mixing (TeleShake 95 platform mixer). The cassette was then washed three times in PBST-P using a 96-well microtiter plate washer (BioTek Instruments, Inc., Winooski, VT). Peptide-bound serum antibodies were detected using either 4.0 nM goat anti-human IgG (H+L) conjugated to AlexaFluor 555 (Invitrogen–Thermo Fisher Scientific, Inc., Carlsbad, CA) or 4.0 nM goat anti-human IgM (H+L) (Novus Biologicals, Centennial, CO), conjugated to DyLight 550 in secondary incubation buffer (0.5% casein in PBST-P) for 1 h with mixing at 37°C. After the final incubation, slides were washed three times with PBST-P followed by distilled water to remove residual salts. Slides were then sprayed with isopropanol and dried by centrifugation.

Dried slides were imaged using an ImageXpress imaging system to detect fluorescently labeled secondary antibodies. The imager used an LED light engine (SemRock) centered at 532-nm wavelength to excite fluorophore-conjugated secondary Ab. Mapix (version 7.2.1; Innopsys, Carbonne, France) was used to place a grid alignment file over the obtained images and extract the median foreground pixel intensities using the central 60% of each feature.

Images were inspected to identify arrays with artifacts and image anomalies. The samples associated with such arrays were re-assayed on arrays from the same production batch as the original assay.

Predictive models were built using machine learning methods based on feed-forward, backpropagating fully connected neural networks, similar to those described previously (36, 42). The peptide sequence was one-hot encoded by transforming each peptide into a vector of length 190. The vector length was derived from a maximum peptide length of 10 residue positions with 19 possible amino acids for each position. Feed-forward neural networks were built individually for each donor using R (version 4.2.2, R Foundation for Statistical Computing, Vienna, Austria) as the programming language and utilizing TensorFlow (version 2.11.0) and Keras (version 2.11.1) as the interface packages. The NN models were constructed using three hidden layers with 100 nodes each with a 10% dropout and no layer bias. Rectified linear unit (RelU) activation was used for each layer. Each NN model was trained 10 times using a random 90:10 split of the dataset each time. The data points were weighted by the frequency of peptides appearing in an intensity interval. To this end, the entire intensity range was subdivided into 100 equal bins, and the number of peptides falling into each bin was calculated. The weight for each peptide was computed using the following formula:

where w_i is the weight of the ith peptide and n_i is the number of peptides in the bin that the ith peptide falls into.

The accuracy of the model to predict Ab binding to the array was evaluated by predicting the binding to the held-out 10% of the data and reported as Pearson’s correlation between the measured and predicted binding intensities. Binding to B. burgdorferi epitopes was accomplished by applying the NN models to the B. burgdorferi B31 reference proteome (UniProt Accession # UP000001807) that had been represented as 10-mers with a sliding window of one amino acid offsets.

Random forest decision trees with XGBoost were used to train classifiers for distinguishing patients from the different cohorts used. Each classifier model was trained 10 times on randomly selected 90% of the patients from the corresponding cohorts, and its performance was assessed on the remaining 10% of patients. All training steps and mean receiver operating characteristic curve calculations were performed in R.

The outlier sum statistics was implemented following the method published by Tibshirani et al. (45). The predicted binding intensity values were first z-score normalized using the following:

where I_i , _j is the predicted binding intensity value of the jth tile in the ith sample, and mean and SD are the mean and standard deviation values of I_i , respectively. In this way, binding intensity values were all normalized to their corresponding mean and standard deviation values. Next, the Z-scores of the binding intensity values of the LD+ and LD− samples (“cases”) were calculated using the means and SD values of the LAD samples (“controls”):

where I_{j ( c )} denotes the predicted binding intensities of the jth tile in the control samples. Afterward, a sliding window smoothing with a window size of 5 was applied to the data. Next, for each protein from the B. burgdorferi proteome, the tile with the maximum I i , j d i f f value for each patient in the case cohort was determined. As a result, each protein is represented as its maximum binding value normalized against the control samples as a reference. These maximum binding values were then used to compute the outlier sum (OS) statistics for each protein:

where N is the number of samples in the case cohort, I i p represents the maximum binding value for protein p in sample i, I   p is the maximum binding value of protein p of the samples in the case cohort, and q _0.75 and IQR are the third quartile and interquartile range of I i p , respectively. The p-values for each O S p were calculated using the t-test, and a null distribution was obtained by randomizing the cohort assignments of the samples 1,000 times. The false discovery rate was calculated using the Benjamini–Hochberg adjustment for multiple comparisons.

The functionalization of Luminex MagPlex (Diasorin, Madison, WI) microspheres (beads) was performed by reacting the carboxylic residues of the microspheres and amine groups of proteins using 1-(3-(dimethylamino)propyl)-3-ethyl-carbodiimide hydrochloride (EDAC)/NHS chemistry. Briefly, 236 μL of each address of stock Luminex MagPlex microspheres (1.27 × 10⁷ beads/mL) was resuspended into 764 μL deionized water (DW; 18.2 MΩ·cm), followed by washing with 1 mL of DW. A total of 46 different regions of microspheres each representing a different spectral region for multiplex detection were used for the binding assays. The carboxylic residues of the microspheres were activated by incubating with 90 μL of 50 mg/mL of NHS (Sigma-Aldrich, St. Louis, MO) and 90 μL of 50 mg/mL of EDAC (Sigma-Aldrich) in 1 mL of 0.1 M sodium phosphate (Sigma-Aldrich) buffer for 20 min at room temperature (RT) under gentle rotation. For each bead region, the reaction using 106 microspheres and 5 μg of protein was carried out in 900 μL of 50 mM of MES (pH 5.0) buffer (Avocado Research Chemicals Ltd., Heysham, Lancashire, UK) for 2 h with gentle rotation at RT. After removing the supernatant, the functionalized microspheres were resuspended into 1 mL of PBS-TBN [PBS buffer (Life Technologies, Burlington, ON, Canada) with 0.02% Tween-20 (Sigma-Aldrich), 0.1% bovine serum albumin (BSA; Sigma-Aldrich), 0.02% sodium azide (Sigma-Aldrich), 150 mM sodium chloride (Life Technologies), and 50 mM sodium phosphate monobasic, pH 7.4]. The surface was then blocked with 1% BSA by rotating for 30 min at RT. The beads were washed three times with PBST. Next, all 46 regions of the functionalized microspheres were combined and mixed with 54 mL of PBST-BSA [PBS buffer with 0.1% Tween-20 (Sigma-Aldrich) and 1% BSA] buffer for further use. Microspheres functionalized with the VlsE protein served as positive control for the LD+ cohort, and “blank” beads that went through the same preparation steps, but were not functionalized with a protein, served as negative control. All washing and supernatant removal steps were performed using a MagJET separation rack (Thermo Fisher Scientific, Carlsbad, CA) to separate the microspheres from the solution.

For assays, 50 μL of functionalized microspheres at a concentration of 40 beads/μL (a total of 2,000 beads) was first dispensed into each well of 96-well a non-binding 96-well plate with a flat bottom (Corning, Corning, NY) using a Bravo automated liquid dispensing system (Agilent, Santa Clara, CA). This step was followed by incubation with 50 μL of serum sample (diluted at 1:500 in PBST) for 1 h at 37°C with shaking at 500 rpm. After washing three times with a magnetic microplate washer (Biotek 405 TS, Agilent), 100 μL of goat anti-human IgG secondary antibody (Jackson ImmunoResearch, West Grove, PA) diluted to 1:125 was added to each well and incubated for 30 min at RT with shaking at 500 rpm. After incubation, the beads were washed three times and resuspended in 100 μL PBST buffer. Binding signal intensities were then measured using a Luminex™ xMAP™ IntelliFlex (Diasorin) system.

The binding intensity values for each sample measured in the bead-based assays were normalized by computing the ratio between the intensity values and the intensity of a reference protein. The reference protein (O51141) was selected as a protein having the lowest CV when measured across all cohorts. The assumption was made that such a protein is least affected by the cohort-specific differences in antibody repertoires and thus can be used as a reference to compare binding intensities across patients.