All recipient and donor blood samples were collected before HCT according to research protocols approved by the Fred Hutchinson Cancer Research Center (FHCRC) Institutional Review Board (IRB) or the National Marrow Donor Program (NMDP). Project-specific IRB approval was obtained for the use of clinical data and research biospecimens.

The FHCRC cohort included 1868 HLA-A, B, C, DRB1, DQA1, DQB1, DPA1, DPB1-matched donor-recipient sibling pairs with mostly European ancestry (~83%) who had a first allogeneic HCT with marrow or growth factor-mobilized blood cells at the FHCRC and Seattle Cancer Care Alliance from 1990 through 2011. Recipients treated with T-cell depleting antibodies or high-dose cyclophosphamide after HCT and those with syngeneic or cord blood donors were excluded. A single prior autologous HCT was allowed. Conditioning regimens were categorized as myeloablative or nonmyeloablative according to the intensity of chemotherapy and total body irradiation. Indications for HCT included hematological malignancy or myelodysplasia. Donors and recipient 4-digit typing of HLA-A, B, C, DRB1, DQB1, DPA1 and DPB1 alleles was determined as described previously (10). We used a cohort from the Center for International Blood and Marrow Transplant Research (CIBMTR), a research collaboration between the NMDP and the Medical College of Wisconsin to test replication of HLA-A*02:01 mHAg GWAS associations discovered in the FHCRC cohort. The CIBMTR cohort consisted of 838 HLA-A, B, C, DRB1, DQB1-matched sibling pairs with self-identified European ancestry and at least 1 HLA-A*02:01 allele. Patients in this cohort had HCT from 2008 through 2018.

Details regarding preparation of genomic DNA samples from donors and recipients for the FHCRC cohort and the related DNA amplification, genotyping platforms, hybridization, genotyping, imputation algorithms, quality control and quality assurance have been described previously (10). Targeted sequencing of variants selected for replication using the CIMBTR sample pairs was done by the Genomics & Bioinformatics shared resource at FHCRC using the AmpliSeq PCR workflow on the MiSeq platform (Illumina). Joint genotyping was done using the GATK pipeline 4.1.8.1 following the best practice workflow (14–16). Briefly, picard 2.18 (17) was used to process Fastq and Bam files to generate analysis-ready reads. BWA 0.7 (18) was used to map paired-end reads to reference genome hg19, and GATK 4.1.8 was used to call variants for each sample (BaseRecalibrator, ApplyBQSR and HaplotypeCaller in GVCF mode) followed by joint-calling (GenomicsDBImport and GenotypeGVCFs). Quality assurance, quality control and variant filtering followed recommendations for targeted sequencing from the GATK website (19–21). We further removed variants with high allelic ratios, ancestry outlier samples identified by clustering in the space of the first two principal components, male samples with high X chromosome heterozygosity and sample pairs whose observed identity-by-descent relationship did not match the expected full-sibling relationship. Whenever possible, HLA-A*02:01 replication testing with CIBMTR samples was done with 2 different variants representing each discovery in case the primary assay failed.

The 6 outcomes tested were acute GVHD categorized as peak grade 2-4, 2b-4 and 3-4 severity (22, 23), stage 2-4 gut GVHD, chronic GVHD and recurrent or progressive malignancy (i.e., relapse). Grade 2b-4 acute GVHD excludes isolated stage 1 gut GVHD, which is frequently recognized at FHCRC (24). Chronic GVHD was diagnosed according to historical criteria (25) because evaluation according to NIH criteria was not available for many patients.

