NK cells were negatively isolated from fresh peripheral blood mononuclear cells (PBMC) derived from 10 different healthy donors and then the most immature CD56^bright/CD16⁻ (referred as CD56^bright) and the mature CD56^dim/CD16^bright (referred as CD56^dim) NK cell subsets were sorted by the expression of CD56 and CD16 and analyzed by flow cytometry for different markers in order to control their purity and phenotypic features (Figures 1A–C). As expected, CD56^bright NK cells were homogeneously CD16⁻, KIR⁻, NKG2A⁺, and CCR7⁺, whereas CD56^dim NK cells were CD16⁺, KIR⁺ and/or NKG2A⁺, and CCR7⁻. This analysis confirmed the high purity of the two sorted cell subpopulations. Importantly, in all instances, both these NK cell subsets did not express the activation marker CD69 (Figures 1B,C).

Then, miRNA expression profiles of human CD56^bright and CD56^dim NK cells subsets were investigated with the human miRNA microarray kit v19.0, which allows for the detection of a total of 2006 different human miRNA sequences (miRBase v19.0). In particular, 10 separate samples for each of CD56^bright and CD56^dim (20 total samples) were hybridized on the arrays. The entire dataset was composed of two separately generated microarray datasets that were analyzed together (Supplementary File 1).

The normalized batch-corrected log₂ intensity values were used to identify miRNAs differentially expressed between the CD56^bright and the CD56^dim NK cell populations (see also Materials and Methods). Considering an arbitrary threshold of 1 light unit (1LU), 251 miRNAs were over the threshold in half or more of the 10 CD56^bright samples, 198 of which were over the threshold in all 10 CD56^bright samples. Similarly, 262 miRNAs were over the threshold in half or more of the 10 CD56^dim samples, 213 of which were over the threshold in all 10 CD56^dim samples. Interestingly, 232 miRNAs were detectable in both NK cell subsets and only 49 were expressed in only CD56^bright (19 miRNAs) or CD56^dim (30 miRNAs).

Figure 2 shows the two-color heatmap plot as result of the unsupervised hierarchical clustering in which we performed a bi-clustering analysis of both miRNA and NK samples. This analysis clearly separates miRNAs differentially expressed and at the same time CD56^bright from CD56^dim NK cell subset. At the same time, we provided evidence that it is possible to separate the two main NK cell subsets by unsupervised hierarchical clustering. In particular, Figure 2 represents a two-color heatmap plot depicting the results of the bi-clustering analysis of both miRNA and NK samples. This analysis clearly separetes CD56^bright from CD56^dim NK cell subset. In particular, the heatmap identify a first level of signature characterized by 14 up-regulated and 23 down-regulated miRNAs in the CD56^dim NK cell population. Importantly, among the most regulated miRNAs, we found 5 miRNAs deeply and homogeneously down-regulated (hsa-miR-31-5p, hsa-miR-130a-3p, hsa-miR-223-3p, hsa-miR-146a-5p, and hsa-miR-92a-3p) and 2 miRNAs deeply and homogeneously up-regulated (hsa-miR-873-5p and hsa-miR-181a-2-3p) in all the CD56^dim samples compared to CD56^bright ones (see arrowheads in Figure 2). Notably hsa-miR-31-5p and hsa-miR-130a-3p expression was below detectable levels on CD56^dim NK cells, whereas hsa-miR-181a-2-3p expression was below detectable levels on CD56^bright NK cells. These miRNAs may be considered specific for CD56^bright or CD56^dim NK cells, respectively.

To better capture the interplay among the miRNAs in the two cell groups (i.e., CD56^dim and CD56^bright), in addition to the univariate analysis shown in Figure 2, we performed a multivariate analysis by using the l1l2-penalized regularization method, a machine learning technique for variable selection. A schematic drawing in Supplementary Figure 2A shows the approach used in such analysis (see also Materials and Methods).

