The study consisted of three HLA-matched sibling cohorts including two populations: Finnish Cohort 1, Spanish Cohort 1, and Finnish Cohort 2. The characteristics of these study cohorts are presented in Table 1.

Finnish Cohort 1 and Spanish Cohort 1 have been described previously in detail (11). Briefly, Finnish Cohort 1 consisted of 239 donor-recipient pairs, 23 individual recipients, and 28 individual donors with available clinical data and imputed genotype. All recipients received related donor (RD)-HSCT from 1993 through 2006 at Helsinki University Hospital, Comprehensive Cancer Center, Stem Cell Transplantation Unit, Finland. Spanish Cohort 1 was composed of 253 donor-recipient pairs, 15 individual recipients, and 30 individual donors with available clinical data and the imputed genotype. The recipients received RD-HSCT from 2002 through 2014 at 13 Spanish transplant centers. HLA-matching was performed at the HLA-A, -B, and -DRB1 loci in both cohorts.

Finnish Cohort 2 included 171 donor-recipient pairs, 3 individual recipients, and 1 individual donor with clinical data and an imputed genotype. RD-HSCT was implemented at two Finnish centers: Helsinki University Hospital, Comprehensive Cancer Center, Stem Cell Transplantation Unit and Turku University Central Hospital during the years 2006–2016. HLA-matching was performed at the HLA-A, HLA-B, HLA-C, and HLA-DRB1 loci (21).

The included clinical outcomes were grades 0, II–IV and III–IV for aGvHD, and grades 0, limited, and extensive for cGvHD. Grading was determined locally according to the European Society for Blood and Marrow Transplantation guidelines (26, 27).

The study conformed to the principles of the Declaration of Helsinki and was approved by the Finnish National Supervisory Authority for Welfare (Dnro V/74832/2017, V/3235/2019) and Health and the ethics committees of Helsinki University Central Hospital (382/13/03/01/2014, HUS/114/2018) and Turku University Central Hospital (ETMK 78/2012). Samples and data from Spanish patients included in this study were provided by the IDIBGI Biobank (Biobanc IDIBGI, B.0000872), integrated in the Spanish National Biobanks Network, and they were processed following standard operating procedures with the appropriate approval of the ethics and scientific committees.

The genotyping and imputation procedures have been previously described in detail (11). Briefly, Finnish Cohort 1 was genotyped using an Illumina Immunochip and Spanish Cohort 1 and Finnish Cohort 2 with Illumina were genotyped with an Immunoarray v2.0. Imputation of autosomal genotype data was performed using IMPUTE2 and the 1000 Genomes Phase 3 reference (28). Quality filtering for variants and samples was performed according to Anderson et al. (29), and an IMPUTE2 INFO-field measure ≥0.5 was used as the cut-off for post-imputation filtering (30). The three datasets were filtered and imputed in separate processes, and the final number of variants included in the analyses was 5041081 for Finnish Cohort 1, 5737173 for Spanish Cohort 1, and 9105726 for Finnish Cohort 2.

The significance of differences in recipient and donor age among the three study cohorts was analyzed using the non-parametric Kruskal–Wallis test. Pearson's chi-square test was used to evaluate the frequencies of diagnosis, gender direction of transplantation, stem cell source, conditioning regimen, and GvHD grades. Due to low frequency, the variation in Hodgkin's lymphoma and aplastic anemia diagnosis were not analyzed. The significance level of tests (alpha level) was set at 0.05.

The population structures were determined using principal component (PC) analysis (PCA). The analysis was implemented with PLINK 1.90b4.1 (www.cog-genomics.org/plink/1.9/) (31). Genotyped variants shared by all three cohorts were pruned to generate a subset of variants in approximate linkage equilibrium (–indep-pairwise 50 5 0.5). Dimensional reduction was executed with the command –pca and the top 20 PCs for each cohort were extracted separately. The combined data were used to generate a scatterplot matrix of the five initial PCAs [Supplementary Figure 1 (Supplementary Data Sheet 1)]. The plot was generated using R version 3.5.0.

Logistic regression analysis (PLINK command –logistic) was employed to investigate the effect of baseline covariates on the GvHD-related outcome. Binary aGvHD or cGvHD status was used as a dependent variable, and recipient gender, recipient age, direction of transplantation, and stem cell source were included as independent variables [Supplementary Table 1 (Supplementary Data Sheet 1)]. The results are presented as odds ratios (ORs) with 95% confidence intervals (CIs), and a P-value < 0.05 was considered statistically significant.

In the preliminary analysis, the associations between recipient and donor genotypes and GvHD clinical outcomes were examined using the PLINK 1.07 (32) 1df chi-square allelic test (command –assoc) and are expressed as odds ratios with 95% confidence intervals. Variants with a minor allele frequency < 0.01, Hardy-Weinberg disequilibrium P-value < 1 × 10⁻⁵, and missing call rate >0.1 were excluded from the analysis. After correction for multiple testing, a genome-wide P-value < 1 × 10⁻⁷ was considered statistically significant.

