INTRODUCTION

Osteoarthritis (OA) is the most common joint disorder worldwide. The knee is the most commonly involved site, with an estimated prevalence between 10% and 13% in elderly people¹. OA is classified as primary when there is no apparent underlying cause, and secondary when it is due to a known cause. OA is a whole joint disease primarily characterized by articular cartilage degeneration; however, synovial inflammation, subchondral bone sclerosis, and osteophyte formation are also observed². Inflammation plays a role in the pathogenesis of OA. Inflammation stimulates chondrocytes and synovium cells to produce cytokines such as interleukin-1 beta (IL-1β) and tumor necrosis factor-alpha (TNF-α), as well as collagenases and aggrecanases such as matrix metalloproteinases (MMPs) and a disintegrin and metalloproteinase with thrombospondin type I motifs (ADAMTS), favoring the catabolism of articular cartilage²,³. Angiogenesis is also involved in OA development. During the osteoarthritic process, blood vessels from the subchondral bone invade the articular cartilage, which is normally avascular, facilitating the progression of OA and the formation of osteophytes. Likewise, angiogenesis occurs in the synovium, where it is strongly linked to macrophage infiltration and inflammation⁴.

Leptin (LEP) belongs to the type 1 cytokine superfamily. LEP is mainly produced by adipocytes and is primarily known for its role as a hypothalamic modulator of food intake; however, it also has a role as a cytokine upregulating the secretion of inflammatory cytokines, including IL-6, IL-12, and TNF-α⁵,⁶. LEP plays a key role in the OA process. LEP levels are increased in synovial fluid and serum of OA patients; furthermore, there is a correlation between the LEP levels of synovial fluid and OA severity⁷,⁸. LEP is also produced by chondrocytes and promote cartilage damage by inducing the production of IL-1β, IL-6, IL-8, MMP-1, MMP-3, MMP-13, ADAMTS-4, ADAMTS-5, and nitric oxide (NO)⁹,¹⁰. Vascular endothelial growth factors A are a family of platelet-derived growth factors being VEGF-A the most prominent¹¹. VEGF-A is a potent angiogenic factor and also plays an active role in OA. In OA patients, VEGF-A levels are increased in the synovial fluid and serum, and there is also a correlation between synovial fluid VEGF-A levels and the grade of OA severity¹²,¹³. Hypoxic, mechanical, and inflammatory stimuli increase VEGF-A production in chondrocytes¹⁴. In chondrocytes and synovial fibroblasts, VEGF-A induces cartilage catabolic mediators such as NO, MMP-1, MMP-3, and MMP-13, as well as IL-1β, IL-6, and TNF-α¹⁵,¹⁶.

Despite the prominent role of LEP and VEGF-A in OA pathogenesis, to date, there are few studies analyzing the association of LEP and VEGFA gene variants with knee OA¹⁷-¹⁹. Genetic association studies in primary OA have typically focused on independent effects of candidate genes. Epistasis is the interaction between two or more genes. This interaction can modify the genetic associations independently of the particular effects of the genes and can be a factor in the susceptibility to developing a disease²⁰. It is, therefore, necessary to explore epistasis in OA, as well as single-gene effects, which prompted us to analyze the association of LEP and VEGFA gene variants with primary knee OA as well as for an interaction between them in the Mexican Mestizo population.

METHODS

RESULTS

A total of 282 participants were included: 103 comprised cases with primary knee OA and 179 were controls. Females were more frequent in cases and controls (71.8% and 67.6%, respectively; p = 0.4). Mean age was 65.5 ± 14.6 years in cases and 60.8 ± 15.6 years in controls (p = 0.01). BMI showed a mean of 25.1 ± 3.1 and 25.0 ± 3.8 in cases and controls, respectively (p = 0.8). Arterial hypertension and Type II diabetes mellitus were the most frequent comorbidities (p = 0.03 and 0.05, respectively) (Table 1).

