Systems vaccinology identifies clinical and immunological correlates of SARS-CoV-2 vaccine response in solid-organ transplant recipients

Systems vaccinology cohort design for identifying correlates with vaccine response in SOT recipients

To identify the pre-vaccination correlates for vaccine response in an immunosuppressed transplantation cohort, 90 patients were recruited to a systems vaccinology study. The study population included 31 lung and 59 kidney transplant recipients. Sex representation within these cohorts was 67.7% and 55.9% male, respectively and the median age of participants was 59 and 63 years. The median time post-transplantation was 5 years in the lung transplant group and 9.6 years in the kidney transplant group. SOT recipients were on different combinations of immunosuppressive treatments during vaccination, including corticosteroids (CS), antimetabolites such as azathioprine (AZA) and mycophenolate mofetil, calcineurin inhibitors such as tacrolimus (FK506) and cyclosporin A (CsA), and the mTOR inhibitor everolimus. Full demographic breakdown of the study cohort is shown in Table 1. Patients were vaccinated with the Pfizer-BioNTech COVID-19 vaccine Comirnaty and were categorized exhibiting no response (<5 binding antibody units (BAU)/mL), detectable response (detectable threshold at >5 BAU/mL) or positive response (positive threshold at >264 BAU/mL as used in previous studies^32–34).

Overall, vaccine response rates were low, as previously reported.^22,24 The transplantation cohort had 0% of positive responses after a single dose (assessed at day 21, before D2), increasing to 3.4% after two doses (assessed at day 49, 28 days after D2) and 24.0% after three doses (assessed at day ∼188, 28 days after D3). Detectable responses were, however, observed in 6.6%, 39.8% and 49.4% of the transplantation cohort, after 1, 2 and 3 doses of the vaccine, respectively (Supplementary Figure 2). This study design allowed for the determination of demographic, clinical and baseline (pre-vaccination) immunological parameters that correlate with detectable or positive responses, after 2 or 3 doses of COVID-19 mRNA vaccine, respectively. As BNT162b2 vaccine has shown >90% efficacy in preventing COVID-19 illness in the general population³⁵, twenty healthy controls were included in our study for comparing their classical immune profile at baseline to that of SOT recipients. Their clinical characteristics are summarized in Table 1.

• Download figure
• Open in new tab

Supplementary Figure 1. Manual gating strategy for T cell phenotyping.

All acquired events were gated on single live cells before excluding CD14⁺ monocytes. CD4⁺ conventional T cells and CD8⁺ cytotoxic T cells were further separated into naive and memory cells based on expression of CCR7 and CD45RA. CXCR5 and ICOS expression was used to identify T follicular helper cells.

• Download figure
• Open in new tab

Supplementary Figure 2. Anti-RBD IgG antibody levels increase after a third vaccine dose in SOT recipients.

Levels of SARS-CoV-2-specific anti-receptor binding domain (RBD) IgG antibodies were measured using ELISA following A. two or B. three doses of mRNA vaccine. C. The change in antibody levels is shown for separate individuals with lung recipients indicated in blue and kidney recipients indicated in red. Statistics in A and B are based on Mann-Whitney-U and in C on Wilcoxon rank comparisons. Of note, binding antibody levels of 57 kidney transplant recipients are represented after dose 2 because of 2 data not available, binding antibody levels of 50 kidney transplant recipients are represented after dose 3 because of 5 data not available and 4 patients with COVID-19 before dose 3, binding antibody levels of 29 lung transplant recipients are represented after dose 3 because of 2 data not available.

Demographic and clinical correlates for vaccine response in SOT recipients