The list obtained from the multivariate analysis includes all the miRNAs found with the univariate analysis (37 miRNAs, Figure 2) and 71 additional miRNAs identified only by this method. Figure 3 shows the heatmap of the 108 miRNAs signature identified (Figure 3A) and the names of the miRNAs listed in the order used in the heatmap (Figure 3B). The result was a miRNA signature associated with a perfect prediction accuracy (100%, corresponding to a zero error score).

We selected 7 of the most regulated miRNAs (indicated in Figures 2, 3) that were differently expressed in CD56^bright and CD56^dim NK cell subsets and we validated them by quantitative Real Time Polymerase Chain Reaction (RT-PCR). The purity of the two NK cell subsets was confirmed by the gene expression analysis of CCR7 (a marker selectively expressed in CD56^bright cells subset but not expressed in CD56^dim NK cell subset as shown in Figure 1C). These analyzes confirmed once again the levels of expression of these miRNAs in the two NK cell subsets (Figure 4A).

Then we wondered if the miRNA signature found in the two main NK cell subsets (CD56^bright and CD56^dim) could be linked to the sequential steps of maturation. Thus, we sorted the additional NK cell subset CD56^bright CD16⁺ (Supplementary Figure 1), that is known to represent an intermediate step of peripheral blood-NK cell maturation, and performed qPCR analysis by looking at the most differentially expressed miRNAs (hsa-miR-130a-3p, hsa-miR-31-5p, hsa-miR-181a-2-3p, hsa-miR-223-3p) (Figure 4B). We showed that the expression levels of these miRNAs on the CD56^bright/CD16⁺ NK cell subset were halfway between those observed on the CD56^bright/CD16⁻ and those on the CD56^dim/CD16⁺ NK cell subsets. In particular, by comparing the qRT-PCR data obtained from CD56^bright/CD16⁻, CD56bright/CD16⁺ and CD56^dim/CD16⁺ NK cells subsets, the existence of a significant progressive change in these miRNAs expression was revealed: this detail represents the trend of NK cell subsets maturation that shows how the expression reaches different levels as the cells goes on maturing. (Figures 4A,B).

In addition, we also sorted the CD56^dim/CD57⁺ and the CD56^dim/CD57⁻ NK cell subsets, recognized as the most mature NK cell subsets. We found that hsa-miR-130 and hsa-miR-31-5p as well as hsa-miR-181a-2-3p were expressed at similar levels in CD56^dim/CD57⁺and CD56^dul/CD57⁻ NK cell subsets (data not shown). This supports the hypothesis that such miRNAs are important molecules for the switch between CD56^bright and CD56^dim phenotypes but not for the terminal differentiation that ends with a CD56^dim/CD57⁺ phenotype.

The miRNA multivariate analysis, as reported in the sections above, identified a signature of 108 predictive miRNAs. In order to validate data derived by multivariate analysis, we selected 12 miRNAs composed of 6 miRNAs validated by univariate analysis (hsa-miR-223-3p, hsa-miR-92a-3p, hsa-miR-31-5p, hsa-miR-181a-2-3p, hsa-miR-130a-3p, hsa-miR-146a-5p) and additional 6 miRNAs identified by multivariate analysis which showed a two-fold change listed in Supplementary Table 1: hsa-miR-215-5p, hsa-miR-513a-5p, hsa-miR-532-5p, hsa-miR-3652, miR-6723-5p, hsa-miR-181a-3p. This filtering procedure was undoubtedly necessary to achieve a more compact and interpretable signature, but it induced an additional verification phase on the multivariate model. In fact, before proceeding with the validation, we had to make sure that a multivariate model based only on the 12 miRNAs of the filtered signature was still capable of a good prediction performance in discriminating CD56^bright and CD56^dim. Therefore, we considered the original dataset limited to the 12 selected miRNAs, obtaining a shortened dataset of 20 samples and 12 variables. Using regularized least squares in a Monte-Carlo Cross Validation setting as in Barbieri et al. (37), we estimated a new predictive model that was able to achieve a 99.7% average prediction accuracy and associated with a p-value as low as 1.020e-20 from the two samples Kolmogorov-Smirnov test (Supplementary Figure 2B). This result showed that the 12 miRNAs model is an excellent classifier of CD56^bright and CD56^dim when considering miRNAs expressions measured with the microarray technology.

