Analysis of Relationships between Immune Checkpoint and Methylase Gene Polymorphisms and Outcomes after Unrelated Bone Marrow Transplantation

Simple Summary Hematopoietic stem-cell transplantation (HSCT) is a curative therapy for blood disorders. Unrelated bone marrow transplantation (uBMT) is a type of allogeneic HSCT that uses the bone marrow of an unrelated donor. While HLA mismatch is a risk factor for poor outcomes in HSCT, such as graft-versus-host disease (GVHD), the importance of non-HLA single-nucleotide polymorphisms (SNPs) remains unclear. The clinical application of immune checkpoint and chromatin methylation inhibitors to cancer has been attracting attention. In the present study, we retrospectively genotyped five SNPs in four immune checkpoint genes, BTLA, PD-1, LAG3, and CTLA4, and two SNPs in methylase genes, DNMT1 and EZH2, in 999 uBMT pairs. Although no correlations were observed between these SNPs and post-uBMT outcomes, recipient EZH2 SNP exhibited a low p-value in the analysis of grade 2–4 acute GVHD (p = 0.010). This SNP may be useful for outcome predictions and needs to be confirmed in a larger-scale study. Abstract Unrelated bone marrow transplantation (uBMT) is performed to treat blood disorders, and it uses bone marrow from an unrelated donor as the transplant source. Although the importance of HLA matching in uBMT has been established, that of other genetic factors, such as single-nucleotide polymorphisms (SNPs), remains unclear. The application of immunoinhibitory receptors as anticancer drugs has recently been attracting attention. This prompted us to examine the importance of immunoinhibitory receptor SNPs in uBMT. We retrospectively genotyped five single-nucleotide polymorphisms (SNPs) in the immune checkpoint genes, BTLA, PD-1, LAG3, and CTLA4, and two SNPs in the methylase genes, DNMT1 and EZH2, in 999 uBMT donor–recipient pairs coordinated through the Japan Marrow Donor Program matched at least at HLA-A, -B, and -DRB1. No correlations were observed between these SNPs and post-uBMT outcomes (p > 0.005). This result questions the usefulness of these immune checkpoint gene polymorphisms for predicting post-BMT outcomes. However, the recipient EZH2 histone methyltransferase gene SNP, which encodes the D185H substitution, exhibited a low p-value in regression analysis of grade 2–4 acute graft-versus-host disease (p = 0.010). Due to a low minor allele frequency, this SNP warrants further investigation in a larger-scale study.


Introduction
Hematopoietic stem-cell transplantation (HSCT) is a curative therapy for blood disorders, such as leukemia. Unrelated bone marrow transplantation (uBMT) is an allogeneic HSCT that uses bone marrow from an unrelated donor as the transplant source. Due to advances in treatment, overall survival (OS) after uBMT, similar to other types of allogeneic HSCTs, has been prolonged; however, further improvements are needed. Acute and chronic graft-versus-host disease (GVHD) pose serious risks to post-transplant recipients. HLA mismatches between a recipient and donor are well-established risk factors for poor HSCT outcomes [1][2][3]. Moreover, even when all 12 HLA-A, -B, -C, -DPB1, -DQB1, and -DRB1 loci are matched at allele levels, severe (grade 3-4) acute GVHD (aGVHD) occurred in 11% of recipients and the 5 year overall survival rate was 53% in uBMT [3]. However, it currently remains unclear whether genetic polymorphisms outside HLA loci affect HSCT outcomes. Non-HLA single-nucleotide polymorphisms (SNPs) are attracting interest because of their potential in outcome predictions [4], as well as donor selection. We previously investigated the relationships between SNPs in innate immune pathway genes and BMT outcomes [5,6].
Chromatin-modifying enzymes, particularly DNA and histone methylases, play important roles in the development of hematological malignancies [14]. Hypomethylating agents, such as azacytidine, are in clinical use for patients with acute myeloid leukemia (AML) and myelodysplastic syndromes (MDS) [15]. The administration of azacytidine causes the demethylation of the PD1 promoter in patients with AML and MDS [16], and it has also been shown to mitigate GVHD in mice [17]. A target of azacytidine, DNA methyltransferase 1 (DNMT1), methylates the cytosine residues of DNA to maintain the DNA methylation pattern [15,18]. The Enhancer of Zeste 2 (EZH2) gene encodes histone H3 lysine 27 methyltransferase and is a component of polycomb repressive complex-2 (PRC-2) [19]. Ezh2 interacts with DNMTs and is required for the DNA methylation of Ezh2 target promoters [20]. Ezh2 appears to play important roles in hematological malignancies [21,22]. Its loss in donor T cells inhibits GVHD in mice [23]; however, this may not be the case in humans [24]. Cancer cells exploit the PRC2-mediated transcriptional silencing of the MHC-I antigen processing pathway, evading T-cell-mediated immunity [25].
Therefore, we herein investigated the relationships between SNPs in the four aforementioned immunoinhibitory receptors genes and two methylases genes and the outcomes of uBMT matched at least at HLA-A, -B, and -DRB1, but not necessarily at HLA-C, which is a typical HLA matching pattern in uBMT [35], coordinated between May 2006 and April 2009 through the Japan Marrow Donor Program (JMDP) [36].