We first investigated the effect of demographic and clinical variables on the vaccine response. Using a linear regression approach, no significant difference was observed between kidney and lung transplant recipients when assessing detectable vaccine responses after 2 or 3 doses (Table 2). However, a smaller proportion of lung transplant recipients showed a positive response after dose 3 compared to kidney transplant recipients. Within all the SOT recipients, no significant effect was observed for BMI or ethnicity. Although not significant, a sex effect was observed, with males more likely to have a positive vaccine response after dose 3. A complex relationship was observed with age – superior vaccine responses were observed in patients aged <50, while vaccine responses were significantly reduced in patients aged 50-70, consistent with a known association of age with reduced vaccine efficacy.³⁶ This effect was however blunted in the oldest patient cohort, aged between 71 and 90 years old, an observation that can potentially be explained by a longer period since transplantation within this group (Supplementary Figure 3A) especially given that years post-transplantation were positively associated with vaccine response (with poorest responses in those <3 years post- transplantation and significant improvements in those >10 years post-transplantation). Moreover, the reduced effect of age on vaccine response in this oldest patient cohort could be associated with a reduced number of cumulative immunosuppressive treatments impacting the immune response too (Supplementary Figure 3B), as well as a selection bias for survivors of infections suggesting a stronger residual immune response in this older group. We were also able to associate comorbidities of SOT recipients with vaccine response, with the only significant association being a detrimental effect of a estimated glomerular filtration rate (eGFR) <30 mL/min/1.73m² after dose 2 as already described in the literature.³⁷ In chronic kidney disease, uremia has notably been related to T cell defects both by decreasing thymic output and increasing T cell senescence.³⁸ In addition, vaccine-induced immune response was impacted by immunosuppressive treatments, particularly in case of combination of multiple treatments for positive response after dose 3. However, as these drugs are ubiquitous within the poorly- responding transplant cohort, effects could only be assessed in comparison to competing treatment regimens. Using this comparison, MMF treatment was the only significant influence over vaccine response, with significantly poorer outcomes after each dose and using both the detectable and positive thresholds.

• Download figure
• Open in new tab

Supplementary Figure 3. Age-association with time since transplant and combined number of treatments.

A. Association between time since transplant in years and age, grouped as 50-70 or >70 years. B. Age of patients receiving either two or three combined treatments. Statistics are based on Mann-Whitney-U. * = p<0.05; *** = p<0.001.

Multiple immune configurations associated with vaccine response in SOT recipients

We next investigated the effect of the baseline immunological profile (immune profile at day 0 of vaccination) on subsequent antibody production. Blood samples from each patient were assessed by flow cytometry for 444 immunological variables, constituting a deep phenotypic analysis of the adaptive immune compartment with extensive innate immune profiling. Two major divisions of the transplant recipients cohort were used: patients capable of some response after two doses of vaccine (“detectable” after second dose versus “undetectable” after second dose), and patients capable of a strong antibody response (>264 BAU/mL) after three doses of vaccine (“positive” after third dose versus “not-positive” after third dose). Of note, these 2 categories represent two different immune profiles at baseline associated with immune response after 2 and 3 doses of mRNA vaccine. However, the majority of patients showing a positive response after dose 3 had a detectable response after dose 2 implying an overlap of patients between these 2 baseline profiles. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, to identify the immunological associates corresponding to vaccine response. Of the immunological parameters associated with vaccine response after demographic adjustment (Supplementary Spreadsheet 3), the 10 parameters with greatest magnitude effects were further analyzed (Figure 1). For some response (“detectable”) after dose 2, the strongest positively associated parameter was CD80 expression, with CD80+ median expression on double negative CD27^-IgD^- B cells DN1 (CD11c^-CD21⁺) and CD21^- B cells including CD21^-CD27⁺ B cells along with the frequency of CD80+ double negative B cells DN2 (CD11c⁺CD21^-) (Figure 1A). For strong vaccine responses (“positive” after dose 3), the key immunological parameters differed: the strongest positive associations were with non-classical monocytes and CD40⁺ plasmacytoid dendritic cells, within their respective parent populations, along with CD40⁺ median expression in CD21^- B cells and double negative B cells DN2 (Figure 1B). In both cases, the immunological dataset was from baseline (day 0) data and showed associations of different immunological parameters between some response after dose 2 and a strong response after dose 3. These associations highlight the different baseline profiles associated with different outcomes. One of the main differences was the baseline T cell activation that was rather predictive of a negative response after dose 3. Thus, the immunological profiles associated with a detectable vaccine response after a two-dose mRNA vaccine schedule were distinct from the immunological profiles associated with a positive vaccine response after a third injection of mRNA vaccine (booster) among SOT recipients.

• Download figure
• Open in new tab

Figure 1. Identification of key immune parameters associated with effective SARS-CoV-2 vaccine response.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, to identify the immunological associates with vaccine response. A. Odds ratio (OR) of the z-score of the 10 highest difference parameters (magnitude cut-off), using the stratification of transplant recipients “detectable” or “undetectable” for SARS-CoV-2 antibodies after two doses of the mRNA vaccine (day 49). B. Odds ratio (OR) of the z-score of the 10 highest difference parameters (magnitude cut-off), using the stratification of transplant recipients “positive” or “non-positive” for SARS-CoV-2 antibodies after three doses of the mRNA vaccine (day ∼188).