We then proceeded with a further validation step with the aim of verifying if such model was still capable of good prediction on new samples measured with a different technology, i.e., qPCR.

This is a very strong requirement on the generalization properties of the model, as the test set was completely independent from the learning set, having been acquired with a different technique and on different samples. To this aim, the expression of the 12 miRNAs was measured on 6 new samples (Supplementary Figure 3). We used the same statistical classifier to predict the status (CD56^bright or CD56^dim) of the 6 samples measured with qPCR obtaining a correct prediction of 4 out of 6 samples. The entire validation procedure is depicted in detail in Supplementary Figure 2B.

Since miRNAs function mainly through the inhibition of target genes, we next tried to identify the mRNA targets regulated by miRNAs differentially expressed in human NK cell subsets that might be involved in their function/maturation.

We considered six prediction tools for mRNA target genes prediction: TargetScan (38), Miranda (39), MirDB (40), Eimmo (41), Pictar (42), and PITA (31). Each tool was used to predict the target genes regulated by miRNAs that belongs to Supplementary Table 1, which contains 37 selected miRNAs up-regulated on CD56^dim and 33 miRNAs up-regulated on CD56^bright, on the basis of multivariate analysis. A complete list of the potential target genes predicted by bioinformatics platforms is shown in Supplementary File 2 both for CD56^dim and CD56^bright.

Since one of the goals of this study was the identification of novel miRNAs that may contribute to the NK cell development and differentiation, we focused our attention on potential target genes encoding for markers that are differently expressed on the two main NK cell subsets. These analyses identified some miRNAs differently expressed on CD56^bright compared to CD56^dim NK cells as potential regulator of targets genes encoding molecules of our interest (Table 1).

In particular, our analysis identified hsa-miR-146a-5p as a possible regulator of targets genes encoding for different members of the most important HLA-I specific NK receptors, the KIR family (12). Notably, hsa-miR-146a-5p seems to be mainly involved in the regulation of inhibitory KIR expression (with the exception of KIR3DL2/3), while other miRNAs appear to be involved in the regulation of activating KIRs. The only exception is represented by KIR3DS1 that again appears to be regulated by hsa-miR-146a-5p (Table 1). In addition, one of the predicted target of hsa-miR-146a-5p is perforin mRNA.

It is important to remember that KIR receptors are selectively expressed by the CD56^dim NK cell population and absent on the CD56^bright cell population [with the exception of KIR2DL4 (15)] and that their expression increases according to the maturation state of NK cells. Furthermore, perforin expression is higher in CD56^dim NK cell compared to CD56^bright cell. In this context, the up-regulation of hsa-miR-146a-5p detected on CD56^bright NK cells may be in line with the down-regulation of KIRs and perforin expression on this subset.

To demonstrate that the KIR genes may represent targets for hsa-miR-146a-5p, we performed additional experiments focusing our attention on genes encoding the inhibitory KIR2DL1 and KIR2DL2 receptors. These KIR receptors exhibit complementary features to enable recognition of both HLA-C1 and HLA-C2 subtypes, the dominant KIR ligands.