For the adjusted analyses, the association between variants and clinical outcomes was determined using PLINK 1.90b4.1 (31) logistic regression analysis (command –logistic). Variant filtering followed the previously described practice. For all recipients, recipient gender, recipient age, graft type, and the first three cohort specific PCs were used as covariates. For the donors, donor gender, donor age, graft type, and the first three PCs were used as covariates. The donor age data were incomplete in Finnish Cohort 1 and Spanish Cohort 1, and the missing data were imputed using the corresponding recipient's age using linear regression.

The Figure 1A presents a schematic diagram of the main steps of the analysis pipeline. The SNP associated P-values of donor aGvHD grades II–IV vs. 0 and cGvHD limited-extensive vs. 0 GWASs from all three study cohorts were mapped into genes using MAGMA (33) v1.07b (https://ctg.cncr.nl/software/magma) with the default 1000 Genomes GRCh37 LD data and gene annotation data provided with MAGMA.

GvHD-related gene expression (GE) data sets (23–25, 34–39) were retrieved from the Gene Expression Omnibus data repository (https://www.ncbi.nlm.nih.gov/geo/). A summary of the selected GE studies is presented in the Supplementary Table 2 (Supplementary Data Sheet 1). Analysis of differential expression in the GE data sets was conducted using the R package limma v3.38.3 functions lmFit for fitting linear models, and eBayes to determine the empirical Bayes moderated standard errors of the models.

Using the gene lists obtained from the GWAS and GE analyses, we conducted three independent meta-analyses: (1) the GWASs together, (2) the GE studies together, and (3) the GWASs and the GE studies combined, for both the aGvHD and cGvHD data sets. These analyses were executed using the robust rank aggregation (RRA) method (40) implemented in the R package RobustRankAggreg v1.1 function aggregateRanks using exacted P-value calculations. Genes with a false discovery rate (FDR) < 0.1 were considered significant and included in the downstream analyses for gene ontology (GO) biological processes (BP) enrichment. GO:BP enrichment was performed using the R package clusterProfiler (41) v3.10.1 function enrichGO with the full list of analyzed genes as the background set. GO terms with an FDR < 0.05 were considered statistically significant.

To control for subtle biases that could potentially cause inflated numbers of false positives, we conducted permutation tests on the gene lists included in the RRA meta-analysis. The gene list permutation was repeated 100 times, and for each iteration, the number of genes with an FDR < 0.1 was calculated. These values were then compared with the numbers of significant genes from the RRA conducted on the original gene lists to estimate whether the observed number of significant genes could be obtained by chance alone.

A similar permutation test was conducted for the GO:BP enrichment analysis, whereby the enrichment at an FDR < 0.05 was re-calculated for each iteration of the permutated RRA for genes with an FDR < 0.1. Each enriched gene combination was accepted only once to limit the number of redundant GO terms.

The clustering and visualization of the GO enrichment results were performed using the REVIGO (42) Web server (http://revigo.irb.hr/).

The visualization of colocalization of GWAS and expression quantitative trait loci (eQTL) events was performed using the R package LocusCompareR (43). The GWAS variants of aGvHD-associated genes were selected using a ± 1Mb range from the gene boundaries. The corresponding blood cis-eQTL SNPs were extracted from the eQTLGen Consortium database (http://www.eqtlgen.org) (44). This database incorporates 31,684 individuals and 37 datasets, resulting in 16,989 cis-eQTL genes. The generated plots visualized the distributions of cis-eQTL and GWAS signals and the linkage disequilibrium to the selected variant.

To visualize the direction and distribution of expression of the significant RRA genes, the gene-wise expression values were extracted from each included GE study. For each study, the expression values for the cases were divided by the mean value of the controls and thereafter scaled by subtracting the mean and dividing by the standard deviation. The processed values were then pooled for each gene from all of the GE studies and plotted as boxplots.

The limitations of ethical permits restrict the public distribution of personal data including individual genetic data. The code for the meta-analyses is publicly available in GitHub (https://github.com/FRCBS/GvHD_meta).

The supplementary information includes Supplementary Data Sheets 1–4.

Table 1 summarizes the key characteristics of the HSCT study cohorts. The three cohorts differed substantially in many aspects. The distributions of recipient and donor ages and the frequencies of diagnoses, stem cell sources, conditioning regimens, and GvHD outcomes significantly varied among the cohorts. Additionally, the GvHD prophylaxis regimens varied markedly. The association between the common GvHD-related risk factors and the disease outcomes diverged as depicted in Supplementary Table 1 (Supplementary Data Sheet 1). Increasing recipient age was significantly associated with all GvHD outcomes in Finnish Cohort 1 (all P-values ≤ 0.036). However, the odds ratios were low (from 1.035 to 1.087) and age had no effect in the other two study cohorts. Female gender was associated with beneficial cGvHD outcomes in Spanish Cohort 1 (P-values = 0.001). Using the bone marrow as a stem cell source reduced the risk for both aGvHD and cGvHD in Finnish Cohort 2 (P-values ≤ 0.045, except for 0.065 in the cGvHD limited-extensive group) and for cGvHD in Finnish Cohort 1 (P-values ≤ 0.007).

In the unadjusted preliminary GWAS in Finnish Cohort 1, we found severe aGvHD-associated loci in the MHC region in the recipient genotype at a genome-wide significance level of P < 5 x 10⁻⁸ and in the donor genotype at a suggestive significance level of P < 5 × 10⁻⁵ [Supplementary Figure 2 (Supplementary Data Sheet 1), panel A for recipients and panel C for donors]. These results were not replicated in the two other HSCT study cohorts, and none of the variants reached a genome-wide significance level (data not shown). Additionally, all genome-wide significance was abolished in Finnish Cohort 1 after adjusting the analyses by recipient age, recipient gender, stem cell source, and the top three PCs [Supplementary Figure 2 (Supplementary Data Sheet 1), panel B for recipients and D for donors].

All data sets included in the RRA meta-analyses are presented in Table 2, and more detailed information on the GE studies performed is listed in Supplementary Table 2 (Supplementary Data Sheet 1). The combined meta-analysis of cGvHD-related data sets revealed only cysteine protease legumain (LGMN) as associated with cGvHD at the FDR < 0.1 level. The analysis of aGvHD-related data sets revealed 51 aGvHD-associated genes at the FDR < 0.1 level, including lymphotoxin beta receptor (LTBR), Janus kinase 1 (JAK1), tumor necrosis factor (TNF) receptor associated factor 6 (TRAF6), signal transducer and activator of transcription 1 (STAT1), vitamin D receptor (VDR), interleukin (IL) 11, IL15, and IL1 receptor 2 (IL1R2) [Figure 1; Supplementary Table 3 (Supplementary Data Sheet 1)]. The GO enrichment analysis of these genes detected 277 aGvHD-associated BPs at the FDR < 0.05 level [Supplementary Table 4 (Supplementary Data Sheet 1)]. The majority of associated GO:BP categories were strongly linked to immune responses and regulation, highlighting T cell function and cytokine-mediated signaling pathways. The top 30 of these detailed GO:BP categories are presented in Figure 2A. The degree of relatedness of all ontological categories defined by their annotation, the sematic similarity, is shown in Figure 2B.

To confirm the validity of the combined meta-analysis including both the GWASs and the GE studies, we performed a meta-analysis of the GWASs and the GE studies separately. Figure 3A presents the comparison of aGvHD-associated genes and unique GO:BP gene combinations in different analysis settings. The meta-analysis of the three GWASs did not result in any associated genes at the FDR < 0.1 level and, consequently, no GO:BPs at the FDR < 0.05 level. We found two genes and three GO:BPs associated with aGvHD in the separate GE meta-analyses, but these numbers were substantially lower than those in the combined meta-analysis (51 genes and 156 unique GO:BP gene combinations).

The mean number of significant (FDR < 0.1) permuted RRA genes was 5.69 with an SD of 0.48, which corresponded with the expected number of 5.1 based on an FDR < 0.1 in the true data. In the GO enrichment analysis of permuted data, the mean number of significant (FDR < 0.05) unique GO terms was 7.95 with an SD of 1.28, which was also close to the expected value of 7.8 obtained from the true data at an FDR < 0.05. The number of significant genes or GO terms did not reach the value obtained with the true data in any of the permutation iterations. Figure 3B depicts the results of the permutation tests.

The colocalization of GWAS and cis-eQTL events of the 51 aGvHD-associated genes were visualized using LocusCompareR and the results for Finnish Cohort 1, Spanish Cohort 1, and Finnish Cohort 2 are presented in Supplementary Data Sheets 2–4, respectively. The majority of events showed no direct positive colocalization signal. However, the cis-eQTL events of TRAF6 and the corresponding variants of aGvHD GWASs were colocalized in Finnish Cohort 1 (Figure 4A). This observation was supported by Spanish Cohort 1, although the lead variants were independent (Figure 4B). We formulated combined effect directions of gene expression for the 51 aGvHD-associated genes from the aGvHD GE data sets (Figure 5). The expression levels of TRAF6, TRAF3, IL1R2, interferon induced protein with tetratricopeptide repeats 5, and bone morphogenetic protein 6 seemed to be reduced among the aGvHD patients compared to the controls.

The limitations of ethical permits restrict the public distribution of personal data including individual genetic data. The code for the meta-analyses is publicly available in GitHub (https://github.com/FRCBS/GvHD_meta).

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.