In identifying immunological parameters associated with effective vaccine response after the second and the third dose, our a priori hypothesis was that the responding phenotype, as opposed to the non-responding phenotype, of SOT recipients would be more similar to the immune status at baseline of healthy subjects. For each of the main immune associates, we compared the values of both responder and non-responder SOT recipients to those of healthy controls. For the ability to have some response after mRNA vaccination (dose 2 “detectable”), of the 10 strongest immune associates in SOT recipients, surprisingly 7 showed a clear trend of the non-detectable SOT population being more similar to the healthy controls, while 1 showed the expected pattern of the detectable SOT cohort being more similar to the healthy controls (Figure 2A). These immune parameters did not seem to be directly impacted by the immunosuppressive treatment regimen. The strongest suppression of vaccine responses was associated with MMF and CS treatment (Table 2), yet segregation of “detectable” and “non-detectable” patients by MMF (Supplementary Figure 4A) or CS (Supplementary Figure 5A) found few consistent patterns. One exception was the frequency of CD38+ monocytes, a rare case where the “detectable” SOT population had a phenotype closer to that of healthy controls. This specific association was driven instead by CS treatment, and assessing only the CS untreated patients, the “detectable” population was actually further away from the healthy control population (Supplementary Figure 5A).

• Download figure
• Open in new tab

Supplementary Figure 4. Influence of MMF treatment on key immune parameters associated with vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, identifying the 10 highest magnitude effects (based on log(OR of the z-score)). A. Raw value graphs, showing healthy controls, dose 2 non-detectable and dose 2 detectable SOT recipients, for the 10 highest magnitude parameters associated with dose 2 antibody detection stratified by treatment with MMF (circle indicates no MMF treatment, cross indicates MMF treatment). Individual values, no imputation. B. Raw value graphs, showing healthy controls, dose 3 non- positive and dose 3 positive SOT recipients, for the 10 highest magnitude parameters associated with dose 3 antibody positive seroconversion stratified by treatment with MMF (circle indicates no MMF treatment, cross indicates MMF treatment). Individual values, no imputation. Statistics are based on two-tailed individual t-tests comparing patients receiving treatment or not within each response group and only significant comparisons are shown. * = p<0.05.

• Download figure
• Open in new tab

Supplementary Figure 5. Influence of CS treatment on key immune parameters associated with vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, identifying the 10 highest magnitude effects (based on log(OR of the z-score)). A. Raw value graphs, showing healthy controls, dose 2 non-detectable and dose 2 detectable SOT recipients, for the 10 highest magnitude parameters associated with dose 2 antibody detection stratified by treatment with corticosteroids (CS; circle indicates no CS treatment, cross indicates CS treatment). Individual values, no imputation. B. Raw value graphs, showing healthy controls, dose 3 non-positive and dose 3 positive SOT recipients, for the 10 highest magnitude parameters associated with dose 3 antibody positive seroconversion stratified by treatment with CS (circle indicates no CS treatment, cross indicates CS treatment). Individual values, no imputation. Statistics are based on two-tailed individual t-tests comparing patients receiving treatment or not within each response group and only significant comparisons are shown. * = p<0.05; ** = p<0.01.

• Download figure
• Open in new tab

Figure 2. Key immune changes associated with effective SARS-CoV-2 vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, identifying the 10 highest magnitude effects (based on log(OR of the z-score)). A. Raw value graphs, showing healthy control, dose 2 non- detectable and dose 2 detectable SOT recipients, for the 10 highest magnitude parameters associated with dose 2 antibody detection. Individual values, no imputation. B. Raw value graphs, showing healthy controls, dose 3 non-positive and dose 3 positive SOT recipients, for the 10 highest magnitude parameters associated with dose 3 antibody positive seroconversion. Individual values, no imputation. Statistics are based on Kruskal-Wallis with Dunn’s post-hoc comparison and only significant comparisons are shown. * = p<0.05; ** = p<0.01.