By using different target prediction software we predicted that hsa-miR-146a-5p has conserved binding sites in the 3′UTR of KIR2DL1 (Figure 5A) and KIR2DL2 (Figure 5B). To functionally verify the putative miRNA-mRNA interaction, we used plasmids containing the entire 3′UTR of both KIR2DL1 and KIR2DL2 and luciferase and Renilla reporter genes. COS-7 cells were transfected with these plasmids with hsa-miR-146a-5p or with negative control (miR-NC). Luciferase results revealed significant inhibition of luciferase activity for hsa-miR-146a-5p with both 3′UTR sequences of KIR2DL1 (Figure 5A) and KIR2DL2 (Figure 5B). Thus, we validated KIR2DL1 and KIR2DL2 as potential hsa-miR-146a-5p targets by luciferase reporter assay.

In order to explore which pathways characterize the predicted gene targets, we performed a in silico functional characterization using the webtoolkit WebGestalt (43) (see Materials and Methods). In particular, we made the enrichment analysis considering the Kyoto Encyclopedia of Genes and Genomes (KEGG) (44).

In details we found 195 common KEGG pathways, 5 pathways specific for CD56^bright and 2 pathways specific for CD56^dim (Supplementary File 3). In Figure 6 we reported a selection of the common KEGG pathways that we considered majorly linked to phenotypic and functional features of NK cells indicating which of the validated miRNAs appear to be directly involved.

Among these KEGG-Pathways, we show, in particular, the “Natural Killer cell mediated cytotoxicity” pathway (Figure 7), where multiple putative targets of hsa-miR-146a-5p are involved, such as CD94, HLA-C, HLA-E, Perforin, and KIRs (including KIR2DL1/L2 that we validated as targets for this miRNA).

Buffy coats from healthy donors were obtained from the Immunohematology and Transfusion Center at the IRCCS San Martino Polyclinical Hospital, Genoa, Italy.

All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Ethical committee of IRCCS San Martino Polyclinical Hospital, Genoa, Italy (39/2012).

Mononuclear cells from heparinized peripheral blood were obtained by density gradient centrifugation over Ficoll (Sigma, St. Louis, MO). Pure populations of NK cells were obtained from PBMC using the NK cell isolation kit (Miltenyi Biotec, Bergisch Gladbach, Germany) according to the manufacturer's instruction. In some experiments, MACS CD15 Micro beads were added to further improve depletion of granulocytes. The purity of NK cells was >98% NK cells (defined as CD56⁺/CD3⁻ CD19⁻/CD14⁻).

Freshly purified NK cells were then analyzed on BD-FACSAria II (BD-Biosciences) and sorted based on the expression of CD56 and CD16 molecules to collect the CD56^bright/CD16⁻, CD56^bright/CD16^dim, and the CD56^dim/CD16⁺ NK cell subsets separately (Figure 1 and Supplementary Figure 1). After each sorting, a complete phenotypic analysis of the populations collected was performed using the BD-FACSVerse (BD-Biosciences). Sorted cells were washed in PBS and the pellet freeze-dried (20 s in liquid nitrogen and then stored at −80°C).

The following purchased mAbs were used in this study: Anti-CD56-PC7 (clone c218), and anti-NKG2A allophycocyanin (Z199 clone) were purchased from Beckman Coulter/Immunotech (Marseille, France); anti-CD3-Viogreen (BW264/56 clone), anti-CD19-VioGreen (LT20 clone), anti-CD14-Viogreen (TÜK4 clone), and anti–KIR2DL1-S1 FITC (11PB6 clone), mAbs were purchased from Miltenyi Biotec (Bergisch Gladbach, Germany); anti-CD16–PerCP-Cy5.5 (3G8 clone), anti–KIR2DL2/L3/S2–FITC (CH-L clone) from BD Bioscience; Anti-hCCR7 (IgG2a, clone 150503) was from R&D Systems Inc (Abingdon, United Kingdom); anti-CD69 (FN50 clone, IgG1) were purchased from BioLegend (San Diego, Calif). The anti-KIR3DL1/L2/S1 mAb (AZ158 clone) was isolated in our laboratory and the specificity has been validated in the corresponding patent (WO2010081890A1). Cytofluorimetric analyses were performed on a FACSVerse (Becton Dickinson, Mountain View, Calif), and data were analyzed with FACSuite software version 1.0.3.

