We searched for association studies on 13 February 2020, the start date for this meta-analysis. Four strings of search terms were used that included combinations of “vascular endothelial growth factor”, “VEGF”, “polymorphism”, “cytokine”, “renal”, “transplant”, “allograft”, and “kidney transplantation” as medical subject heading and text in MEDLINE using PubMed, Google Scholar, Science Direct and Mednar, unrestricted by language. Details of the search strategies for each of these four databases are shown in Table S1 ( Extended data ²⁵ ).

References cited in the retrieved articles were also hand-screened to identify additional eligible studies. In case of duplicate articles, we selected the one with a later date of publication.

The following PICO elements were applied in the meta-analysis: (i) Population: renal allograft patients; (ii) Intervention: VEGF gene polymorphisms; (iii) Comparators: rejectors (RJ) versus non-rejectors (NRJ); and (iv) Outcome: allograft rejection post-KT.

Inclusion criteria were: (i) case–control design evaluating the association between VEGF SNPs and risk of allograft rejection; (ii) available VEGF genotype frequencies in the presence and absence of allograft rejection and (iii) sufficient genotype frequency data to enable calculation of the odds ratios (ORs) and 95% confidence intervals (CIs). Exclusion criteria were studies that: (i) did not involve renal allografts; (ii) were review articles; (iii) were functional studies; (iv) did not involve VEGF SNPs and with genotype or allele frequencies that were unusable/absent or, when available, combined with SNPs in other genes, preventing proper data extraction.

We examined four SNPs ( Table 1; Extended data: S2 Table ²⁵ ). Observed phenotypic associations have been attributed to the proximity of SNPs in the VEGF gene ^{26– 28} , termed linkage disequilibrium (LD). LD is the correlation between alleles located near each other ²⁹ and is measured in terms of D′ and r ² with a value of 1 indicating complete LD ^{30, 31} . LD values were based on the European (CEU), and the Indian populations (Gujarati: GIH and Telugu: ITU) from LDlink. Complete LD between rs699947 (-2578C/A) and rs144854329 (-2549 insertion/deletion) merited combination, labeled VEGF1. -1154G/A (rs1570360), and 938C/T (rs3025039) were not in complete LD, thus analyzed separately, notated as VEGF2 and VEGF3, respectively ( Table 2).

Two investigators (TE and NP) independently extracted data and arrived at a consensus. Authors of the component articles were contacted is cases of missing data. The following information were obtained from each publication: first author’s name, year of the study, country of origin, ethnicity, age of the subjects, comparators, VEGF SNPs (rs number), including transplant outcome in term of type of allograft rejection and values needed to tally the Clark-Baudouin score ( Table 1). Sample sizes as well as genotype data in RJ and NRJ were also extracted along with calculated outcomes of the minor allele frequency. HWE was assessed using the application in https://ihg.gsf.de/cgi-bin/hw/hwa1.pl, HWE was reported as P-values of the controls from the Pearson's goodness-of-fit χ ²-square test.

Using the G*Power program ³² , we evaluated statistical power. Assuming an OR of 1.5 at a genotypic risk of α = 0.05, power was considered adequate at ≥80%. Methodological quality of the included studies was assessed with the Clark-Baudouin scale ³³ . In this scale, scores of <5, 5–6 and ≥7 represent low, moderate and high quality, respectively.

Given the hypothesis of association between VEGF SNPs and risk of allograft rejection following KT, we estimated the ORs with 95% CIs for each study by comparing RJ with NRJ among transplant recipients. Table 2 shows the frequencies of the variant ( var) and wild-type alleles, as well as wt- var or heterozygous genotype ( het). Non-uniformity of the variant ( var) allele in VEGF1 and VEGF2 warranted the use of the allele-genotype model for VEGF1 and VEGF2. On the other hand, the var alleles in VEGF3 (rs3025039) were uniform (all < 0.50), so the standard genetic models were suitable: (i) homozygous: var– var and wt– wt genotypes compared with wt– wt; (ii) recessive: var– var versus het + wt– wt; (iii) dominant: var– var + het versus wt– wt; and (iv) codominant: var versus wt. Using raw data for frequencies, study specific risks (ORs) of allograft rejection were estimated and pooled ORs were calculated by comparing the effects on the same baseline. Multiple comparisons were corrected with the Bonferroni test. Subgrouping was based on ethnicity (Indians/Caucasians) and type of rejection (AR/CR). High significance (P ^a < 0.0001) indicated strong evidence for association.

Heterogeneity in meta-analysis ³⁴ was addressed with the following: (i) its presence warranted use of the random-effects model ³⁵ , otherwise fixed-effects model ³⁶ was used; (ii) estimated with the χ ²-based Q test ³⁷ ; (iii) quantified with the I ²statistic ³⁸ ; and (iv) sources were outlier treated. Outlier treatment divided the comparisons into pre-outlier and post-outlier.

