Skip to main content

Genome-wide pathway analysis identifies VEGF pathway association with oral ulceration in systemic lupus erythematosus



Systemic lupus erythematosus (SLE) is a genetically complex rheumatic disease characterized by heterogeneous clinical manifestations of unknown etiology. Recent studies have suggested the existence of a genetic basis for SLE heterogeneity. The objective of the present study was to identify new genetic variation associated with the clinically relevant phenotypes in SLE.


A two-stage pathway-based approach was used to identify the genetic variation associated with the main clinical phenotypes in SLE. In the discovery stage, 482 SLE patients were genotyped using Illumina Human Quad610 microarrays. Association between 798 reference genetic pathways from the Molecular Signatures Database and 11 SLE phenotypes was tested using the set-based method implemented in PLINK software. Pathways significantly associated after multiple test correction were subsequently tested for replication in an independent cohort of 425 SLE patients. Using an in silico approach, we analyzed the functional effects of common SLE therapies on the replicated genetic pathways. The association of known SLE risk variants with the development of the clinical phenotypes was also analyzed.


In the discovery stage, we found a significant association between the vascular endothelial growth factor (VEGF) pathway and oral ulceration (P value for false discovery rate (P FDR) < 0.05), and between the negative regulation signaling pathway of retinoic acid inducible gene-I/melanoma differentiation associated gene 5 and the production of antinuclear antibodies (P FDR < 0.05). In the replication stage, we validated the association between the VEGF pathway and oral ulceration. Therapies commonly used to treat mucocutaneous phenotypes in SLE were found to strongly influence VEGF pathway gene expression (P = 4.60e-4 to 5.38e-14). Analysis of known SLE risk loci identified a strong association between PTPN22 and the risk of hematologic disorder and with the development of antinuclear antibodies.


The present study has identified VEGF genetic pathway association with the risk of oral ulceration in SLE. New therapies targeting the VEGF pathway could be more effective in reducing the severity of this phenotype. These findings represent a first step towards the understanding of the genetic basis of phenotype heterogeneity in SLE.


Systemic lupus erythematosus (SLE) is a disabling multisystem rheumatic disease with substantial epidemiological variation [1]. SLE is characterized by the dysregulation of the immune system and high phenotypical diversity [2]. This phenotypical heterogeneity includes a wide range of clinical manifestations that are exemplified by the current use of multiple clinical phenotypes as criteria to diagnose the disease [3]. So far, however, little is known about the causes of this phenotypic variation. Understanding the molecular mechanisms associated with the pathogenesis of SLE phenotypes could therefore be of high relevance to develop more efficient therapeutic approaches and preventive strategies.

SLE is characterized by a strong genetic component, with a sibling recurrence rate (λs) of 8–29 and estimated heritability of approximately 66% [4]. To date, nine genome-wide association studies (GWAS) of SLE risk have been performed in European and Asian populations [5]. Together these studies have led to the identification of >40 loci associated with SLE susceptibility. Despite this extraordinary success, there is still a lack of understanding of the genetic variation that is relevant for the development of specific phenotypes within the disease. There is evidence, however, that the main clinical phenotypes in SLE aggregate in families [6], suggesting a genetic basis underlying disease heterogeneity.

To date, only a few candidate gene studies have been performed to uncover the genetics of clinical heterogeneity in SLE [7]. These studies have identified immunity-related genes associated with clinically relevant SLE phenotypes [8]. From these, the most significant associations have been detected between renal disorder and genetic variation in the ITGAM and STAT4 genes, which have been also associated with discoid rash and oral ulceration, respectively [9, 10]. Other significant findings include the association between renal disorder and TNFSF4, malar rash and FCGR2A and hematological disorder and variation in the IL21 gene. So far, however, the genetic component for most SLE phenotypes has been only partially explained. Therefore, the analysis of genetic variation at a genome-wide scale is needed to identify additional variation associated with SLE clinical heterogeneity.

Complex traits like disease risk or clinical phenotypes have been shown to be caused by multiple genes of small effect size [11]. The identification of these small-effect genes is currently one of the major challenges in the characterization of the genetic background for disease phenotypes [12]. Importantly, single-marker GWAS do not allow the identification of genetic variants with small effect sizes, unless extremely large sample sizes are used [13]. In this common type of GWAS, a large number of markers are tested for association and, consequently, stringent significance thresholds are applied, which makes the identification of small-effect variants very difficult [14]. In addition, single-marker GWAS ignore the joint contribution of multiple genes that act coordinately in the same biological process [15]. The characterization of the genetic basis of many complex traits will therefore require the development of new powerful methods that are able to leverage biological knowledge and efficiently integrate the evidence from multiple loci with moderate to small effect sizes.

Recently, novel statistical methodologies that are able to test genetic risk associations at the pathway level have been developed [16]. Pathway-based approaches test whether sets of functionally related genes are jointly associated with a particular phenotype [17]. This methodology strongly reduces the number of association tests and, therefore, it can substantially increase the power to identify new genetic variation compared to single-marker GWAS [18]. The genome-wide pathway analysis (GWPA) has been recently used to characterize the genetic basis of several complex diseases like cancer [19]. Very recently, using the GWPA approach we have identified new genetic variation associated with psoriasis, an autoimmune disease of the skin [20]. This result confirms the utility of GWPA in the study of the genetics of autoimmune diseases.

To gain a better understanding of the genetic basis underlying phenotype heterogeneity in SLE we have performed, for the first time, a GWAS of clinical phenotypes using the GWPA approach. In this study we have analyzed a discovery cohort of 482 SLE patients of European ancestry to determine the association between 798 reference biological pathways and the main clinical phenotypes of SLE represented that are used as diagnosis criteria. Using an independent cohort of 425 SLE patients from the same ancestry, we have then performed a validation study of the most significant genetic pathways. Based on these results, we have performed an in silico validation analysis to evaluate the functional impact of drugs commonly used to treat the associated phenotype. Our findings provide new insights into the biological mechanisms associated with clinical phenotypes of SLE.


