We consecutively selected all patients who had received an HLA-identical sibling-donor allo-HSCT in our center between 2000 and 2015 and for whom we had DNA samples available for NGS analysis. GVHD prophylaxis included cyclosporine A and methotrexate for all patients.

Those patients who died or relapsed without having developed GVHD before day +100 for aGVHD (n=2) or day +180 for cGVHD (n=28) were censored. Those patients who died due to causes not related to transplant toxicity were also censored (n=20), and in our analysis of relapse, patients who died during the first year after transplantation, without having relapsed, were also excluded (n=13).

Clinical variables considered were donor and recipient sex, donor and recipient age, underlying disease, stem-cell source, conditioning regimen, prior radiation therapy and previous HSCT ( Table 1 ).

The local ethics committee approved the study, and all recipients and donors provided written informed consent in accordance with the Declaration of Helsinki.

Post-transplant complications analyzed were grade II-IV aGVHD, III-IV aGVHD, cGVHD, moderate-severe cGVHD, TRM, relapse and OS.

GVHD classification and clinical data collection were performed at the moment of GVHD diagnosis by the attending physician following the 1994 Consensus Conference on aGVHD grading and the National Institutes of Health criteria for diagnosis and staging of cGVHD.

In the case of TRM, we have considered TRM every death that occurs while the patient is in remission, any death not attributable to relapse.

Pre-transplant samples were selected, mostly peripheral blood. The selected samples were in complete remission as the objective was to identify germline variants (polymorphism). Subsequently, genomic DNA was purified automatically (Maxwell^® 16 Blood DNA Purification Kit; Promega, Madison, Wisconsin) following manufacturer instructions. The DNA extracted was frozen at -80°C in our biobank (ISCIII N°C.0000915).

We designed a custom panel of 132 cytokine genes that included 73 interleukin and 59 chemokine genes ( Supplementary Table 1 ). The probes were designed to detect coding and non-coding regions, namely untranslated regions (UTR), splicing (± 1,2 base pairs), and upstream and downstream regions (± 200 base pairs).