In contrast to the analysis for some ability to respond, for the ability to have a strong response after a booster dose of mRNA vaccine (“positive” after dose 3), most of the parameters showed a clear trend for the positive SOT population being more similar to the healthy controls, while parameters on double negative CD4^-CD8^- T cells that had the non-positive SOT population were more similar to the healthy controls (Figure 2B). Again, these immune parameters were largely not directly impacted by treatment regimen (Supplementary Figure 4B, Supplementary Figure 5B). These results suggest two distinct functional immune profiles in SOT recipients: a specific non-classical immune profile, capable to develop a detectable humoral response after 2 doses of mRNA vaccine, and a distinct immune profile closer to the healthy population, which is related to a better response to vaccination after a third dose. The non-classical immune pattern might be driven by “unconventional” immune cell subsets, such as double negative (DN) B cells, a cellular subset associated with antibody-mediated response which is less abundant in the general population but has been involved in autoimmune diseases and in vaccine response.^39,40 In our cohort, all responders displayed a higher amount of these DN B cells (Supplementary Figure 6A). In contrast, an immune profile of SOT recipients closer to the one found in the healthy population was able to elicit higher antibody levels after additional vaccine administration (dose 3). Recently, identification of an adaptive immune signature at baseline prior to vaccination in a healthy population was associated with a weaker antibody response²⁹, and interestingly, our data show that T cell activation (such as HLA-DR expression on CD4⁺ T cells) reduced the chance to develop high antibody levels after the third vaccine dose (Figure 2B).

• Download figure
• Open in new tab

Supplementary Figure 6. Changes in DN and transitional B cell subsets associated with effective SARS-CoV-2 vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables. A. Raw value graphs, showing healthy control, dose 2 non-detectable and dose 2 detectable SOT recipients (left panel) and healthy controls, dose 3 non-positive and dose 3 positive SOT recipients (right panel), for double negative (DN) B cells. Individual values, no imputation. B. Raw value graphs, showing healthy controls, dose 2 non-detectable and dose 2 detectable SOT recipients (left panel) and healthy controls, dose 3 non-positive and dose 3 positive SOT recipients (right panel), for double negative (DN) B cells. Individual values, no imputation. Statistics are based on Kruskal-Wallis with Dunn’s post-hoc comparison and only significant comparisons are shown. *** = p<0.001; **** = p<0.0001.

To assess this immunological configuration model, we performed two analyses of the global immune profile. First, we used a Uniform Manifold Approximation and Projection (UMAP) approach of the baseline immunological profile of SOT recipients, responders and non- responders, and healthy controls; considering the baseline immune profile of this latter group associated with a high immune response to 2 doses of mRNA vaccination as reported previously.³⁵ For some level of response (“detectable” after dose 2), a slight separation was observed between “detectable” and healthy controls when displaying all 444 immune parameters; a stronger separation when only the 10 most highly associated immune parameters were displayed (Figure 3A). By contrast, when assessing the immunological profile of strong responders (“positive” after dose 3), the responders clustered closer to the healthy controls (Figure 3B), consistent with the individual parameter assessment described above. Second, we looked at the relationship between immune parameters using Spearman’s correlations. This allowed us to look at the relationship of the key immune parameters in healthy controls, and whether that relationship was preserved or disrupted in responding or non-responding SOT recipients. For the top 10 associated immune parameters after dose 2, we observed strong correlations between different immune characteristics within the healthy controls, for example a strong co-correlation of the CD80-associated phenotypes, and a strong negative correlation between CD57 expression on double negative T cells and the CD80-associated phenotypes (Figure 4A). In the SOT recipients capable of some response (“detectable”) after dose 2, large disturbances to these correlations were observed, while SOT recipients that did not respond (“not detectable”) after dose 2 showed more similar correlations as the healthy population (Figure 4A,B). By contrast, the disruptions in co-correlations in the non-responding and responding SOT recipients groups after dose 3 were similar (Figure 4C,D). Together, these analyses support a model where some SOT recipients exhibited a non-classical immune profile, capable of sub-optimal response after 2 vaccine doses, while some SOT recipients showed a more “normalized” immune profile in certain immune subsets, allowing higher levels of antibodies if an additional vaccine dose was given.

• Download figure
• Open in new tab

Figure 3. Immune configuration associated with effective SARS-CoV-2 vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables, and SOT recipients were assessed for ability to generate binding antibodies after vaccination. Immune configurations displayed as UMAP plots for listed parameters for each individual, annotated based on patient characteristics. Distances between the groups were calculated using the Calinski and Harabasz score and are shown on the right of each panel. A. UMAP for all 444 all immune parameters, with individuals annotated as healthy controls, or SOT recipients with non-detectable or detectable antibodies after dose 2. B. UMAP for 10 highest magnitude parameters associated with “detectable” responses after dose 2, with the same patient annotation. C. UMAP for all 444 all immune parameters, with individuals annotated as healthy controls, or SOT recipients with non-positive or positive antibodies responses after dose 3. D. UMAP for 10 highest magnitude parameters associated with “positive” responses after dose 3, with the same patient annotation.

• Download figure
• Open in new tab

