Relationship of common variants in VEGFA gene with osteonecrosis of the femoral head: A Han Chinese population based association study

The pathology of non-traumatic osteonecrosis of the femoral head (ONFH) is complex. Several studies have linked some polymorphisms of vascular endothelial growth factors A (VEGFA) with ONFH, but the results are not consistent and are even conflicting. In the study, 22 single nucleotide polymorphisms (SNPs) in VEGFA were genotyped in 1,762 subjects (489 cases and 1,273 controls). Genetic association analyses were performed in single markers and haplotype levels. Stratification analysis was conducted for ONFH patients. Gene by environment interactions were tested between VEGFA and the smoking status of the subjects. Gene expression and eQTL data of significant SNPs were extracted from GTEx to examine their potential biological function. The SNP, rs2010963, was identified to be significantly associated with ONFH (χ2 = 11.66, P = 0.0006, OR = 1.29). Haplotypes including rs2010963 were also identified to be correlated with ONFH in the haplotype-based analyses. After stratifying by the causes of ONFH, a significant signal from rs2010963 could only be identified in alcohol-induced patients (Pallelic = 0.0009) but not in steroid-induced patients (Pallelic = 0.055). No significant results were obtained from the gene by environmental interaction analyses. Significant expression differences of VEGFA were identified in multiple human tissues for different genotypes of rs2010963. Our findings indicate that SNP rs2010963 is significantly associated with ONFH.

is difficult to tell whether these SNPs are surrogates of some underlying susceptible variants or genetic markers with biologically functional significance.
In this study, we aimed to investigate the genetic association between VEGFA and ONFH in a large Chinese study sample. Through genotyping several pre-selected genetic markers covering VEGFA in our study subjects, we examined the statistical association between genetic polymorphisms and ONFH in both single-marker and haplotype-based methods. In addition, combined with relevant bioinformatics tools, we aimed to examine the potential biological function of the significant SNPs identified in the association analysis.

Methods
Study Subjects. In the study, a total of 489 unrelated male patients with non-traumatic ONFH and 1,273 unrelated control subjects were consecutively recruited at the Luoyang Orthopedic Hospital of Henan Province (Luoyang, China) from 2013 to 2016. Patients were diagnosed according to assessment by X-rays, magnetic resonance imaging (MRI), and bone scans. Based on the etiological factors of ONFH, patients were divided into a steroid-induced group (254 cases) and an alcohol-induced group (235 cases). Steroid-induced ONFH was defined by a history of a mean daily dose of ≥16.6 mg or a highest daily dose of 80 mg of prednisolone equivalent within 1 year before the development of symptoms or radiological diagnosis in asymptomatic cases. When steroid and other factors were excluded, patients with a history of ethanol consumption of at least 400 ml per week for at least 1 year were categorized under alcohol-induced ONFH. Patients with a demonstrable history of direct trauma or with possible combined causes were excluded. Those who had a chronic metabolic disorder of the heart, kidney, or liver were also excluded. Control subjects were matched with patients for age and BMI and were enrolled from subjects attending routine medical checkups. The controls had no hip pain, and anteroposterior and frog-leg lateral pelvic radiographs did not show any lesions. The controls had a history of ethanol consumption of at least 400 ml per week for at least 1 year; however, they had no alcohol-induced ONFH or other related diseases, no history of thromboembolic events and no symptoms of hip disease. All participants were restricted to the Han Chinese population who lived in Luoyang city and surrounding areas. Informed consent was obtained from all groups. The study protocol conformed to the ethical guidelines of the 1975 Declaration of Helsinki and was approved by the Ethics Committee of Luoyang Orthopedic Hospital of Henan Province.