Subjects
The characteristics and transplantation outcomes of the donors and patients in the present study were previously described (Table S1, of Takahashi et al. [6]). Briefly, 999 donor-recipient pairs who satisfied all of the following criteria were retrospectively genotyped:  the pair underwent uBMT between May 2006 and April 2009 through the JMDP; Japanese  ethnicity; donor age of at least 20 years; days of survival were available; HLA-A, -B, -C,  -DPB1, -DQB1, and -DRB1 alleles were retyped; HLA-A, -B, and -DRB1 were confirmed to  be matched (HLA-C, -DPB1, and -DQB1 were not necessarily matched) [6]. Among them, 822 malignant-disease patients without a previous transplantation history and their uBMT donors were the subjects of the main analysis (Table S1), similar to the previous study [6]. Clinical data were collected by the Japan Society for Hematopoietic Cell Transplantation using the Transplant Registry Unified Management Program (TRUMP) [37].

SNP Selection
We selected SNPs that encode nonsynonymous amino-acid substitutions and have a minor allele frequency (MAF) >0.05 in 104 Japanese residents of Tokyo (JPT104) from the 1000 Genomes Project [38]. We used rs9288952 and rs76844316 for BTLA, rs2227982 for PDCD1, rs870849 for LAG3, rs231775 for CTLA4, rs2302427 for EZH2, and rs2228612 for DNMT1 (Tables S2 and S3, Supplementary Materials). Among these SNPs, rs76844316 and rs231775 are functional SNPs, with one of the two alleles of each SNP exhibiting stronger activity than the other allele in cultured cells (Table S2) [39,40].

SNP Genotyping
Genomic DNA preparation, TaqMan SNP Genotyping Assays (Applied Biosystems, Waltham, MA, USA), and direct Sanger sequencing following PCR were performed as previously described [6]. In TaqMan assays, the default threshold set by the software (quality = 95) was used for all SNPs, except rs2227982, for which quality = 90 was used as the threshold. To confirm the results of these TaqMan assays, approximately 10 samples for each of the three genotypes, two homozygous and one heterozygous, were re-genotyped by direct Sanger sequencing following PCR using the primers shown in Table S2. DNA samples unsuccessful in TaqMan genotyping were then genotyped by the aforementioned direct Sanger sequencing following PCR.

Outcomes
The primary outcome of the present study was grade 2-4 aGVHD within 100 days of transplantation. Secondary outcomes were grade 3-4 aGVHD within 100 days, chronic GVHD (cGVHD), extensive cGVHD (ecGVHD), OS, non-relapse mortality (NRM), and relapse. The detailed definitions and competing events of these outcomes used in the present study were identical to those previously described [6].

Covariates
Recipient characteristics used as covariates were sex, age, disease stage, body mass index, the cytomegalovirus (CMV) serostatus, and performance status. The donor characteristics used were sex, age, and the CMV serostatus. Transplantation characteristics used as covariates were myeloablative conditioning, cyclosporine A, arabinofuranosyl cytidine (Ara-C), cyclophosphamide, the number of nucleated cells infused, and days from diagnosis to BMT. Mismatches between a recipient and the corresponding donor of the ABO blood type, HLA-C, HLA-DPB1, and HLA-DQB1 were also used as covariates. Total HLA mismatches denote the sum of the numbers of mismatches at HLA-C, -DPB1, and -DQB1 because HLA-A, -B, and -DRB1 were matched. These covariates are shown in Table  S1 and were previously described [6].