Study population

In the discovery stage, a total of 482 SLE patients were recruited. SLE patients were collected from the outpatient clinics of the rheumatology departments of 15 Spanish University Hospitals belonging to the Immune-Mediated Inflammatory Disease (IMID) Consortium [21]. All patients were diagnosed by a rheumatologist. Only those patients with SLE that fulfilled ≥4 of the 1982 revised American College of Rheumatology (ACR) diagnosis criteria were included in the present study [3]. All patients included in this study were >16 years old at the time of sample collection and had >3 years of evolution from the diagnosis date. SLE patients with psoriasis, inflammatory bowel disease (Crohn’s disease or ulcerative colitis) or other rheumatic diseases like rheumatoid arthritis, or multiple sclerosis were excluded from the study. All SLE patients were Caucasian European with all four grandparents born in Spain.

In the validation stage, an independent cohort of 425 SLE patients was used to replicate the genetic pathways that were significantly associated with the SLE phenotypes in the discovery stage. All patients from the validation cohort fulfilled the ACR diagnostic criteria for SLE and were also collected from the IMID Consortium, following the same inclusion and exclusion criteria as for the discovery cohort. All the procedures were followed in compliance with the principles of the Declaration of Helsinki.

The main epidemiological and clinical variables of the discovery and validation cohorts are summarized in Table 1. The distribution of each variable was compared between the discovery and validation cohorts using Fisher’s exact test or Student’s t test for categorical and quantitative variables, respectively.

Table 1 Main epidemiological and clinical features of the discovery and validation patient cohorts

SLE phenotypes

The diagnosis of SLE is of major importance to guide both the disease classification and the patient therapy [22]. Given the high phenotypic heterogeneity of SLE, in order to analyze the most relevant clinical manifestations for disease diagnosis we defined the SLE phenotypes according to the established ACR diagnostic criteria for SLE [3]. Consequently, the 11 SLE phenotypes represented by the ACR diagnostic criteria were analyzed using the GWPA approach. These criteria include malar rash, discoid rash, photosensitivity, oral ulcers, arthritis, serositis, renal disorder, neurologic disorder, hematologic disorder, immunologic disorder and antinuclear antibodies. The distribution of each clinical phenotype in the discovery and replication cohorts is shown in Table 1.

DNA extraction and genome-wide genotyping in the discovery and validation patients

In the discovery stage, the genome-wide genotyping of the 482 SLE patients was performed using the Illumina Quad610 Beadchips (Illumina, San Diego, CA, USA) at the Centro Nacional de Genotipado (CeGen, Madrid, Spain). The genotyping quality control analysis was performed using PLINK software (Additional file 1: Figure S1) [23]. To evaluate the presence of potential population stratification in the SLE patient cohorts, we used the principal component analysis (PCA) implemented in the EIGENSOFT (v4.2) software [24]. Using the first 10 PCs of variation over 10 iterations we identified 14 samples showing an outlier genetic background and were excluded from downstream analysis (Additional file 1: Figure S1). After the quality control analysis, a final dataset of 507,051 single nucleotide polymorphisms (SNPs) and 395 SLE patients was available for the GWPA.

The validation of the two genetic pathways associated with SLE in the discovery stage required the genotyping and analysis of a total 1347 SNPs. Given the large number of variants to be tested and the utility of genome-wide data for accurate genetic ancestry identification, the 425 SLE patients in the validation cohort were genotyped using the same microarray platform. Genotyping for the validation stage was performed at the HudsonAlpha Institute for Biotechnology (Huntsville, AL, USA). The same quality control analysis as in the discovery stage was performed (Additional file 1). A total of 394 SLE patients and all 1347 SNPs from the two genetic pathways passed the quality control and were available for the pathway-based analysis of the validation stage.

Analysis of association between established SLE risk SNPs and SLE phenotypes

Genetic variants associated with SLE risk

The list of established genetic variants (P < 5e-8) for SLE risk was obtained from a recent GWAS meta-analysis in a case-control cohort of European ancestry [5]. A total of 43 genetic variants associated with SLE risk were identified and selected for the analysis of association with SLE clinical phenotypes (Additional file 1: Table S6).

Imputation of genetic variants associated with SLE risk

From the established autosomal SLE risk SNPs (N = 41 SNPs), the genetic variants that were not directly genotyped by the GWAS Quad610 genotyping array (N = 17 SNPs) were imputed (Additional file 1). Those SNPs that did not pass the stringent imputation quality control filter (N = 1 SNP, information quality metric <0.8) were excluded from the study. Therefore, after excluding two non-autosomal variants and a low-quality imputed SNP, a total of 40 from the initial 43 established SLE risk SNPs were finally available for analysis of association with SLE phenotypes in the discovery cohort. In the validation cohort, the same procedure was followed to obtain the genotypes of the risk variants to be tested for replication.

Statistical association analysis

The statistical association analysis between the allele dosage of the established SLE risk SNPs and the SLE clinical phenotypes was performed using the logistic regression model implemented in the SNPTEST v2 software (Oxford, UK) [25]. In this model, the allele dosage was defined as follows:

$$ \mathrm{Allele}\ {\mathrm{Dosage}}_{\mathrm{i}}={\displaystyle \sum_{\mathrm{g}=0}^2} \Pr \left(\mathrm{G}=\mathrm{g}\right)*\mathrm{g} $$

Where g represents each genotype of a particular genetic variant i and Pr(G = i) is the marginal posterior probability obtained by imputation. The allele dosage takes values between 0 and 2. SLE patients without phenotypical data available for the phenotype analyzed were excluded from the association analysis. Finally, the P values obtained from the discovery and replication stages were combined using the METAL software [26].