Variants were analyzed in 2 categories according to whether they are predicted to encode differences in the amino acid sequence of an annotated protein in the recipient when compared to the donor, hereafter termed “coding” variants, or not (hereafter termed “noncoding” variants). The coding category included variants whose alleles were predicted to have missense, inframe insertion or deletion, frameshift, start lost, stop gained or lost, stop retained, incomplete terminal codon, or protein altering effect on the overlapping transcript as predicted by the Ensembl Variant Effect Predictor (VEP) (26). The annotations were extracted from the Annotation Explorer application hosted on Biodata Catalyst (https://biodatacatalyst.nhlbi.nih.gov/) (27).

For each variant, recipient mismatching was evaluated separately for the major allele and the minor allele. For a variant with major and minor alleles “a” and “b”, pairs with “bb” donors and “ab” or “aa” recipients were categorized as mismatched for the major allele, and pairs with “aa” donors and “ab” or “bb” recipients were excluded to avoid confounding by mismatching for the minor allele. Similarly, pairs with “aa” donors and “ab” or “bb” recipients were categorized as mismatched for the minor allele, and pairs with “bb” donors and “aa” or “ab” recipients were excluded to avoid confounding by mismatching for the major allele. Recipient allele mismatch associations (RAMAs) for each outcome were based on cause-specific hazard ratio analysis using Cox regression comparing pairs with recipient mismatching versus those without mismatching for a given major or minor allele, treating death as a competing risk for all endpoints and relapse as a competing risk for acute and chronic GVHD.

In the discovery phase of the HLA-A*02:01 analysis, RAMAs were tested by proportional hazards analysis in the FHCRC cohort using 824 HLA-A*02:01-positive donor-recipient pairs with at least 1 HLA-A*02:01 allele and in 929 HLA-A02 supertype-negative donor-recipient pairs who did not have HLA-A*02:01 or other HLA-A02 supertype alleles (A*02:XX, 68:02, 68:15, 68:57, 68:28 or 69:01) (28). Coding RAMAs with likelihood ratio p-values ≤.01 in HLA-A*02:01-positive pairs and p-values ≤.01 for A*02:01 interaction in a combined HLA-A*02:01-positive and HLA-A02-negative analysis were identified as discoveries that could be attributed to HLA-A*02:01-specific mHAgs. The number of noncoding RAMAs was much larger than the number of coding RAMAs. Therefore, the threshold p-values for likelihood ratio tests and interaction tests used to identify noncoding discovery candidates were set at ≤10^-4. Discovery RAMAs were tested in the CIBMTR replication cohort with Bonferroni adjustments for multiple comparisons. For Bonferroni adjustment, acute GVHD, chronic GVHD and relapse were treated as separate analyses, with grade 2-4, 2b-4, 3-4 acute GVHD and stage 2-4 gut GVHD considered as a single category. In addition, discoveries with expected hazard ratios >1.0 for acute and chronic GVHD and <1.0 for relapse were analyzed separately from discoveries with unexpected hazard ratios <1.0 for acute and chronic GVHD and >1.0 for relapse.

In a separate MHC-agnostic analysis, we screened for RAMAs independent of whether the donor and recipient pairs had any specific HLA allotype. For this purpose, the FHCRC cohort was randomized at a 3:2 ratio into discovery and replication cohorts with 1125 and 743 pairs, respectively. We initially identified coding RAMAs with likelihood ratio p-values <1.0 x 10^-3 and noncoding RAMA with p-values <1.0 x 10^-5 as discovery candidates. Because the number of candidates was too large to tolerate Bonferroni adjustment for multiple comparisons, we used p-values <1.0 x 10^-4 for coding RAMAs and 1.0 x 10^-6 for noncoding RAMAs as more stringent thresholds to prespecify discovery candidates for replication testing.

European ancestry linkage disequilibrium (LD) groups of genomic variants with pairwise r² ≥0.70 were identified by the LDmatrix tool (ldlink.nci.nih.gov) (29). The LDproxy tool (ldlink.nci.nih.gov) was used to identify external missense proxy variants with r² ≥0.70 in the European population to be included in linkage groups for analysis of HLA-A*02:01 RAMAs and to identify external proxy variants predicted to encode non-annotated open reading frames for analysis of MHC-agnostic RAMAs. When necessary, the LDproxy tool was also used to identify proxies with r² ≥0.95 that could be used as backup assays for the primary variant being tested for replication in the CIBMTR cohort. Peptides up to 13 residues upstream and downstream from the amino acid encoded by genomic variants (i.e., source peptides) were recovered from Haplosaurus (ensembl.org) (30), the dbSNP database (ncbi.nlm.nih.gov/snp) (31) or from manual reading of the UCSC genome browser (genome-euro.ucsc.edu) (32). Some source peptides were shorter than 27 residues due to start or stop codons. The source peptides were then analyzed for predicted binding of 8 to 14-mer processed peptides to selected MHC class I allotypes by NetMHCpan4.1 at the default setting (33). The binding predictions are based on both the estimated binding affinity percentile rank and on publicly available data from mass spectrometry analysis of peptides eluted from HLA molecules. Processed peptides predicted to bind HLA-A*02:01 but not containing the variant residue were culled from further analysis.

The initial discovery analysis identified 604 variants associated with an outcome in the HLA-A*02:01-positive cohort and having a statistical interaction with the absence of HLA-A*02:01 in the combined HLA-A*02:01-positive and negative cohorts, indicating that these associations could be attributed to HLA-A*02:01-specific mHAgs ( Figure 2 and Table S1 ). Among the 604 discoveries, 197 were coding variants and 407 were non-coding variants. To characterize noncoding variants more accurately according to the presence or absence of any associated coding variants, we identified 82 external missense proxies with r² ≥0.70 as additional RAMA discoveries. Of the 686 variants, 527 were combined into 93 linkage groups with pairwise r² ≥0.70, and 159 did not have r² ≥0.70 for association with any other variant in the data set. The 93 linkage groups contain between 2 and 105 members. Among the 686 RAMAs, 279 involved a variant that encodes at least 1 peptide, yielding a total of 284 unique source peptides ( Table S2 ).

To validate the use of NetMHCpan4.1 (33) to identify source peptides that generate processed peptides predicted to bind to HLA-A*02:01, we tested its performance with a set of 49 previously identified genomic variants reported to encode mHAgs presented by HLA-A*02:01 (4, 8). The NetMHCpan4.1 algorithm predicted HLA-A*02:01 binding for all but 3 of the source peptides tested ( Table S3 ). In 4 cases, the predicted peptide differed from the reported peptide, but the NetMHCpan4.1 algorithm predicted alternatives that approximated the reported peptide. Although this analysis could not determine the false-positive rate, the results show that the NetMHCpan4.1 algorithm has a low false-negative rate in predicting the binding of processed peptides to HLA-A*02:01.

Of the 284 unique source peptides identified in our discovery analysis, 102 have at least 1 processed peptide that i) contains 8 to 14 residues, ii) includes the amino acid encoded by the genomic variant, and iii) is predicted to bind to HLA-A*02:01 ( Figure 3 and Table S4 ). These 102 source peptides produce a total of 188 unique processed 8 to 14-mer peptides that are predicted to bind to HLA-A*02:01. Of these, 47 have a 9-mer HLA-A*02:01-binding pattern or sequence logo characterized by having leucine, isoleucine, or methionine at position 2 and valine, leucine, isoleucine, or methionine at position 9 (34, 35).