Statistical Analysis
Statistical analyses were conducted as previously described [6]. The additive genetic model was defined as the risk per number of the minor allele. The recessive model was defined as the risk of the minor homozygous genotype relative to the two other genotypes, and the dominant model was similarly defined. OS was examined using Cox propor-tional hazard regression. Other outcomes were analyzed by the Fine-Gray proportional sub-distribution hazard (SH) regression. Reported covariates fixed in multivariable regression of aGVHD, cGVHD, and OS were previously described [6]. Fixed covariates in multivariable regression of NRM and relapse were the same as those used in multivariable regression of OS. The Wald test was used to calculate p-values in the regression analysis. All p-values are two-tailed unadjusted values, and p < 0.005 was considered to be significant. The R software (ver. 3.2), the EZR software [41], VCFtools (ver. 0.1.11), and Microsoft Excel were used.

Genotyping
We genotyped all SNPs in 999 subjects, except for the LAG3 SNP, rs870849. It was not possible to genotype one recipient for rs870849. The TaqMan assay produced a very weak signal and three independent trials of PCR amplification for the Sanger sequencing failed for an unknown reason; therefore, the rs870849 genotype of this recipient was finalized as undetermined (Table S3). However, this did not affect the statistical analysis because this recipient underwent a previous transplantation and was not a subject in the main analysis.
Allele frequencies were similar among the donors, recipients, and 104 Japanese individuals in Tokyo (JPT104) from the 1000 Genomes Project (Table S4). The null hypothesis for the Hardy-Weinberg equilibrium (HWE) was not rejected for any SNPs (p > 0.005; Table S4, Supplementary Materials). Linkage disequilibrium (LD) between two SNPs in a single gene, namely, rs9288952 and rs76844316 in BTLA, was similar among the recipients, donors, and JPT104 (Table S5). Therefore, there was no strong selection for or against a certain genotype at these SNPs in uBMT donors or recipients, and, as a consequence, all SNPs were included in the statistical analysis.

Grade 2-4 Acute GVHD
We conducted a univariable regression of grade 2-4 aGVHD, the primary outcome, on these seven SNPs under additive, dominant, and recessive genetic models. None of these SNPs correlated with grade 2-4 aGVHD (p > 0.005; Table S6). We then performed the multivariable regression analysis previously described, which is based on BIC-based covariate selection and takes interaction terms into account. However, no SNP variable or SNP-covariate interaction terms correlated with grade 2-4 aGVHD (data not shown). To make the results comparable with previous findings, we also performed a multivariable regression adjusted with reported covariates for grade 2-4 aGVHD. The results obtained were essentially unchanged, and no SNP correlated with grade 2-4 aGHVD (p > 0.005; Table 1). Among them, the recipient EZH2 SNP in the additive model exhibited low p-values (p = 0.011 and 0.010 in univariable and multivariable analyses, respectively).