Total RNA including miRNAs was isolated using the miRNeasy Mini Kit (Qiagen). RNA samples were quality-checked via the Agilent 2100 Bioanalyzer platform (Agilent Technologies). RNA integrity (RIN) was determined by means of capillary electrophoresis by using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA). The threshold for RNA quality suitable for microarray analysis was considered to be a RIN>6.0 (67). All RNA samples used in the microarray assays and real-time PCR revealed RIN values between 7.6 and 9.4.

Sample labeling and hybridization was performed according to the Agilent miRNA Complete Labeling and Hyb Kit (Agilent Technologies Santa Clara, USA) following the manufacturer's protocol. Briefly, 100 ng of total RNA was labeled using T4 RNA ligase, incorporating Cyanine 3-Cytidine biphosphate (pCp). The Cyanine-3-labeled miRNA samples were prepared for One-Color based hybridization with complete miRNA Labeling and Hyb Kit (Agilent Technologies, Santa Clara, USA) according to the manufacturer's instructions. Labeled miRNA samples were hybridized overnight (20 h, 55°C) to Agilent Human microRNA Microarrays 8x60K v19 using Agilent's recommended hybridization chamber and oven. Afterwards, microarrays were washed with increasing stringency using Gene Expression Wash Buffers (Agilent Technologies, Santa Clara, USA). Fluorescent signal intensities were detected with Agilent's Microarray Scanner System (Agilent Technologies) and extracted from the images using Feature Extraction 10.7.3.1 Software (Agilent Technologies, Santa Clara, USA). For determination of differential miRNA expression FES derived output data files were further analyzed using the GeneSpringGX software (Agilent Technologies).

The experimental design of this study contains groups of matched sample pairs [CD56^bright and CD56^dim] obtained from 10 individuals. Therefore, statistical tests were performed to identify miRNAs differentially expressed between the sample groups.

The intensity data of the 10 individual microarrays were subjected to probe summarization, thresholding (1.0), and log₂-transformation using GeneSpring v12.6 (Agilent Technologies). Subsequently, the arrays were normalized between each other by quantile normalization.

The miRNA data sets were normalized by quantile normalization and base-2 logarithms (log₂) intensities were calculated. The variability between samples was visualized by centering the normalized data relative to the median of all samples for each miRNA. These centered data were calculated by subtracting the median of all samples from expression value of each individual. The recentered values are then represented in a two-color heatmap format on a log₂ scale (see also the heatmap-coded columns median centered log₂ intensities in the Supplementary File 1). The transformed log₂ intensity values were further adjusted to reduce a batch effect. The corresponding batch-adjusted data are provided in the table “normalized log₂ data bc” within the Supplementary File 1.

These normalized batch-corrected log₂ intensity values were used in subsequent analyses to identify miRNAs differentially expressed between the CD56^bright and the CD56^dim NK cell populations. For the identification of differentially expressed miRNAs, the robustness of detection and the statistical significance were taken into account. With regard to statistical significance, we applied a standard selection criteria based on the following criteria: in the comparison of two groups, a miRNA was classified as induced if its Adj. p-value was ≤ 0.05 with difference in the median expression level (fold-change) of at least two-fold; a miRNA was considered repressed if its Adj. p-value was ≤ 0.05 with a fold-change value ≤ - 2.0 (Supplementary File 2).

The unsupervised hierarchical clustering was performed by using the online tool MORPHEUS (https://software.broadinstitute.org/morpheus/). Specifically, Figure 2 shows the results of the bi-clustering analysis of both miRNAs and NK cell samples.

The data presented here have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) with the GEO accession number GSE116743 [Internet address: (68)].