Libraries were prepared using a capture gene panel according to the manufacturer’s protocol (Agilent, Santa Clara, California). Paired-end sequencing (2×101 bp) was performed using the Illumina HiSeq platform (Illumina, San Diego, California). FASTQ files were aligned against the human reference genome (GRCh38/hg38 version) using the Burrows Wheeler Alignment tool v0.7.15-r1140. Variant calling and indel-realignment were performed using a combination of two different algorithms: GATK and VarScan. GeneSystems software (Sistemas Genómicos, Valencia, Spain) was used for variant annotation in order to provide the infrastructure and interface for bioinformatic analysis. Identified variants were annotated using the Ensembl database, population databases (the Exome Aggregation Consortium and 1000 Genomes), and specific variant databases (ClinVar, Catalogue of Somatic Mutations in Cancer (COSMIC). The population bases GenomAD and 1000 Genomes were used to determine population frequency (Minor Allele Frequency, MAF).

All polymorphisms located in coding sequences, splicing (± 1,2 base pairs), UTR and upstream and downstream regions (± 200 base pairs) were analyzed. Non-synonymous variants with ≥30X depth in the canonical isoform, a variant allele frequency (VAF) greater than 40%, a MAF ≥1% and representing at least 5% of our cohort were selected. From the total number of polymorphisms (genetic variants) analyzed, we selected those in which the presence of the minor variant (allele) is associated with the development of a complication.

Given the high dimensionality of the data and the limited sample size, conventional methods were not flexible enough to find discriminative features. To address this, we proposed Bayesian probabilistic methods to manage uncertainty while performing dimensionality reduction.

Data analysis consisted of two main steps: feature selection and classification or prediction. To address the first step, we developed a novel probabilistic variant of the Logistic Regressor, incorporating L1 regularization, referred to as the Bayesian Logistic Regression (BLR) model. The proposed BLR model follows a Variational Autoencoder (VAE) approach, in which the weight vector is approximated through a multivariate Gaussian distribution. Additional details of the BLR model and VAE approach can be found in the Supplementary Material .

In the classification step, we wanted to account for the potential non-linear relationships between the selected input variables and the complications. As the BLR is limited in its ability to capture non-linear relationships, it was only utilized for the feature selection stage. As an alternative, a probabilistic non-parametric kernel model, specifically a Gaussian Process (GP) classifier, was used to generate predictions. Additional details can be found in the Supplementary Material .

Taking into account the number of variables, BLR models were used to select the most important variables involved in each complication. Initially, the 15 most relevant genetic/clinical variables were selected. The strategy implemented was to select those predictive models that obtained the highest area under the curve (AUC), i.e. the greatest predictive power, with the lowest number of variables ( Supplementary Table 3 ). Depending on the presence of the variables selected for each model, a score was calculated for each patient. Subsequently, taking into account sensitivity, specificity, false positives, and false negatives (in both cases setting the limit at 15%), the best cutoff point was selected for the subsequent stratification of patients into high or low risk according to the model selected for each post-transplant complication ( Supplementary Table 4 ).

The selected variants were represented in Tables 2 – 8 , in these tables display the following characteristics:

An enrichment analysis of the genes included in the models ultimately selected was conducted using the Enrichr web-based software application in order to establish common functions and establish a functional relationship between them.

Of the global cohort, 51/88 (58%) of patients presented grade II-IV and 27/88 (30.6%) III-IV aGVHD. Additionally, 39/63 (62%) and 25/61 (41%) developed cGVHD and moderate-severe GVHD respectively. Regarding TRM, 22/70 (31,4%) patients presented. Finally, 17/81 (23%) and 44/90 (49%) patients relapse or died respectively.

The cumulative incidence (CI) of grade II-IV and III-IV aGVHD at 100 days post-transplant was 56.7% and 23.3%, respectively. Cumulative incidence of cGVHD and moderate-severe cGVHD (mod-sev cGVHD) at three years post-transplant was 43.3% and 27.8%, respectively. Three-year CI for TRM was 24.4%, and 2-year OS was 60%.

The median times and the survival plot of each complication are shown in Supplementary Figure 1 :

Using the filters previously defined in Material and Methods, 737 polymorphisms in the 119 genes were detected in 90 donor-recipient pairs ( Supplementary Table 2 , Figure 1 ). Of these polymorphisms, 93 (12.6%) were located in coding regions, and 644 (87.4%) in non-coding regions: 3’ and 5’ UTR, splicing, upstream and downstream regions.

Different predictive models, including clinical variables and 737 genetic variants, were built to anticipate the development of post-transplant complications such as GVHD, TRM, relapse and OS following allo-HSCT.

Initially, the 15 most relevant genetic/clinical variables were selected to reduce the dimensionality of the data ( Supplementary Table 3 ). Of these initial 15 variables, the final variables for predictive models were selected taking into account the AUC and the number of variables. The models ultimately obtained for each of the post-transplant complications (aGVHD II-IV, aGVHD III-IV, cGVHD, mod-sev cGVHD, TRM, relapse and OS) are shown in Tables 2 – 8 . Finally, 41 polymorphisms in 30 cytokine genes were selected (17 interleukins and 13 chemokines). Of these polymorphisms, 5 (12.2%) were located in coding regions, and 36 (87.8%) in non-coding regions.

Based on the sensitivity and specificity data, as well as the ratio of false positives and false negatives (a limit of 15% was established), cut-off points were selected to classify the patients at high or low-risk for the development of each post-transplant complication ( Supplementary Table 4 ). Finally, the cumulative incidence of each complication according to patient risk stratification ( Figure 2 ).

Considering the genes included in the predictive models for post-transplant complications, 41 polymorphisms in 30 cytokines. Four signaling pathways functionally associated with the genes were obtained by enrichment analysis using web-based Enrichr software ( Table 9 ) with a p-value of 2.5*10^-56, 1.48*10^-27, 4.04*10^-18 and 4.97*10^-9, respectively.

All genes selected were included in the cytokine signaling pathway. Another signaling pathway selected by Enrichr was JAK/STAT, in which 13/30 cytokines, specifically interleukins (43.3%), were included. Finally, the IL-17/TH17 signaling pathway consisted of 6/30 cytokines (20%), including five interleukins and one chemokine.

GVHD models Interleukins related to the activation, proliferation and differentiation of immune system cells involved in the development of GVHD were included in acute and chronic GVHD models. In addition, different chemokines responsible for driving these cells to the sites of inflammation were included ( Figure 3A ). IL-12RB1 is a CD4⁺ T cell receptor that, depending on the binding of the IL-12 or IL-23 ligand, leads to a differentiation to a Th1 or Th17 phenotype (17), characteristic of aGVHD and cGVHD, respectively, that was selected for the aGVHD and cGVHD models. The IL-2RA receptor, which is responsible for the proliferation of T cells involved in GVHD and which in soluble form has been reported to be increased in stages of immune activation (18), was also included in these models. Specifically, the aGVHD model included interleukins such as IL-17A, IL-17RA and IL-26 related to the synthesis of proinflammatory cytokines (19, 20), which may be involved in the development of this complication. CCL25, although described in inflammatory diseases, is also associated with the proliferation of regulatory T cells (Tregs) that inhibit GVHD (21, 22). IL-21R has been reported as responsible for the proliferation and differentiation of lymphoid cells, but there are already studies that demonstrate that blocking IL-21/IL-21R in the gut increases Tregs, thereby inhibiting GVHD (23). Regarding cGVHD, it is worth noting that the chemokines selected in the predictive models have already been described in the literature as related to the development of GVHD. The interleukins IL-7/IL-7R, related to survival, proliferation and activity of T and B cells (24); as well as IL-25, which is related to TH2 activation but inhibits TH17 (25), T cells related to cGVHD, were also included ( Figure 3B ).

Notably, a greater number of genes were selected in the more severe GVHD models. This may be explained by the fact these could be expected to involve a more inflammatory context affecting various organs and thus producing a greater severity of GVHD. This is consistent with the selected models for severe GVHD obtained in a previous study (8).

Regarding the genetic model for TRM, all interleukins selected signaled through the JAK/STAT pathway, leading to the synthesis of proinflammatory cytokines. These cytokines could lead to an inflammatory environment that can affect TRM due to tissue damage. This is further enhanced by the chemokines selected for this model, which are related to the migration of immune system cells to inflammatory sites.

For the relapse model, three polymorphisms were selected in a single interleukin, IL-21R, which through its signaling leads to the proliferation and differentiation of lymphoid cells (26, 27). Thus, it could positively or negatively regulate immune responses. It worth noting that mainly chemokines were selected in this model (7/10). This could be explained as proposed in a study of pediatric patients diagnosed with acute lymphoblastic leukemia where the expression of chemokines and their receptors is increased in relapsed patients. In this study the authors propose that the cells secrete chemokines that recruit those leukocytes that express the specific receptor. Thus, tumor cells that secrete and/or respond to chemokines would have a selective advantage and may show resistance to chemokine therapy (28).

Finally, regarding the OS model, IL-10RB signaling exerts anti-inflammatory functions (29); hence we must identify the effect of the polymorphism selected on the expression of this gene. In turn, signaling of the other interleukin included, IL-12RB1, leads to the synthesis of IFN-γ, which may have an antitumor role in the activation of effector T cells, which will destroy tumor cells. However, it may also have a protumor role, involving proliferative and antiapoptotic signals (30). The chemokines selected are related to the migration and differentiation processes of different cells in the immune system. It is likely that these chemokines are able to recruit different types of immune cells in the tumor which, in turn, can modulate tumor growth and metastasis, as has been demonstrated in the literature (31).

An analysis of enrichment with the Enrichr application (32) found that in addition to the cytokine and chemokine signaling pathway, as expected, two signaling pathways represented in the inflammatory context were also selected: the IL17/TH17 and JAK/STAT signaling pathways. TH17 cells are CD4⁺ T cells that produce proinflammatory cytokines such as IL-17A, IL-17F, IL-21, IL-22, TNF-α, G-CSF and some chemokines. IL-17 is related to innate immunity and inflammation and relates T cell activation to neutrophil mobilization and activation (33). Several studies provide evidence of the role of IL-17 in GVHD, demonstrating that activated TH1 and TH17 lymphocytes secrete proinflammatory cytokines that lead to apoptosis of cells in target tissues in GVHD, primarily in the intestine, liver and skin (34, 35). The JAK/STAT signaling pathway regulates the activation of immune cells related to GVHD, including APCs, T cells, neutrophils and B cells (36, 37). It is therefore involved in regulating cell activation, proliferation, migration and cytokine production, increasing the severity of GVHD. In this context, drugs such as ruxolitinib that inhibit this signaling pathway are already available and used for the treatment of GVHD (38).

Based on the results obtained, it should be taken into account that most of the genetic variables selected for the predictive models of each of the post-transplant complications were located in non-coding regions, and that most of these polymorphisms identified for their clinical relevance (association with post-transplant complications) are not described in the literature. Therefore, it will be important in the future to conduct functional studies to determine whether these genetic variables produce changes in protein expression. In addition, these models should be validated with a larger sample size and in other HSCT regimens.

In conclusion, the incorporation of these predictive models in the management of transplanted patients may contribute to optimize the treatment and improve the outcomes of these patients in the context of a personalized precision medicine.