GWPA method

The gene-set analysis, also referred to as pathway analysis, is a very powerful methodology to analyze the genetic architecture of complex diseases using GWAS data [27, 28]. An important advantage of this approach is that the hypothesis space is significantly reduced compared to single-marker GWAS. While extremely large number of markers (>500,000 to several millions) are independently tested for association in single-marker GWAS, the number of simultaneous tests is several orders of magnitude lower in pathway-based GWAS (typical range 500–2000 pathways). Consequently, the threshold for significance is much less stringent than the consensus threshold used for single-marker GWAS (P < 5e-8) [18, 27]. In addition, pathway-based studies integrate the effects of multiple genetic risk variants that participate in the same biological processes. For these reasons, pathway-based studies can have high statistical power to discover new susceptibility genetic variants provided that they operate within the analyzed pathways.

In the present study, the GWPA was performed using the set-based test implemented in the PLINK software as described previously [20, 23]. Compared to other methods, this set-based test uses genotype data to estimate pathway association instead of P values for significance. Importantly, this approach accounts for the linkage disequilibrium between SNPs and therefore avoids an increase in false positive results due to genes with multiple, highly correlated markers. For each pathway, independent SNPs are first identified (linkage disequilibrium of r 2 < 0.2 here), and from these an average statistic is calculated. Finally, the statistical significance of the pathway is computed using permutation, thereby efficiently correcting by the number of SNPs within the pathway (Additional file 1). As described for the analysis of association with established SLE risk SNPs, those patients with missing data for the phenotype tested for association were excluded from the analysis. In order to account for multiple testing, the false discovery rate (FDR) method was used. The corrected empirical P values from the discovery and validation stages were combined using Fisher’s method.

Gene set definition

Reference biological pathway annotation databases BioCarta, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome were used for the present study [29]. A total of 217, 186 and 674 curated biological pathways from the Biocarta, KEGG and Reactome databases, respectively, were included, respectively (5th October 2015). Very small uninformative pathways (i.e. <=15 genes) were excluded from the analysis. As described previously, we also excluded large genetic pathways (i.e. >300 genes) [17]. The SNP-gene mapping was performed using the NCBI RefSeq database release 63 (12th October 2015) and an SNP-gene distance window of 20 Kb [30]. The final gene set used for the present GWPA was composed of 211,724 SNPs mapping to 798 different pathways. The list of genetic pathways included in the GWPA is shown in Additional file 1: Table S7.

In silico analysis of VEGF pathway genes after treatment with topical immunotherapies for cutaneous SLE