In supervised multivariate statistics, the goal is to model the relationship between input data (samples) and their corresponding labels (cell type, bright, or dim) by looking for a function depending on several variables (miRNAs) at a time. Such a function should be able to: (i) have good prediction performance on the given data, (ii) generalize on previously unseen data, and (iii) effectively describe the interplay between the measured variables. The final outcome of a classifier is the prediction of labels associated to a set of input examples. To assess a prediction performance, one should define appropriate metric, such as the accuracy score. Accuracy is defined as the ratio of correctly predicted labels.

In this context, regularization methods are a popular class of machine learning techniques that can be expressed as the minimization problem of one loss function V that measures the adherence of the objective function to the data and one or more regularization penalites R that introduce additional information used to solve the problem: min_f V(f(X),y) +; λR(f).

Different choices of V and R lead to different algorithms. The parameter λ controls the trade-off between the adherence of the model to the data and the regularity of the function f. Choosing the regularization parameter appropriately leads to unbiased statistical models (69). In the remainder of this section we will illustrate the two regularization methods used in the present work. The software library we used is implemented in Python and publicly available at https://github.com/slipguru/palladio.

The microarray data was validated by real-time PCR according to the manufacturer's instructions (Qiagen; miScript Reverse Transcription kit, catalog no. 218060; miScript Primer Assays for the following miRNAs: hsa-miR-92a-3p, catalog no. MS00006594, hsa-miR-146a-5p catalog no. MS00003535, hsa-miR-133b catalog no. MS00031430, hsa-miR-130a-3p catalog no. MS00003444, has-miR-31-5p catalog no. MS00003290, has-miR-223-3p catalog no. MS00003871, hsa-miR-181a-2-3p catalog no. MS00008834; hsa-miR-873-5p catalog no. MS00010605; hsa-miR-181a-3p catalog no. MS00006692; hsa-miR-215-5p catalog no. MS00003829; hsa-miR-513a-5p catalog no. MS00009912; hsa-miR-532-5p catalog no. MS00004571; hsa-miR-3652 catalog no. MS00023128; hsa-miR-6723-5p catalog no. MS00045794; miScript SYBR Green PCR kit, catalog no. 218073). Then, miRNA was normalized to the small nucleolar RNA RNU6B (primer assay, catalog no. MS00014000). The real-time PCR reactions were performed for target miRNAs and for the internal reference (U6), with each sample analyzed in triplicate. For qPCR data analysis, the ΔΔCt method, as previously described by Schmittgen and Livak (73), was used to calculate the relative expression levels of the target miRNAs. For cDNA synthesis we used Qiagen RT² First Strand Kit (catalog no. 330401). RT² PCR Assay Gene Expression was performed for CCR7 chemokine receptor 7 (NM_001838). The GAPDH was used for normalization in the fold change expression data calculations. After qPCR relative expression is determined with the ΔΔCt method as previously described.