Figure 4. Immunological relationships associated with effective SARS-CoV-2 vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables, and SOT recipients were assessed for ability to generate binding antibodies after vaccination. Multivariate logistic regression was performed adjusting by type of transplant, sex, age and the number of years since transplantation, identifying the 10 highest magnitude effects (based on log(OR of the z-score)). Spearman’s correlations were calculated between each of the 10 highest magnitude effect variables. A. For those immune parameters associated with dose 2 “detectable”, Spearman’s correlations between each parameter pair are shown for healthy controls (left). For dose 2 “undetectable” (center) and “detectable” (right) SOT recipients, the upper-right triangular matrix shows the Spearman’s correlations within the respective patient population, while the lower-left triangular matrix visualizes by color and size the deviations in the correlations observed between the patient population and the healthy population with green indicating high correlation within the healthy population and low correlation within the patient population and brown low correlation within the healthy population and high correlation within the patient population. B. Values of the differential correlation for each cell population associated with a detectable response after dose 2, comparing the healthy population and the respective patient population as individual values (left) and violin plot (right). C. For those immune parameters associated with dose 3 “positive”, Spearman’s correlations between each parameter pair are shown for healthy controls (left). For dose 3 “non-positive” (center) and “positive” (right) SOT recipients, the upper-right triangular matrix shows the Spearman’s correlations within the respective patient population, while the lower-left triangular matrix visualizes by color and size the deviations in the correlations observed between the patient population and the healthy population. D. Values of the differential correlation for each cell population associated with a positive response after dose 3, comparing the healthy population and the respective patient population as individual values (left) and violin plot (right).

Machine learning identifies immunological parameters predicting vaccine response

Our results demonstrate the influence of baseline immune configuration on vaccine response in SOT recipients. To identify those immune parameters with the highest potential utility in diagnosis, we used a machine learning approach. Multiple machine learning models were used and compared by recursive feature elimination, with 10-fold cross-validation to assess the explanatory power of the model, with the area under the ROC curve (AUC). For both the dose 2 “detectable” and dose 3 “positive” classifications, the support vector machine approach gave the highest accuracy (Table 3). For some level of response (“detectable” after dose 2), the optimal machine learning model had a predictive capacity (AUC) of 95% (Figure 5C), while for strong responses (“positive” after dose 3), the optimal machine learning model had a predictive capacity (AUC) of 93.3% (Figure 5D). To identify the immunological parameters contributing to the predictive model, we used a cluster permutation approach. For some response (“detectable” after dose 2), the strongest single parameter was the frequency of CD40⁺ B cells, with NK and NKT phenotypes also strongly contributing, and CD80⁺ median expression in switched memory B cells as a strong influence (Figure 5A). For strong response (“positive” after dose 3), the strongest single parameter was CD40⁺ median expression in non-switched memory B cells, with CD40 expression in naive B cells also ranking highly, along with NK phenotypes (Figure 5B). Notably, while there is some overlap with the linear regression multivariate analysis used above, the lists of features are distinct, as highly co-correlated parameters (such as CD80 expression across multiple B cell subtypes) provide decreasing amounts of additional information, and therefore do not dominate the influential parameter list in machine learning approach. Together, this analysis suggests the key markers aiming to stratify SOT recipients into different functional groups regarding vaccine response outcome.

• Download figure
• Open in new tab

Figure 5. Machine learning predicts SARS-CoV-2 vaccine response in SOT recipients.

110 individuals were recruited to a systems vaccinology study, including 20 healthy controls, 31 lung transplant recipients and 59 kidney transplant recipients. Immunological profiles were assessed at baseline (day 0 of vaccination) for 444 immunological variables, and SOT recipients were assessed for ability to generate binding antibodies after vaccination. The support vector machine approach was selected to model vaccine response, based on the highest performance. A. Best features according to the machine learning model to identify the detection of antibodies against COVID-19 after the second dose of the vaccine. The best features were chosen using the feature selection algorithm RFE-SVM and the evaluation of the importance of each feature in the model was done using clustered permutation. B. Best features according to the machine learning model to identify the positive seroconversion against COVID-19 after the third dose of the vaccine. The best features were chosen using the feature selection algorithm RFE-SVM and the evaluation of the importance of each feature in the model was done using clustered permutation. C. ROC curve of the machine learning algorithm developed for predicting “detectable” after dose 2, as evaluated by nested cross-validation. Curve is mean performance in test groups during 10-fold cross-validation, with fill area the standard deviation. D. ROC curve of the machine learning algorithm developed for predicting “positive” after dose 3, as evaluated by nested cross-validation. Curve is mean performance in test groups during 10-fold cross- validation, with fill area the standard deviation.