For purposes of replication, 372 RAMAs were removed from the list of 686 discoveries and external proxies, so that no linkage group contained more than 2 variants, and 137 assay proxies with r² >0.95 were added to enable backup testing of singleton RAMAs that did not fit within any linkage group. A total of 49 RAMAs were not evaluated due to lack of an appropriate assay or because the results did not pass QC after testing ( Figure 4 and Table S5 ). Grade 2b GVHD could not be ascertained in the CIBMTR cohort because gut staging was not available. Due to differences in the cumulative incidence frequencies of acute GVHD in the discovery and replication cohorts, we tested grades 2-4 GVHD in the CIBMTR cohort for replication of grades 2b-4 GVHD discoveries in the FHCRC cohort. We also tested replication for grades 3-4 GVHD, chronic GVHD and relapse, and excluded 117 RAMAs for grade 2-4 GVHD and stage 2-4 gut GVHD from replication testing.

Of the 285 RAMAs tested for replication, 121 represented backup assays within a linkage group, leaving 164 independent genomic signals (IGS). Among the 164 IGS, 83 (51%) had hazard ratios >1.0 as expected for acute and chronic GVHD or <1.0 as expected for relapse, while the remainder had unexpected hazard ratios <1.0 for acute or chronic GVHD or >1.0 for relapse. Each IGS was characterized according to whether any member of its linkage group encoded a source peptide, whether any processed peptides were predicted to bind to HLA-A*02:01, and whether any processed peptides had an HLA-A*02:01 sequence logo as defined above. Of the 164 IGS, 133 contained at least one coding RAMA. Among the 133, 57 contained at least one RAMA that was predicted to produce processed peptides that bind to HLA-A*02:01, and 30 of the 57 contained at least one RAMA that was predicted to produce a processed peptide with an HLA-A*02:01-binding sequence logo.

Among the 30 IGS associated with at least 1 HLA-A*02:01 sequence logo peptide, 16 had hazard ratios >1.0 as expected for acute and chronic GVHD and <1.0 as expected for relapse. These IGS with expected hazard ratios were prespecified to be evaluated formally for replication with statistical correction for multiple comparisons. These pre-specified IGS all had p-values >.05 for replication ( Table 2 ). We also evaluated the 14 IGS that were associated with at least 1 HLA-A*02:01 sequence logo peptide but had unexpected hazard ratios. Likewise, these IGS all had p-values >.05 for replication ( Table 2 ).

