ORIGINAL RESEARCH article
Sec. Alloimmunity and Transplantation
Meta-Analysis of Genome-Wide Association and Gene Expression Studies Implicates Donor T Cell Function and Cytokine Pathways in Acute GvHD
- 1Finnish Red Cross Blood Service, Helsinki, Finland
- 2Stem Cell Transplantation Unit, Comprehensive Cancer Center, Helsinki University Hospital, Helsinki, Finland
- 3Turku University Hospital, Turku, Finland
- 4Department of Hematology, Genomics Unit, Hospital General Universitario Gregorio Marañón, Instituto de Investigación Sanitaria Gregorio Marañón (IiSGM), Madrid, Spain
- 5Department of Hematology, Institut Català d'Oncologia, Institut d'Investigació Biomèdica de Girona (IDIBGI), Girona, Spain
Graft-vs.-host disease (GvHD) is a major complication after allogeneic hematopoietic stem cell transplantation that causes mortality and severe morbidity. Genetic disparities in human leukocyte antigens between the recipient and donor are known contributors to the risk of the disease. However, the overall impact of genetic component is complex, and consistent findings across different populations and studies remain sparse. To gain a comprehensive understanding of the genes responsible for GvHD, we combined genome-wide association studies (GWAS) from two distinct populations with previously published gene expression studies on GvHD in a single gene-level meta-analysis. We hypothesized that genes driving GvHD should be associated in both data modalities and therefore could be detected more readily through their combined effects in the integrated analysis rather than in separate analyses. The meta-analysis yielded a total of 51 acute GvHD-associated genes (false detection rate [FDR] <0.1). In support of our hypothesis, this number was significantly higher than that in a permutation meta-analysis involving the whole data set, as well as in separate meta-analyses on the GWAS and gene expression data sets. The genes indicated by the meta-analysis were significantly enriched in 277 Gene Ontology terms (FDR < 0.05), such as T cell function and cytokine-mediated signaling pathways, and the results highlighted several established immune mediators, such as interleukins and JAK-STAT signaling, and presented TRAF6 and TERT as potential effector candidates. Altogether, the results support the chosen methodological approach, implicate a role of gene-level variation in donors' key immunological regulators predisposing patients to acute GVHD, and present potential targets for therapeutic intervention.
Allogeneic hematopoietic stem cell transplantation (HSCT) is utilized as a curative treatment for life-threatening hematological malignancies. Despite improvements in survival and reduced morbidity over time, graft- vs.-host disease (GvHD) remains a major severe complication of HSCT (1). In the early stages of the treatment, donor T cell-mediated alloimmune reactions against the recipient's tissues give rise to acute GvHD (aGvHD), and consequently the majority of GvHD prophylaxis is directed at reducing aGvHD. Conditioning leads to gastrointestinal track-derived diffusion of microbial antigens, which, in turn, induces extensive immune activation and cytokine storm (2). However, the immune pathology of chronic GvHD (cGvHD) more closely resembles the symptoms of common autoimmune diseases involving dysregulation of immune tolerance and chronic tissue damage (3).
The success of HSCT treatment is strongly influenced by genetic disparities between the recipient and donor (4). The best characterized are classical human leukocyte antigen (HLA) genes encoded within the human major histocompatibility complex (MHC) on chromosome 6. Mismatching recipient and donor for HLA-molecules is a well-established risk factor for GvHD and overall survival (5, 6). However, both HLA and minor histocompatibility antigen disparities have been shown to result in the beneficial graft- vs.-leukemia effect reducing the risk of relapse (7, 8). Emphasizing the complex genetic background of the condition, the susceptibility to GvHD is also potentially altered by donor or recipient polymorphisms in genes involved in immune responses (9–11) and drug metabolism (12, 13), common deletions (14), and regulatory elements (15). Nevertheless, owing to the lack of replication of many of the results in large independent cohorts, a full understanding of relevant genetic factors remains elusive.
The first genome-wide association studies (GWASs) on GvHD were published by Sato-Otsubo et al. (16) and Bari et al. (17) in 2015. Sato-Otsubo et al. studied unrelated HSCT in a Japanese population and identified an association between HLA-DPB1 allele disparity and aGvHD. Additionally, they discovered three novel loci, including one in HLA-DP region, linked to severe aGvHD. Bari et al. identified recipient SUFU rs17114808 as a susceptibility locus for aGvHD in a cohort of US children (17). However, this result was not replicated in German pediatric and adult cohorts (18). In 2017, Goyal et al. reported three HLA-DP region loci in recipient genomes associated with severe aGvHD in a population with European-American ancestry (19). A genome-wide approach was also employed by Martin et al. (20) and Ritari et al. (21, 22), who both studied the effect of genome-wide recipient-donor mismatching. They concluded that an increase in genome-wide recipient mismatching is associated with an increased risk for GvHD.
The incomplete understanding of molecular mechanisms of GvHD-pathology and the shortage of solid biomarkers have prompted studies on GvHD-related gene expression. Donor gene expression profiling was utilized to detect “stronger alloresponders” by Baron et al. (23) in 2007, and in 2008, Buzzeo et al. (24) reported the first preliminary molecular signature of aGvHD. In 2015, Furlan et al. suggested that the recipient T cell transcriptional profile could be employed in identification of novel therapeutics, and aurora kinase A was presented as a potential target (25).
We hypothesized that combining different data sets and types could enhance the ability to detect weak but systematic common gene-level effects underlying GvHD-pathogenesis. In our present study, we explore the cooperative relationship between GWAS signals mapped into genes and published gene expression profiles in aGvHD and cGvHD. We carried out GWAS on Finnish and Spanish sibling HSCT cohorts and integrated these data with donor gene expression data sets by a meta-analysis based on gene rankings. We sought to determine whether the two data modalities lend support to each other and to understand the functional categories that may explain the outcome of HSCT in terms of GvHD pathogenesis.
Materials and Methods
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.
Genotyping and Imputation
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.
Statistical Analysis of the Characteristics of the Study Cohorts
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.
Population Structure Inference and Logistic Regression Analysis of the Covariates
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−5, and missing call rate >0.1 were excluded from the analysis. After correction for multiple testing, a genome-wide P-value < 1 × 10−7 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.
Figure 1. Overview of the study setup and associated genes. (A) Schematic diagram showing the main steps of the analysis pipeline. (B) Manhattan plot of donor aGvHD meta-analysis results. The meta-analysis includes the donor GWAS II–IV vs. 0 results from Finnish Cohort 1, Spanish Cohort 1, Finnish Cohort 2, and six previously published aGvHD gene expression studies. The analysis was conducted as depicted in the Method section. The red line indicates a false detection rate of <0.1.
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.
Visualization of Colocalization of GWAS and eQTL Events
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.
Gene Expression Directions
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.
Data and Code Availability Statement
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).
Characteristics of the HSCT Study Cohorts and the GvHD-Related Risk Factors
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).
GWASs of the HSCT Study Cohorts
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−8 and in the donor genotype at a suggestive significance level of P < 5 × 10−5 [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.
Figure 2. Visualization of enrichment and sematic similarity of the aGvHD-associated gene ontologies. The aGvHD-associated genes discovered in the meta-analysis were analyzed for enriched Gene Ontology Biological Process (GO:BP) categories. (A) T cell activation and cytokine response-focused GO:BP categories. The X-axis shows the enrichment level, and the size of the circle depicts the false detection rate (FDR). (B) Sematic similarity analysis of all the GO:BP categories. The bubbles represent the individual GO:BP categories, and more related terms are closer in the plot. The uniqueness from the total mean is shown by a color scale, with blue indicating a less unique and red indicating a more unique category. The size of the circle indicates the number of detected genes within the underlying GO term. The labels present some of the cluster representatives.
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).
Figure 3. Validation of the aGvHD meta-analysis. (A) Results of meta-analyses on the GWAS and GE studies separately and the meta-analysis on the GWAS and GE data together. The left side panel shows the number of associated genes at a false detection rate (FDR) level of <0.1, and the right side panel shows the number of Gene Ontology biological process (GO:BP) categories at an FDR < 0.05. (B) Results of 100 meta-analyses with permuted gene order. The left side panel shows the numbers of aGvHD-associated genes at an FDR level of <0.1, and the right side panel shows the numbers of GO:BP categories at an FDR level of <0.05. The vertical red line depicts the corresponding values from the original meta-analysis. The meta-analyses include the donor GWAS II–IV vs. 0 results from Finnish Cohort 1, Spanish Cohort 1, Finnish Cohort 2, and six previously published aGvHD gene expression (GE) studies. The analyses were conducted as described in the Methods section.
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.
Colocalization of GWAS and eQTL Events and the Combined Effect Directions of Gene Expression Studies
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.
Figure 4. Colocalization of aGvHD GWAS and eQTL events of TRAF6 in Finnish Cohort 1 and Spanish Cohort 1. Visualization of colocalization events was performed using LocusCompareR as described in the Methods section. The cis-eQTL P-values of the TRAF6 locus were extracted from the eQTLGen Consortium database (http://www.eqtlgen.org) and the aGvHD-related P-values were obtained from the donor GWAS II–IV vs. 0 results. The dots in the scatter plots are colored according to their linkage disequilibrium to the colocalization lead variant. (A) Shows the colocalization of the eQTL and GWAS distributions of TRAF6 in Finnish Cohort 1, and (B) shows these distributions in Spanish Cohort 1.
Figure 5. Direction and distribution of expression of aGvHD-associated genes. The effect directions were derived from all aGvHD-linked gene expression studies as described in the Methods section. The box plots represent the interquartile range (IQR) with the median, and the whiskers represent a maximum IQR of 1.5. The X-axis depicts the normalized distribution of gene expression values of cases relative to the mean of controls. SD, standard deviation.
Data and Code Availability
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).
In this study, we performed a novel gene-level meta-analysis on GvHD by integrating GWAS data and gene expression results from different populations and study settings. We discovered pathways associated with aGvHD pathogenesis, implicating immunological responses in processes such as T cell function, cytokines, JAK-STAT signaling, and regulation of the TRAF6 gene.
The risk for GvHD has mainly been investigated during past decades by concentrating on the MHC region and specific candidate genes (10, 11), and only in recent years has the use of genome-wide approaches emerged in the field (16, 17, 19–22). These studies have shown the importance of genetic component in HSCT complications, but the results remain diverse, and their replication is incomplete. Even though individual gene expression studies can yield several differentially expressed genes, the dynamic nature of gene expression and heterogeneity in the treatment settings and genetic and environmental backgrounds of the subjects may make generalization from a single study difficult. Additionally, interpreting the functionality and downstream effects of disease-associated GWAS variants has generally been problematic because the majority of detected polymorphisms fall within non-coding regions of the genome. Our results highlight the strength of extensive meta-analysis; neither the independent adjusted GWASs nor the separate meta-analyses of gene expression or GWAS yielded statistically significant results alone. In contrast, the full data combination produced biologically meaningful genes and biological processes that were confirmed by permutation analysis. Thus, these results may indicate genes that have an impact on GvHD regardless of the various differences between the cohorts.
In the present study, the 51 aGvHD-associated genes included several prominent immunological effectors, supporting their role in GvHD pathogenesis. Many of these genes have been linked with GvHD in previous studies [Supplementary Table 5 (Supplementary Data Sheet 1)]. Connection with aGvHD has been demonstrated for telomerase reverse transcriptase (TERT), IL11, IL15, Kruppel like factor 2 (KLF2), STAT1, interferon regulatory factor 5 (IRF5), histone deacetylase 4 (HDAC4), JAK1, receptor interacting serine/threonine kinase 1 (RIPK1), TRAF3, and TRAF6 (45–58). As most of these genes have been studied mainly in murine models, their detection by our meta-analysis supports the role of these immune effectors in humans, and therefore provides candidates for potential drug targets. For instance, HDAC inhibitor Vorinostat has been shown to be promising in phase I/II trials (28784598).
Many of the implicated genes are parts of key immunoregulatory pathways such as JAK-STAT, TNF, and nuclear factor kappa B (NFκB). JAK1, an essential kinase for cytokine receptor signaling molecules, and STAT1, an activating transcription factor in response to pathogens, were involved in 33% of the significantly enriched biological processes. Importantly, JAK-STAT signaling assembles a connective molecular mechanism for multiple observations involving IL11, IL15 and other cytokines and links these results to the aGvHD-predisposing cytokine storm. Currently, the JAK1/2 inhibitor ruxolitinib is considered a potent salvage therapy for corticosteroid-refractory GvDH (59) and is being studied in a prospective randomized phase 2 trial (NCT02396628). The role of TNF signaling was also evident in the results; the herpesvirus entry mediator (HVEM), a member of the TNF receptor superfamily, is expressed on T-cells, and soluble LTBR protein has been shown to compete with HVEM receptor to inhibit T-cell activity and prevent GvHD in murine models. The RIPK1 receptor interacting serine/threonine kinase 1 is a downstream effector of TNFR2 (58), the activation of which has been shown to protect from aGvHD in a murine model and act concordantly in human cells in vitro (57). TNF down-stream signaling and transcription are regulated by TRAFs (50, 60). Three genes identified by the present study are closely involved with the NFκB pathway, a central regulator of cytokine production: NFκB inhibitor β, inhibitor of NFκB2 kinase subunit β, and NFκB2. Even though NFκB inhibition has not been successful in human trials, the important role of this pathway in immunological functions may warrant further studies in members of its signaling network.
The biologically active form of vitamin D, 1,25-dihydxoyvitamin D3, has been reported to exert beneficial immunosuppressive modulation via binding to the nuclear VDR (61). In an HSCT setting, vitamin D deficiency is considered a common complaint that may affect the treatment outcome (62). Genetic polymorphisms of the VDR have been shown to alter the immunomodulatory effect of vitamin D (63). Consistent with these observations, our results suggest a role for VDR in aGvHD-pathogenesis and present 24 VDR-associated biological processes. A possible mechanism could involve a positive effect of vitamin D on telomere length of blood cells directly or via increasing sex hormone levels, supporting proliferative capacity, reconstitution, and long-term immunological functionality of the graft (64–66). This hypothesis is supported by our results identifying TERT as significantly decreased in severe aGvHD, and enrichment of GOs involving VDR in hormone metabolism.
eQTL databases may be utilized to detect downstream effects of variants on gene expression. However, due to the abundance of cis-eQTL data, generating causal hypotheses for functional mechanisms of variants has been strongly affected by the high false-positive rate (43). To overcome this limitation, the colocalization analysis visualizes the distributions of eQTL and GWAS signals instead of focusing solely on the lead variants. Our results suggest colocalization of cis-eQTL events of TRAF6 and the corresponding variants of aGvHD in GWASs, indicating potential functional importance. TRAF6 has been identified as an important nuclear factor κB activator and an essential mediator for the GvHD-suppressing ability of thymic-derived regulatory T cells (49). Our results showing reduced TRAF6 expression among aGvHD patients are in consistent with this result. In contrast, high levels of TRAF6 have been associated with increased GvHD severity in mice (50). TRAF6 is regulated by miR-146a (67) and together with TRAF3, which was also implicated by our analysis, is a downstream effector of TLR molecules (47) involved in innate immunity-driven pathogenesis of GvHD (48).
There were several limitations of the study. The HSCT study cohorts were treated and collected over a long period of time in different centers, rendering treatment protocols heterogeneous. The numbers of cases were low in the separate GWASs and the case-control balance was therefore hardly optimal. Moreover, the gene expression studies were very heterogeneous in their study designs and sample collection procedures, and the overall number of samples was relatively low. This is to some extent mitigated by the employed meta-analysis method which was designed for noisy and heterogeneous data (40). Our approach was also not appropriate for analyzing the complex relationship between genetic variants implicated by GWAS and expression of the corresponding genes. To determine the effect directions of gene expression for GWASs, a transcriptome-wide association study method (68) would be required, but this could not be pursued here due to lack of appropriate reference material. This approach could be improved in future studies by using larger gene expression data set with more shared genetic variants. While the data we used were from Caucasian populations, it would be important to extend these analyses into other populations. Single nucleotide polymorphisms may vary among the populations which may impact the gene-level summary of the GWASs or regulation of gene expression.
We obtained no results for cGvHD, which may be due to the smaller number of available studies or lack of strong genetic component in this condition. We recently reported that the severity of cGvHD is associated with the overall immunogenetic differences between HSCT pairs (22); hence, its genetic background may differ from that of aGvHD.
We conclude that the meta-analysis of varied data types and populations identified common gene-level effects underlying the development of acute GvHD. As our approach was able to discover several previously known genes and immunologically relevant functional categories, it may be applied in different setting and to other data types to help identify novel genes and pathways.
Data Availability Statement
The datasets generated for this study will not be made publicly available because the limitations of ethical permits restrict the public distribution of personal data including individual genetic data.
The studies involving human participants were reviewed and 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. Written informed consent from the participants' legal guardian/next of kin was not required to participate in this study in accordance with the national legislation and the institutional requirements.
KH, SK, JP, and JR designed the study. LV, RN, and AN evaluated the clinical end points of Finnish Cohort 1, and MI-R, US, and MP evaluated Finnish Cohort 2. DG and IB provided the DNA samples and clinical data for Spanish Cohort 1. SK managed all DNA samples, verified HLA assignations and performed the preliminary GWAS of Finnish Cohort 1. KH performed the genomic imputation, PCA, GWASs, and logistic regression analyses. JR performed the meta-analyses, detection of gene expression effect directions, and visualization of co-localization. KH and JR interpreted the results and drafted the manuscript. All authors critically revised the final version of the manuscript and approved its submission for the publication.
This research was supported by the Finnish Funding Agency for Innovation Tekes to the SalWe GET IT DONE (GID) Personalized Diagnostics and Care program (ID 3982/31/2013) and the Academy of Finland (grant 288393).
Conflict of Interest
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.
We want to acknowledge the patients and the IDBGI Biobank (Biobanc IDIBGI, B.0000872), integrated in the Spanish National Biobanks Network, for their collaboration. We wish to acknowledge Sisko Lehmonen for the skillful and precise technical assistance, the FIMM Technology Center of the University of Helsinki, especially Drs. Janna Saarela and Kati Donner for the high-quality genotyping, and the CSC – IT Center for Science, Finland, for the computational resources.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2020.00019/full#supplementary-material
4. Warren EH, Zhang XC, Li S, Fan W, Storer BE, Chien JW, et al. Effect of MHC and non-MHC donor/recipient genetic disparity on the outcome of allogeneic HCT . Blood. (2012) 120:2796–806. doi: 10.1182/blood-2012-04-347286
5. Morishima Y, Kashiwase K, Matsuo K, Azuma F, Morishima S, Onizuka M, et al. Biological significance of HLA locus matching in unrelated donor bone marrow transplantation. Blood. (2014) 125:1189–97. doi: 10.1182/blood-2014-10-604785
6. Fürst D, Müller C, Vucinic V, Bunjes D, Herr W, Gramatzki M, et al. High-resolution HLA matching in hematopoietic stem cell transplantation: a retrospective collaborative analysis. Blood. (2013) 122:3220–9. doi: 10.1182/blood-2013-02-482547
7. Spierings E, Kim Y-H, Hendriks M, Borst E, Sergeant R, Canossi A, et al. Multicenter analyses demonstrate significant clinical effects of minor histocompatibility antigens on GvHD and GvL after HLA-matched related and unrelated hematopoietic stem cell transplantation. Biol Blood Marrow Transplant. (2013) 19:1244–53. doi: 10.1016/j.bbmt.2013.06.001
8. Hobo W, Broen K, Van der Velden WJFM, Greupink-Draaisma A, Adisty N, Wouters Y, et al. Association of disparities in known minor histocompatibility antigens with relapse-free survival and graft-versus-host disease after allogeneic stem cell transplantation. Biol Blood Marrow Transplant. (2013) 19:274–82. doi: 10.1016/j.bbmt.2012.09.008
9. Harkensee C, Oka A, Onizuka M, Middleton PG, Inoko H, Hirayasu K, et al. Single nucleotide polymorphisms and outcome risk in unrelated mismatched hematopoietic stem cell transplantation: an exploration study. Blood. (2012) 119:6365–72. doi: 10.1182/blood-2012-01-406785
10. Chien JW, Zhang XC, Fan W, Wang H, Zhao LP, Martin PJ, et al. Evaluation of published single nucleotide polymorphisms associated with acute GVHD. Blood. (2012) 119:5311–9. doi: 10.1182/blood-2011-09-371153
11. Hyvärinen K, Ritari J, Koskela S, Niittyvuopio R, Nihtinen A, Volin L, et al. Genetic polymorphism related to monocyte-macrophage function is associated with graft-versus-host disease. Sci Rep. (2017) 7:15666. doi: 10.1038/s41598-017-15915-3
12. Khaled SK, Palmer JM, Herzog J, Stiller T, Tsai N-C, Senitzer D, et al. Influence of absorption, distribution, metabolism, and excretion genomic variants on tacrolimus/sirolimus blood levels and graft-versus-host disease after allogeneic hematopoietic cell transplantation. Biol Blood Marrow Transplant. (2016) 22:268–76. doi: 10.1016/j.bbmt.2015.08.027
13. Byun JM, Kim H-L, Shin D-Y, Koh Y, Yoon S-S, Seong M-W, et al. The Impact of Methylenetetrahydrofolate Reductase C677T polymorphism on patients undergoing allogeneic hematopoietic stem cell transplantation with methotrexate prophylaxis. PLoS ONE. (2016) 11:e0163998. doi: 10.1371/journal.pone.0163998
14. McCarroll SA, Bradner JE, Turpeinen H, Volin L, Martin PJ, Chilewski SD, et al. Donor-recipient mismatch for common gene deletion polymorphisms in graft-versus-host disease. Nat Genet. (2009) 41:1341–4. doi: 10.1038/ng.490
15. Gam R, Shah P, Crossland RE, Norden J, Dickinson AM, Dressel R. Genetic association of hematopoietic stem cell transplantation outcome beyond histocompatibility genes. Front Immunol. (2017) 8:380. doi: 10.3389/fimmu.2017.00380
16. Sato-Otsubo A, Nannya Y, Kashiwase K, Onizuka M, Azuma F, Akatsuka Y, et al. Genome-wide surveillance of mismatched alleles for graft-versus-host disease in stem cell transplantation. Blood. (2015) 126:2752–63. doi: 10.1182/blood-2015-03-630707
17. Bari R, Hartford C, Keung Chan W, Vong Q, Li Y, Gan K, et al. Genome-wide single-nucleotide polymorphism analysis revealed SUFU suppression of acute graft-versus-host disease through downregulation of HLA-DR expression in recipient dendritic cells. Sci Rep. (2015) 5:11098. doi: 10.1038/srep11098
18. Katz MC, Michaelis S, Siegmund DM, Koch R, Bethge W, Handgretinger R, et al. Association analysis between SUFU polymorphism rs17114808 and acute graft versus host disease after hematopoietic stem cell transplantation. Bone Marrow Transplant. (2018) 53:377–82. doi: 10.1038/s41409-017-0059-3
19. Goyal RK, Lee SJ, Wang T, Trucco M, Haagenson M, Spellman SR, et al. Novel HLA-DP region susceptibility loci associated with severe acute GvHD. Bone Marrow Transplant. (2017) 52:95–100. doi: 10.1038/bmt.2016.210
20. Martin PJ, Levine DM, Storer BE, Warren EH, Zheng X, Nelson SC, et al. Genome-wide minor histocompatibility matching as related to the risk of graft-versus-host disease. Blood. (2017) 129:791–8. doi: 10.1182/blood-2016-09-737700
21. Ritari J, Hyvärinen K, Koskela S, Itälä-Remes M, Niittyvuopio R, Nihtinen A, et al. Genomic prediction of relapse in recipients of allogeneic haematopoietic stem cell transplantation. Leukemia. (2019) 33:240–8. doi: 10.1038/s41375-018-0229-3
22. Ritari J, Hyvärinen K, Koskela S, Niittyvuopio R, Nihtinen A, Salmenniemi U, et al. Computational analysis of HLA-presentation of non-synonymous recipient mismatches indicates effect on the risk of chronic graft-versus-host disease after allogeneic HSCT. Front Immunol. (2019) 10:1625. doi: 10.3389/fimmu.2019.01625
23. Baron C, Somogyi R, Greller LD, Rineau V, Wilkinson P, Cho CR, et al. Prediction of Graft-versus-host disease in humans by donor gene-expression profiling. PLoS Med. (2007) 4:e23. doi: 10.1371/journal.pmed.0040023
25. Furlan SN, Watkins B, Tkachev V, Flynn R, Cooley S, Ramakrishnan S, et al. Transcriptome analysis of GVHD reveals aurora kinase a as a targetable pathway for disease prevention. Sci Transl Med. (2015) 7:315ra191. doi: 10.1126/scitranslmed.aad3231
27. Shulman HM, Sullivan KM, Weiden PL, McDonald GB, Striker GE, Sale GE, et al. Chronic graft-versus-host syndrome in man. A long-term clinicopathologic study of 20 Seattle patients. Am J Med. (1980) 69:204–17. doi: 10.1016/0002-9343(80)90380-0
28. Howie BN, Donnelly P, Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. (2009) 5:e1000529. doi: 10.1371/journal.pgen.1000529
29. Anderson CA, Pettersson FH, Clarke GM, Cardon LR, Morris AP, Zondervan KT. Data quality control in genetic case-control association studies. Nat Protoc. (2010) 5:1564–73. doi: 10.1038/nprot.2010.116
31. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets . Gigascience. (2015) 4:7. eCollection 2015. doi: 10.1186/s13742-015-0047-8
32. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. (2007) 81:559–75. doi: 10.1086/519795
34. Glauzy S, Peffault de Latour R, André-Schmutz I, Lachuer J, Servais S, Socié G, et al. Alterations of circulating lymphoid committed progenitor cellular metabolism after allogeneic stem cell transplantation in humans. Exp Hematol. (2016) 44:811–6.e3. doi: 10.1016/j.exphem.2016.05.008
35. Lupsa N, Érsek B, Horváth A, Bencsik A, Lajkó E, Silló P, et al. Skin-homing CD8+ T cells preferentially express GPI-anchored peptidase inhibitor 16, an inhibitor of cathepsin K. Eur J Immunol. (2018) 48:1944–57. doi: 10.1002/eji.201847552
36. Takahashi N, Sato N, Takahashi S, Tojo A. Gene-expression profiles of peripheral blood mononuclear cell subpopulations in acute graft-vs-host disease following cord blood transplantation. Exp Hematol. (2008) 36:1760–74. doi: 10.1016/j.exphem.2008.07.007
37. Kohrt HE, Tian L, Li L, Alizadeh AA, Hsieh S, Tibshirani RJ, et al. Identification of gene microarray expression profiles in patients with chronic graft-versus-host disease following allogeneic hematopoietic cell transplantation. Clin Immunol. (2013) 148:124–35. doi: 10.1016/j.clim.2013.04.013
38. Hakim FT, Memon S, Jin P, Imanguli MM, Wang H, Rehman N, et al. Upregulation of IFN-inducible and damage-response pathways in chronic graft-versus-host disease. J Immunol. (2016) 197:3490–503. doi: 10.4049/jimmunol.1601054
39. Miller AM, Schneider R, Vance EA, Agura E, Powell S, Pascual V, et al. Differential Gene Expression Analysis Predicts For Response To Bortezomib In The Treatment Of Steroid Refractory Chronic Graft-Versus-Host Disease (cGVHD). Miller AM, editor. Blood. (2013) 122:4587. doi: 10.1182/blood.V122.21.4587.4587
44. Võsa U, Claringbould A, Westra H-J, Jan Bonder M, Deelen P, Zeng B, et al. Unraveling the polygenic architecture of complex traits using blood eQTL meta-analysis. bioRxiv. (2018) 18:10. doi: 10.1101/447367
45. Brestoff JR, Vessoni AT, Brenner KA, Uy GL, DiPersio JF, Blinder M, et al. Acute graft-versus-host disease following lung transplantation in a patient with a novel TERT mutation. Thorax. (2018) 73:489–92. doi: 10.1136/thoraxjnl-2017-211121
46. Hill GR, Cooke KR, Teshima T, Crawford JM, Keith JC, Brinson YS, et al. Interleukin-11 promotes T cell polarization and prevents acute graft-versus-host disease after allogeneic bone marrow transplantation. J Clin Invest. (1998) 102:115–23. doi: 10.1172/JCI3132
47. Häcker H, Redecke V, Blagoev B, Kratchmarova I, Hsu L-C, Wang GG, et al. Specificity in Toll-like receptor signalling through distinct effector functions of TRAF3 and TRAF6. Nature. (2006) 439:204–7. doi: 10.1038/nature04369
48. Xiao HW, Luo Y, Lai XY, Shi JM, Tan YM, He JS, et al. Donor TLR9 gene tagSNPs influence susceptibility to aGVHD and CMV reactivation in the allo-HSCT setting without polymorphisms in the TLR4 and NOD2 genes. Bone Marrow Transplant. (2014) 49:241–7. doi: 10.1038/bmt.2013.160
49. Lu Y, Hippen KL, Lemire AL, Gu J, Wang W, Ni X, et al. MiR-146b antagomir-treated human Tregs acquire increased GVHD inhibitory potency. Blood. (2016) 128:1424–35. doi: 10.1182/blood-2016-05-714535
50. Stickel N, Prinz G, Pfeifer D, Hasselblatt P, Schmitt-Graeff A, Follo M, et al. MiR-146a regulates the TRAF6 / TNF-axis in donor T cells during GVHD. Blood. (2014) 124:2586–96. doi: 10.1182/blood-2014-04-569046
51. Blaser BW, Roychowdhury S, Kim DJ, Schwind NR, Bhatt D, Yuan W, et al. Donor-derived IL-15 is critical for acute allogeneic graft-versus-host disease. Blood. (2005) 105:894–901. doi: 10.1182/blood-2004-05-1687
52. Pabbisetty SK, Rabacal W, Volanakis EJ, Parekh VV, Olivares-Villagómez D, Cendron D, et al. Peripheral tolerance can be modified by altering KLF2-regulated Treg migration. Proc Natl Acad Sci USA. (2016) 113:E4662–70. doi: 10.1073/pnas.1605849113
53. Ma H-H, Ziegler J, Li C, Sepulveda A, Bedeir A, Grandis J, et al. Sequential activation of inflammatory signaling pathways during graft-versus-host disease (GVHD): Early role for STAT1 and STAT3. Cell Immunol. (2011) 268:37–46. doi: 10.1016/j.cellimm.2011.01.008
54. Atarod S, Ahmed MM, Lendrem C, Pearce KF, Cope W, Norden J, et al. miR-146a and miR-155 expression levels in acute graft-versus-host disease incidence. Front Immunol. (2016) 7:56. doi: 10.3389/fimmu.2016.00056
55. Reddy P, Maeda Y, Hotary K, Liu C, Reznikov LL, Dinarello CA, et al. Histone deacetylase inhibitor suberoylanilide hydroxamic acid reduces acute graft-versus-host disease and preserves graft-versus-leukemia effect. Proc Natl Acad Sci USA. (2004) 101:3921–6. doi: 10.1073/pnas.0400380101
56. Carniti C, Gimondi S, Vendramin A, Recordati C, Confalonieri D, Bermema A, et al. Pharmacologic inhibition of JAK1/JAK2 signaling reduces experimental murine acute GVHD while preserving GVT effects. Clin Cancer Res. (2015) 21:3740–9. doi: 10.1158/1078-0432.CCR-14-2758
57. Chopra M, Biehl M, Steinfatt T, Brandl A, Kums J, Amich J, et al. Exogenous TNFR2 activation protects from acute GvHD via host T reg cell expansion. J Exp Med. (2016) 213:1881–900. doi: 10.1084/jem.20151563
59. Zeiser R, Burchert A, Lengerke C, Verbeek M, Maas-Bauer K, Metzelder SK, et al. Ruxolitinib in corticosteroid-refractory graft-versus-host disease after allogeneic stem cell transplantation: a multicenter survey. Leukemia. (2015) 29:2062–8. doi: 10.1038/leu.2015.212
62. Ros-Soto J, Anthias C, Madrigal A, Snowden JA. Vitamin D: is it important in haematopoietic stem cell transplantation? A review. Bone Marrow Transplant. (2018) 54:810–20. doi: 10.1038/s41409-018-0377-0
63. Carrillo-Cruz E, García-Lozano JR, Márquez-Malaver FJ, Sánchez-Guijo FM, Montero Cuadrado I, Ferra-i-Coll C, et al. Vitamin D modifies the incidence of graft-versus-host disease after allogeneic stem cell transplantation depending on the vitamin D receptor (VDR) polymorphisms. Clin Cancer Res. (2019) 25:4616–23. doi: 10.1158/1078-0432.CCR-18-3875
64. Richards JB, Valdes AM, Gardner JP, Paximadas D, Kimura M, Nessa A, et al. Higher serum vitamin D concentrations are associated with longer leukocyte telomere length in women. Am J Clin Nutr. (2007) 86:1420–5. doi: 10.1093/ajcn/86.5.1420
65. Calado RT, Yewdell WT, Wilkerson KL, Regal JA, Kajigaya S, Stratakis CA, et al. Sex hormones, acting on the TERT gene, increase telomerase activity in human primary hematopoietic cells. Blood. (2009) 114:2236–43. doi: 10.1182/blood-2008-09-178871
66. Törlén J, Gaballa A, Remberger M, Mörk L-M, Sundberg B, Mattsson J, et al. Effect of graft-versus-host disease prophylaxis regimens on T and B cell reconstitution after allogeneic hematopoietic stem cell transplantation. Biol Blood Marrow Transplant. (2019) 25:1260–8. doi: 10.1016/j.bbmt.2019.01.029
67. Tang Y, Luo X, Cui H, Ni X, Yuan M, Guo Y, et al. MicroRNA-146A contributes to abnormal activation of the type I interferon pathway in human lupus by targeting the key signaling proteins. Arthritis Rheum. (2009) 60:1065–75. doi: 10.1002/art.24436
Keywords: GvHD, GWAS, gene expression, meta-analysis, HSCT
Citation: Hyvärinen K, Koskela S, Niittyvuopio R, Nihtinen A, Volin L, Salmenniemi U, Putkonen M, Buño I, Gallardo D, Itälä-Remes M, Partanen J and Ritari J (2020) Meta-Analysis of Genome-Wide Association and Gene Expression Studies Implicates Donor T Cell Function and Cytokine Pathways in Acute GvHD. Front. Immunol. 11:19. doi: 10.3389/fimmu.2020.00019
Received: 11 November 2019; Accepted: 07 January 2020;
Published: 03 February 2020.
Edited by:Hermann Einsele, Julius Maximilian University of Würzburg, Germany
Reviewed by:Shigeo Fuji, Osaka International Cancer Institute, Japan
Robert Zeiser, University of Freiburg, Germany
Copyright © 2020 Hyvärinen, Koskela, Niittyvuopio, Nihtinen, Volin, Salmenniemi, Putkonen, Buño, Gallardo, Itälä-Remes, Partanen and Ritari. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kati Hyvärinen, firstname.lastname@example.org