Study design and patient’s characteristics

The study was designed to assess clinical and baseline immunological parameters predicting the humoral response to BioNTech/Pfizer BNT162b2 mRNA (Comirnaty®) vaccination in lung transplant and kidney transplant recipients, naive for SARS-CoV-2 infection. A cohort of health- care workers from Hôpital Saint-Pierre in Brussels was used as a control group regarding baseline immune parameters.

The ethics committee of Hôpital Erasme, Brussels, Belgium (references P2021/182 and A2021/156 for lung transplant recipients; P2020/284 and A2021/131 for kidney transplant recipients) and the Belgian Federal Agency for Medicines and Health Products (FAMHP, EudraCT 2021-004992-16 for lung transplant recipients; EudraCT 2021-000-412-28 for kidney transplant recipients) approved the monocentric prospective phase IV investigator-initiated study of the immunogenicity of the BNT161b2 vaccine (Pfizer-BioNTech). The ethics committee of CHU Saint-Pierre approved a 6-month study on the carriage and seroprevalence of SARS-CoV- 2 in health-care workers (CE/20-04-17).⁸⁰

All participants of at least 18 years of age provided written informed consent for each study, respectively. Lung transplant and kidney transplant recipients from Chest department and Nephrology, dialysis, and transplantation department in Hôpital Erasme, Belgium, respectively and health-care workers from Hôpital Saint-Pierre were enrolled before administration of the first dose of BNT162b2 vaccine. Participants with a previous laboratory-confirmed SARS-CoV-2 infection or with a baseline anti-RBD IgG titer >5 BAU/mL were considered previously infected and excluded from this study. Demographics of the healthy and patient populations are given in Table 1. SOT recipients were vaccinated with the Pfizer-BioNTech COVID-19 Vaccine Comirnaty on the day of blood sampling, with second and third doses given on days 21 and between day 150 and 170, respectively. Participants had blood drawn to collect peripheral blood mononuclear cells (PBMCs) just prior to the first vaccine dose (baseline), and serum three weeks after D1 prior to the second dose (post-D1), four weeks after the second dose (post-D2) and four weeks after the third dose (post-D3).

Sample processing and isolation of PBMCs

Serum tubes were centrifuged at 3000 rpm for 10 min to separate serum. Serum was stored at - 80°C for antibody analysis. PBMCs were isolated from heparinized whole blood by density gradient centrifugation as recommended (Leucosep tubes, Greiner). Whole blood was diluted with an equal amount of Dulbecco’s Phosphate-Buffered Saline with 2% Fetal Bovine Serum and centrifuged. Upper plasma layer was discarded without disturbing the interface. The PBMC layer was retained and washed three times with medium. Cell pellets were diluted with HBSS and aliquoted. PBMCs were frozen progressively with 10% DMSO in FBS and stored in liquid nitrogen.

SARS-CoV-2-specific binding antibodies by enzyme-linked immunosorbent assay