One IGS result that was not pre-specified drew scrutiny because the replication p-values for the primary and backup assays were 0.001 and 0.004, respectively. The replication hazard ratios for the association of recipient minor allele mismatching of rs62572859 and rs12554984 with the risk of grade 3-4 acute GVHD were 0.12 [95% confidence interval (CI), 0.02-0.84)] and 0.14 (95% CI, 0.02-0.99), respectively. In addition to the implausibility of the direction, the hazard ratios reflect a result based on an expectation of only 7 events in a few dozen patients (10% incidence of grade 3-4 acute GVHD, 9% probability of recipient allele mismatching for these variants, and ~800 genotyped pairs).

With the test for statistical interaction as a criterion, the 604 discovery RAMAs and 82 external missense proxies in Table S1 included only 2 of the 49 previously identified variants that produce HLA-A*02:01-associated mHAgs (rs9051-G and rs9051-C) in Table S3 . Very few of these 49 HLA-A*02:01-associated mHAgs have been tested for association with HCT outcomes. Therefore, we evaluated these alleles as candidates for association with all 6 HCT outcomes. For this analysis, the FHCRC HLA-A*02:01-positive cohort was randomized at a 3:2 ratio into discovery and replication cohorts with 486 and 338 donor and recipient pairs, respectively. Associations with p-values ≤.05 in the discovery cohort were tested in the replication cohort. The discovery analysis identified 17 RAMAs, 6 with expected hazard ratios and 11 with unexpected hazard ratios ( Tables S6 ). These RAMAs all had p-values >.05 for replication ( Table 3 ). Table S7 summarizes result for the association of the 49 RAMAs with the 6 outcomes in the combined FHCRC discovery and replication cohorts to enable future meta-analysis with results from other cohorts.

Our study was designed to identify mHAgs that can affect outcomes when presented by HLA-A*02:01 but not when presented by other HLA allotypes. Peptides that bind HLA class-1 molecules, however, show promiscuous binding to multiple allotypes both within and across supertypes (36–39). The extent to which the peptides presented by HLA-A*02:01 could bind to other HLA supertypes has not been determined. To address this question, we used the PromPDD algorithm (40) to evaluate the 42 9-mer and 10-mer peptides representing the 49 known HLA-A*02:01-associated mHAgs for binding to other HLA-A and B allotypes with allele frequencies >.05 in the USA NMDP European Caucasian population (A*01:01, A*03:01, A*11:01, A*24:02, B*07:02, B*08:01, B*15:01, B*35:01, B*40:01, B*44:02), none of which is included in the A02 supertype (41, 42). Approximately 80% of the individuals in this data set have at least 1 of the 4 HLA-A alleles or HLA-A*02:01, and approximately 66% have at least 1 of the 6 HLA-B alleles. As a positive control, the algorithm predicted that 34 (79%) of the 42 tested peptides bind to HLA-A*02:01 ( Figure 5A and Table S8 ). Ten (24%) of the 42 peptides were predicted to bind to HLA allotypes outside the A02 supertype, including HLA-A*24:02 (n =1), HLA-*B08:01 alone (n = 4) and both HLA-B*08:01 and HLA*B44:02 (n = 5).

Minor histocompatibility antigens originate from source proteins or peptides that are processed through the MHC class I and class II antigen processing pathways for binding to MHC molecules (1, 2). To address the possibility that the source peptides that produce HLA-A*02:01-associated mHAgs could also produce processed peptides predicted to bind to other HLA allotypes, we used the NetMHCpan4.1 algorithm at the default setting to identify 9-mer peptides that include the variant amino acid residue and were identified as strong binders having a predicted affinity <500 nM for binding to the same HLA-allotypes listed above. Twenty-six (53%) of the 49 tested source peptides produced processed peptides that met these stringent specifications ( Figure 5B and Table S9 ). Elution scores ranged from 0.33 to 0.96, indicating a high likelihood that these peptides have been detected in mass spectrometry studies. Binding score percentile ranks ranged from 0.01 to 1.03, indicating strong binding to the respective HLA allotype, and predicted binding affinities ranged from to 6.5 to 383 nM ( Table S9 ). At least one processed 9-mer peptide was predicted to have strong binding to HLA-A*03:01 (n =5), HLA-A*11:01 (n = 6), HLA-B*07:02 (n = 6), HLA-B*08:01 (n = 4), HLA-B*15:01 (n = 4), HLA-B*35:01 (n = 1), HLA-B*44:02 (n = 9), and HLA-B*44:02 (n = 2) ( Figure 5B ). Four processed peptides were predicted to bind to 2 of the 10 tested HLA allotypes, 2 were predicted to bind to 3 of the 10 tested allotypes, and 1 was predicted to bind 4 of the 10 tested allotypes. No processed 9-mer peptides that include the variant amino acid were predicted to have strong binding to HLA-A*01:01 or HLA-A*24:02.

