54 patients from Cruces University Hospital and Donostia University Hospital, who suffered from MM and received autoHSCT were included in the study. Clinical data of the patients can be found in Table 1 . Blood samples were taken from patients at six different time points. Sample 1 (S1) was collected the day before the start of the conditioning treatment; sample 2 (S2) was collected after leukocyte recovery (>1000 leukocytes/µL) (median: 13 days after autoHSCT, range: 10-25 days); sample 3 (S3) was collected 30 days after cell infusion; sample 4 (S4) was collected 100 days after cell infusion; sample 5 (S5) was collected 180 days after cell infusion and sample 6 (S6) was collected 365 days after cell infusion. Sample collection was carried out through the Basque Biobank for Research (http://www.biobancovasco.org), which complies with the quality management, traceability and biosecurity, set out in the Spanish Law 14/2007 of Biomedical Research and in the Royal Decree 1716/2011. The study was approved by the Basque Ethics Committee for Clinical Research (BIO14/TP/003 and PI+CES+INC-BIOEF 2017-03). All subjects provided written and signed informed consent in accordance with the Declaration of Helsinki.

Plasma and peripheral blood mononuclear cells (PBMCs) were obtained from blood samples. PBMCs were isolated by Ficoll (GE Healthcare) density gradient centrifugation, cryopreserved in heat-inactivated Fetal Bovine Serum (FBS) (GE Healthcare Hyclone) containing 10% Dimethylsulfoxide (DMSO) (Thermo Scientific Scientific) and stored in liquid nitrogen until they were used. Plasma samples were stored at -80⁰C until they were used.

Prior to flow cytometry experiments, cryopreserved PBMCs were thawed at 37⁰C in a water bath and washed twice with RPMI 1640 medium supplemented with L-Glutamine (Lonza). Then, cells were incubated for 1 hour at 37⁰C and 5% CO₂ with 10U DNase (Roche) in R10 medium (RPMI 1640 medium containing GlutaMAX from Thermo Fisher Scientific, 10% FBS and 1% Penicillin-Streptomycin (P-S) from Thermo Fisher Scientific). Afterwards, cells were washed once, resuspended in NK cell medium (RPMI 1640 medium containing GlutaMAX, 10% FBS, 1% P-S, 1% MEM Non-Essential Amino Acids Solution and 1% Sodium Pyruvate, both from Thermo Fisher Scientific), filtered using 70μm cell strainers and counted just before further use.

The absolute number of NK cells was measured in each blood sample. NK cells were defined as CD45+, CD3-, CD56+ and/or CD16+. The following clinical grade fluorochrome conjugated monoclonal antibodies (mAbs) were used: FITC anti-CD16 (CLB/FcGran1), PE anti-CD56 (MY31), PerCP-Cy5.5 anti-CD3 (SK7) and V450 anti-CD45 (2D1) from BD Biosciences. The absolute number per µL of blood was calculated based on the following formula: (percentage of NK cells within lymphocyte gate X absolute number of lymphocytes)/100. The absolute number of lymphocytes was obtained from the hemogram.

For phenotypic analysis of NK cells, PBMCs were first stained with LIVE/DEAD Fixable Aqua Dead Cell Stain Kit (Invitrogen) reagent to exclude dead cells, following manufacturer´s recommendations. Afterwards, cells were washed with PBS containing 2.5% bovine serum albumin (BSA) (Sigma-Aldrich) and extracellular staining was performed. For that, cells were incubated for 30 minutes on ice in the dark with the following fluorochrome conjugated mAbs: BV421 anti-CD56 (NCAM 16.2), BV510 anti-CD3 (UCHT1), BV510 anti-CD14 (MφP9), BV510 anti-CD19 (SJ25C1) from BD Bioscience; PE anti-NKG2C (134591) from R&D Systems; PE-Vio770 anti-NKp80 (4A4.D10), APC-Vio770 CD57 (REA769) and APC-Vio770 CD69 (REA824) from Miltenyi Biotec and APC anti-NKG2A (Z199) from Beckman Coulter. Cells were then washed again with PBS containing 2.5% BSA and fixed and permeabilizated prior to intracellular staining. Two different intracellular staining protocols were used. For “phenotype panel 1” (see Supplementary Table 1 ), PBMCs were fixed by incubating them with 4% paraformaldehyde (PFA) (Sigma-Aldrich/Merck) for 15 minutes on ice. Next, cells were washed twice with PBS containing 2.5% BSA and permeabilized by incubating with 1x BD Perm/Wash Buffer (BD Bioscience) for 15 minutes at room temperature (RT). Then, cells were incubated for 30 minutes on ice with FITC anti-FcϵRIγ (Merck) antibody. On the other hand, for “phenotype panel 2” (see Supplementary Table 1 ), cells were fixed and permeabilized with Foxp3/Transcription Factor Staining Buffer Set (eBioscience) following manufacturer´s recommendations. Then, cells were incubated for 30 minutes at RT with the following mAbs: FITC anti-Ki67 (20Raj1, Invitrogen) and PE anti-Granzyme B (GB11, BD Bioscience). In both intracellular staining protocols, after the incubation with the corresponding antibodies, cells were washed with the corresponding permeabilization buffer and resuspended in PBS. Samples were then acquired in a MACSQuant Analyzer 10 flow cytometer (Miltenyi Biotec). Flow cytometry panels used to study NK cell phenotype are shown in Supplementary Table 1 .

For functional assays, after thawing and counting the cryopreserved PBMCs, cells in NK cell culture medium were plated at 0.5 x 10⁶ cells/well in 48-well plates. PBMCs were then primed with 10 ng/mL of recombinant human (rh) interleukin (IL)-15 (Miltenyi Biotec) and cultured for 20-21 hours. Then, for cytokine stimulation, 10 ng/mL of rhIL-12 (Miltenyi Biotec) and 50 ng/mL of rhIL-18 (MBL International) were also added in addition to rhIL-15. For target cell stimulation, 721.221 cells were added after the 20-21 hours of culture with rhIL-15 at Effector:Target (E:T) 1:1 ratio (0.5x10⁶ PBMCs and 0.5x10⁶ 721.221 cells). Then, IL-12+IL-15+IL-18 and 721.221 stimulated PBMCs were cultured for 5 additional hours. PE anti-CD107a (REA792, Miltenyi Biotec) was added at the start of the co-culture period to all the conditions and GolgiStop (monensin) and GolgiPlug (brefeldin A) protein transport inhibitors (BD Biosciences) were added after 1 hour for the rest of the incubation time following manufacturer’s protocol. Afterwards, viability and extracellular staining was performed as explained in the Phenotypic Analysis section. Next, cells were fixed with 4% PFA and permeabilized using 1x BD Perm/Wash Buffer. Then, cells were incubated for 30 minutes on ice with PerCP-Cy5.5 anti-interferon (IFN)γ (B27, BD Biosciences) and APC anti-tumor necrosis factor (TNF) (MAb11, BioLegend). Samples were acquired in a MACSQuant Analyzer 10 flow cytometer. The flow cytometry panel used to study NK cell functionality is shown in Supplementary Table 1 . The percentage of NK cells positive for CD107a, IFNγ and TNF was calculated after subtracting the non-stimulus condition. 721.221 cell line was cultured in NK cell medium supplemented with 5 µg/mL Plasmocin (InvivoGen). 721.221 cell line was routinely tested for mycoplasma infection with Venor GeM Classic detection kit (Minerva Biolabs).

For the measurement of IL-15 plasma levels, the Human IL-15 Quantikine ELISA Kit (R&D Systems) was used, following manufacturer’s recommendations. IL-15 plasma levels were measured in S1, S2, S3 and S4 patients´ samples. The optical density was determined using Varioskan Flash fluorimeter (Thermo Fisher Scientific) set to 450nm.

Flow cytometry data were analyzed using FlowLogic v.7.3 (Inivai Technologies) software. GraphPad Prism v.9 was used for graphical representation and statistical analysis. Data were tested for normal distribution with Shapiro-Wilk normality test. If data were normally distributed, Student´s t test for paired values was used to determine significant differences. Non-normal distributed data were compared with Wilcoxon matched-pairs signed rank test. Statistical analysis were done by comparing S2 with the rest of the samples. Data were represented as boxplots with the median and 25–75th percentiles, as violin plots with the median and quartiles or as bar plots as mentioned in each figure legend.

For bivariate analysis, normality of continuous variables was tested with Shapiro-Wilks test. Mean and standard deviation (SD) was presented when the variable followed a normal distribution or median interquartile range (IQR) otherwise. The Student´s t test or Mann-Whitney U test were used to compare continuous variables between groups. Qualitative variables were compared using the Chi-square or Fisher exact test. Correlation plots between variables were calculated and visualized as correlograms. Spearman´s Rank Correlation coefficient was indicated by square size and heat scale. Bivariate analysis was performed using Cox proportional hazard regression models. Survival curves were estimated with Kaplan-Meier method and compared by log-rank test. All the analysis were performed with R (version 4.0.4): A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria

Data from ELISA experiments was processed following manufacturer´s recommendations. GraphPad Prism software was used to fit the standard curve with nonlinear regression and log-log line model.

We first studied NK cell recovery of 54 MM patients after autoHSCT. Although not significant, a decrease in the absolute number of NK cells was noticed from S1 to S2 [(median S1: 108.2 cell/µl (range: 1 – 406 cell/µl) vs S2: 69.00 cell/µl (range: 4 - 806.5 cell/µl)]. Afterward, a significant increase was observed at S3, followed by a decrease at S4. The absolute number of NK cells was maintained at the same levels after that; although, the values were significantly higher than what was observed in S1 and S2 ( Figure 1A ). Thus, these results suggest that the recovery of NK cells is not completely achieved until 30 days after autoHSCT (S3).

We next studied the distribution of different NK cell subsets. NK cells were identified as viable CD3-/CD14-/CD19-/NKp80+ cells (11, 12). Then, based on the expression of the CD56 marker, two major NK cell subsets were defined: CD56^bright and CD56^dim ( Supplementary Figure 1A ). CD56^dim NK cells were the main subset at all time points. However, their frequency significantly decreased at S2, while the frequency of CD56^bright greatly increased. Nevertheless, this redistribution of NK cell subsets was not maintained over time ( Figure 1B ). We observed no differences in terms of the frequency of CD56^bright and CD56^dim NK cells at S2 between patients having a leukocyte recovery time period of <13 days vs ≥13 days ( Supplementary Figures 2A, B ). In addition, we studied the distribution of a third NK cell subset known as CD56^neg NK cells. We have recently described an expansion of this subset in MM patients (11). In the case of autoHSCT setting, the frequency of CD56^neg NK cells significantly dropped at S2 and showed variations over the studied period ( Figure 1B ). These results suggest that NK cell subset distribution is altered at leukocyte recovery, but resembles pre-transplant distribution 30 days after autoHSCT.

We next used CD57, NKG2C and FcϵRγ markers to study the different adaptive subpopulations within CD56^dim NK cells and within CD56^dimCD57+ NK cells ( Supplementary Figure 1B ). First, we determine the frequencies of FcϵRγ-, CD57+FcϵRγ-, NKG2C+ and CD57+NKG2C+ adaptive NK cell subpopulations. Results showed that adaptive NK cells lacking the expression of FcϵRγ were more abundant than the NKG2C+ ones ( Figures 1C, D ). In addition, the percentage of FcϵRγ- adaptive cells was reduced at S2 and maintained after that ( Figure 1C ). In contrast, no significant differences in the frequency of NKG2C+ NK cells were noticed from S1 to S3. However, an expansion of these subsets was observed at S4 and their frequency continued elevated thereafter ( Figure 1C ). Next, we analyzed the co-expression of NKG2C and FcϵRγ markers identifying four different populations. NKG2C-FcϵRγ+ population represents conventional NK cells and was therefore not included in the analysis. Regardless of CD57 expression, the main population at all time points were those with the NKG2C-FcϵRγ- phenotype. Moreover, the frequency of both CD57+NKG2C-FcϵRγ- and NKG2C-FcϵRγ- populations was reduced at S2 and maintained low after that ( Figure 1E ). In contrast, the frequency of NKG2C+FcϵRγ+ and CD57+NKG2C+FcϵRγ+ was not significantly altered immediately after autoHSCT, while the frequency of NKG2C+FcϵRγ- and CD57+NKG2C+FcϵRγ- was reduced at S2 and expanded at S4 ( Figure 1E ). These results suggest a different reconstitution pattern among adaptive NK cell populations after autoHSCT.

As the frequency of CD56^bright NK cells notably increases at leukocyte recovery, and these cells are considered more immature than CD56^dim NK cells (13), we next studied the NK cell differentiation and maturation status using NKG2A and CD57 markers. While CD57 has been described to be expressed by terminally differentiated NK cells, NKG2A is expressed in earlier differentiation stages (14–16). As others have previously reported (17), the frequency of NKG2A+ NK cells increased from S1 to S2 and, although decreased at S3, remained elevated until S4 ( Figure 2A , upper panel). In contrast, CD57+ NK cells behave the opposite way ( Figure 2A , lower panel). The variations observed in the frequency of NKG2A+ and CD57+ cells were mainly due to an increment in the frequency of CD56^dimNKG2A+ and CD56^negNKG2A+ NK cells, and a decrease in the frequency of CD56^dimCD57+ NK cells respectively ( Figure 2A ). Surprisingly, an expansion of CD56^bright NK cells expressing CD57 was observed at leukocyte recovery (S2) ( Figure 2A ). However, as they are immature NK cells, CD56^bright NK cells are not expected to express CD57 (16).

During the final steps of differentiation, NK cells lose the expression of NKG2A while acquire CD57 (15). Therefore, we analyzed the co-expression of both markers after the autoHSCT. We found that NKG2A+CD57- NK cells were the main subset at S2, but over time, the frequency of NKG2A+CD57+ NK cells gradually increased, turning out to be the predominant subset after day 180 post-autoHSCT (S5). Moreover, the percentage of NKG2A-CD57+ NK cells also gradually augmented after the initial decrease at S2 ( Figure 2B ). When CD56^bright, CD56^dim and CD56^neg NK cell subsets were separately analyzed we found that while within CD56^bright NK cells, NKG2A+CD57- cells were the most abundant subset at all time points, CD56^dim NK cells were mostly NKG2A+CD57- at leukocyte recovery and NKG2A+CD57+ after day 100 post-autoHSCT. In contrast, the majority of CD56^neg NK cells were NKG2A-CD57- at all time points ( Figure 2B ). Together, these data reveal that NK cells have a more immature phenotype at leukocyte recovery, as demonstrated by the increment in the frequency of NKG2A+CD57- cells, and that this phenotype last for more than 30 days after autoHSCT (% of NKG2A+ cells in S1 vs S3: p<0,0001; S1 vs S4: ns; S1 vs S5: ns; and S1 vs S6: ns; % of CD57+ cells in S1 vs S3: p = 0,0468; S1 vs S4: ns,; S1 vs S5: ns; and S1 vs S6: p = 0,0139; % of NKG2A+CD57- cells in S1 vs S3: p<0,0001; S1 vs S4: ns; S1 vs S5: ns; and S1 vs S6: p = 0,0072; % of NKG2A+CD57+ cells in S1 vs S3: ns; S1 vs S4: ns; S1 vs S5: ns; and S1 vs S6: p = 0,0088; % of NKG2A-CD57+ cells in S1 vs S3: p <0,0001; S1 vs S4: ns; S1 vs S5: ns; and S1 vs S6: ns; % of NKG2A-CD57- cells in S1 vs S3: p <0,0001; S1 vs S4: p = 0,0487; S1 vs S5: ns; and S1 vs S6: ns).

Next, we further analyzed the phenotype of NK cells during cell reconstitution after autoHSCT. For that, we studied the activation status and the cytotoxic potential of NK cells. The analysis of the expression of CD69, an activation marker (18), showed that the frequency of CD69+ NK cells slightly decreased at leukocyte recovery, and not significant differences were noticed after that. This decrease might be principally due to a lower percentage of CD56^bright NK cells expressing CD69 at S2 ( Figure 3A ). On the other hand, we evaluated the cytotoxic potential of NK cells after autoHSCT by assessing the expression of granzyme B. The expression of this marker did not differ on total NK cells from S1 to S2, but it was significantly incremented at S3. When the three NK cell subsets were separately analyzed, we observed that the levels of granzyme B on CD56^bright NK cells were greatly increased at S2, maintained elevated at S3 and gradually reduced thereafter. Conversely, CD56^dim NK cells expressed less granzyme B at S2 than at S1, but the expression increased again at S3 and it was higher than at pre-autoHSCT. In the case of CD56^neg NK cells, a reduced expression was also observed at S2, and the expression was gradually incremented after that ( Figure 3B ). These data suggest that the activation status and the cytotoxic potential is different among the three studied NK cell subsets during NK cell reconstitution. Moreover, these data showed that the cytotoxic potential of CD56^bright NK cells is incremented after autoHSCT.

Finally, we used the Ki67 nuclear protein as a marker to study NK cell proliferation (19). Results showed that there was a dramatic increase in the frequency of Ki67-expressing cells in total NK cells and in the three different NK cell subsets at leukocyte recovery that was not found at other time point ( Figure 3C ), indicating that NK cells exhibit a high proliferation rate immediately after autoHSCT while the majority of them enter into a quiescent state 30 days after autoHSCT.

Once phenotypic features of NK cells during their reconstitution after autoHSCT were analyzed, we next studied their functional capacities. For that, we examined degranulation and cytokine production by measuring CD107a and TNF, respectively, after stimulation with the 721.221 cell line and IFNγ production after IL-12+IL-15+IL-18 stimulation. Results revealed that, in comparison to CD56^dim NK cells, CD56^bright NK cells exhibited a trend towards an increase in their effector functions ( Supplementary Figure 3A ). No statistically significant differences were observed in the effector functions during NK cell reconstitution in any subset ( Supplementary Figure 3B ). Thus, as others have previously described (17, 20), these data suggest that the effector functions of NK cells were completely recovered early after autoHSCT.

IL-15 is a cytokine that plays an important role not only in proliferation, cytotoxic activity and cytokine production of NK cells, but also in their development and differentiation (21, 22). In addition, Porrata et al. showed that IL-15 levels at day 15 after autoHSCT affect survival of non-Hodking lymphoma (NHL) patients through NK cell recovery (23). Thus, we determined plasma levels of IL-15 at different time points (S1, S2, S3 and S4). Results showed a large increment of IL-15 levels at S2 ( Figure 4A ). Moreover, patients with ≤13 days between autoHSCT and S2 showed significantly higher levels of IL-15 at S2 than those with >13 days ( Figure 4B ). Nevertheless, no association was found between IL-15 levels at S2 and OS (overall survival) and PFS of these patients (data not shown).

Next, we performed correlation analysis and we found that IL-15 levels at S2 negatively correlated with absolute number of NK cells at S2. However, a positive association was observed at S2 between IL-15 levels and the percentage of Ki67+ NK cells and CD69+ NK cells. On the other hand, we found a negative correlation between IL-15 levels and the expression of granzyme B on NK cells at S2 ( Figure 4C and Supplementary Figures 4A–D ). Then, we categorized the patients in two groups based on their IL-15 plasma levels at S2 (median: 17.45pg/mL, range: 3.96-52.83pg/mL). We observed that patients that had higher IL-15 plasma levels showed significantly lower NK cell absolute numbers, significantly higher frequency of Ki67+ NK cells and expressed less granzyme B in NK cells at S2 ( Supplementary Figure 4E ). However, no significant differences were noticed in the percentage of CD69+ cells between groups ( Supplementary Figure 4E ). Therefore, these results suggest that IL-15 play a major role in NK cell proliferation capacity and also affect the cytotoxic potential and NK cell numbers after autoHSCT.

As NK cell count at 15 days post-autoHSCT has been defined as an independent predictor for survival after autoHSCT, we finally investigated whether a specific NK cell phenotype was associated with the clinical outcome of MM patients after autoHSCT. As shown in Table 1 , from the 54 patients included in this study, 29 patients experience disease progression (53.7%) and 9 patients died (16.7%) during the studied period (60 months). Thus, we performed bivariate analysis for PFS, OS and TTNT, and the phenotypic data generated in the present study, but we did not find any association (data not shown). Each phenotypic variable was then dichotomized into upper (>66^th percentile) and middle/low (≤66^th percentile) groups for further analysis. However, we did not obtain any significant result when bivariate analysis was done with dichotomized variables (data not shown).

Then, we carried out survival analysis with the dichotomized phenotypic variables and we observed statistically significant differences in terms of PFS when the frequency of terminally differentiated NKG2A-CD57+ NK cells at S3 and S4 was taken into account. In addition, we also observed differences in terms of TTNT when the frequency of NKG2A-CD57+ NK cells at S4 was analyzed. Specifically, we noticed higher PFS in patients in the middle/low percentile group of NKG2A-CD57+ NK cells at S3 and at S4: median PFS 44.4 vs 26.6 months in both S3 (HR=0.45; 95% CI=0.20-0.99; p=0.042), ( Figure 5A ) and S4 (HR=0.36; 95% CI=0.15-0.86; p=0.016) ( Figure 5B ). Furthermore, patients in the middle/low percentile group of NKG2A-CD57+ NK cells also showed a longer TTNT at S4; median TTNT 49.2 vs 37.2 months (HR=0.41; 95% CI=0.17-0.98; p=0.039) ( Figure 5C ). Therefore, these data suggest that the differentiation and maturation status of NK cells at +30 and +100 days after autoHSCT affects the clinical outcome of MM patients and that a lower frequency of mature and terminally differentiated NK cells is associated with protection against disease progression after autoHSCT in these patients.