Binding antibodies at baseline and after vaccination were assessed manually using an enzyme- linked immunosorbent assay (ELISA) for the quantitative detection of IgG-class antibodies to RBD (Receptor Binding Domain, Wuhan strain) [Wantai SARS-CoV-2 IgG ELISA (Quantitative); CE-marked; WS-1396; Beijing Wantai Biological Pharmacy Enterprise Co., Ltd, China] as previously described.^81,82 For quantification of antibodies, diluted serum samples (1/10, 1/100, 1/400, 1/1600 and 1/6400) were tested with an internal standard, calibrated against NIBSC 20/136 (First WHO International Standard Anti-SARS-CoV-2 Immunoglobulin), and an external positive control sample included on each plate. Diluted samples were incubated (37°C, 30 min) with pre-coated microwells and washed five times. Next, plates were incubated with horseradish peroxidase (HRP)-conjugated anti-human IgG antibodies (37°C, 30 min) and washed five times before adding a TMB and urea peroxide solution for 15 min (37°C°, dark). After incubation, a stop solution (0.5 M H2SO4) was added, and optical density (OD) was measured at 450nm using a Bio-Rad iMark Microplate Absorbance reader (cat#1681130). Net OD values were converted to arbitrary IgG units per mL by interpolation from a point-by-point plot fitted with the standard concentrations and net OD values (correlation coefficient R2≥0.9801), using GraphPad Prism version 9.0.0 for Windows (GraphPad Software, San Diego, California, USA) and exported to Microsoft Excel. Antibody measurements were adjusted for any sample dilution, converted to binding antibody units per mL (BAU/mL) and reported as such. Lower limit of quantification (LLQ) was 5 BAU/mL according to the WHO standard. Clinical performance characteristics of the assay, evaluated in 69 PCR-confirmed COVID-19 patients (comprising mild and severe clinical outcomes, ≥ 15 days post onset of symptoms) and 167 pre-pandemic sera, resulted in a specificity of 100% (95% CI 97,75-100) at a sensitivity of 100% (95% CI 94,73-100) for a cut-off of 6 BAU/mL. Data lower than the limit of detection were attributed the value 2.5 BAU/mL before transformation. Responses reaching the threshold of >264 BAU/mL, as previously used^32,34, were considered “positive”, while responses >5 BAU/mL were “detectable”.

Flow cytometry

Frozen PBMCs were thawed in complete RPMI, washed twice with PBS, and stained with live/dead marker (fixable viability dye eFluor780; eBioscience) as well as fluorochrome- conjugated antibodies against surface markers: anti-CD3 (REA613, Miltenyi Biotec), anti-CD123 (7G3), anti-CD80 (L307.4), anti-CD21 (B-ly4), anti-CD27 (L128), anti-BAFF-R (11C1), anti-CD94 (HP-3D9), anti-CD8 (SK1), anti-CD86 (2331 (FUN-1)), anti-CD141 (1A4), anti-CD56 (NCAM16.2), anti-CD4 (SK3), anti-CD16 (3G8), anti-CD40 (5C3, all BD Biosciences), anti-CD57 (HNK-1), anti-CD24 (ML5), anti-HLA-DR (L243), anti-CD19 (HIB19), anti-IgM (MHM-88), anti-CD11c (3.9), anti-IgD (IA6-2), anti-CD10 (HI10a), anti-CD38 (HB-7, all BioLegend), anti-CD14 (TuK4, eBioscience). Cells were then fixed using the Foxp3/Transcription factor buffer staining set (eBioscience) and data were collected on a BD FACSymphony A3 Cell Analyzer (BD Biosciences). A maximum of 5×10⁵ events were acquired for each sample. The data of 5 kidney transplant recipients were excluded from downstream analysis due to insufficient sample quality.

A subset of samples was also used to perform an additional staining using markers focused on the T cell compartment. Cells were thawed and plated as described above, before staining with a live/dead marker (fixable viability dye eFluor780, eBioscience) was performed. After washing, cells were stained with fluorochrome-conjugated antibodies against surface markers: CD14 (61D3), CD4 (RPA-T4), CD45RA (HI100), CCR7 (3D12), CXCR5 (MU5UBEE), ICOS (ISA-3, all eBioscience), CD8 (RPA-T8, Biolegend). Cells were washed again and fixed using the Foxp3/Transcription factor buffer staining set (eBioscience) before acquisition on a BD Canto II (BD Biosciences). A maximum of 2.5x10⁵ events were acquired for each sample. Raw .fcs data were compensated using AutoSpill.⁸³ The complete set of FCS files used has been deposited on FlowRepository and annotated in accordance with the MIFlowCyt standard. These files may be downloaded for further analysis from http://flowrepository.org/id/FR-FCM-Z662.

Flow cytometry gating

Compensated .fcs files were manually pre-processed in FlowJo (LLC, Ashland, Oregon; v10.8) to only contain events with stable fluorescent signals. For samples acquired on the FACSymphony, gating on singlets, leukocytes, live cells, and subpopulations was performed in R (version 4.0.2) by concatenating 7 individual samples and automatically calculating adequate boundaries between populations using the gating strategy outlined in Supplementary Spreadsheet 1. These boundaries were then applied to every individual sample to determine cell frequencies and median fluorescence intensity (MFI) values of activation markers on relevant populations (Neumann et al., manuscript in preparation). For samples acquired on the BD Canto II, manual gating was performed selecting singlets, lymphocytes and living cells before excluding CD14⁺ monocytes. CD4⁺ and CD8⁺ T cells were further divided based on their CCR7/CD45RA expression profile and expression of ICOS and CXCR5 was evaluated on non- naive cells (excluding CCR7⁺CD45RA⁺ cells). For the gating strategy see Supplementary Figure 1.

Association between demographic and medical variables with vaccine response

To assess the association of demographic characteristics, comorbidities and treatments with the vaccine response induced by mRNA vaccination, sex and age were considered a priori confounders, i.e., always considered confounders regardless of their effect on the relationship. The other demographic variables and treatment regimens were considered potential confounders. They were only considered confounding if a significant impact on the model between exposure evaluated and the outcome was observed. The relationship was evaluated using multivariate logistic regression. Stepwise regression bidirectional elimination was used to select among potential confounders. The stepwise criterion of permanence was a significant p- value (<0.05) of the potential confounder, thus characterizing its importance for the model.

Batch correction

For flow cytometry, a batch variation in channel intensity is inherent to the measurement process. Here, we first verified whether there were batch differences in the cytometry data. We found flow cytometry readouts with a significant difference between batches (t Student/Mann- Whitney test), while other demographic/treatment variables of the study did not show significant differences between batches (t Student/Mann-Whitney test for continuous variable and Pearson’s chi-squared for categorical variable). This indicated that the batch effect differences were due to technical variation in the flow cytometer and staining procedure. To correct measurement distortions due to day-to-day variation we located a cell population that had the smallest proportional variation between samples for each immune marker (considering populations of cells with more than 30 cells) and used their mean to correct the respective marker channel.

Identification of immune variables associated with vaccine response

For each of the 444 immunological variables, we assessed their association with vaccine response by considering only SOT recipients divided based on: i) the presence of antibodies against SARS-CoV-2 detectable/undetectable after the second dose; ii) the presence of antibodies against SARS-CoV-2 detectable/undetectable after the third dose; and iii) the presence of antibodies against SARS-CoV-2 above the positive threshold/below the positive threshold after the third dose. This stratification of patients was based on the very low number of SOT recipients reaching the positive threshold, particularly after only two doses, requiring analysis to focus on achieving the detectable antibody reading, or on reaching the positive threshold after the third dose only. As each immune marker may have different levels of measurements, to compare the magnitude of the effect of each marker on the vaccine response, a z-score transformation was applied. The evaluations of the associations of the immune markers in vaccine response were performed through multivariate logistic regression. Sex and age were considered a priori confounders. Stepwise regression bidirectional selection⁸⁴ was used to select additional potential confounders. The stepwise criterion for selecting features was Akaike information criterion⁸⁵, keeping only potential confounders that improve the model. This resulted in adjustment for sex, age, type of transplant, and the time since transplant. Among the immunological markers that obtained p-values lower than 0.05, the 10 that had the greatest magnitudes of effect were chosen for further analysis. Logistic regression models were run in Python using the package statsmodels version 0.13.2.