The test for statistical interaction with HLA-A*02:01 as a criterion for discovery eliminated many possible mHAgs from consideration. To address this limitation, we used an MHC-agnostic approach to identify mHAgs associated with outcomes in the FHCRC cohort (10). The initial screen identified 802 variants, 403 with expected hazard ratios and 399 with unexpected hazard ratios ( Table 4 and Table S10 ). Pruning for linkage disequilibrium (LD) yielded 285 IGS, 142 with expected hazard ratios and 143 with unexpected hazard ratios. By using 10-fold smaller p-value thresholds, we prespecified 28 signals to be tested for replication in categories defined according to whether the discovery hazard ratio was expected or not and whether the IGS included coding variants or not. None of the 17 prespecified signals with unexpected hazard ratios met criteria for replication ( Table S11 ). Of the 11 prespecified signals with expected hazard ratios, 1 signal met criteria for replication with Bonferroni adjustment for multiple comparisons, and 1 other signal came close to meeting criteria ( Table S11 ). Table S12 summarizes result for the association of the 802 MHC-agnostic RAMAs with the 6 outcomes in the combined FH discovery and replication cohorts to enable future metanalysis with results from other cohorts.

The prespecified signal that met criteria for replication contained 2 variants, rs56040842 and rs61851500. The replication hazard ratios for the association of recipient minor allele mismatching of rs56040842 and rs61851500 with the risk of grade 2b-4 acute GVHD were 3.38 (95% CI, 1.9-6.1; p = .0006) and 2.71 (95% CI, 1.5-4.8; p = .004), respectively. These hazard ratios reflect results based on an expectation of only 9 events in a score of patients (40% incidence of grade 2b-4 acute GVHD, 3% probability of recipient allele mismatching for these variants, and ~750 genotyped pairs). According to LDproxy, the signal defined by recipient minor allele mismatching of rs56040842 and rs61851500 contains only 1 other variant in LD with pairwise r² ≥ 0.7 (rs61851533), and all 3 variants are in a PRKG1 intron. Because intron sequences can encode non-annotated proteins (43–45), we used the ExPASy Translate tool (web.expasy.org/translate) to identify open reading frames that include the variant allele and encode an amino acid sequence that differs from the major allele. With this process we identified 4 such open reading frames ( Table S13 ), none of which map to an annotated coding region as determined by the BLATP tool (ensembl.org). We used NetMHCpan4.1 to determine whether any of the potential source peptides are predicted to produce processed 8-11mer peptides with elution scores >0.05 and predicted percentile ranks ≤2.0 for binding to highly prevalent HLA-A and B allotypes as described above. This screening identified 27 peptides that are collectively predicted to bind to a range of HLA-A and B allotypes, including HLA-A*02:01, HLA-A24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*15:01, HLA-B*39:01 and HLA-B*58:01 ( Figure 6A and Table S14 ).

One other variant that came close to meeting criteria for replication was rs17811627. The signal containing this variant included 4 other variants in the database, all having pairwise r² >0.82 with the others. For rs17811627, the replication p-value for the association of mismatching for the minor allele with grade 3-4 GVHD was 0.013, whereas the Bonferroni adjusted threshold was 0.125. Replication p-values for the other 4 variants in the signal ranged from.02 to.04. The replication hazard ratio for rs17811627 was 1.84 (95% CI, 1.2-2.9).

According to LDproxy, the signal containing rs17811627 includes 38 variants having pairwise r² ≥0.70 with the others, all located in an intergenic region spanning 46 Kb on chromosome 15. Because intergenic regions can encode non-annotated proteins (45), we used the process described above to identify 30 open reading frames that include the minor allele and encode an amino acid sequence that differs from the major allele ( Table S15 ), none of which map to an annotated coding region as determined by the BLATP tool. NetMHCpan4.1 screening identified a total of 151 8-11mer peptides each with predicted binding to at least 1 of the 12 HLA-A and B allotypes tested ( Figure 6B and Table S16 ). Each tested HLA-allotype was predicted to bind between 13 and 28 peptides, and 11 of the 151 peptides were predicted to bind between 3 and 5 of the HLA-allotypes tested.