Both gene variants, rs2167270 of LEP and rs2010963 of VEGFA, were in agreement with HWE (p = 0.3 and 0.4, respectively). The allelic distribution of both SNPs did not show differences among the groups (p = 0.02 and 0.4; respectively) (Table 2). In the genotypic association analysis, a tendency toward an increased risk under a codominant model and a recessive model was observed for AA homozygous of rs2167270. In the case of rs2010963, a tendency toward a lower risk for CC homozygous was observed under a codominant and a recessive model. However, these tendencies were not statistically significant (Table 2).

Through MDR analysis, an interaction was observed between the LEP and VEGFA. This model had a testing accuracy of 0.5199 and a maximum CVC of 10/10. The model was significant at the level of 0.02, indicating that all testing accuracies were greater than 0.5 during cross validation (Table 3). In addition, this model was associated with an increased risk to knee OA [OR (95% CI) = 1.8 (1.1-2.9)]. Figure 1A summarizes the distribution of the two-locus genotype combination of LEP and VEGFA in cases and controls and for the new MDR-constructed variable. As can be observed, the pattern of high- and low-risk cells for the LEP variant differed depending on the value of VEGFA variant, which provides evidence of epistasis. The interaction graph showed in figure 1B indicates the amount of information gained about case–control status by putting the two gene variants together. A low percentage of information gain of case–control status explained by LEP or VEGFA considered independently was observed (0.53% and 0.36%; respectively). However, a large percentage of information gain explained by the interaction between these two loci was observed (1.53%), which is interpreted as a synergistic or non-additive relationship between the VEGFA and LEP.

Figure 1 (A) Graphical model showing the distribution of knee OA cases and controls for each genotype combination of vascular endothelial growth factor A (VEGFA) and LEP gene variants. The cells shaded dark are high risk for disease and cells shaded light are low risk. Left bars represent the number of cases and right bars represent the number of controls. The new variable constructed by multifactor dimensionality reduction is shown on the right. The pattern of high and low risk for the VEGFA polymorphism differs depending on the value of the LEP polymorphism, which is evidence of epistasis or gene-gene interaction. (B) Interaction graph comprises a node of each gene variant and a pairwise connection between them. The percentage of entropy removed (i.e., information gain) by each variable is visualized for each node shown in a box. The pairwise Cartesian product of variables is visualized for the connection depicted as a line accompanied by a percent of entropy explained by that interaction. The bold line indicates a strong synergistic interaction. 

DISCUSSION