Several published target prediction algorithms have been used for the prediction of miRNA targets. Generally, such software mainly use sequence complementarity, evolutionary conservation among different species, and thermodynamic criteria to estimate the probability of a miRNA:mRNA duplex formation. For further reading, the review by Bartel published in 2009 (30) is recommended (http://www.ncbi.nlm.nih.gov/pubmed/19167326). For gene target prediction, we considered six prediction tools: miRanda [http://www.microrna.org (39)], mirDB [http://mirdb.org/miRDB/index.html (40)], Pictar [http://pictar.mdc-berlin.de/ [42]], PITA [http://genie.weizmann.ac.il/pubs/mir07 (31)]. TargetScan [http://targetscan.org (74)], ElMMo [http://www.mirz.unibas.ch/ElMMo2 (41)]. Concerning TargetScan we used the Predicted Targets (default predictions) file belonging to the 7.0 release. Concerning PITA we considered the PITA_sites_hg18_0_0_ALL file of the last release (number 6). Concerning PicTar we downloaded and used all the targets from doRiNA database (75), as suggested in the home page of this prediction tool. Concerning miRanda we considered the predictions from the August 2010 release. Specifically we used the most completed file that includes the “Non-good mirSVR Score and Non-conserved miRNAs.” Concerning mirDB we considered the last version (5.0). Concerning EImmo we used the last release, number 3, updated in January 2009.

For this purpose, we used the pEZX-MT06 target reporter vectors containing the entire 3'UTR sequence for KIR2DL1 and KIR2DL2 (GenBank Accession: NM_014218.2, HmiT010083-MT06 from GeneCopoeia; GenBank Accession: NM_014219.2, HmiT010084-MT06 from GeneCopoeia) containing putative miR-146a-5p binding sites and transiently transfected them into COS-7 cells by using the Lipofectamin 2000 reagents according to the manufacturer recommendations (Invitrogen). COS-7 cells were grown up in 24-well plates to a confluency of about 80%, followed by a co-transfection of hsa-miR-146a-5p mimic (two different concentrations: 36.4 nM and 54 nM, Qiagen) together with 150 ng of each pEZX-MT06 constructs (HmiT010083-MT06 and HmiT010084-MT06). The luciferase activity was compared to values obtained after we co-transfected the negative control (miR-NC, Allstars negative control, Qiagen) and the plasmids (HmiT010083-MT06 and HmiT010084-MT06). Forty-eight hours after transfection, cells were lysed in 100 μl of passive lysis buffer according to the Luc-Pair Luciferase Assay Kit 2.0 (GeneCopoeia); 20 μl of the lysate was used for the luciferase activity measurements following the instructions of the Dual-Luciferase Reporter Assay System (Promega). Luciferase assays were run on a LUMIstar Luminometer (BMG Labtech) in three independent biological replicates.

For the functional analysis of the miRNA signature we used the on-line gene set analysis toolkit WebGestalt (43)¹WEB-based GEne SeT AnaLysis Toolkit (WebGestalt): update 2013. Nucleic Acids Res, 41 (Web Server issue, W77-83). The toolkit performs the functional characterization by a gene set enrichment analysis in several databases including KEGG (44). Given a KEGG pathway and a reference set (such as the entire human genome) the enrichment is based on the comparison between the fraction of signature genes in the pathway and the fraction of pathway genes in the reference set. The signature is enriched in the KEGG pathway if the former is larger than the latter fraction. To perform the enrichment analysis in KEGG, we selected the WebGestalt human genome as reference set, p ≤ 0.05 as level of significance, 3 as the minimum number of genes, the Bonferroni correction to correct for multiple hypotheses and the default Hypergeometric test as statistical method.

One of the specific KEGG pathways we found enriched (NK cell mediated cytotoxicity, hsa04650) was colored using the “Color Pathway” tool available in KEGG Mapper that is a suite of mapping tools provided by KEGG (KEGG Copyright Permission 180201).

Statistical analysis was performed with Graphpad Prism (La Jolla, CA) software. For statistical analysis of cytofluorimetric experiments (Figure 1C) and of RT-PCR data on CD56^bright/CD16⁻, and CD56^dim/CD16⁺ NK cells subsets (Figure 4A and Supplementary Figure 3) was used non-parametric t-test (Mann Whitney test). P-value of < 0.05 (^*), < 0.01 (^**), < 0.001 (^***), and < 0.0001 (^****) was considered statistically significant, when not indicated, data were not statistically significant. Multiple comparisons between RT-PCR data on CD56^bright/CD16⁻, CD56^bright/CD16⁺, and CD56^dimCD16⁺ NK cells subsets were performed using univariate analysis of variance (One way Anova with Bonferroni's post-test; Figure 4B). Cytofluorimetric experiments are reported as the sample mean ± the standard deviation (SD) (Figure 1C) while RT-PCR data are reported as the sample mean ± the standard error of the mean (SEM) (Figures 4A,B and Supplementary Figure 3).

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