A total of 177 patients from nine centres were included utilising a cohort of the Cooperative European Paediatric Renal Transplant Initiative (CERTAIN) registry ( Table 1 ). Patients were representative of a paediatric transplant population: median age 10.8 years (IQR 4.9-14.7), of which 60% were male and 40% had underlying congenital abnormalities of the kidney and urinary tract. Patients represented a low HLA sensitisation group: 82% of patients had panel reactive antibody reactivity of 0%, and 92% were first transplants. 64 (36%) patients had induction treatment with basiliximab and the majority had steroid and calcineurin inhibitor-based maintenance immunosuppression (172, 97%). Median follow-up time was 4.5 (IQR, 3.0-5.0) years. There were two patient deaths at 4.0 and 4.5 years post-transplant due to sepsis, and five allograft failures at 1.5, 2, 3.5, 4.5 and 5 years post-transplant.

As shown in Figure 1 , there was wide variability in the degree of score correlation among molMM methods. Correlation with amino acid sequence polymorphism (as reflected in the AAMS) was highest with the netMHC1k score across all HLA loci examined (r=0.67-0.87, p<0.001), reflecting the low affinity threshold used in the latter method to select polymorphic donor HLA peptides as putative CD4⁺ T-cell epitopes. In contrast, correlation between AAMS and netMHC was poor (r=0.42-0.59, p<0.001) due to the stringent criteria used in the latter method to select HLA class II peptide binders. There was poor to moderate correlation between AAMS and EMS3D for HLA class I (r=0.48-0.53, p<0.001) and class II (r=0.61-0.71, p<0.001) loci respectively, consistent with the structural nature of the EMS3D score. Overall, the correlation among different molMM scores was highest for HLA class II comparisons.

De novo DSA were detected in 56/177 (32%) patients: 9/56 (14%) class I only, 29/56 (52%) class II only, and 18/56 (32%) to both class I and class II. DSA incidence was highest against HLA-DQ (30/56, 54% of DSA positive patients) ( Table S1 ).

For each patient, we first considered each donor mismatched HLA allele individually and used the single HLA molecule method to assess the ability of molMM scores to predict DSA development (10). Model discrimination was assessed using area under receiver operating characteristic (ROC) curves (AUC, Figure S1 ). For transplant pairs with two allelic mismatches at a particular locus, this method utilizes the highest scoring mismatch for risk stratification at the patient level. Table S2 shows the comparison between antMM and molMM scores in predicting locus-specific DSA development at the patient level. Other than for HLA-A DSA, molMM methods had significantly better discrimination in their association with locus-specific DSA compared to antMM. Overall, EMS3D was the most consistent method for DSA prediction with AUCs >0.7 for all loci examined.

Our principal aim was to stratify patients according to risk of DSA development at a specific HLA locus, rather than simply establish the association between HLA incompatibility and alloantibody responses. We therefore used the Youden index from the single molecule analysis ( Table S3 ) and categorised patients as HLA matched, those with low risk HLA mismatches (molMM score <Youden index) and patients with high risk HLA mismatches (molMM ≥Youden index) at a particular locus. Figure 2 presents the results of analyses with DSA as a time dependent variable. For simplicity, and because of its discriminative ability in ROC curve analyses above, patient risk stratification according to EMS3D versus antMM is depicted (similar analyses for other molMM methods are shown in Figure S2 ). Overall, there was a significant association between HLA risk strata, as determined by EMS3D, and risk of DSA development, whereas this was not the case using traditional antigen mismatching, other than for risk strata within the HLA-A locus. More importantly, the incidence of DSA increased significantly in high compared to low EMS3D risk groups for all HLA loci (HLA-A: HR 3.2, 95% CI: 1.1-9.7, high vs. low risk, p<0.05 log rank test; HLA-B: 8.0, 95% CI 1.1-61, p<0.05; HLA-DQ: 5.0, 95% CI 1.2-11.7, p<0.001 and HLA-DR: 2.8, 95% CI 1.1-7.4, p<0.05). In contrast, the cumulative incidence of DSA development was not significantly different in patients with one versus two antigen mismatches for all HLA loci examined. Stratification using AAMS, netMHC and netMHC1k was more variable ( Figure S2 ). Only AAMS stratification showed a significant difference in cumulative DSA incidence against HLA-B mismatches between low and high risk groups (but essentially stratified the majority of patients in the high risk group, n=121), whereas high versus low risk strata were not significantly different for DSA against HLA-A or -B using the other methods. For HLA-DQ mismatches, all molMM methods showed a significant difference in DSA free survival between high versus low risk patient groups. Finally, for HLA-DR mismatches, high versus low risk group comparisons were significantly different for AAMS and netMHC molMM scores, but not for netMHC1k. To exclude a confounding effect of immunosuppression levels on the supp associations, we examined how often patients with and without DSA had CNI levels above 5ng/ml ( Figure S3A ) and found no significant difference between the two groups (other than for high HLA-B risk patients where those with DSA were more likely to have CNI levels above 5ng/ml). We also compared CNI levels in the 6 months preceding DSA development and during the entire duration of follow up in patients without DSA and found no difference; this was the case for patients in the low and high risk groups (other than for low HLA-DQ risk patients where patients with DSA had higher mean CNI levels) ( Figure S3B ). Finally, the intra-patient variability (IPV) of CNI levels was similar in patients who did and did not develop DSA (data not shown).