Sensitivity analysis was used to test for robustness of the summary effects. Publication bias was considered for significant (P ^a < 0.05) comparisons with ≥ 10 studies ⁴⁴ . Significance was set at a two-sided P-value of < 0.05, except for heterogeneity estimation, which was set at P _het < 0.10) ³⁷ . Data for the meta-analysis were analyzed using Review Manager 5.3 (Cochrane Collaboration, Oxford, England), SIGMASTAT 2.03, and SIGMAPLOT 11.0 (Systat Software, San Jose, CA).

Figure 1 outlines the study selection process in a flowchart following guidelines form the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA; Reporting guidelines). Table S1 ( Extended data ²⁵ ) shows the initial search using combinations of four search strings applied to four databases resulted in 1,949 citations, followed by a series of omissions that mostly involved duplications (n = 1,924), The gray literature database (Mednar) yielded no additional papers for inclusion. Thus, the final number of included articles for this meta-analysis was seven ^{23, 24, 39– 43} .

Of the seven articles, five ^{23, 24, 39– 41} examined more than one VEGF SNP ( Table 1). The number of studies VEGF1 (rs699947 and rs144854329), VEGF2 (rs1570360) and VEGF3 (rs3025039) were 10, five and three, respectively ( Table 2). Of the 10 VEGF1 studies, seven and three were in Indian ^{39– 41} and Caucasian ^{23, 24, 43} populations, respectively. Of the five VEGF2 studies, three and two were in Indian ^{39– 41} and Caucasian ^{23, 24} populations, respectively. One Caucasian ⁴² and two Indian ^{40, 41} studies comprised VEGF3. Table 1 shows two publications ^{41, 43} that investigated CR, which translated to three studies for VEGF1 ( Table 2), otherwise, the rest focused on AR ( Table 1 and Table 2).

Table 2 shows an aggregate total sample size (663 RJ/956 NRJ) and a statistical power of 97.7% for VEGF1. In contrast, both VEGF2 (105 RJ/254 NRJ) and VEGF3 (265 RJ/290 NRJ) were underpowered (40.5% and 65.2%). Mean age of the subjects was 39.96±6.6 years (± standard deviation) indicating a near to middle-age demographic profile of the KT subjects. The Clark-Baudouin scores (median 6.0, interquartile range 6.0–6.75) indicated that the methodological quality of the component studies was moderate. Control frequencies deviated from the HWE in three studies (from two articles) for VEGF1 ^{40, 41} , two studies ^{39, 40} for VEGF2, and one study for VEGF3 ⁴⁰ .

VEGF1 associations with KT. Table S2 ( Extended data ²⁵ ) shows 32 comparisons, six of which were significant (P ^a = 0.0009–0.04). Of the six, five were post-outlier derived and four survived the Bonferroni correction ( Table 3). Of the four, three were in wt indicating increased risk (overall: 1.41, 95% CI 1.14-1.75, P ^a = 0.002 [ Figure 2], Indian: OR 1.44, 95% CI 1.13-1.84, P ^a = 0.004, CR: OR 2.10, 95% CI, P ^a = 0.0009) and one in var, indicating reduced risk (Indian: OR 0.61, 95% CI 0.45-0.820, P ^a = 0.001). Only the CR outcome had zero heterogeneity (I ² = 0%).

VEGF2 associations with KT. Table S2 ( Extended data ²⁵ ) shows 18 comparisons, four of which were significant (P ^a = 0.001–0.04), were in the wt model and had moderate heterogeneity (I ² = 30%-51%). Three of the four were products of pre-outlier analysis, where the HWC outcome (OR 1.39, 95% CI 1.01-1.91, P ^a = 0.04) confirmed the overall outcome (OR 1.48, 95% CI 1.01-2.15, P ^a = 0.04). The other overall outcome was post-outlier derived and survived the Bonferroni correction (OR 1.58, 95% CI 1.19-2.09, P ^a = 0.0001). The significant Caucasian outcome (OR 1.55, 95% CI 1.06-2.28, P ^a = 0.02) contrasted with the non-significant Indian outcome (OR 1.36, 95% CI 0.72-2.58, P ^a = 0.34).

VEGF3 associations with KT. Table S3 ( Extended data ²⁵ ) shows eight comparisons, two of which were significant (P ^a = 0.008–0.01) but did not withstand Bonferroni correction. These two homogeneous (I ² = 0%) pooled ORs indicated reduced risk in the dominant and codominant models (ORs 0.66–0.69, 95% CIs 0.47-0.92).

Summary of significant VEGF associations with KT. Table 3 summarizes the information on the 12 significant outcomes, five of which survived the Bonferroni correction, four in VEGF1 and one in VEGF2, all deemed robust. These outcomes identified three VEGF polymorphisms (rs699947, rs144854329 and rs1570360) that were associated with allograft rejection post-KT. VEGF1 subgroup outcomes identified CR associations and Indians to be at risk. Depending on the genetic model, the Indian population were both susceptible ( wt: OR 1.44, 95% CI 1.13-1.84) and protected ( var: OR 0.61, 95% CI 0.45-0.82).