SNP Selection and Genotyping.
We searched for all SNPs with minor allele frequencies (MAF) ≥ 0.05 within the region of the VEGFA gene in the 1000 Genomes Chinese Han Beijing population (CHB). Then, MAF ≥ 0.05 with pair-wise tagging and r 2 ≥ 0.8 were used as the cut-off criteria during tag SNP selection, which generated 22 tag SNPs covering the region of the VEGFA gene for our study. Basic information on the 22 selected SNPs is summarized in Supplemental Table S1. Genomic DNA was isolated from peripheral blood using a Tiangen DNA extraction kit (Tiangen Biotech Co. Ltd, Beijing, China) according to the manufacturer's protocol. SNP genotyping was performed using a Sequenom MassARRAY platform with iPLEX GOLD chemistry (Sequenom, San Diego, CA, USA) based on the manufacturer's protocols. The results were processed using Sequenom Typer 4.0 software, and genotype data were generated from the samples 15 . Genotyping was conducted by laboratory personnel blinded to the case-control status, and the genotyping results, data entry and statistical analyses were independently reviewed by two authors. We randomly re-performed the analysis on 5% of the sample, with a concordance of 100%.
Statistical and Bioinformatic Methods. χ 2 tests were performed using Plink 16 for each marker to examine the potential association between SNPs and ONFH disease status. Genomic control (GC) was conducted to identify and correct potential false positive signals due to population stratification. Linkage disequilibrium (LD) blocks were constructed for the 22 selected SNPs, and haplotype-based association tests were performed using Haploview 17 . Significant SNPs identified by single-marker-based tests were re-analyzed in stratification analysis. In this analysis, our patients were stratified by the clinical type of ONFH (alcohol-induced or steroid-induced). To investigate the potential gene by environment interactions between selected SNPs of VEGFA and the smoking status of our study subjects, we performed G-by-E interaction analysis by fitting logistic models with a multiplying term. R project 18 was utilized for general statistical computing and G-by-E interaction analysis. Bonferroni corrections were applied to address multiple comparisons. For single-marker-based analysis, our P value threshold was 0.05/22 ≈ 0.002.
RegulomeDB was utilized to examine the potential biological functions of selected SNPs 19 . We investigated the potential effects of significant SNPs on the gene expression of VEGFA using the database of GETx (https://www. gtexportal.org/home/) 20 . Data for gene expression of VEGFA in 47 human tissues were extracted and compared among different genotypes of significant SNPs identified in association tests.

Results
Genetic association between polymorphisms of VEGFA and ONFH. The clinical characteristics of all subjects are presented in Table 1. There were no differences in age, body mass index (BMI) or smoking status between the patients and controls. All 22 selected SNPs (Supplemental Table S1) passed the Hardy-Weinberg Equilibrium (HWE) test. We identified one significant SNP, rs2010963 (χ 2 = 11.66, P = 0.0006, OR = 1.29, RegulomeDB Score = 4), as being associated with the disease status of ONFH ( Table 2). The C allele of rs2010963 is related to a higher risk of ONFH. Two other SNPs showed nominal significance but failed to persist after multiple comparison corrections. The median of the χ 2 statistics for our selected SNPs was 0.15, which is far smaller than the expected value of 0.456. Therefore, no significant population stratification could be detected in our data (the Q-Q plot is shown in Supplemental Fig. S1).
Six LD blocks were constructed based on our data (Fig. 1). The P value threshold used here was 0.05/6 = 0.008. One significant LD block was identified as being associated with the disease status of ONFH (χ 2 = 11.66,  (Table 3). Another LD block, rs699947-rs1570360, showed nominal significance (P = 0.023). Stratification analyses for rs2010963 were performed at both the genotypic and allelic level (Table 4). Interestingly, after we stratified patients by their clinical type, significant signal from rs2010963 was only identified in alcohol-induced patients (P allelic = 0.0009) but not in steroid-induced patients (P allelic = 0.055). This discordance was identified in both genotypic and allelic analysis.
No significant results were obtained through gene-by-environment-interaction analyses. The most significant signal identified was from rs3025020 (P = 0.0184, OR = 0.69). However, it failed to persist following Bonferroni correction (Supplemental Table S2).

Effects of significant SNPs on gene expression of VEGFA.
We investigated the effects of SNP rs2010963 on the expression of VEGFA by examining the eQTL data from 47 normal human tissues extracted from GTEx (Supplemental Table S3). Significant expression differences were identified in four human tissues: the adrenal gland, esophagus muscularis, pancreas and thyroid (Fig. 2). The most significant difference was in the thyroid, with P < 1 × 10 −6 . The C allele of rs2010963 was related to higher expression of VEGFA. These findings indicate that rs2010963 is an eQTL for VEGFA. In addition, SNPs capturing eQTL signals with genome-wide significance were summarized in Supplemental Table S4. As shown, six out of these 22 selected SNPs had eQTL signals with genome-wide significance.
Earlier studies have shown that this growth factor promotes the proliferation and migration of vascular endothelial cells and plays an important role in the physiology of angiogenesis. Knock-out of this gene in mice results in abnormal blood vessel formation in the embryonic stage. In addition to its role in angiogenesis, this growth factor has also been shown to be essential for the formation of endochondral bone. Gerber et al. investigated the potential role of VEGFA in the formation of endochondral bone by inactivating VEGF in 24-day-old mice 26 .
The results showed that proliferation, differentiation and maturation of chondrocytes were basically normal, but resorption of terminal chondrocytes was significantly inhibited. These findings indicate that VEGF is an essential coordinator of bone formation in the growth plate 26 . Compared to previous studies based on Asian populations, which resulted in inconsistent findings about the significant SNPs of VEGFA 12,13 , in this study, we identified rs2010963 (−634C/G or +450 C/G as indicated in some other studies) to be associated with the disease status of ONFH. Some SNPs, such as rs699947 (−2578A/C) and rs1570360 (−1154A/G), which were shown to be significantly associated with this bone disorder in other studies, were shown to be only surrogates of rs2010963 in our study. These SNPs had nominal significance with ONFH in our data only because they are in medium LD with rs2010963. Given that limited SNPs analyses were difficult to draw reliable and stable conclusions 27-31 ,our further haplotype-based analyses provided more evidence about the role of rs2010963 and its surrounding genomic regions in the susceptibility to ONFH. In our study, the C allele would increase the risk of ONFH by approximately 20% compared to the G allele. The eQTL analyses showed that rs2010963 is a potential eQTL of VEGFA and can affect the gene expression of VEGFA in multiple normal human tissues. The C allele of rs2010963 is related to higher expression of VEGFA. Combined with the results from gene expression analysis, we believe that SNP rs2010963 (−634C/G) might be more than just a surrogate of some underlying ungenotyped genetic variants but a variant with specific biologically functional significance. However, because this was a genetic association study, it is impossible for us to unravel the potential link between increased expression of VEGFA and disease risk of ONFH from this study alone, and more research is needed in the future.
One interesting result of this study is that by stratifying our ONFH patients, we found that a significant signal for rs2010963 could only be identified among alcohol-induced ONFH patients and not among steroid-induced patients. Because the alcohol-induced group had a smaller sample size compared to the steroid-induced group, this discordance cannot be explained by decreased statistical power due to a smaller sample size caused by   stratification. One potential explanation is that this difference might be caused by some undetected selection bias in our patient recruitment process. On the other hand, an alternative explanation, which might be more informative, is that this difference might indicate some deeper difference in the biological mechanisms for the two types of ONFH patients. Still, more studies with larger sample sizes are needed in the future to replicate and validate these findings. This study has several limitations. Firstly, in this study, we only recruited males as the study subjects, and this may impair the generalization of the study results. Most previous similar studies enrolled both males and females, although the ratios of the genders are always imbalanced [12][13][14] . Secondly, due to the limitations of our study design, although we extracted gene expression data from GTEx for normal human tissues, it might be more meaningful to quantify the gene expression of VEGFA and estimate the differences of gene expression between ONFH patients and controls in our own study subjects. In addition, another potential limitation is that we only included common polymorphisms (SNPs with MAF > 0.05). Low-frequency and rare variants were not included in this study. In this sense, the present study is incomplete and is unable to systematically unravel the potential genetic architecture for osteonecrosis, and our results should be considered to be preliminary and confirmed in the future research. Thirdly, we did not control the exposure to steroids in our control samples and this would be a potential flaw in study design. This flaw might be, at least partly, responsible for the non-significant findings in stratification analysis for steroid-induced patients, but the main significant results in the study would not be affected by the limitation.
In summary, our study systematically examined the genetic association between SNPs in the VEGFA gene based on study subjects of Chinese ethnic groups. Our findings indicate that SNP rs2010963 is significantly associated with the risk of ONFH, which provide additional supportive evidence of the relationship between VEGFA gene and ONFH. Comprehensive investigation with more SNPs, different populations, larger sample sizes and functional experiments are prospected to validate our results, understand the effects of VEGFA on the risk of ONFH, and elucidate the potential biological mechanisms of ONFH.