Notably, in patients with two allelic mismatches within a locus, the DSA targeted the allele with the highest molMM score in the majority of cases (HLA-A: 5/5, HLA-B: 5/7, HLA-DQ: 9/10, HLA-DR: 7/7, denominator denotes the number of patients with two allelic HLA mismatches in the respective locus and one DSA; results were the same for all molMM methods, Figure 3 depicts molMM scores and DSA target using the EMS3D score). In two out of the total of three mis-classified cases the DSA targeted HLA with molMM score above the Youden index.

For TCMR, we focused on late acute rejection, defined as occurring >6 months post-transplant, due to its importance in allograft failure (28). Overall, the incidence of late TCMR (23 patients, 13%) and of ABMR (11 patients, 6%), as diagnosed on ‘for-cause’ biopsies, in this cohort was low. Notwithstanding this limitation (especially for ABMR models), we used multivariable Cox models to investigate the relationship between donor-recipient HLA incompatibility and rejection risk ( Table 2 ). Neither the number of antMM nor any of the molMM scores for individual HLA loci, HLA class I or HLA class II were associated with TCMR risk in this cohort.

In ABMR models, only HLA-DQ mismatching, assessed using the total molMM scores for mismatches in this locus, was significantly associated with ABMR-free survival for AAMS (HR: 1.02, 95% CI: 1.00-1.05 per amino acid mismatch), EMS3D (HR: 25, 95% CI: 1.01-640 per unit increase), and netMHC (HR: 1.08, 95% CI: 1.00-1.16 per donor peptide) ( Table 2 ). However, fitting multivariable models using HLA-DQ low versus high EMS3D risk group molMM stratification did not show a significant association with ABMR (adjusted HR: 2.3, 95% CI: 0.58-9.6 for EMS3D; Figure S4 ). The same was the case after HLA-DQ risk stratification using the other molMM methods (data not shown). Nevertheless, the small number of ABMR cases in this cohort did not enable a robust assessment of risk stratification at the patient level.

Taken together, the above analyses indicate that molecular HLA mismatching was superior to serological HLA mismatching at predicting the risk of DSA development and of ABMR. No single molMM clearly outperformed the others, although a robust comparison of molecular mismatch methods was not possible in this relatively small patient cohort. Overall, EMS3D was the most consistent method for stratifying risk of primary humoral alloimmunity in this patient cohort.

During the study timeframe, only five patients lost their grafts – two were associated with ABMR, one each due to BK nephropathy and recurrence of primary disease, and one not stated. As a surrogate outcome for long-term graft loss, we used eGFR50 defined as persistent 50% reduction in graft function from the baseline eGFR at month 3. This outcome was met in 27 (15%) patients. Based on the above results, we stratified patients using the EMS3D risk threshold from the DSA analysis in a ‘time-zero’ Cox model (i.e. utilising information available at the time of transplantation) adjusted for donor and recipient factors ( Table 3 ). As expected, donor and recipient age, baseline eGFR and re-transplantation were significant predictors of eGFR50 ( Table S4 ). Patients with high risk HLA-DQ mismatches, as assessed using EMS3D, were more likely to reach the eGFR50 endpoint (high versus low risk: adjHR 3.5, 95% CI 1.6-8.2, p<0.005) ( Figure 4A ), whereas molecular mismatching at HLA-A, -B and -DR was not associated with this outcome. Importantly, analysis based on serological HLA mismatches did not show any significant associations with eGFR50-free survival for any HLA-locus ( Figure 4B shows the adjusted survival curve for HLA-DQ). Of note, there were more patients stratified into the low risk EMS3D group (110 patients) than patients with zero split antigen mismatches (66 patients). Accounting for IPV of CNI levels did not significantly change the results. To assess whether the effect of high risk HLA-DQ mismatching on the eGFR50 outcome is no longer apparent when graft rejection after transplantation is accounted for, we next included rejection episodes as time-varying events in the above ‘time-zero’ multivariable model. As might be expected, TCMR and ABMR were significantly associated with graft function deterioration. Nevertheless, the association of high risk EMS3D HLA-DQ mismatching with eGFR50 remained significant (high versus low risk: adjHR 2.9, 95% CI 1.2-7.0, p<0.05) in this model, suggesting an independent effect beyond rejection episodes detected on ‘for-cause’ graft biopsies ( Table S5 ).