The five Bonferroni-filtered findings ( wt and var alleles) were either products of outlier treatment and/or subgrouping. Subgrouping identified the ethnicity and rejection type that was significant, thus specifying associations of the VEGF polymorphisms with allograft rejection post-KT. Subgrouping provided contrasts regarding significant outcomes: (i) In VEGF1, significant in Indians (P ^a = 0.001–0.004), non-significant in Caucasians (P ^a = 0.78–1.00); (ii) in VEGF2, significant in Caucasians (P ^a = 0.02) and non-significant in Indians (P ^a = 0.34); (iii) in VEFG1, significant in CR (P ^a = 0.0009), non-significant in AR (P ^a = 0.12). Subjecting these P ^a-values to Bonferroni correction and sensitivity treatment raised the level of evidence that facilitated interpretation with greater confidence. We have shown that meta-analytical tools such as subgrouping, outlier and sensitivity treatments are instrumental in generating evidence for association. By design, such features are not present in the component single-study outcomes. This underpins the value of meta-analysis in systematically synthesizing primary study results and providing insight into associations of VEGF SNPs with allograft rejection post-KT. Conflicting outcomes between primary studies may be due to small sample sizes, hence, lack of power. Underpowered outcomes appear to be common in candidate gene studies ⁴⁵ and are prone to the risk of Type 1 error. In spite of the evidence for associations, the complexity of allograft rejection involves interactions between genetic and non-genetic factors allowing for the likelihood of environmental involvement. Gene-gene and gene-environment interactions have been reported to have roles in associations of other SNPs with post-KT allograft rejection. Two articles ^{39, 43} examined polymorphisms in other genes that included interleukin 18 ( IL18), transforming growth factor beta 1 ( TGFB1) and angiotensin II receptor type 1 ( AGTR1) . None of the seven articles acknowledged gene-environment interaction. Four ^{23, 39– 41} of the included articles mentioned haplotype analysis with three presenting haplotype data ^{23, 39, 40} . Additional well-designed studies exploring other parameters would confirm or modify our results in this study and add to the extant knowledge about the association of the VEGF SNPs and renal allograft outcome.

VEGF plays a crucial role in kidney physiology with its involvement in maintaining the integrity and permeability of the glomerular capillary basement membrane ¹⁷ . Adaptive response of VEGF toward renal allograft tissue may be related to its angiogenic property on endothelial cells since VEGF contributes to tissue repair response of damaged capillaries ²³ . After KT, the recipient’s neutrophils and macrophages infiltrate the allograft after reperfusion of the transplanted tissue leading to the production of VEGF ²⁴ . Shahbazi et al. showed that genetically directed variations in VEGF production with increased frequency of VEGF producing alleles seemed to influence susceptibility to acute allograft rejection ²⁴ . However, Lemos et al. also suggested that renal allograft recipients with genetic potential for high VEGF production had significantly better graft survival compared to recipients with low VEGF production ²³ . Our results along the timeline of post-KT outcomes indicated increased risks, both for AR and CR in the wt allele, which agreed with Shahbazi et al. ²⁴ but contrasted with Lemos et al. ²³ . However, the significance of our increased risk CR finding may require caution in its interpretation given the low number of studies (n = 2) and low statistical power (64.4%). More studies may be needed to clarify our CR outcome. In terms of ethnicity, Indians carriers of the wt CC genotype in rs699947 (-2578C/A), were afforded better graft survival than the CA and AA genotypes ⁴¹ . In contrast, Shahbazi et al. found that the -2578 C allele (rs699947) and the -1154 G allele (rs1570360) were associated with increased risk of acute renal allograft rejection in Caucasians conferring greater risk among wt homozygotes (-2578C/C and -1154 G/G) compared to -2578C/A and -1154G/A heterozygous genotypes ²⁴ . These inconsistent associations among previous studies may be due to the variations in genetic background influenced by differential ethnicities of the patients.

Interpreting our findings should consider its limitations and strengths. Strengths include: (i) VEGF1 combined sample sizes translated to high aggregate statistical power (97.7%); (ii) significant HWC outcomes validated the overall pooled effects in wt. These validations served to reduce the risk of genotyping errors and minimize methodological weaknesses in our study; (iii) subgroup outcomes in CR and Indians point to potential clinical utility in the genetics of renal transplantation; (iv) efficiency of outlier treatment was the key to generating associative significance and eliminating or reducing heterogeneity and (v) stability of the core overall outcomes are underpinned by surviving the Bonferroni correction (minimizing Type 1 error risk) and robustness (determined with sensitivity treatment). On the other hand, limitations include: (i) all the component studies were underpowered; (ii) most of the moderately significant outcomes (67%) were non-robust.

All data underlying the results are available as part of the article and no additional source data are required.

Dryad: Influence of polymorphisms in the vascular endothelial growth factor gene on allograft rejection after kidney transplantation: a meta-analysis, https://doi.org/10.5061/dryad.gqnk98skz ²⁵ .

This project contains the following extended data:

Dryad: PRISMA checklist for ‘Influence of polymorphisms in the vascular endothelial growth factor gene on allograft rejection after kidney transplantation: a meta-analysis’, https://doi.org/10.5061/dryad.gqnk98skz ⁴⁶ .

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).