Our results showed no association between LEP and VEGFA gene variants with primary knee OA when their independent effects were analyzed. Nevertheless, MDR analysis provided evidence of epistasis between LEP and VEGFA, which could be implicated in the development of primary knee OA. To the best of our knowledge, there are no studies analyzing the gene-gene interaction between those genes in the context of OA. Moreover, to date, only few studies have analyzed the association of LEP and VEGFA gene variants with primary knee OA. In the Han Chinese population, three variants of LEP were assessed for their contribution to knee OA: the rs11761556 and rs12706832, located in the intron region, and the rs2071045 of the 3´UTR region. In general, the genotypes of those variants were significantly associated with a lower risk to knee OA in normal weight and overweight groups, but not in the obese group¹⁹. Two independent studies analyzed the association of VEGFA variants with knee OA in the Mexican population. In a population sample settled in West Mexico, the −460 T/C (rs833061) and +405 C/G (rs2010963), located in intron and 5´UTR regions of VEGFA, respectively, were assessed using polymerase chain reaction-restriction fragment length polymorphisms. The variants did not show significant association with knee OA. Notably, the cases had a significantly high BMI (29.9 [17.5-46.9])¹⁷. The other study analyzed the variants rs699947 and rs3025039 (located in intron and 3´UTR regions of VEGFA, respectively) in cases with knee OA who had a BMI of 29.0 ± 4.19; however, there was no association. This study also analyzed epistasis among several genes, and a non-statistically significant interaction of VEGFA with other genes such as MMP3, COL3A1, HF1AN, and EGF was found (LEP was not analyzed)¹⁸. The rs2167270 of LEP has not been previously investigated in OA. This variant was not associated with knee OA, and we only observed a non-statistically significant tendency toward an increased risk for AA homozygous. Regarding rs2010963 of VEGFA, there was no association, and we merely observed a tendency toward a lower risk. This variant was previously analyzed in the population of West Mexico, and the authors did not find an association with OA, which seems to be in line with our findings. However, according to the genotype frequencies reported in that study, there is an HWE deviation indeed (p = 0.0006), which could affect their results¹⁷. Of notice, in our study, the LEP and VEGFA variants follow HWE, and the minor allele frequencies are not importantly different to those reported in population with Mexican ancestry in The 1,000 Genomes Project (https://www.ensembl.org/index.html/). Therefore, we consider that our results are reliable.

Single-nucleotide variants (SNVs) are the most common variant type in the human genome. Genome-wide association studies have uncovered robust associations between SNVs and numerous complex diseases. However, only ~4% of the associated variants identified result in differences in the protein product. As a result, an SNV can be the “best available marker,” rather than the “probable causal variant”. SNVs in the UTRs, generally do not change the protein sequence, but may regulate gene expression. The 5´UTR possesses a consensus sequence containing the translation initiation codon, as well as regulatory elements (e.g., RNA-binding proteins sites). Thus, variants at 5´UTR can have an impact on the production of the protein²⁵,²⁶. For example, the 90 A/G variant at the 5´UTR of the tryptophan hydroxylase 2 gene alters the mRNA structure and/or RNA-protein interaction, affecting gene expression²⁷. However, the functional consequences of variants within UTRs are not always immediately obvious and often not characterized. To date, there are no functional consequences described for the rs2167270 and rs211093, both located in the 5´UTR region; nevertheless, they can be a useful genetic marker to analyze the association of LEP and VEGFA with knee OA.

LEP and VEGF-A are implicated in inflammation and angiogenesis, respectively; nevertheless, they may act synergistically, favoring both processes as demonstrated in different types of cancer²⁸-³⁰. LEP and VEGF-A are involved in the pathophysiology of OA; it is, therefore, possible that, as in cancer, they could play a synergistic role¹². During osteoarthritic process, LEP and VEGF-A can modulate the production of IL-1β, IL-6, IL-8, and different MMPs¹⁵,³¹,³². Strikingly, VEGF-A per se may cause characteristic features of OA since intra-articular injection of VEGF-A in the knee joints of mice induces OA-like changes (i.e., proteoglycan loss, cartilage calcification and degradation, bone sclerosis, osteophyte formation, and synovial hyperplasia)³³.

We are aware that our study could have some limitations. Being a hospital-based case–control study, there is an inherent possibility of incurring in selection bias; therefore, to avoid bias, we assessed incident cases with primary knee OA with their proper controls during the same period of time through specific inclusion criteria, performing the same procedures on both study groups, performing a blinded radiologic evaluation, and rigorously controlling potential confounders through exclusion criteria or through multivariate statistical analysis. Sample size could also be a disadvantage. In this regard, it may be a challenge to achieve a large sample size in hospital-based case–control studies; to compensate for this, we attempted to increase our internal validity through the strategies mentioned above. Obesity is one of the most relevant risk factors associated with OA, mainly with secondary OA. The overload conferred by obesity is an important factor, but there is also a chronic inflammatory state observed in overweight people34. Genetic associations in OA could be influenced by differences in BMI, as observed in the study in the Han Chinese population¹⁹. Subsequently, we controlled for BMI, and the participants were cases and controls which were not overweight. In consequence, we believed that our results are not influenced by overweight. In fact, our patients have lower BMI in comparison to other studies. In summary, while cognizant of our limitations, we are confident that our main strength is that our study groups were adequately selected and that we are properly including primary knee OA.

In conclusion, our results suggest that the independent effects of LEP and VEGFA gene variants are not associated with primary knee OA; nevertheless, the interaction between them is associated with the genetic susceptibility to develop primary knee OA.