To delineate how molMM might enhance patient risk stratification, we next examined the relationship between traditional antigen HLA-DQ mismatching versus molecular HLA-DQ mismatching and graft function deterioration. Within each level of antigen mismatch, there was a wide range of EMS3D scores ( Figure 4C ). We stratified each antigen mismatch into low or high risk, based on the EMS3D Youden index derived from the DSA analysis above, to generate six strata, i.e. two strata (low- or high-EMS3D) per HLA-DQ antigen mismatch (0–2) ( Figure 4D ). Patients with low risk HLA-DQ mismatches had equivalent eGFR50-free survival regardless of the number of antigen mismatches present (0-MM/low-EMS3D, 1-MM/low-EMS3D and 2-MM/low-EMS3D clustered together, p>0.1), whereas patients with high-risk antigenic mismatches had 3-10 fold higher odds of graft deterioration (using the 0-MM/low-EMS3D group as reference) ( Figure 4D ). Patients with 0-MM/high-EMS3D had worse eGFR50-free survival compared to patients with 2-MM/low-EMS3D. Taken together, the above results indicate that patient risk for graft function deterioration can be stratified more accurately using molecular HLA-DQ mismatching rather than antigenic mismatching.

We further compared the performance indices of risk stratification models for the eGFR50 outcome using EMS3D (high and low risk) versus antMM. Model discrimination (the performance of the multivariable ‘time-zero’ model in correctly ranking patients according to risk of eGFR50 outcome) was assessed using time-dependent AUC. Stratification with EMS3D had consistently higher AUCs throughout the study period (81% at 2 years, and 70-73% after 2 years post-transplant) compared to antMM (77% at 2 years, reducing to 64-65% after 3 years post-transplant) ( Figure 5A ). Figure 5B shows the model calibration plots, derived using leave-one-out cross-validation, for prediction of the eGFR50 outcome at 4-years post-transplant. The 45-degree line represents a match between the predicted risk and actual risk. Overall, both models had similar calibrated risk prediction for patients at low risk but the EMS3D model had better calibration performance for patients at higher risk of graft function deterioration.

