Putative functional variants of lncRNA identified by RegulomeDB were associated with knee osteoarthritis susceptibility

Background Knee osteoarthritis (KOA) is the most common form of chronic degenerative joint disease worldwide. Its incidence has increased in recent years. Aberrant expression profile of lncRNAs in damaged bone and cartilage of KOA patients has been reported recently, indicating its potential contributions in KOA development and a promising target for disease diagnosis and treatment. The aim of this study was to identify the association between genetic variation in lncRNA and KOA. Methods We retrieved relevant articles from the PubMed, Medline and Embase databases up to Jul 2017 investigating the association between lncRNA and the risk of osteoarthritis. There are 15 lncRNAs which show connection with osteoarthritis. We selected potential functional polymorphisms identified by RegulomeDB database in these lncRNAs. A case-control study was conducted which contained 278 KOA patients and 289 OA-free controls. Results Logistic regression analyses revealed that H19 rs2067051 T allele was significantly associated with decreased risk of KOA after adjusted for age, gender and BMI in recessive genetic model (OR = 0.63, P = 0.03) and additive genetic model (OR = 0.79, P = 0.03). MEG3 rs4378559 T allele was significantly associated with increased risk of KOA in additive genetic model (OR = 1.32, P = 0.04). Heterogeneity tests proved that H19 rs2067051, MEG3 rs4378559 and HOTTIP rs202384’s risk effects on KOA were more remarkable for female, BMI ≥ 25 and younger age (age < 60), respectively. Conclusion The results indicate that potential functional genetic variation in lncRNA plays an important role in the pathogenesis of KOA.


Background
Knee osteoarthritis (KOA) is the most common joint disease worldwide, characterized by progressive degeneration of articular cartilage, synovitis, osteophyte formation, and subchondral bone sclerosis [1]. It can cause severe pain and physical disability thus substantially reduce elder people's quality of life [2,3]. Although the incidence of osteoarthritis is still increasing, affecting approximately 10% of men and 18% of women over 60 years of age, its pathophysiology is still evolving and undetermined [4]. Twin-pair studies and family-based segregation analyses have showed clear evidence of a heritable component in the susceptibility of OA [5]. Researches have shown the evidence that genetic factors may play a vital role in the development of OA, although little of them have been identified so far.
Long noncoding RNAs (lncRNA) have been tentatively defined as a type of RNA molecule more than 200 nucleotides in length (while microRNAs are of 20-22 nucleotides in size) and are characterized by their complexity and diversity of sequences and mechanisms of action [6]. Several lncRNAs have been reported play an important role in the pathogenesis of osteoarthritis [7,8]. Aberrant expression profile of lncRNAs in damaged bone and cartilage of OA patients has been reported recently, indicating its potential contributions in OA development and a promising target for disease diagnosis and treatment [9]. However, the role of lncRNAs played in cartilage metabolism and their overall contributions to the degradation of chondrocyte extracellular matrix and the pathogenesis of OA are still not fully understood. Considering the important roles lncRNA played in cartilage anabolism and catabolism, we hypothesized that single nucleotide polymorphisms (SNPs) in lncRNA gene may individually or jointly contribute to the risk of osteoarthritis.
RegulomeDB, a database integrates a big collection of regulatory information from ENCODE and other data sources, is a powerful tool to choose functional SNPs in a specified chromosome region [10]. RegulomeDB presents a scoring system with categories ranging from 1 to 6 by the way of integrated annotations data on methylation, chromatin structure, protein motifs and binding. The lower RegulomeDB score indicates the stronger evidence for a variant to be located in a functional region.
Recently, many studies have used RegulomeDB to identify causal polymorphisms associated with human diseases such as cancer, polycystic ovary syndrome and so on [11,12]. Thus in this case-control study, we explored KOA risk with lncRNA polymorphisms which are annotated in RegulomeDB.

Patients
In this study, Han Chinese KOA patients were recruited between June 2013 and August 2016 at the orthopaedics department, Changzhou No.1 people's hospital according to the American College of Rheumatology (ACR) diagnostic criteria for KOA [13]. Each patient had taken the affected knee's anteroposterior weight-bearing radiographs. Two orthopedics doctors provide a Kellgren-Lawrence (KL) score ranging from 0 to 4 assessed by the radiographic data [14]. Other knee joint diseases were excluded such as inflammatory arthritis (septic arthritis, rheumatoid, or autoimmune disease), posttraumatic arthritis, or knee joint developmental dysplasia. All healthy volunteer controls were frequency-matched to the case subjects by gender and age. The healthy control subjects reported no history of OA or other rheumatic diseases, which were included from the same hospital during the same period. Each participant was interviewed face-to-face to gather demographic data by two trained investigators. After signing informed consent at recruitment, peripheral blood of each participant was collected. All subjects were weighed with a calibrated beam balance to 0.1 kg, wearing the least possible clothes. Also their height was measured in centimeters using a stadiometer. The body mass index (BMI) was then calculated using the above data. This study was approved by the Human Research Ethics Committees of the Changzhou No.1 people's hospital.
Variant with lower RegulomeDB score indicates the variant is more strongly located in a functional region. Consequently, we selected SNPs with RegulomeDB scores ranging from 1 to 2b. Based on the data from UCSC database (GRCh37/hg19), regulomeDB database, the criteria of minor allele frequency (MAF) > 0.05 and linkage disequilibrium (LD) < 0.8 in Han Chinese, we found 10 potentially functional SNPs in these 15 lncRNAs.

DNA extraction and genotyping
Genomic DNA was extracted from leukocyte pellets by proteinase K digestion and phenol chloroform as described in previous study [20]. The genotype detecting was performed by Sequenom MassARRAY assay without knowing the status of case or control according to instructions of the manufacturer. Three SNPs genotyping were failure due to probe design. Finally, the left 7 SNPs successfully genotyped call rate were all above 95%. Over 10% specimens were randomly slected to be repeated with over 99% consistency and two blank (water) controls were tested in each 384-well plate for quality control.

Statistical analyses
The χ 2 test and Student's t-tests were used to analyze distribution differences of demographic characteristics and genotypes between cases and controls, for categorical variables and continuous variables respectively. Hardy-

3-4 143
Bold font presents P < 0.05 Weinberg equilibrium (HWE) was tested using a goodness-of-fit χ 2 test to compare the expected and observed genotype frequencies in controls for the distribution of each SNP. We use logistic regression to estimate odd ratios (ORs) and 95% confidence intervals (CIs) to evaluate the association with KOA susceptibility after adjusted for age, sex and BMI. We use the χ 2 -based Q-test to measure the heterogeneity (ORs and 95% CIs) derived from corresponding subgroups to examine the differences between different subgroups. We also assessed the cumulative effects of the genotyped 6 SNPs using a risk score analysis with a linear of the SNP genotypes (coded as 0, 1, and 2) weighted by the regression coefficient. All of the statistical analyses were

Results
The distribution of selected clinical variables between OA cases and controls was summarized in Table 1. In brief, the age and gender between two groups were comparable (P > 0.05). Rs2067079 of GAS5 was excluded due to violation of Hardy-Weinberg equilibrium. The genotype distributions of the left 6 SNPs and their associations with KOA risk were presented in Table 2. The success rates of genotyping for all these SNPs were above 98%. Logistic regression analyses revealed that H19 rs2067051 T allele was significantly associated with decreased risk of KOA adjusted for age, gender and BMI in recessive genetic model (OR = 0.63, 95% CI: 0.42-0.97, P = 0.03) and additive genetic model (OR = 0.79, 95% CI: 0.64-0.98, P = 0.03). MEG3 rs4378559 T allele was significantly associated with the increased risk of KOA in additive genetic model without adjustment (OR = 1.32, 95% CI: 1.01-1.74, P = 0.04), but the association no longer existed after adjustment (P = 0.08). HOTTIP rs2023843 C allele showed boundary positive in additive genetic model after adjustment (OR = 1.23, 95% CI: 0.96-1.58, P = 0.10).
Further stratified analyses were conducted in four positive models as above mentioned (Table 3). Results showed that H19 rs2067051 T allele on KOA was more evidenced for female in recessive model (TT vs.TC+ CC) (P for heterogeneity was 0.01). MEG3 rs4378559 T allele on KOA risk was more obvious in BMI ≥ 25 in additive model (P for heterogeneity was 0.02). HOTTIP rs2023843 C allele on KOA was more distinct for younger age (age < 60) in additive model (P for heterogeneity was 0.03).
Cumulative effects of the 6 SNPs were conducted using a risk score analysis weighted by the regression coefficient ( Table 4). The mean of cumulative risk score among KOA cases (0.38 ± 0.40) was higher than that among controls (0.29 ± 0.34) (P value for T test < 0.01). Logistic regression analyses revealed that risk score was significantly associated with increased risk of KOA adjusted for age, gender and BMI (OR = 2.21, 95% CI: 1.36-3.60, P < 0.01).

Discussion
In the current study, we conducted a case-control study to explore the association between genetic variants identified by RegulomeDB in candidate lncRNA genes and KOA. Our results confirm that putative functional variants H19 rs2067051 and MEG3 rs4378559 were associated with KOA susceptibility. Heterogeneity existed between different gender, BMI and age group for H19 rs2067051, MEG3 rs4378559 and HOTTIP rs202384 respectively, which means gene-environment interactions between genetic variants in lncRNA genes and these clinical data for KOA risk.
All RegulomeDB score of H19 rs2067051, MEG3 rs4378559 and HOTTIP rs2023843 were 2b. ChIP-seq data showed these polymorphisms binding to many proteins including EZH2, E2F6, REST and IKZF1 (http:// regulome.stanford.edu/). The development of osteoarthritis was ameliorated by inhibition of EZH2 through the Wnt/β-catenin pathway [21,22]. E2F6 encodes a member of transcription factors which play an important role in the cell cycle controling [23]. A transcriptional repressor which represses neuronal genes in non-neuronal tissues was encoded by REST [24]. IKZF1 encodes a transcription factor which belongs to the zinc-finger DNA-binding protein family, which correlates with chromatin remodeling [25]. IKZF1 participated in inflammation suggested that IKZF1 may contribute to osteoarthritis [26]. A 2.5 kb RNA Polymerase II dependent transcript was expressed from H19 and it is spliced, capped, polyadenylated, and finally exported into the cytoplasm [27]. MiR-675, which suppresses growth, is developmentally reversed in H19 [28]. It also participates in the development of OA through affecting the pathological processes including inflammatory response, extracellular matrix (ECM) disruption, angiogenesis and apoptosis. Research showed that H19 was up-regulated in OA compared with normal tissue and was a metabolic marker in damaged cartilage of OA patients or in cultured chondrocytes under hypoxic signaling condition [29,30]. Rs2067051 of H19 has been verified connected with birth weight and coronary artery disease [31,32]. Both high birth weight and coronary artery disease have relationship with higher BMI, which are crucial risk factors for KOA. In this study, we also found heterogeneity was existed for H19 rs2067051 in different gender group. It is already found that H19 represents a key factor which causes sex differences in the incidence of cholestatic liver injury in mice whose multidrug resistance 2 gene was knockout [33].
A maternally expressed lncRNA, which was named MEG3, was closely connected to inflammation-related diseases, for example knee osteoarthritis [34]. It revealed that MEG3 demonstrated low expression level in OA cartilage samples in rat model. It was also shown that MEG3 was down regulated and negatively correlated with VEGF expression levels in cartilage samples from knee osteoarthritis patients [34]. MEG3 knockdown could promoted chondrocytes proliferation and inhibited chondrocytes apoptosis induced by IL-1β possibly through miR-16/SMAD7 system [35].
HOTTIP, which participate in the process of osteoarthritis advance and endochondral ossification, showed higher expression level in OA cartilage in comparison to normal cartilage [36]. Through modulating intergrin-α1 either transcriptionally by HOXA13 or epigenetically by DNMT-3B, HOTTIP control cartilage destruction and development. Research data also indicated that HOTTIP could be potential diagnostic biomarker or therapeutic target for OA and other joint cartilage related disease.

Conclusions
In conclusion, our study suggested that potential functional genetic variation in lncRNA plays a crucial role in the pathogenesis of KOA. Further independent studies with different races and larger sample sizes are necessary to elucidate the molecular mechanisms underlying these findings.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Authors' contributions KJW and MJC drafted the first manuscript. KJW, MJC, WGD and QJ participated in study design, data acquisition, analysis and interpretation of data, critical review and final approval. KJW, WGD and QJ contributed with clinical expertise, conduct of patient visit procedures, data collection, data interpretation and final approval of the manuscript. KJW and MJC performed all statistical analyses. KJW, MJC and QJ contributed with data analysis and interpretation, critical review, and final approval of the manuscript.
Ethics approval and consent to participate This study was approved by the Human Research Ethics Committees of the Changzhou No.1 people's hospital. All participating subjects signed a written, informed consent to participate in the study.

Consent for publication
Not applicable.

Competing interests
Qing Jiang is a member of the editorial board (Section Editor) of BMC Musculoskeletal Disorders. The authors declare that they have no competing of interests.

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