Similarities between immune markers in healthy and transplant recipients

To understand the similarity of the immune response between healthy controls and SOT recipients, UMAP dimension reduction technique⁸⁶ was performed. To do so, markers with more than 30% of missing values were first removed and subsequently the remaining missing values were imputed using the mean. After that, a standardization transformation was performed. To assess the distances between individuals, the Calinski and Harabasz score⁸⁷ was used. In addition to understanding the immune similarity between SOT recipients and healthy controls, it is also important to understand whether the relationships between immune markers change within individuals. For that we used a heatmap for Spearman correlation and compared these correlations in SOT recipients with healthy controls.

Machine learning

Data was cleaned and standardized as described above with standardization applied within each group. For selection of important markers contributing to the outcome by machine learning, we used recursive feature elimination⁸⁸, with 10-fold cross-validation using the area under Receiver Operating Characteristic (ROC) curve as readout. The model used for feature selection was support vector machine (SVM). The choice of model hyperparameters was performed by grid search with 3-fold cross validation using balanced accuracy as a selection criterion.

Due to the multicollinearity shared between the markers, the potential importance of the marker can be underestimated as there are other similar markers contributing similar information. Therefore, to limit this potential loss of important markers, we applied hierarchical clustering of markers based on the dissociation of the markers given 1 minus modulus of the Spearman correlation. After the identification of these clusters, the assessment of the importance of each marker was performed using permutation importance⁸⁹ using all markers and excluding markers contained in the same group of the evaluated marker. Thus, the marker’s importance in the phenotype is its impact on the model containing only itself in its group and a marker from every other group. The model used in the permutation presented the same hyperparameters selected in the feature selection. Permutation was performed with 100 repetitions and evaluated by area under ROC curve (AUC).

The machine learning models were created using the best features to reduce dimensionality. The best models were selected and evaluated by nested cross-validation, applying 10 external and 5 internal cross-validations. The selection of the hyperparameters of the models was done by grid search using balanced accuracy as a criterion.

Data and code availability

Raw data for the study is available from https://flowrepository.org/experiments/6338 (FR-FCM- Z662). Source data for the analysis presented in the main and supplementary text is available as Supplementary Spreadsheet 2. Source code for the analysis is available on GitHub at: https://github.com/rafael-veiga/immunological_covid.