Patients were recruited from the web-based CERTAIN Registry, which is an international research registry with funding for patient entry (52). To facilitate recruitment, data is exchanged bi-directionally with Eurotransplant and ERA-EDTA Registry for centres within Europe. The registry collects clinical and laboratory data of paediatric kidney transplant recipients in a detailed and comprehensive manner every 3 months in the first year and six-monthly thereafter. This allows a thorough patient characterization in general and also of specific patient subcohorts (52). Written informed consent was obtained from all parents/guardians to participate in the registry, with assent from patients when appropriate for their age. The CERTAIN Registry has been approved by the respective ethics committee of each contributing centre and is kept in full accordance with the principles of the Declaration of Helsinki and Good Clinical Practice guidelines. The study was designed, analysed and reported according to the STROBE guidelines (https://www.strobe-statement.org).

We included patients from centres which routinely performed prospective monitoring of post-transplant HLA antibodies, at minimum within the first year post-transplant and yearly thereafter. No exclusion criteria were applied. HLA antibody testing was performed following local protocols, using screening panels followed by single antigen HLA specific beads. HLA antibody positivity was defined using local thresholds ( Table S6 ). Only de novo donor-specific HLA antibodies (DSA) were included in this study, and HLA specificity was assigned at the allelic level by comparing donor and recipient HLA allelic types and referencing results from individual allelic beads within single antigen panels. Histology reports were verified by an independent histopathologist based on the Banff 2017 criteria (53). The CERTAIN Registry was accessed at 15/05/2020. The data were not sufficient to perform DSA analysis for HLA-C, but otherwise had no missing data.

HLA typing in the registry was provided for HLA-A, -B, -DR and -DQ at the split serotype level along with donor/recipient ethnicity. Allelic imputation was performed using population specific haplotype frequency data from the National Marrow Donor Program for HLA-A, -B, -C, -DQA1, -DQB1, -DRB1 and DRB3/4/5 loci (54, 55). The probable haplotype pair was assigned for each donor/recipient genotypes.

Amino Acid Mismatch Score (AAMS) and Electrostatic Mismatch Score 3D (EMS3D) calculations: The number of amino acid mismatches present on donor HLA was calculated using the Cambridge HLA Immunogenicity algorithm, as previously described (7, 8). Similarly, electrostatic potential differences on the surface of donor compared to recipient HLA were calculated using the EMS3D algorithm, as previously described (7, 8, 20–22, 56).

NetMHCIIpan (version 4.0) was used to predict the affinity and percentage rank of all potential 15-mer peptides present in donor HLA extracellular sequences (which were absent from recipient sequences), using recipient HLA-DRB1, HLA-DRB3/4/5 and HLA-DQ as presenting molecules. The resultant 9-mer core, affinity and percentage rank was used to enumerate the number of unique 9-mer cores from all potential peptides per donor HLA that either met a 1000 nM threshold (abbreviated netMHC1k) or both of a 500 nM threshold and were in the top 2% percentage rank (abbreviated netMHC). The former (1000 nM) peptide threshold is that employed by the PIRCHE algorithm (6, 24).

We analysed the ability of the molecular HLA mismatch scores to predict de novo locus-specific DSA development using the ‘single molecule’ method (10). Briefly, each donor HLA mismatch was scored and correlated to DSA development against that molecule specifically, independent of patient specific effects. Thresholds were derived from the Youden index using receiver operating characteristic (ROC) curve analysis for each HLA locus and high/low molMM categories were defined.

We then applied the molecular mismatch scores to predict DSA development at the patient level. In patients with two HLA mismatches within a locus, the highest scoring mismatch was used to represent the patient’s molMM score. We assessed the ability of molMM scores to predict locus-specific DSA development at the patient level and used the area under ROC curve (AUC) to compare the discrimination ability of each molMM score. Pair-wise comparisons between models were performed using the “DeLong” method in pROC package (version 1.18.0) in R (57).

We used cumulative event curves and Cox proportional hazards models for time to event analyses. Patients were right-censored at time of last follow-up or patient death. We performed univariate analysis for each locus-specific DSA outcome using the molecular risk thresholds and compared groups using log-rank tests. For TCMR and ABMR, we performed multi-variate analysis with molMM scores as the main predictor and adjusted for recipient age, recipient ethnicity, graft number, cold ischaemic time and panel reactive antibodies (%), defined a priori. Patient adherence was assessed using the surrogate marker, intra-patient variation (IPV) in immunosuppression levels (28). All patients were prescribed tacrolimus except four patients using ciclosporin. CNI levels were taken from month 6 post-transplant and were available every 6 months. Due to the limited number of measurements, CNI levels for the whole duration of follow-up was used to calculate the IPV using the mean absolute deviation. A threshold of 25% was used to categorise patients into high/low IPV as previously published (28). We also analysed CNI levels in the months preceding DSA formation. For the four patients on ciclosporin, the level was divided by 10 to approximate tacrolimus levels. We compared CNI levels between positive and negative patients for each molMM risk group (high/low) using the Mann Whitney test. We also assessed the number of times the CNI level was above 5ng/ml for patients with and without DSA, as described by Wiebe et al. (31).

Graft function deterioration was modelled using the eGFR50 outcome, defined as persistent 50% reduction in graft function (on two consecutive results six months apart) from the baseline eGFR taken at month 3. As the aim of the study was to assess the predictive ability of HLA mismatching, we first constructed a base model using donor and recipient characteristics ( Table S4 ). Four variables were kept in the step-wise multiple regression analysis: baseline eGFR, recipient age, donor age and transplant number, which were used as covariates for further analyses of eGFR50 outcome. Statistical analyses were performed using ‘survival’ package (version 3.2-13) in R (version 4.1.2) (58, 59) Model accuracy and discrimination were analysed using time-dependent AUC and calibration plots. Results were validated using leave-one-out internal cross-validation with 200 bootstrap iterations (30) P-values of <0.05 were considered statistically significant and results were considered exploratory as multiple models were fitted.