Association between EN1 rs4144782 and susceptibility of knee osteoarthritis: A case-control study

Osteoarthritis (OA) is a complex disease that affects the whole joint, resulting from the combined influence of biomechanical factors and genetic factors. The heritable component for primary OA accounts for about 60% of variation in population liability to the disease. So far, genome-wide association studies (GWAS) and candidate gene studies have established many OA-related loci. However, these findings account for only a rather small fraction of the genetic component. To further reveal the genetic architecture of OA, we conducted this case-control study to explore the association of locus EN1 rs4144782 and knee OA susceptibility in a Chinese population. EN1 rs4144782 was significantly associated with increased risk of knee OA (OR=1.26; 95% CI: 1.05-1.50, P value=0.012). In dominant model, compared with carriers of GG genotype, those with AG or AA genotype have an 1.44-fold increased risk of OA (OR: 1.44; 95% CI: 1.10-1.88; P value=0.008). Subgroup analyses didn't change the results materially. This should be the first association study of EN1 locus on risk of OA, and our finding suggested that the EN1 rs4144782 might contribute to the susceptibility of knee OA.

Despite of numerous efforts to identify genetic factors associated with OA, the identified loci could still explain a small part of the hereditary of OA, due to the disease heterogeneity and limitation of large sample size [4,13,17,27,28,[30][31][32][33][34][35][36].Recently, whole-genome sequencing identifies locus Homeobox protein engrailed-1 (EN1) gene as a determinant of bone mineral density (BMD) [37].Higher BMD was reported to be associated with OA in many cross-sectional and longitudinal epidemiological studies [38][39][40][41].One important mechanism that could underpin the observed association between BMD and OA is genetic pleiotropy-the existence of genetic variants that contribute to both BMD variation and OA risk [42].For example, Yerges-Armstrong et al [43] evaluated the associations between GWAS-identified genetic loci for variance in BMD and knee OA, and identified four SNPs significantly associated with prevalent radiographic knee OA, including rs2016266, rs10226308, rs3736228, and rs10835187.This results further supported the hypothesis that BMD, or its determinants, might be a risk factor for OA development.Thus, to clarify the potential relationship between locus EN1 rs4144782 and knee OA susceptibility, we conducted this case-control study in a Chinese population.

Clinical characteristics between OA cases and healthy controls
Totally, we included 500 OA patients and 500 matched healthy controls in this study.As shown in Table 1, the proportion of the matching variable which included age, gender, smoking status, alcohol status, and regular physical activity in cases and healthy controls was comparable (all P value > 0.05).However, BMI in OA cases was larger than that in healthy controls (P < 0.001), and it's more likely to be diabetes patients for OA cases (P value = 0.027).

Stratified analyses of associations between EN1 rs4144782 and knee OA risk
Table 3 presents stratified analyses of association between EN1 rs4144782 and knee OA risk by age, gender, smoking and drinking status, physical activity, diabetes, and BMI.Generally, the results were not changed materially.The significant ORs ranged from 1.23 to 1.28 in allelic model (P value: 0.013-0.049).While in dominant model, the significant ORs ranged from 1.34 to 1.53 (P value: 0.007-0.049).We also detected an significant association for non-smokers in recessive model (OR:1.43;95% CI: 1.00-2.05;P value = 0.048).

DISCUSSION
In this case-control study, we investigated the association between EN1 rs4144782 and susceptibility to OA in a Chinese population.Our results showed that EN1 rs4144782 was significantly associated with increased risk of knee OA.The results reconfirmed the close connection of BMD and OA.To the best of our knowledge, this should be the first association study of EN1 locus on risk of OA.
The emerging findings of susceptibility loci of many disease have provided insight into the early diagnosis in the past years [44,45].OA is considered as a multifactorial disease, resulting from the combined influence of biomechanical factors and genetic factors [46].COL2A1 was the first gene identified to be associated with OA by genetic linkage study in 1990 [47].Following, SNPs in gene VDR, ER, TGF-beta1, GDF5, and VNTR were found in sequence [48][49][50][51].Recently, Gok et al [52] first reported the CA repeat length of >/= 20 in the promoter region of ADAMTS9 gene has a six-fold increase in probability for having severe OA.The most intriguing hypothesis that high systemic BMD increased OA risk conclude that genetic variations associated with high BMD will lead to a corresponding increase in OA risk, which was confirmed by the detection of SNP rs2016266, rs10226308, rs3736228, and rs10835187 [43].
The engrailed gene plays an important role in the development of segmentation course, where it is required for the formation of posterior compartments [53][54][55].En-1 is in the central portion of the mouse chromosome 1, while En-2 is in the proximal portion of mouse chromosome 5 [56].Zhang et al [57] also indicate that EN1 is a positive regulator of brown adipogenesis.It was shown that polymorphisms (rs1861972, rs1861973) of homeobox transcription factor gene ENGRAILED 2 and its haplotype A-C was associated with autistic disorder in Chinese children [58].In current study, Our results showed that EN1 rs4144782 was significantly associated with increased risk of OA.EN1 is a homeobox gene that contributes to the development in the dorsal midbrain and anterior hindbrain of human beings [59].
SNP rs4144782 was located in the intron region of EN1 gene.Using HaploReg v4 software, we found that a G to A transform will result in the expression decrement of motifs CAC-binding-protein, EWSR1-FLI1, PU.1_disc3, SP1_disc3, SP1, STAT_disc7, TATA_disc7, and TFII-I, as well as increased expression of motifs HDAC2_disc6, TCF12_disc5, and VDR_2 [60].RegulomeDB also shows that allele exchange of EN1 rs4144782 will affect the binding efficiency of many transcription factors [61].These will then affect the expression of En1 gene in sequence.
There are several limitations that can be derived from the present study.Although the sample size of the current study was moderate, the number of subjects on the interaction analyses was relatively small.In addition, all the patients and controls involved in this study were enrolled from the hospital, thus inherent bias may affect the results.Also, this study just considered common variation (minor allele frequency > 0.05) like other candidate gene studies.Rare variants, which owns larger effect sizes than the average effect size of common variants, should be explored for its contribution to the development of knee OA.Considering the different genetic architecture (different minor allele) of EN1 rs4144782 in different ethnicities (Asian, European and American), results should be interpreted discreetly and Conclusively, in this case-control study, we identified that EN1 rs4144782 was significantly associated with increased risk of knee OA.This conclusion contributes to the hypothesis for an association of high BMD-associated SNPs with higher knee OA risk, which will help to elucidate the connection of BMD and OA, and to perform risk prediction, early prevention, diagnosis and individualized treatment of OA.Further work is needed to determine the precise role of EN1 locus on the occurrence of OA, and how rs4144782 regulates the function of EN1 gene and regulates its downstream genes, related pathways in vitro and in vivo.

Study population
The case-control study was conducted in accordance with STROBE (strengthening the reporting of observational studies in epidemiology) guidelines (STROBE checklist in Supplementary File 1).The sample size was calculated using Quanto software [62].Finally, OA cases (n=500) were recruited at Renmin Hospital of Wuhan University (Wuhan, Hubei province, China) between January 2008 and April 2015 in current study.The diagnosis of knee OA was based on having a K/L grade of 2 or higher in at least one knee [63].The representative X-ray images were presented in Supplementary File 2. Healthy controls (n=500)

DNA extraction and genotyping
Peripheral venous blood (2ml) from all participants was collected by vacuum tube with ethylene diamine tetraacetic acid (EDTA).Genomic DNA was extracted from the whole blood samples of all participants using QIAamp DNA Blood Mini kit (Qiagen, Berlin, Germany) according to the manufacturer's instructions.The DNA concentration and purity were estimated by measuring ratio of the absorbance at 260 and 280 nm.The tagSNP rs4144782 was selected using the web-based software SNP info (https://snpinfo.niehs.nih.gov/).Genotyping was done with Sequenom iPLEX matrix assisted laser desorption/ionization-time-of-flight mass spectrometry technology using 15 ng DNA.The MassARRAY Typer 4.0 software was used for proper data acquisition and analysis.A random 10% of quality control samples were checked for the concordance.Laboratory personnel were blinded to the case-control status of the specimens, and to the quality control samples.

Statistical analyses
All analyses were performed using SAS version 9.2 and P values were two-sided in current study.The Hardy-Weinberg equilibrium (HWE) was assessed using a goodness-of-fit χ2 test performed to identify possible genotyping errors among the controls of each study.Categorical variables were compared between patients with knee OA and healthy controls using the χ2 test or Fisher's exact test where appropriate.Logistic regression was used to estimate ORs and 95% confidence intervals (CIs) as a measure of the association with the susceptibility of knee OA, adjusted for the potential confounding factors, including age, gender, smoking and drinking status, physical activity, diabetes, and BMI.We also conducted subgroup analyses by age, gender, smoking and drinking status, physical activity, diabetes, and BMI.