Discussion
In the present study, we initially examined the relationships between polymorphisms in the immune checkpoint genes, BTLA, PD-1, LAG3, and CTLA4, and BMT outcomes, and we found no correlations.
In a murine model, BTLA played distinct roles in GVHD; the administration of an anti-BTLA monoclonal antibody inhibited donor anti-host T-cell responses, whereas BTLA also served as a ligand that sent a prosurvival signal in donor T cells [42]. The two BTLA SNPs tested in the present study were not associated with aGVHD. One explanation for this result is a difference between mice and humans. For example, many murine BMT models are MHC-mismatched and humans do not have a genetically uniform background, which is in contrast to experimental animals. Another explanation is the nature of these two SNPs. Both SNPs are non-synonymous in terms of amino-acid substitutions. One of them, rs9288952, is associated with malignant breast cancer in Chinese women [43] and has an MAF > 0.25; however, no molecular function has been reported. Therefore, this SNP may not be a functional SNP or may be associated with a functional SNP. The other BTLA SNP, rs76844316, was shown to be functional, with the T, but not G, allele being inhibitory against concanavalin A-and anti-CD3 Ab-induced IL-2 production when retrovirally introduced into Jurkat T cells [39]. Although this SNP is associated with susceptibility to rheumatoid arthritis (RA), the MAF is <0.1. Therefore, the lack of a correlation between this SNP and aGVHD may be attributed to its low MAF. The location of a variant within the BTLA gene may also be important. These two SNPs encode amino-acid substitutions within the cytoplasmic domain of the BTLA protein (https://www.uniprot.org/uniprot/Q7Z6A9 (last accessed in 1 March 2021)). There was no SNP within the extracellular domain affecting ligand binding that met our SNP selection criteria. Collectively, the present results do not necessarily indicate the true lack of a relationship between the BTLA gene and aGVHD.
PD-1 is another immune checkpoint protein, and its blockade is of widespread interest as a treatment for cancer, including Hodgkin's lymphoma [8,44]. Blockade of the PD-1/PD-L pathway has been shown to enhance aGVHD lethality in mice [45]. In a Spanish study, two PD-1 SNPs in donors were associated with grade 2-4 aGVHD in HLA-identical sibling HSCT [34]. In the present study, another PD-1 SNP, rs2227982, was not associated with any post-uBMT outcomes. Although rs2227982 may not be a functional PD-1 SNP, previous studies on east Asians reported correlations with the risk of type 1 diabetes, with glucose and insulin levels after the oral glucose tolerance test [46], and with ankylosing spondylitis [47]. Therefore, even if PD-1 is associated with BMT outcomes, this relationship may depend on specific clinical characteristics. A similar explanation may be applied to the LAG3 SNP used in the present study, rs870849, which was previously shown to be associated with multiple sclerosis [48].
CTLA4 SNPs correlated with transplant outcomes in some [27,28,[30][31][32], but not all HSCT studies [33,49], and discrepancies were noted in the findings obtained in the studies that reported correlations [26,29]. The CTLA4 SNP used in the present study, rs231775, encodes an Ala-to-Thr substitution at the 17th amino-acid position; the CTLA4 protein with threonine 17 more strongly inhibits T-cell activation [40]. Despite this reported functionality, no correlation was found between this SNP and transplant outcomes in the present study. Therefore, CTLA4 may not be associated with BMT outcomes, or, even if it is, the relationship may depend on some factors, such as patient characteristics.
We also examined the relationships between polymorphisms in the methylase genes, DNMT1 and EZH2, and BMT outcomes. The DNMT1 SNP used, rs2228612, was identified as a cancer risk factor in a meta-analysis [50]. However, no correlations were observed between this SNP and post-BMT outcomes in the present study. Although DNA methyltransferase inhibitors are used to treat AML and MDS, azacitidine inhibits not only DNMT1, but also the de novo DNA methylases, DNMT3a and DNMT3b [15]. This is one possible explanation for the lack of a relationship between the DNMT1 SNP and post-BMT outcomes. DNMT3a and DNMT3b SNPs were not included in the present study because we did not detect any SNPs in these genes with sufficiently high MAFs in the Japanese population.
As already discussed, EZH2 is involved in the silencing of MHC-I [25], as well as in GVHD in mice [23]. The EZH2 SNP, rs2302427, encodes the Ezh2 D185H substitution, which is located in the domain that interacts with DNA methyltransferases (Table S2). This SNP has been identified as a risk factor for colorectal cancer in a Chinese population [51]. This SNP was also detected in an AML patient with a morphology resembling acute promyelocytic leukemia without the PARA gene rearrangement [52]. In the present study, recipient rs2302427 exhibited a low p-value in the analysis of grade 2-4 aGVHD (p = 0.010) in the additive model (Table 1) despite its low MAF of 0.09 (Table S3). Therefore, recipient rs2302427 warrants further study with a larger subject size, particularly in a recessive model that was not examined in the present study due to the low MAF ( Table 1). The rs2302427 recessive model will test the GG genotype, which encodes the homozygous Ezh2 H185 protein (Table S2).
The present study had a number of limitations. The study design was retrospective, and clinical decisions should be based on prospective studies. Furthermore, SNP-SNP interactions were not considered. Several of the genes analyzed may be redundant for HSCT outcomes in humans. Therefore, future studies on a larger sample size are needed. In addition, the molecular functions of five of the SNPs are unknown (Table S2).

Conclusions
Although we expected immune checkpoint gene polymorphisms to be useful for predicting BMT outcomes, no correlations were observed. The recipient EZH2 SNP, rs2302427, which encodes the D185H substitution, exhibited a low p-value in analysis of grade 2-4 aGVHD despite a low MAF and, thus, warrants further investigation in a large-scale study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.