Abstract
Kidney transplantation (KT) clinical outcomes are highly variable across patients and would benefit from predictive biomarkers to achieve personalized/precision medicine. The B cell activating factor (BAFF) system signaling plays an essential role in B lymphocytes' homeostasis, and is implicated in activation and survival of B lymphocytes. Single nucleotide polymorphisms (SNPs) in BAFF system genes are therefore strong candidates to identify the genetic mechanisms underpinning variable clinical outcomes in KT. We report here new findings on BAFF system genetic polymorphisms in KT patients in relation to two key phenotypes of clinical interest: graft survival and acute rejection (AR). A total of 168 KT patients, of which 29 suffered AR, participated in this study. The BAFF system polymorphisms in five genes TNFSF13B, TNFSF13, TNFRSF13C, TNFRSF13B, and TNFRSF17 were characterized using TaqMan SNP genotyping. Patients with KT who had an AA genotype in polymorphism rs3803800 of the TNFSF13 gene had a higher risk of suffering AR (p = 0.046; odds ratios = 3.38, 95% CI: 1.02–11.2). Moreover, patients with AA genotype (rs3803800) in the TNFSF13 gene had a significantly lower AR-free time than the GG/GA genotypes (69.2% vs. 85.7%; p = 0.037). Of importance, bioinformatics analysis showed that the polymorphism rs3803800 could alter splicing regulation and affect the proliferation-inducing ligand (APRIL) expression levels. The analysis of graft survival did not show a significant association with the polymorphisms analyzed in this study. In conclusion, the rs3803800 genetic polymorphism from this study of BAFF system genes appears to display importance in AR-free time for KT patients, and thus, warrants further research in independent populations as a putative predictive biomarker of AR. These findings also inform future personalized/precision medicine efforts and functional genomic studies in KT patients.
Introduction
Over the last decade, biomarker discovery science has focused on advancing our ability to forecast the highly variable outcomes in renal allograft monitoring. New technologies in omics systems science have contributed to the breakthroughs in diagnostic and biomarker innovation (Stapleton et al., 2018) and are often regarded as the cornerstone of personalized/precision medicine. The latter is a proactive strategy to detect and prevent pathological processes by identifying and enacting on earlier and comprehensive information that forecasts clinical outcomes (Naesens and Anglicheau, 2018).
Kidney transplantation (KT) clinical outcomes are highly variable across patients and would benefit from predictive biomarkers to achieve personalized/precision medicine. Of note, B cells play an important part in the biological processes that mediate renal graft rejection and tolerance through a variety of effector mechanisms, including antibody production, antigen presentation to T lymphocytes, and cytokine release (Luu et al., 2014). B lymphocytes have considerable functional and phenotypic diversity, embodying several subpopulations with specialized effector functions.
During the development and proliferation of B lymphocytes, the availability of cytokines that provide survival signals, such as IL-7 or B cell activating factor (BAFF), is crucial for the correct immune response mediated by lymphocytes B (Harris et al., 2000). The components of the BAFF system have been proposed as biomarkers of kidney graft dysfunction and alloantibody production. B cell depletion by agents such as alemtuzumab or rituximab has been associated with increased antibody production and antibody-mediated rejection (AMR) (Chhabra et al., 2013). The latter patients had higher post-transplantation levels of BAFF than the patients with a control intervention, hypothesizing that high levels of BAFF may produce a dysregulation in the microenvironment that promotes the expansion and activation of alloreactive B cells (Chhabra et al., 2013; Hoffmann et al., 2015).
Therefore, the BAFF system is a strong candidate pathway to forecast post-transplantation outcomes and inform immunosuppression changes (Galián et al., 2016; Thibault-Espitia et al., 2012).
BAFF and the APRIL (a proliferation-inducing ligand) are cytokines belonging to the tumor necrosis factor (TNF) family that possess a homotrimeric type II transmembrane structure (Luu et al., 2014), which are proteolytically processed by furins proteases in consensus sequences generating the soluble forms. The BAFF system has three receptors: the BAFF receptor (R-BAFF), the transmembrane activator and calcium modulator, and the transmembrane activator and calcium modulator and cyclophilin ligand interactor (TACI), as well as the B cell maturation antigen (BCMA).
BCMA, also known as TNFRSF17, and the TACI receptor, also known as TNFRSF13B, can bind to BAFF and APRIL. R-BAFF, also known as TNFRSF13C, is a protein that BAFF can bind to (Bossen and Schneider, 2006). BAFF and APRIL receptors are mainly expressed on B lymphocytes (Luu et al., 2014).
Numerous studies have investigated single nucleotide polymorphisms (SNPs) of genes related to the functions of B cells, such as the BAFF gene, which is related to the activation and survival of B cells. The rs9514828 polymorphism present in the BAFF gene's promoter region is reportedly associated with BAFF expression levels in an autoimmune disease model (Kawasaki et al., 2007). Although it has not been tested in KT, this polymorphism can have important repercussions as BAFF levels are associated with the development of donor-specific antibodies (DSAs) de novo and AMR (Banham et al., 2013).
We report here new findings on BAFF system genetic polymorphisms in KT patients in relation to two key phenotypes of clinical interest: graft survival and acute rejection (AR).
Materials and Methods
Demographic information, clinical features, and the study design
The study followed the Declaration of Helsinki, and the HCUVA protocol (PI15/01370 and PI19/01194) was approved by the Ethics Committee. Before taking part in the trial, all patients gave written informed consent.
At the University Clinic Hospital “Virgen de la Arrixaca” (Spain), a total of 168 consecutive adult KT patients were retrospectively studied. DSA Luminex findings, DNA storage, follow-up over 5 years post-transplantation, and complete clinical data were only included in this investigation for patients whose kidney graft had been in function for at least 1 month after transplantation. Return to dialysis was considered as allograft loss.
Table 1 shows demographic information as well as clinical features. All transplant patients had their estimated glomerular filtration rate (eGFR) and creatinine measured (normal values between brackets): creatine (0.7–1.2 mg/dL) and eGFR (>90 mL/min/1.73 m2). Before the transplant, the patients had the following values: creatinine (mg/dL; 2.92.1; mean SD) and an eGFR <60 mL/min/1.73 m2 for more than 3 months, which indicates chronic kidney disease.
Demographic and Clinical Characteristics of the Kidney Transplant Patients Included in This Study
Boldface = values of p <0.05 considered statistically significant.
Total differences in HLA-A. HLA-B and HLA-DR between donor and recipient. Comparisons made using Fisher's exact test or χ2 as appropriate.
The results were expressed as mean ± SD.
AR, acute rejection; HLA, human leukocyte antigen; KTRs, kidney transplants recipients; SD, standard deviation.
The main KT indications were glomerulonephritis (34.2%), polycystic kidney disease (20.3%), type I diabetes mellitus (11.9%), chronic obstructive pyelonephritis (8.4%), unknown renal insufficiency (6.1%), lupic nephritis (3.6%), reflux nephropathy (2.4%), and others (1%).
Twenty-nine patients suffered acute rejection (AR), 15 suffered acute cellular rejection (ACR), and 14 suffered AMR. The mean age of the group recipients without acute rejection (NAR) was significantly higher than that in the AR group (57.2 ± 10.9 and 52.5 ± 13.0 years, respectively; p = 0.043). Regarding the human leukocyte antigen (HLA) mismatches between donor and recipient, the NAR group presented significantly fewer differences than those of the AR group (4.05 ± 1.1 vs. 4.67 ± 0.8; respectively; p = 0.004). No significant differences were observed in the variables sex, cold ischemia time, and primary disease.
Immunosuppressive treatment
Oral tacrolimus (Prograf, Astellas, Ireland), mycophenolate mofetil (MMF, CellCept; Roche, Switzerland), and prednisolone were given to all participants (Dacortin, Merck, Spain). The tacrolimus (FK)-based protocol was started at a dose of (0.10–0.15 mg/kg/day), and the dose was adjusted to maintain a trough level of FK in whole blood between 8 and 12 ng/mL during the first month after transplant, between 7 and 10 ng/mL during the second and third months after transplant, and between 5 and 8 ng/mL after that. MMF was begun at a dose of 2000 mg/day, and then gradually reduced to 1000–1500 mg/day during the first month after surgery, depending on white blood cell counts. MMF was begun at a dose of 2000 mg/day, and then gradually reduced to 1000–1500 mg/day during the first month after surgery, depending on white blood cell level.
On the day of transplantation, day 12, and day 34 after the procedure, methylprednisolone was given intravenously at doses of 500, 250, and 125 mg/day, respectively. Oral prednisolone was started at 20 mg on day 5 after the operation and gradually decreased to 510 mg/day after 2–3 months. Some instances were treated with induction treatment based on thymoglobulin or basiliximab based on the immunological risk before transplantation.
Kidney rejection diagnosis
Allograft ACR was defined as a 20% increase in serum creatinine over baseline and biopsy-proven rejection [specimens were evaluated by light microscopy and immunofluorescence staining with a marker of classical complement activation (C4d) and classified using the Banff classification as updated in 2017] (Wallin, 2018). As previously documented, a diagnosis of AMR for the kidney needs the presence of DSAs, distinct histological abnormalities, and C4d deposition in peritubular capillaries all at the same time (Cohen et al., 2012; Legaz et al., 2021).
Mild ACR (Banff grade I) was treated with pulse steroids (500 mg methylprednisolone boluses) and increased maintenance immunosuppression. All other ACR were treated with anti-thymocyte globulin. AMR was also treated with pulse steroids and intravenous immunoglobulin (0.25 g/kg and the last session 1 g/kg, maximum 140 g), divided into two doses associated with plasmapheresis (three sessions a day, every 5 days). Later, 500 mg anti-CD20 (Rituximab; Roche Pharmaceuticals) was administered intravenously. Anti-AMR treatment was also administered in two patients receiving anti-proteasome inhibitor Bortezomic (Velcade®, formerly PS-341).
DNA extraction
Genomic DNA from peripheral blood leukocytes was extracted using the Maxwell16 Blood DNA Purification Kit (Promega, MA, USA), according to the manufacturer's instructions on a Maxwell16 kit (Promega). The concentration of the DNA and its purity according to the ratio 260:280 nm was measured with the NanoDrop 2000 (ThermoScientific, MA, USA). About 260:280 nm ratios between 1.8 and 2.2 were considered samples with the appropriate purity, as previously published (Boix et al., 2021).
Genotyping of BAFF system polymorphisms
The polymorphisms' genotyping was performed using TaqMan SNP Genotyping type hydrolysis probes (ThermoFisher) using an ABI-7500 Fast Real-Time PCR System (Applied Biosystems, Singapore). The characteristics of the SNPs studied and of the TaqMan probes used are summarized in Table 2.
Characteristics of the Studied Polymorphisms of the B Cell Activating Factor System
According to the reference genome.
Codigo Probes Taqman (ThermoFisher) with FAM and VIC as reporter fluorophores. All wild-type alleles are tagged with fluorochrome VIC and mutants with FAM.
MAF, minor allele frequency by ALFA Allele Frequency (https://www.ncbi.nlm.nih.gov/snp/rs9514828#frequency_tab).
BAFF, B cell activating factor; BCMA, B cell maturation antigen; L, Leu; N, Asn; N/A, not applicable; P, Pro; R, Arg; R-BAFF, BAFF receptor; rs, reference SNP; S, Ser; SNP, single nucleotide polymorphism; TACI, Transmembrane Activator and Calcium Modulator and Cyclophilin Ligand Interactor.
Approximately 10 ng of DNA were mixed with 0.25 μL of a mix of specific primers for the gene of interest and specific polymorphism probes together with 2.5 μL of polymerase chain reaction (PCR) master mix for a final volume of 5 μL. The thermal cycler parameters were a step of initial activation of Taq polymerase 10 min at 95°C and 40 cycles of denaturalization during 15 sec at 95°C and primers binding and extension, for 60 sec at 60°C. As a negative control, pure RNase-free water was used. For the genotypes assignment, the ABI-7500 Fast Real-Time PCR System software calculates the Rn values of the VIC and FAM reporter fluorophores for each sample. The Rn value is the fluorescence signal normalized by the fluorescence of Rox dye (ThermoFisher Scientific), a passive reference fluorophore.
Five genes of the BAFF system were analyzed in this study: TNFSF13B, TNFSF13, TNFRSF13C, TNFRSF13B, and TNFRSF17 (Table 2). A polymorphism was analyzed for each of the genes. For the study of polymorphisms, three types of possible genotypes were established; Homozygous Allele 1 [Rn(Allele1)/Rn (Allele 2)>1.5], Homozygous Allele 2 [Rn(Allele 1)/Rn (Allele 2)<0.66] and heterozygous [Rn(Allele 1)/Rn (Allele 2)<0.66].
Analysis of the sequences of the TNFSF13 gene and the APRIL protein
The transcript of the TNFSF13 gene was obtained with Ensembl ID ENST00000396545 and was used for bioinformatic analyses at the genomic level.
For protein level analyses, the APRIL sequence was obtained from the UniProt with code O75888 (“UniProt,” 2021). The amino acid sequence of APRIL and the exon where the rs3803800 polymorphism is located are given in Supplementary Table S1.
Prediction of protein alterations by substitution in the DNA sequence of TNFSF13
A total of two different algorithms were used to estimate the prediction of protein alteration owing to polymorphism analyzed in this study (rs3803800).
PolyPhen two predicts the effect of substituting one amino acid for another on the structure and function of human proteins (Adzhubei et al., 2013; “PolyPhen-2: prediction of functional effects of human nsSNPs,” 2021). PolyPhen 2 returns a value between 0 and 1, with which it classifies the substitution as benign, possibly harmful, and probably harmful.
MutPredes predicts the probability that an amino acid change will be harmful and cause disease. The result of the tool indicates the probability that the amino acid substitution is detrimental. These obtained probabilities greater than 0.75 are potentially considered deleterious substitutions (Li et al., 2009).
Prediction of protein stability
I-Mutant suite allows predicting changes in the stability of a protein in mutations by substituting a single amino acid from the protein sequence (“Prediction of Protein Stability Changes upon Mutations,” 2021). MUpro is a predictor of Protein Stability Changes for Single-Site Mutations from Sequences (“Prediction of Protein Stability Changes upon Mutations,” 2021). Both vector machine-based tools for predicting protein stability changes by nonsynonymous substitutions.
The value of the free energy change is predicted, and a score between −1 and 1 is calculated to measure the confidence of the prediction. A score less than zero means that the variant decreases the stability of the protein; conversely, a score greater than zero means that the variant increases the stability of the protein (Cheng et al., 2006).
The energy difference between the native and the mutant protein was calculated based on the Gibbs free energy value (G) for both algorithms. The Gibbs free energy change predicted by the substitution (ΔΔG), expressed in Kcal/mol, was calculated using the following formula: ΔΔG = ΔG (New Protein) − ΔG (Native Protein).
Prediction of post-transductional modifications
The ModPred webserver (“ModPred—Predictor of Post-translational Modification Sites in Proteins—My Biosoftware—Bioinformatics Softwares Blog,” 2021) was used to predict post-translational modification (PTM) sites of APRIL. ModPred is based on a set of logistic regression models for each type of PTM (Pejaver et al., 2014). The prediction is considered to have acceptable confidence when the results are >0.5.
NetOGlyc 4.0 and NetNGlyc 1.0 are neural network-based web tools that predict O- and N-glycosylation sites in mammalian proteins (Steentoft et al., 2013).
Prediction of transcription factor splice and binding sites
SNP2TFBS is a web program (“Prediction of Protein Stability Changes upon Mutations,” 2021) that examines nucleotide variations in transcription factor binding sites in the human genome (Kumar et al., 2017). The algorithm is based on the PWM (position weight matrix) mathematical model, the most widely used to describe a transcription factor's DNA-binding specificity. The difference between the PWM score for the reference and alternative nucleotide is calculated.
The HSF (Human Splicing Finder) tool (“Prediction of Protein Stability Changes upon Mutations,” 2021), predicting splicing donor and acceptor sites, points of branching and splicing regulatory elements such as ISE (Intronic Splicing Enhancers) and ISS (Intronic Splicing Silencers), and whether some of these elements is altered by the variant form of the polymorphism (Desmet et al., 2009).
Statistical analysis
Previously, each SNP was evaluated for the Hardy–Weinberg equilibrium in a control cohort, comparing the frequencies of the observed genotypes and the expected genotypes using a p-value >0.05 were indicative of Hardy–Weinberg equilibrium (“MAF | 1000 Genomes,” 2021). Using univariate logistic regression, the correlation between SNPs and the chance of developing AR or producing DSAs was assessed using odds ratios (ORs).
Four inheritance models were used to assess the association between polymorphisms and AR: codominant, dominant, recessive, and additive. The models' definition assumes one SNP with two alleles (A and a), where A modifies the risk of rejection or producing DSAs. The dominant model assumes that a single copy of A is sufficient to modify the risk and that two copies modify it by the same magnitude. Recessive model, two copies of A are necessary to modify the risk. Therefore, heterozygous Aa and homozygous aa have the same risk. The additive model assumes that each copy of A modifies the risk by an additive amount. Therefore, AA homozygotes have a higher risk than Aa heterozygotes. Codominant model, each genotype provides a different and nonadditive risk. Values of p < 0.05 was considered significant for all statistical tests.
The statistical analyses and graphs were created using the Statistical Package for the Social Sciences (version 22; SPSS, Chicago, IL, USA) software and GraphPad Prism (version 6; San Diego, CA, USA) software, as well as the R programming language, using the Integrated Development R Studio version 3.4 environment for the latter.
Results
Genotypic frequencies of BAFF system polymorphisms and Hardy–Weinberg equilibrium
The observed frequency in KT patients was compared with the expected frequency obtained in a healthy population (“MAF | 1000 Genomes,” 2021). The Hardy–Weinberg distribution for each polymorphism was also analyzed in Table 3.
Analysis of the Frequencies Observed for Each of the Selected Polymorphisms in the B Cell Activating Factor System Genes in Patients with Kidney Transplantation
Hardy–Weinberg distributions according to: https://www.internationalgenome.org/category/maf/ Values of p < 0.05 are considered significant.
KT, kidney transplantation; n, number of individuals; w/m, wild/mutant.
The observed frequencies of the alleles were consistent with the expected frequencies according to the Hardy–Weinberg equilibrium (p > 0.05) for the polymorphisms of the gene TNFSF13B (BAFF, rs9514828), TNFSF13 (APRIL, rs3803800), TNFRSF13C (R-BAFF, rs77874543), and TNFRSF17 (BCMA, rs11570136). Conversely, for the studied polymorphism of the TNFRSF13B gene (TACI, rs34562254), an imbalance concerning the expected frequencies (p < 0.001) is observed, with a frequency for homozygous wild genotypes AA and GG higher than expected (7.2% vs. 1.7% and 81.3% vs. 75.8%, respectively), and frequencies lower than expected of the heterozygous AG genotype (11.5% vs. 22.5%).
Association of polymorphisms with AR and the generation of DSAs
The association of genetic polymorphisms with the incidence of AR was analyzed using different genetic models. The results, summarized in Table 4, show that using the recessive model, the AA genotype of the polymorphism rs3803800 of the TNFSF13 gene is associated with AR apparition (OR = 3.38, 95% CI: 1.02–11.2; p = 0.046) compared with NAR group. In patients with AR, the presence of the AA genotype (the minor allele) amounts to 17.2% (5/29), whereas in patients in the NAR group, it is only observed in 5.8% (8/139).
Analysis of the Risk for Suffering Acute Rejection in Relation to Genetic Polymorphisms and Based on the Inheritance Model
Assuming that M represents the most frequent allele (wild allele) and m the minor frequent, the inheritance models are described below: Codominant: M/M versus M/m versus m/m; dominant: M/M versus (M/m+m/m); recessive: (M/M+M/m) versus m/m; additive: M/m and m/m were weighted respectively as 1 and 2 with respect to the M/M genotype. Underlined the genotypes used as reference in each inheritance model for the calculation of the odd ratios.
OR and p-value were obtained by logistic regression. Values of p < 0.05 were considered significant. Significant results are highlighted in bold.
OR of genotype M/m compared with M/M.
OR of genotype m/m compared with M/M.
OR could not be calculated as there were no m/m genotypes in the acute rejection group.
CI, confidence interval; NAR, without acute rejection; OR, odds ratios.
No significant differences were observed in the rest of the polymorphisms studied in any of the inheritance models used for AR and NAR groups.
On the contrary, the results of the association of the polymorphisms with the generation of DSA during the follow-up period post-transplant were analyzed for the different inheritance models (Table 5). No significant differences were observed for any of the polymorphisms studied, only to highlight a slight increase in the frequency of the AA genotype in the polymorphism rs3803800 of TNFSF13 in transplants with the presence of DSAs (20%) compared with transplants without DSAs (6.8%). However, the differences were borderline near to signification (OR = 3.42, 95% CI: 0.96–12.2; p = 0.057).
Analysis of the Risk for Producing Donor-Specific Antibodies Based on Genetic Polymorphisms and the Inheritance Pattern
Assuming that M represents the most frequent allele and m the least frequent. the inheritance models are described hereunder: Codominant: M/M versus M/m versus m/m; dominant: M/M versus (M/m+m/m); recessive: (M/M+M/m) versus m/m; additive: M/m and m/m were weighted, respectively, as 1 and 2 with respect to the M/M genotype. Underlined the genotypes used as reference in each inheritance model for the calculation of the odd ratios.
OR and p-value were obtained by logistic regression. Values of p < 0.05 were considered significant.
OR of genotype M/m compared with M/M.
OR of genotype m/m compared with M/M.
OR could not be calculated as there were no m/m genotypes in the acute rejection group.
DSA, donor-specific antibodies.
Survival analysis of gene polymorphisms of the BAFF system
The results of the graft survival and rejection-free time analyses are given in Table 6. The follow-up time for graft survival was 51 ± 20.9 months, and during this period, 16 transplant recipients (9.4%) lost the kidney graft and returned to dialysis.
Kaplan–Meier Survival Analysis of B Cell Activating Factor Polymorphisms in Relation to the Time Free from Acute Graft Rejection and Graft Survival at 5 Years
The comparison between the groups was carried out using the log-rank test. Values of p < 0.05 were considered significant. Significant results are highlighted in bold. Assuming M the most frequent allele and m the least frequent allele. The inheritance models that were used for the survival analyses using Kaplan–Meier curves are described hereunder: Codominant: MM versus Mm versus mm. Dominant: MM versus (Mm+mm). Recessive: (MM+Mm) versus mm.
Codom, codominante; Dom, dominant; Rec, recessive.
Regarding the AR-free time, the AA genotype of rs3803800 of the TNFSF13 gene had a lower rejection-free time than the GG and GA genotypes (Table 6 and Fig. 1). At 24 months post-transplantation, 69.2% of those patients transplanted with the AA genotype were AR free compared with 85.7% of those transplanted and who had the GG or GA genotype.

Survival curve for acute rejection-free time according to the rs3803800 genotype of the TNFSF13 gene (APRIL). Comparisons were made using the log-rank test. Values of p < 0.05 were significant. AR, acute rejection.
The analysis of graft survival did not show a significant association with the polymorphisms analyzed in this study.
In silico analysis of the rs3803800 polymorphism of the TNFSF13 gene
Owing to the association of rs3803800 of the TNFSF13 gene with the AR in KT patients, different predictive analyses were carried out to know the effect of this polymorphism at the genomic and protein levels (Table 7).
Analysis of the Impact of the APRIL Protein Polymorphism According to Different Algorithms
GPI, glycosylphosphatidylinositol; PTM, post-translational modification.
The rs3803800 polymorphism generates a nonsynonymous change in position 96 of the APRIL protein, resulting in a substitution of asparagine for a serine (N96S) (Table 2). Analysis of the impact of the APRIL protein polymorphism according to different algorithms was performed (Table 7).
Two predictions of protein alterations algorithms were used, PolyPhen-2 and MutPred, to predict whether the N96S change produced an effect on the coding protein. Both the PolyPhen-2 and MutPred algorithms classify the substitution as benign (p < 0.01 and 0.007, respectively).
To study the changes in the stability of the protein, I-Mutant 2.0 and MUpro algorithms were used. These algorithms predict protein stability by examining the Gibbs free energy difference between the native protein and the protein with the amino acid substitution (ΔΔG), considering relevant stability alterations from | ΔΔG | > 0.5.
The ΔΔG value obtained by I-Mutant 2.0 shows that the N96S substitution caused by the polymorphism produces a decrease in protein stability of −0.30 kcal/mol, whereas MUpro gives a value of ΔΔG = −0.54 kcal/mol, confirming the decrease in protein stability when the change is produced.
Next, prediction of post-transductional modifications to identify possible changes to APRIL's PTM was performed.
For the identification of O- and N-glycosylation sites, we used NetOGlyc 3.1 and NetNGlyc, respectively. Both algorithms indicate that the N96S change does not generate potential changes in the protein's glycosylation pattern.
For the analysis of the prediction of transcription factor splice and binding sites, the ModPred web tool was used. ModPred showed that the N96S substitution removes a possible glycosylphosphatidylinositol (GPI) amidation modification site (score = 0.53). The algorithm also identifies a potential proteolytic cleavage site (score = 0.73) not altered by the N96S change.
At the genomic level, the rs3803800 polymorphism generates a change from an A to a G at position 287 of the coding sequence of the TNFSF13 gene (c.287A>G) that codes for APRIL. To evaluate if this nucleotide change could affect any transcription factor's binding site, we analyzed it using the SNP2TFBS bioinformatics tools. The results show that the transcription factor MZF1 1–4 binding site is found in the polymorphism area, but it is not affected by the nucleotide change.
On the contrary, we use the HSF web tool, which consists of a set of algorithms capable of predicting splicing sites and splicing motifs [Exonic Splicing Enhancer (ESE)] or silencers [Exonic Splicing Silencer (ESS)] of splicing, and if a nucleotide change alters them. The results of the HSF web tool for rs3803800 polymorphism are summarized in Table 8.
Results with the Human Splicing Finder Application
Algorithms included in the HSF application version 3.1.
The site of the rs3803800 polymorphism is highlighted in bold. The analyses were performed with exon 2 of the TNFSF13 transcript with Ensembl ID: ENST00000396545.
Consensus scores vary from 0 to 100 and the threshold is defined by each algorithm. Each score above the threshold is considered a potential ESS/ESE. For enhancers, if the SNP causes the score to fall below the threshold, the SNP is considered to break the ESE. For silencers, if the SNP causes it to exceed the threshold, it is considered that a new ESS is created.
If the SNP exceeds the threshold and produces an alteration of ±10%, it is considered that the variant alters the splicing site.
For the RESCUE EXE algorithm, if the reason is found in the database, it is considered a potential ESE. If the SNP causes formation of an absent reason in the database, it is considered to break the ESE.
ESE, Exonic Splicing Enhancer; ESS, Exonic Splicing Silencer; HSF, Human Splicing Finder.
The HSF matrix algorithm predicts that the G allele can generate new splicing donor and acceptor sites lost in individuals with the A allele. On the contrary, the ESE and ESS motif prediction algorithms infer different effects between the A alleles and G. Although the prediction algorithms infer that the G allele would be part of ESS motifs, which would contribute to the silencing of the expression of exon 2 of the TNFSF13 gene, the A allele, on the contrary, would be part of ESE motifs, promoting its expression. These data indicate that the polymorphism rs3803800 is part of an essential regulatory zone of splicing. Therefore, these results suggest that this polymorphism could alter the splicing efficiency, which has biological consequences.
Discussion
During the last decades, research efforts have focused on the search for new biomarkers that allow in KT to predict and improve the monitoring of renal allograft to achieve the goal of personalized/precision medicine. Thus, our objective was to study the gene polymorphisms of the BAFF system in patients with KT to find predictive biomarkers of AR that can inform future adoption of preventive measures, establish a personalized follow-up of the patient at risk of rejection, and use bioinformatics applications to search crucial functions of these polymorphisms.
The B lymphocyte activating factor (BAFF) plays a crucial role in the survival, maturation, and activation of B lymphocytes (Bossen and Schneider, 2006). Dysregulation of BAFF production has been associated with autoimmune diseases (Davidson, 2010; Wei et al., 2015) and in the field of KT with humoral rejection and the production of DSAs (Min et al., 2016; Zhao et al., 2007). For this reason, the study of genetic polymorphisms that may affect the expression or functionality of the BAFF system molecules may help estimate the risk of the transplanted patient suffering AR or generating DSAs.
Therefore, in this study, we evaluated the association of five polymorphisms of the BAFF system (rs9514828, rs3803800, rs7787543, rs34562254, and rs11570136) with AR development and the production of DSAs in kidney transplant recipients. Our results show that the AA genotype of the polymorphism rs3803800 of the TNFSF13 gene (APRIL) is associated with AR. The in silico analyses suggest that there may be a change in the regulation of splicing of the TNFSF13 gene that could explain the predisposition to graft rejection of RTRs carrying the AA genotype. Other studies found an association with IgA nephropathy and systemic lupus erythematosus susceptibility (Han et al., 2016; Lee et al., 2007; Zhong et al., 2017).
Furuya et al. (2014) hypothesized that this SNP could affect splicing efficiency (which agrees with our bioinformatics results) and favor the expression of one of the two isoforms of the APRIL protein [although the article by Furuya et al. (2012) fails to prove this hypothesis]. It would have important consequences at the biological level because the beta isoform (lacking exon 3) makes the protein have 16 less amino acids and is not secreted into the extracellular environment. The causes of why it is not secreted are not yet known, but it is speculated that it may be owing to a defect when forming the homotrimers of APRIL (which is the way it binds to its receptors and therefore with biological activity).
Therefore, this SNP could alter the efficiency of splicing between exon 2 (where the SNP is) and exon 3 and favor the expression of one of the two isoforms, but this has not been proven yet and remains as hypotheses. Furthermore, the study by Furuya et al. (2021) suggested that what really affects the expression of one isoform or another is not this SNP individually but the haplotype formed by this SNP and others close to it.
R-BAFF is weakly expressed on B lymphocytes, beginning to be expressed in the late transition stage, and its expression progressively increasing as maturation progresses. Its expression decreases when B lymphocytes enter the germinal center and are re-expressed in memory B lymphocytes, absent in plasmablasts and plasma cells.
The polymorphism rs77875443 of the TNFRSF13C gene, which codes for R-BAFF, produces an alteration in the multimerization of R-BAFF on the surface of BL, increasing the risk of suffering from diseases such as common variable immunodeficiency (Pieper et al., 2014). Besides, the rs77874543 polymorphism has been associated with a lower risk of sepsis and multiple sclerosis (Kompoti et al., 2015; Ntellas et al., 2020) and innate and adaptive immunity as predictors of outcome in critically ill patients.
Therefore, patients who carry the homozygous CC genotype have defects in the development of BL because of the structural alteration of the R-BAFF (Lougaris et al., 2016), which, in the context of a transplant, be a carrier of at least one C allele, could have hypothetical and essential clinical consequences. However, in our research, we did not obtain any association, although it must be taken into account that the low frequencies of the C allele, which represents only 9.6% in frequency and with a single homozygous patient (0.1%) in our cohort, make it difficult to draw consistent conclusions.
On the contrary, TACI is an inducible receptor in humans and is expressed by memory B lymphocytes, plasma cells, and a small population of activated CD27-negative B lymphocytes (Sakai and Akkoyunlu, 2017). The polymorphism rs34562254 of the TNFRSF13B gene generates a variation in a proline's amino acid, leucine. The rs34562254 has been associated with multiple myeloma, coronary artery damage in Kawasaki patients (Rand et al., 2016; Zhu et al., 2017), and primary antibody deficiency syndromes and tonsillar hypertrophy and sarcoidosis. Our study did not find an association, although it must be considered that the distribution of the TNFRSF13B genotypes was not in Hardy-Weinberg equilibrium (HWE) in our control cohort.
The data show a bias toward homozygous genotypes in our cohort, their frequency being higher than expected. The lack of HWE may indicate errors in the genotyping method, an evolutionary advantage of homozygosity, or an unexpected population structure. Therefore, the results warrant future replication and validation.
The polymorphism rs11570136 of the TNFRSF17 gene, which codes for BCMA, is located 2 kb upstream of the first exon. Although this region does not constitute a coding zone, it is located in the vicinity of the promoter zone so that polymorphic changes could affect gene expression levels because of the possible alteration in the anchoring of the polymerases in the promoter zone. This polymorphism was previously evaluated in a cohort of patients with chronic B lymphocytic leukemia associated with diagnosis (Meinl et al., 2018). In our study, we did not find an association with the rejection or production of anti-HLA antibodies. However, given the importance of the BCMA receptor in plasma cells' survival, it suggests that other functional polymorphisms may affect humoral responses in kidney recipients.
On the contrary, the polymorphism rs9514828 of the TNFSF13B gene, which codes for BAFF, is located in the binding site of the transcription factor MZF1 (Jasek et al., 2015). The CC genotype of rs9514828 has been associated with higher expression of BAFF (Novak et al., 2006) and a higher incidence of diseases such as Sjögren's syndrome or T cell lymphomas (De Almeida and Petzl-Erler, 2013; Nezos et al., 2014). Different polymorphisms of the BAFF gene have been associated with graft-versus-host disease (Zhai et al., 2012). This polymorphism has not been studied in KT, but others such as rs12583006 have been associated with an increased risk of producing class II DSAs (Clark et al., 2011).
Our study has obtained an association neither with AR development nor with DSA production. Elevated BAFF levels have been mainly associated with AMR (Chang et al., 2017; Pongpirul et al., 2018). In our cohort, mixing rejections of various types and a limited number in terms of the presence of humoral rejections may have caused us not to detect a possible association. Future studies will be necessary to evaluate an association with a cohort exclusively of patients with AMR in greater populations.
The polymorphism rs3803800 of the TNFSF13 gene, which codes for APRIL, has been associated with IgA nephropathy in several studies, in addition to serum immunoglobulin levels (Sango et al., 2016; Yu et al., 2012; Zhong et al., 2017). Our results show that using the recessive model, the AA genotype is associated with a higher risk of developing AR. There are no previous studies of the rs3803800 polymorphism in KT; so, future studies will be necessary to confirm these results. Jasek et al. (2020), through bioinformatics analysis, suggested that this SNP is part of a regulatory element of gene expression and that the A allele is associated with an increase in APRIL gene expression in most tissues. However, they did not conduct laboratory studies to confirm this. Various studies link the A allele of this SNP to pathologies of immunological origin.
Given the linkage of this polymorphism with serum immunoglobulin levels and the development of B lymphocytes, it is likely that it may be involved in allospecific responses to graft, favoring the activation process, recognition of nonshared HLA antigens between donor and recipient production of DSA alloantibodies, and eventual tissue damage of the kidney graft. However, some APRIL haplotypes, which include this polymorphism, have been shown to have a protective effect against lupus, associated with lower soluble APRIL levels (Osman et al., 2012). However, another study did not find an association of rs3803800, neither with the soluble levels nor with the levels of APRIL transcripts (Kawasaki et al., 2007). That is why we decided to carry out a bioinformatics study to establish the possible biological effects resulting from the genetic variant.
The rs3803800 polymorphism generates a change from asparagine to serine at position 96 of the protein's amino acid sequence (N96S). The algorithms for predicting deleterious changes in the protein do not indicate that the N96S substitution generates significant changes in the structure and/or function of APRIL. On the contrary, the algorithms to see the influence on the protein's stability show that N96S generates a decrease of approximately −40 kcal/mol, insufficient to explain essential alterations in the protein structure. The ModPred tool predicts that the amino acid change can eliminate a potential GPI anchoring site, although structural studies of APRIL have not reported so far in the literature that this post-transductional modification occurs in vivo.
At the genomic level, the rs3803800 polymorphism generates a change from an A to a G at position 287 of the coding sequence of the TNFSF13 gene (c.287A>G). The genomic position of rs3803800 is the binding site of the transcription factor MZF1_1–4, but the SNP2TFBS algorithm does not indicate that its binding changes occur because of the nucleotide change. The HSF web tool showed that the nucleotide change could alter the splicing enhancer or silencing motifs. The data suggest that the G allele is crucial in regulating splicing that could be affected in carriers of allele A.
Therefore, these results suggest that rs3803800 could produce an alteration in the expression of the TNFSF13 gene, although according to the results of Furuya et al. (2012), there were no differences, neither in the levels of transcripts nor in the expression of the different soluble isoforms of APRIL (Osman et al., 2012; Zhong et al., 2017).
It should be taken into account that another study (Furuya et al., 2012) was carried out with a low number of samples, and they used APRIL monomers, when in vivo, this molecule acts in the form of trimers, so future functional studies will be necessary to clarify the role of rs3803800 in the expression of APRIL.
Despite the novel results obtained, some important limitations should be mentioned with regard to this study. In the first place, the cohort of RTRs with AR, even considering a time greater than 2 years, is too small to draw definitive conclusions, and it will be necessary to expand it to confirm the results. On the contrary, it would be interesting to repeat the analysis dividing the cohort into subgroups by type of rejection. The dysregulation of the BAFF system molecules has been associated mostly with rejection mediated by antibodies (Banham et al., 2013; Furuya et al., 2012; Legaz et al., 2020; Schuster et al., 2017), so it is possible that in this study, we are losing some relevant association.
On the contrary, when analyzing the results, different inheritance models were compared without using correction method for multiple comparisons such as the Bonferroni or Benjamini–Hochberg method, so that the results obtained could be the result of chance, the result of multiple comparisons made, and should be taken with caution until later confirmation in future studies.
In addition, the multivariate analysis could not be conducted because the potential covariates from clinical records of the patients such as ischemia times, type of therapy of induction, and immunosuppression dose were not available in each patient.
Conclusions
In summary, the rs3803800 genetic polymorphism from this study of BAFF system genes appears to display importance in AR-free time for KT patients, and thus, warrants further research in independent populations as a putative predictive biomarker of AR. These findings also inform future personalized/precision medicine efforts and functional genomic studies in KT patients. Also, studies on gene expression and proteomics should be carried out to analyze the possible biochemical and clinical implications of the genetic variants in the context of transplant prognosis in the early post-transplant years.
Footnotes
Authors' Contributions
M.M. and I.L.: substantial contributions to the conception and/or the design of the article, R.A., E.K.E.B.J., S.L., V.J.-C., H.M.-B., J.A.G., C.B., M.R.M.-Q., J.d.l.P.-M., A.M., M.M. and I.L. to acquisition, analysis and interpretation of the data. All authors have met the ICMJE criteria for authorship and made a significant intellectual contribution, read and approved the final version of the article.
Acknowledgments
The authors thank all study participants and the reviewers who have dedicated their time to improve our work. The authors also thank Instituto de Salud Carlos III (ISCIII) and the Spanish Ministry of Economic and Competitiveness for making our work possible.
Author Disclosure Statement
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.
Funding Information
Grant number PI15/01370 and P19/01194 and co-funding from the European Union with European Fund of Regional Development (FEDER) with the principle of “A manner to build Europe” are gratefully acknowledged.
Abbreviations Used
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