In the GWPA, we identified significant genetic association between oral ulceration and the VEGF pathway. The VEGF pathway plays a crucial role in angiogenesis, and there is increasing evidence supporting the implication of this biological process in the pathogenesis of SLE cutaneous phenotypes [31]. These disease phenotypes are commonly treated with steroid and non-steroid topical immunotherapies in SLE [32, 33]. Consequently, we hypothesized that the topical immunotherapies prescribed for cutaneous SLE mediate their therapeutic effect in this tissue through the VEGF pathway and, therefore, should induce significant transcriptional changes in the pathway genes. In order to test this hypothesis, we used transcriptional data from microarray experiments in the NCBI Gene Expression Omnibus microarray database (GEO, In this database, we searched for whole genome expression profiling datasets generated from cutaneous/mucocutaneous human samples or cell cultures (5 November 2015). From these, we looked for tissue or cell cultures treated with any of the common steroid and non-steroid topical immunotherapies most widely used in SLE (Additional file 1). We found a total of three datasets analyzing the transcriptional variation after treatment with four common immunotherapies: betamethasone valerate and pimecrolimus (GSE32473), diphencyprone (GSE52360) and imiquimod (GSE68182). The first two transcriptional datasets (i.e. GSE32473 and GSE52360) were obtained from skin biopsies and the latter (GSE68182) from an in vitro study on vaginal mucosal cells (i.e. cell line Vk2/E6E7). For each gene expression dataset, we performed quality control analysis and subsequent normalization on the log2 scale using the quantile normalization method. The analysis of differential expression of the VEGF pathway genes between treated and non-treated samples was performed using Student’s t test. The statistical significance of the global perturbation of the VEGF pathway was assessed using the binomial test. All analyses were performed using the R statistical software [34].


Phenotypic characterization of the studied cohorts

The epidemiological and phenotypical characteristics of the discovery and replication SLE populations are shown in Table 1. There were no significant differences between the discovery and replication cohorts in the distribution of the epidemiological and phenotypical variables (P > 0.05, Table 1).

Identification of SLE risk genetic variants associated with SLE phenotypes

In the discovery stage, we found that 19 out of the 43 SNPs previously associated with SLE risk were also significantly associated with one or more clinical phenotypes (P Discovery < 0.05, Table 2). The association results between each established genetic variant and SLE phenotype are summarized in Additional file 1: Table S1. However, in the independent validation cohort, only the association between PTPN22 and hematologic disorder (P Replication = 0.043 , P Combined = 8.25e-4) and between PTPN22 and the production of antinuclear antibodies (P Replication = 0.028, P Combined = 0.001) were significantly replicated. Combining the statistical evidence from the two cohorts, an additional seven loci were found to be associated with SLE phenotypes at the nominal level (P Combined < 0.05, Table 2).

Table 2 SLE risk SNPs association with clinical phenotypes

Identification of genetic pathways associated with SLE phenotypes

In the GWPA, two genetic pathways were found to be significantly associated with an SLE phenotype after multiple test correction (P value for false discovery rate (P FDR) < 0.05, Table 3). The VEGF pathway was associated with the presence of oral ulcers (P FDR = 0.044) and the RIG-I/MDA5 negative regulation signaling pathway was associated with the production of antinuclear antibodies (P FDR = 0.016). The results for analysis of association between each genetic pathway and SLE phenotype are shown in Additional file 1: Table S2.

Table 3 Genetic pathways associated with the SLE phenotypes in the discovery stage

Using the independent validation cohort, the association between the VEGF genetic pathway and oral ulcers in SLE was significantly replicated (P FDR = 0.026, Table 3). The details of association between the VEGF pathway and oral ulcers are shown in Additional file 1: Table S3.

Perturbation of the VEGF genetic pathway by current therapies for cutaneous SLE

Topical immunotherapies are drugs commonly used to treat cutaneous phenotypes of SLE like oral ulceration. Given the observed genetic association between oral ulcers and VEGF pathway, we performed an in silico analysis to evaluate the effect of four current topical immunotherapies on the expression of its constituent genes. Using whole genome expression datasets from patients and relevant cell types treated with these therapies, we found that three out of the four analyzed drugs significantly perturb the expression of VEGF pathway genes (Fig. 1, Additional file 1: Table S4). A total of 16, 12 and 7 genes out of the 29 genes from the VEGF pathway were significantly differentially expressed after imiquimod (P = 5.38e-14), betamethasone valerate (P = 5.69e-9) and diphencyprone (P = 4.59e-4) treatment, respectively (Additional file 1: Table S5).

Fig. 1
figure 1

Vascular endothelial growth factor (VEGF) pathway perturbation after topical immunotherapy. Network representation of VEGF genes according to the differential gene expression after treatment with common topical immunotherapies: imiquimod (a), betamethasone valerate (b), diphencyprone (c) and pimecrolimus (d). Genes are represented as nodes and are connected by edges according to experimental or computational evidence of interaction between their encoded proteins. The diameter of each node is proportional to the significance of differential expression, with significant genes (P < 0.05) in red and non-significant genes (P ≥ 0.05) in blue


One of the major challenges in the pathogenesis of SLE is to understand the biological mechanisms responsible for its phenotypic heterogeneity. Although single-marker GWAS have successfully identified a large number of genetic variants associated with SLE risk, the genetic basis of SLE phenotypes has so far been analyzed only in a few candidate gene studies. In order to identify new genetic variation, we have here performed the first GWAS on SLE phenotypes using a pathway-based approach. Using a discovery cohort of individuals with European ancestry and a validation cohort with the same ancestry, we have identified and validated the association between VEGF genetic pathway and oral ulcers, a common manifestation of SLE. The results of this study provide new insights into the genetic basis and the biological mechanisms associated with the clinical heterogeneity of SLE.

The VEGF pathway is a network of genes that are involved in the transduction of different intracellular signals and act coordinately to modulate inflammatory and angiogenic processes [35, 36]. The dysregulation of angiogenesis has been described as an important biological mechanism in the pathogenesis of SLE [37]. In particular, there is growing evidence that angiogenesis is also involved in the development of skin manifestations in SLE patients [38, 39]. The serum levels of VEGFA protein itself have been suggested as a useful marker for disease activity monitoring in SLE patients [40]. Importantly, the serum levels of VEGFA have also been found to be significantly elevated in SLE patients with cutaneous manifestations [41]. Despite this clinical evidence, the VEGF pathway had not been previously associated with SLE at the genetic level. Our study provides the first evidence of a genetic association between the VEGF pathway and oral ulceration in SLE.

Oral ulcers are frequently chronic mucocutaneous lesions that affect up to 54% of patients with SLE [42, 43]. This clinical manifestation is characterized by high angiogenic activity and a loss of the epithelial and connective tissue in the oral mucosa [44, 45]. Accordingly, anti-angiogenic therapies like thalidomide have been suggested to promote oral ulcer healing and to control ulcer recurrence [46, 47]. From a clinical perspective, oral ulceration has been associated with an increase in the disease activity and a worse prognosis in SLE [48, 49]. The early detection of oral ulcers is therefore highly valuable as it contributes to an earlier diagnosis of SLE and, consequently, to faster initiation of treatment.

The genetic association identified in this study is consistent with previous evidence from other ulcer-related diseases like Behçet disease (BD), recurrent aphtosus ulceration (RAU) and gastroduodenal ulcers. BD is an inflammatory disorder characterized by an extremely high frequency of oral ulcers (>95%) [50]. Clinical evidence suggests that VEGF cytokine could be directly implicated in the formation of oral ulcers in BD [51]. In RAU, the most common oral mucosal disease, the salivary levels of VEGF have been also associated with oral ulceration [52]. Finally, genetic variation in the VEGF gene has been also associated with the risk of developing gastroduodenal ulcers [53]. Evidence from these studies implicates angiogenesis in the pathophysiological development of oral ulcers, at both the genetic and at the functional level. Accordingly, the genes in the VEGF pathway are strong candidates for susceptibility in diseases with a high prevalence of oral ulcers like BD or RAU. Future studies aimed at testing the association between VEGF pathway genes and these diseases are therefore warranted.

Topical steroid and non-steroid immunotherapies have been successfully used for the treatment of the cutaneous manifestations in SLE. However, topical immunotherapies are not exempt from side effects, including skin atrophy or telangiectasias with corticosteroids and the exacerbation of the inflammatory processes with non-steroid immunotherapies like imiquimod [32]. In this study, we have demonstrated that topical steroid and non-steroid immunomodulators significantly affect the expression of the VEGF pathway genes. The results of this in silico analysis indicate that the VEGF genetic pathway could be a key mediator of the benefits of topical immunotherapy to reduce oral ulceration. Therefore, our findings suggest that the VEGF pathway is a relevant source of new drug targets for oral ulceration that could allow more specific treatment while reducing the undesirable side effects of current therapies. In order to confirm the VEGF pathway as new source for drug discovery in the treatment of SLE oral ulceration, further prospective studies using oral ulcer samples from SLE patients are clearly needed.

In the present study, we have also identified and validated the association between SLE risk locus PTPN22 and the production of antinuclear antibodies and with the presence of a hematologic disorder. It is the first time that this coding SNP (rs2476601) has been associated with the development of hematologic disease in SLE. A previous study reported a non-significant trend for association between PTPN22 and antinuclear antibody positivity [54]. The results from our study provide the evidence to confirm the association between this susceptibility gene and the production of a common autoantibody in SLE. Also, the previously identified association between genetic variation in the TNFSF4 gene and renal disorder was replicated in the discovery cohort [9]. Conversely, the reported association between renal disorder and ITGAM and STAT4 genes was not replicated. The renal disorder phenotype encompasses different clinical manifestations (e.g. persistent proteinuria or cellular casts) and, therefore, differences in frequencies in any of these sub-phenotypes could have prevented the replication. Additional studies performing specific analysis of association for each renal disease subtype could help to further refine this genetic association. Oral ulceration was also previously found to be associated with variation in STAT4. This association was not replicated in the present single-marker analysis. The lack of replication could be explained by the comparatively smaller sample size of our study cohorts and by the small effect size reported for the association (OR ~ 1.12). For example, in the European cohort, STAT4 association was not statistically significant despite analyzing >4,000 individuals [9]. Finally, after combining both patient cohorts, we also found seven other SLE risk loci - IL10, IKZF2, MHC class III, UHRF1BP1, ETS1-FLI1, SH2B3 and SLC15A4 - to be nominally associated with different phenotypes. These also represent new genetic associations with SLE heterogeneity. Further studies using independent cohorts of phenotypically well-characterized SLE patients like the present one are needed to corroborate these associations.

GWPA represents a new and powerful approach to identify genetic variation associated with complex phenotypes. However, this study has limitations. First, the sample size of the discovery and replication cohorts is moderate compared to recent case-control GWAS and this could have led to missed pathway associations with SLE clinical phenotypes. Second, the statistical power to detect significant associations is lower for those clinical phenotypes with more extreme frequencies (e.g. 7% of SLE patients have a neurologic disorder). Therefore, larger cohorts of well-characterized SLE patients will be needed to identify additional genetic variation associated with clinical phenotypes. Finally, the GWPA methodology also has intrinsic limitations, mainly related to the current knowledge of biological processes and subsequent definition of the pathways [55]. For many human genes, the functional annotation is far from complete, which precludes their mapping to reference biological pathways. Also, genomic variation located far from the transcribed region itself could be relevant for the regulation of the gene expression. Conversely, variants within genes could be influencing the activity of other distant genes or genes in other chromosomes (i.e. trans-regulation). With the increase in knowledge of genomic regulation [56], the mapping of SNPs to their functionally related genes will clearly improve. Consequently, the integration of this knowledge into GWPA will likely increase the power of this approach to identify new pathways associated with human diseases.


In conclusion, we have performed, for the first time, a genome-wide association analysis to identify genetic risk factors for the main phenotypes in SLE. To do this, we have used a pathway-based analysis approach. Using this approach, we have identified and validated the association between the VEGF genetic pathway and the presence of oral ulcers in SLE. These findings show the existence of a genetic basis underlying SLE heterogeneity that is independent from the genetic component associated with disease risk. The results of the present study could contribute to the development of more efficient therapies to treat cutaneous manifestations of SLE in the near future.



American College of Rheumatology


Behçet disease


Caucasian European




Confidence interval

ETS1 :

ETS Proto-Oncogene 1 gene


False discovery rate


False discovery rate discovery cohort


False discovery rate validation cohort

FLI1 :

Fli-1 Proto-Oncogene gene


Gene Expression Omnibus


Genome-wide association study


Genome-wide pathway analysis


Zinc finger DNA binding protein helios, subfamily 1A, 2 gene

IL10 :

Interleukin-10 gene


Immune-mediated inflammatory disease


Kyoto Encyclopedia of Genes and Genomes



MDA5 :

Melanoma differentiation-associated protein 5 gene


Major histocompatibility complex


Sample size

N- :

Sample size of negative individuals for a particular clinical variable

N+ :

Sample size of positive individuals for a particular clinical variable


Odds ratio

P :

P value

P C :

P value combined

P D :

P value discovery cohort


single nucleotide polymorphism base pair in build GRCh37/hg19

PTPN22 :

Protein tyrosine phosphatase, non-receptor type 22 gene

P V :

P value validation cohort


Disease risk allele


Recurrent aphtosus ulceration


Retinoic acid-inducible gene 1 gene


Standard deviation

SH2B3 :

Lymphocyte-specific adapter protein Lnk 3 gene

SLC15A4 :

Solute carrier family 15 member 4 gene


Systemic lupus erythematosus


Single nucleotide polymorphism


Ubiquitin-like containing PHD and RING finger domains 1-binding protein 1 gene


Vascular endothelial growth factor

λs :

Sibling recurrence rate


  1. Danchenko N, Satia JA, Anthony MS. Epidemiology of systemic lupus erythematosus: a comparison of worldwide disease burden. Lupus. 2006;15(5):308–18.

    Article  CAS  PubMed  Google Scholar 

  2. Lisnevskaia L, Murphy G, Isenberg D. Systemic lupus erythematosus. Lancet. 2014;384(9957):1878–88. doi:10.1016/S0140-6736(14)60128-8. Epub 62014 May 60131.

    Article  PubMed  Google Scholar 

  3. Tan EM, Cohen AS, Fries JF, Masi AT, McShane DJ, Rothfield NF, Schaller JG, Talal N, Winchester RJ. The 1982 revised criteria for the classification of systemic lupus erythematosus. Arthritis Rheum. 1982;25(11):1271–7.

    Article  CAS  PubMed  Google Scholar 

  4. Harley JB, Alarcon-Riquelme ME, Criswell LA, Jacob CO, Kimberly RP, Moser KL, Tsao BP, Vyse TJ, Langefeld CD, Nath SK, et al. Genome-wide association scan in women with systemic lupus erythematosus identifies susceptibility variants in ITGAM, PXK, KIAA1542 and other loci. Nat Genet. 2008;40(2):204–10. doi:10.1038/ng.81.Epub2008Jan1020.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Bentham J, Morris DL, Cunninghame Graham DS, Pinder CL, Tombleson P, Behrens TW, Martin J, Fairfax BP, Knight JC, Chen L, et al. Genetic association analyses implicate aberrant regulation of innate and adaptive immunity genes in the pathogenesis of systemic lupus erythematosus. Nat Genet. 2015;47(12):1457–64. doi:10.1038/ng.3434.Epub2015Oct1426.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Olson JM, Song Y, Dudek DM, Moser KL, Kelly JA, Bruner GR, Downing KJ, Berry CK, James JA, Harley JB. A genome screen of systemic lupus erythematosus using affected-relative-pair linkage analysis with covariates demonstrates genetic heterogeneity. Genes Immun. 2002;3 Suppl 1:S5–S12.

    Article  CAS  PubMed  Google Scholar 

  7. Ceccarelli F, Perricone C, Borgiani P, Ciccacci C, Rufini S, Cipriano E, Alessandri C, Spinelli FR, Sili Scavalli A, Novelli G, et al. Genetic factors in systemic lupus erythematosus: contribution to disease phenotype. J Immunol Res. 2015;2015:745647. doi:10.1155/2015/745647.Epub742015Dec745621.

    Article  PubMed  PubMed Central  Google Scholar 

  8. De Azevêdo SJ, Addobbati C, Sandrin-Garcia P, Crovella S. Systemic lupus erythematosus: old and new susceptibility genes versus clinical manifestations. Curr Genomics. 2014;15(1):52–65. doi:10.2174/138920291501140306113715.

    Article  Google Scholar 

  9. Sanchez E, Nadig A, Richardson BC, Freedman BI, Kaufman KM, Kelly JA, Niewold TB, Kamen DL, Gilkeson GS, Ziegler JT, et al. Phenotypic associations of genetic susceptibility loci in systemic lupus erythematosus. Ann Rheum Dis. 2011;70(10):1752–7. doi:10.1136/ard.2011.154104.Epub152011Jun154130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Taylor KE, Remmers EF, Lee AT, Ortmann WA, Plenge RM, Tian C, Chung SA, Nititham J, Hom G, Kao AH, et al. Specificity of the STAT4 genetic association for severe disease manifestations of systemic lupus erythematosus. PLoS Genet. 2008;4(5):e1000084. doi:10.1371/journal.pgen.1000084.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Robinson MR, Wray NR, Visscher PM. Explaining additional genetic variation in complex traits. Trends Genet. 2014;30(4):124–32. doi:10.1016/j.tig.2014.02.003.Epub2014Mar1011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, McCarthy MI, Ramos EM, Cardon LR, Chakravarti A, et al. Finding the missing heritability of complex diseases. Nature. 2009;461(7265):747–53. doi:10.1038/nature08494.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Du Y, Xie J, Chang W, Han Y, Cao G. Genome-wide association studies: inherent limitations and future challenges. Front Med. 2012;6(4):444–50. doi:10.1007/s11684-012-0225-3.Epub12012Nov11683.

    Article  PubMed  Google Scholar 

  14. Johnson RC, Nelson GW, Troyer JL, Lautenberger JA, Kessing BD, Winkler CA, O'Brien SJ. Accounting for multiple comparisons in a genome-wide association study (GWAS). BMC Genomics. 2010;11:724. doi:10.1186/1471-2164-1111-1724.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Zhang K, Cui S, Chang S, Zhang L, Wang J. i-GSEA4GWAS: a web server for identification of pathways/gene sets associated with traits by applying an improved gene set enrichment analysis to genome-wide association study. Nucleic Acids Res. 2010;38(Web Server issue):W90–95. doi:10.1093/nar/gkq1324.Epub2010Apr1030.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Gui H, Li M, Sham PC, Cherny SS. Comparisons of seven algorithms for pathway analysis using the WTCCC Crohn's disease dataset. BMC Res Notes. 2011;4:386. doi:10.1186/1756-0500-1184-1386.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Ramanan VK, Shen L, Moore JH, Saykin AJ. Pathway analysis of genomic data: concepts, methods, and prospects for future development. Trends Genet. 2012;28(7):323–32. doi:10.1016/j.tig.2012.03.004.Epub2012Apr1013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, Lin X. Powerful SNP-set analysis for case-control genome-wide association studies. Am J Hum Genet. 2010;86(6):929–42. doi:10.1016/j.ajhg.2010.05.002.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Chen D, Enroth S, Ivansson E, Gyllensten U. Pathway analysis of cervical cancer genome-wide association study highlights the MHC region and pathways involved in response to infection. Hum Mol Genet. 2014;23(22):6047–60. doi:10.1093/hmg/ddu304.Epub2014Jun6016.

    Article  CAS  PubMed  Google Scholar 

  20. Aterido A, Julia A, Ferrandiz C, Puig L, Fonseca E, Fernandez-Lopez E, Dauden E, Sanchez-Carazo JL, Lopez-Estebaranz JL, Moreno-Ramirez D, et al. Genome-wide pathway analysis identifies genetic pathways associated with psoriasis. J Invest Dermatol. 2016;136(3):593–602. doi:10.1016/j.jid.2015.11.026.Epub2015Dec1029.

    Article  CAS  PubMed  Google Scholar 

  21. Julia A, Tortosa R, Hernanz JM, Canete JD, Fonseca E, Ferrandiz C, Unamuno P, Puig L, Fernandez-Sueiro JL, Sanmarti R, et al. Risk variants for psoriasis vulgaris in a large case-control collection and association with clinical subphenotypes. Hum Mol Genet. 2012;21(20):4549–57. Epub 2012 Jul 4519.

    Article  CAS  PubMed  Google Scholar 

  22. Yu C, Gershwin ME, Chang C. Diagnostic criteria for systemic lupus erythematosus: a critical review. J Autoimmun. 2014;48-49:10–3. doi:10.1016/j.jaut.2014.1001.1004. Epub 2014 Jan 1021.

    Article  CAS  PubMed  Google Scholar 

  23. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75. Epub 2007 Jul 2025.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38(8):904–9. Epub 2006 Jul 2023.

    Article  CAS  PubMed  Google Scholar 

  25. Marchini J, Howie B, Myers S, McVean G, Donnelly P. A new multipoint method for genome-wide association studies by imputation of genotypes. Nat Genet. 2007;39(7):906–13.

    Article  CAS  PubMed  Google Scholar 

  26. Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1. doi:10.1093/bioinformatics/btq340. Epub 2010 Jul 2198.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Wang K, Li M, Hakonarson H. Analysing biological pathways in genome-wide association studies. Nat Rev Genet. 2010;11(12):843–54. doi:10.1038/nrg2884.

    Article  CAS  PubMed  Google Scholar 

  28. de Leeuw CA, Neale BM, Heskes T, Posthuma D. The statistical properties of gene-set analysis. Nat Rev Genet. 2016;17(6):353–64. doi:10.1038/nrg.2016.29.

    Article  PubMed  Google Scholar 

  29. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739–40. doi:10.1093/bioinformatics/btr260. Epub 2011 May 1735.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Edwards YJ, Beecham GW, Scott WK, Khuri S, Bademci G, Tekin D, Martin ER, Jiang Z, Mash DC, ffrench-Mullen J, et al. Identifying consensus disease pathways in Parkinson's disease using an integrative systems biology approach. PLoS One. 2011;6(2):e16917. doi:10.1371/journal.pone.0016917.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Lesiak A, Narbutt J, Kobos J, Kordek R, Sysa-Jedrzejowska A, Norval M, Wozniacka A. Systematic administration of chloroquine in discoid lupus erythematosus reduces skin lesions via inhibition of angiogenesis. Clin Exp Dermatol. 2009;34(5):570–5. doi:10.1111/j.1365-2230.2008.03006.x. Epub 02008 Dec 03010.

    Article  CAS  PubMed  Google Scholar 

  32. Chang AY, Werth VP. Treatment of cutaneous lupus. Curr Rheumatol Rep. 2011;13(4):300–7. doi:10.1007/s11926-011-0180-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Khandpur S, Sharma VK, Sumanth K. Topical immunomodulators in dermatology. J Postgrad Med. 2004;50(2):131–9.

    CAS  PubMed  Google Scholar 

  34. Ihaka R, Gentleman G. R: A language for data analysis and graphics. J Comput Graph Stat. 1996;5(5):299–314.

    Google Scholar 

  35. Shibuya M, Claesson-Welsh L. Signal transduction by VEGF receptors in regulation of angiogenesis and lymphangiogenesis. Exp Cell Res. 2006;312(5):549–60. Epub 2005 Dec 2005.

    Article  CAS  PubMed  Google Scholar 

  36. Ferrara N. Vascular endothelial growth factor: basic science and clinical progress. Endocr Rev. 2004;25(4):581–611.

    Article  CAS  PubMed  Google Scholar 

  37. Sakly N, Mirshahi P, Ducros E, Soria J, Ghedira I, Mirshahi M. Angiogenic activity in sera of patients with systemic lupus erythematosus. Lupus. 2009;18(8):705–12. doi:10.1177/0961203309103087.

    Article  CAS  PubMed  Google Scholar 

  38. Walchner M, Meurer M, Plewig G, Messer G. Clinical and immunologic parameters during thalidomide treatment of lupus erythematosus. Int J Dermatol. 2000;39(5):383–8.

    Article  CAS  PubMed  Google Scholar 

  39. Cuadrado MJ, Karim Y, Sanna G, Smith E, Khamashta MA, Hughes GR. Thalidomide for the treatment of resistant cutaneous lupus: efficacy and safety of different therapeutic regimens. Am J Med. 2005;118(3):246–50.

    Article  CAS  PubMed  Google Scholar 

  40. Kuryliszyn-Moskal A, Klimiuk PA, Sierakowski S, Ciolkiewicz M. Vascular endothelial growth factor in systemic lupus erythematosus: relationship to disease activity, systemic organ manifestation, and nailfold capillaroscopic abnormalities. Arch Immunol Ther Exp (Warsz). 2007;55:179–85. Epub 2007 Jun 2008.

    Article  CAS  Google Scholar 

  41. Moneib HA, Salem SA, Aly DG, Khedr HT, Wafaey HA, Hassan HE. Assessment of serum vascular endothelial growth factor and nail fold capillaroscopy changes in systemic lupus erythematosus with and without cutaneous manifestations. J Dermatol. 2012;39(1):52. doi:10.1111/j.1346-8138.2011.01322.x. Epub 02011 Sep 01328.

    Article  CAS  PubMed  Google Scholar 

  42. Mustafa MB, Porter SR, Smoller BR, Sitaru C. Oral mucosal manifestations of autoimmune skin diseases. Autoimmun Rev. 2015;14(10):930–51. doi:10.1016/j.autrev.2015.06.005.Epub2015Jun1024.

    Article  PubMed  Google Scholar 

  43. Khatibi M, Shakoorpour AH, Jahromi ZM, Ahmadzadeh A. The prevalence of oral mucosal lesions and related factors in 188 patients with systemic lupus erythematosus. Lupus. 2012;21(12):1312–5. Epub 2012 Jul 1325.

    Article  CAS  PubMed  Google Scholar 

  44. Munoz-Corcuera M, Esparza-Gomez G, Gonzalez-Moles MA, Bascones-Martinez A. Oral ulcers: clinical aspects. A tool for dermatologists. Part II. Chronic ulcers. Clin Exp Dermatol. 2009;34:456–61. doi:10.1111/j.1365-2230.2009.03219.x. Epub 02009 Apr 03214.

    Article  CAS  PubMed  Google Scholar 

  45. Arbiser JL, Johnson D, Cohen C, Brown LF. High-level expression of vascular endothelial growth factor and its receptors in an aphthous ulcer. J Cutan Med Surg. 2003;7(3):225–8.

    Article  PubMed  Google Scholar 

  46. Jacobson JM, Greenspan JS, Spritzler J, Ketter N, Fahey JL, Jackson JB, Fox L, Chernoff M, Wu AW, MacPhail LA, et al. Thalidomide for the treatment of oral aphthous ulcers in patients with human immunodeficiency virus infection. National Institute of Allergy and Infectious Diseases AIDS Clinical Trials Group. N Engl J Med. 1997;336(21):1487–93.

    Article  CAS  PubMed  Google Scholar 

  47. Grinspan D. Significant response of oral aphthosis to thalidomide treatment. J Am Acad Dermatol. 1985;12(1 Pt 1):85–90.

    Article  CAS  PubMed  Google Scholar 

  48. Urman JD, Lowenstein MB, Abeles M, Weinstein A. Oral mucosal ulceration in systemic lupus erythematosus. Arthritis Rheum. 1978;21(1):58–61.

    Article  CAS  PubMed  Google Scholar 

  49. Parodi A, Massone C, Cacciapuoti M, Aragone MG, Bondavalli P, Cattarini G, Rebora A. Measuring the activity of the disease in patients with cutaneous lupus erythematosus. Br J Dermatol. 2000;142(3):457–60.

    Article  CAS  PubMed  Google Scholar 

  50. Taylor J, Glenny AM, Walsh T, Brocklehurst P, Riley P, Gorodkin R, Pemberton MN. Interventions for the management of oral ulcers in Behcet's disease. Cochrane Database Syst Rev. 2014;9:CD011018. doi:10.1002/14651858.CD14011018.pub14651852.

    Google Scholar 

  51. Yalcin B, Arda N, Tezel GG, Erman M, Alli N. Expressions of vascular endothelial growth factor and CD34 in oral aphthous lesions of Behcet's disease. Anal Quant Cytol Histol. 2006;28(6):303–6.

    PubMed  Google Scholar 

  52. Brozovic S, Vucicevic-Boras V, Mravak-Stipetic M, Jukic S, Kleinheinz J, Lukac J. Salivary levels of vascular endothelial growth factor (VEGF) in recurrent aphthous ulceration. J Oral Pathol Med. 2002;31(2):106–8.

    Article  CAS  PubMed  Google Scholar 

  53. Kim YS, Park SW, Kim MH, Jang EJ, Park JS, Park SJ, Baik HW, Chung G, Hahm KB. Novel single nucleotide polymorphism of the VEGF gene as a risk predictor for gastroduodenal ulcers. J Gastroenterol Hepatol. 2008;23 Suppl 2:S131–9. doi:10.1111/j.1440-1746.2008.05404.x.

    Article  CAS  PubMed  Google Scholar 

  54. Namjou B, Kim-Howard X, Sun C, Adler A, Chung SA, Kaufman KM, Kelly JA, Glenn SB, Guthridge JM, Scofield RH, et al. PTPN22 association in systemic lupus erythematosus (SLE) with respect to individual ancestry and clinical sub-phenotypes. PLoS One. 2013;8(8):e69404. doi:10.1371/journal.pone.0069404. eCollection 0062013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Khatri P, Sirota M, Butte AJ. Ten years of pathway analysis: current approaches and outstanding challenges. PLoS Comput Biol. 2012;8(2):e1002375. doi:10.1371/journal.pcbi.1002375. Epub 1002012 Feb 1002323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. ENCODE:Project:Consortium. The ENCODE (ENCyclopedia Of DNA Elements) Project. Science. 2004;306(5696):636–40.

    Article  Google Scholar 

Download references


The authors thank the patients and clinical specialists collaborating in the IMID Consortium for participation.


This study was funded by the Spanish Ministry of Economy and Competitiveness (grant numbers: PSE-010000-2006-6 and IPT-010000-2010-36) and by the “Agència de Gestió d’Ajuts Universitaris i de Recerca” (AGAUR, Generalitat de Catalunya, FI-DGR 2016, grant number: 00587). The study sponsor had no role in the collection, analysis or interpretation of the data.

Availability of data and materials

Not applicable.

Authors’ contributions

Study design: AA, AJ, AFN and SM. Data acquisition: PC, RB, FJLL, JJPV, AO, JLA, MAAZ, PV, JMN, JLMF, AZ, JMP, MF, ED, MLL, MLC, JLG, AFN and SM. Statistical analysis and data interpretation: AA, AJ, AFN and SM. Sample processing: NP and RT. Genotyping: RMM and DA. Manuscript drafting: all authors. All authors read and approved the final manuscript.

Authors’ information

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

The study and the consent procedure were reviewed and approved by the Vall d’Hebron Hospital Ethics Committee. Informed consent was obtained from all participants.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Antonio Julià or Antonio Fernández-Nebro.

Additional file

Additional file 1:

Supplementary information is available in the online Supplementary Material file. (PDF 13965 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Aterido, A., Julià, A., Carreira, P. et al. Genome-wide pathway analysis identifies VEGF pathway association with oral ulceration in systemic lupus erythematosus. Arthritis Res Ther 19, 138 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: