Genome-Wide Differentially Methylated Region Analysis to Reveal Epigenetic Differences of Articular Cartilage in Kashin–Beck Disease and Osteoarthritis

Introduction

Kashin–Beck disease (KBD) is an endemic and chronic osteoarthropathy, widely affects the population in certain areas of Russia, North Korea, and China (Stone, 2009). According to China Health and Family Planning Statistical Yearbook 2016, over 60,000 people are currently affected by KBD, and over 1.1 million people are potentially affected by KBD (Commission, 2016). In contrast, osteoarthritis (OA) is also a chronic and degenerative joint disease, and affects the population in the worldwide, over 250 million people are currently affected by OA (Mora et al., 2018). Although KBD and OA share several common characteristics in manifestation and pathological of articular cartilage, extensive studies have shown the significant differences in the etiology and molecular mechanism (Schepman et al., 2011).

KBD is diagnosed as an environmental-affected disease, and more than 50 potential environmental risk factors have been investigated, including iodine and selenium deficiency (Moreno-Reyes et al., 1998), the fungal mycotoxins (Peng et al., 1992; Yang et al., 2017), and the humic acids in drinking water (Guo et al., 2014), etc. The degenerative changes in KBD first appear in the deep zone (approximately 30–40% of articular cartilage thickness) of the cartilage and starts early in children as young as 2- or 3-year-old, while OA is diagnosed as the age-related disease and frequently compounded by various risk factors and biological age-related changes, including high-impact sports (Saxon et al., 1999), cardiovascular disease (Visman et al., 2019), and diabetes (Piva et al., 2015), etc. The degenerative changes in OA first appear in the superficial (i.e., tangential) zone (approximately 10–20% of articular cartilage thickness), and increases rapidly between the 50- and 75-year olds. Therefore, clarifying the pathogenetic similarities and differences between KBD and OA may provide novel clues for understanding the distinct pathogenesis of the both.

DNA methylation is associated with transcriptional regulation, a major epigenetic factor often involving in gene regulation and cell differentiation during chondrocyte differentiation (Compere and Palmiter, 1981; Miranda-Duarte, 2018). For example, methylcytosine dioxygenase TET1 is an important enzyme in catalyzing cytosine demethylation in chondrocytes (Rice et al., 2020); the KBD-specific differentially methylated site cg25325949, located in the HAPLN1 gene body, is a key component of cartilage extracellular matrix (Wang et al., 2017), etc. However, most studies of DNA methylation profiling have shown that functional changes occur in particular regions, such as promoters, CpG islands (Jaenisch and Bird, 2003) or CpG island shores, and replication of differentially methylated regions (DMRs) are often more successful than that changes at individual sites (Ventham et al., 2016). In addition, statistically, we found that the statistical power in testing a single CpG site might be too small to detect, especially for small sample size studies (e.g., five cases vs. five controls). In contrast, interrogating the CpG sites within a particular region rather than individual CpG site would be able to improve the sensitivity and specificity for testing epigenetic alterations in tissue differentiation (Yan et al., 2020).

In this paper, we performed a comprehensive differential methylation region analysis on genome-wide DNA methylation data from five KBD patients, five OA patients, and five healthy subjects. We first illustrated that the DNA methylation levels of OA is highly correlated with age, while KBD is not. Next, we performed the paired differentially methylated region analyses to show the distinct epigenetic alternations in KBD and OA. Finally, we annotated the DMRs to the associated differentially methylated genes (DMGs) to examine the exact role of DNA methylation in gene expression, and further performed the enrichment analysis to uncover the distinct functional pathways in KBD and OA.

Materials and Methods

Sample Collections

The knee cartilage specimens of five KBD patients, five OA patients, and five healthy controls were collected from the femoral condyles of knee joints undergoing total knee joint arthroplasty, and the normal knee cartilage specimens were collected from subjects undergoing amputation caused by traffic accidents. All study subjects were Chinese Han and matched with age and sex, and with no other bone or joint disease history. All subjects were dissected and rapidly frozen in liquid nitrogen. DNA was then extracted from the cartilage specimens using QIAamp DNA Mini Kit (QIAGEN, Germany) following the standard protocol. Agarose gel electrophoresis was employed to evaluate the quality of extracted DNA.

Genome-Wide DNA Methylation Profiling

Genome-wide DNA methylation profiling was performed using the Illumina Infinium HumanMethylation450 BeadChip, resulting in 485,577 methylation sites across the whole genome. According to the standard protocol for the EZ DNA Methylation Kit (Zymo Research, USA), 500 ng DNA was used for bisulfite conversion, amplified, hybridized to the HumanMethylation 450 array, stained, and washed. The raw image intensities were scanned using iScan SQ scanner (iScan System, Illumina USA). Raw DNA methylation data were processed using Minfi package (version 1.32.0) (Aryee et al., 2014). Specifically, methylation probes with intensities indistinguishable from that of the background (detection p-value > 0.01) in more than one sample were excluded. Besides, SNP probes, cross-reactive probes, probes with a bead count <3 in at least 5% of samples and the probes located on the X and Y chromosomes were also removed, resulting in 405,389 total probes, with each probe corresponding to a CpG site. Those 405,389 probes could be further annotated to a total of 19,553 genes.

DNA Methylation Level Measurements

To perform the methylation data using existing methods (e.g., limma), DNA methylation level measurement is a necessary step prior to downstream analysis. There are two commonly used DNA methylation measurement methods, β − value and M − value (Du et al., 2010), to measure the methylation levels from microarray-based Infinium methylation assay. M − value is a logit-transformed of β − value, i.e., $M-value=log(\frac{\beta }{1-\beta })$. The range of β − value is from 0 to 1 and shows the more intuitive biological interpretations, while the range of M − value is from −∞ to +∞ and shows more statistical stabilities. The previous study utilized the β − value to perform the differentially methylated analysis, showing β − value severely compressed at the low and high when compared with M − value, as well as its bimodal distribution was less obvious. Therefore, M − value may provide more insights into the methylation level difference between case and control across the genome, but it is difficult to visualize such difference using β − value measurement (Supplementary Figure 1). Therefore, we used M − value as the measurements for methylation level for downstream analysis.

Differentially Methylated Regions

To identify DMRs that were related to KBD and OA patients, we first performed the empirical Bayes moderated t-test implemented in limma package (version 3.42.2) and controlled the effects of age and sex as covariates. We then performed ProbeLasso function (Butcher and Beck, 2015) implemented in the R/Bioconductor package ChAMP (version 2.16.2) (Tian et al., 2017) to detect DMRs. The parameter settings of ProbeLasso were set as default, i.e., the lasso radius is equal to 375, the minimum separation (bp) between neighboring DMRs is equal to 1,000, and the minimum DMR size (bp) is equal to 50. For significant DMSs (differential methylated sites), we used a threshold of false-discovery rate (FDR) < 0.05 after correcting by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995), while the DMRs were defined as significant regions if the DMR contained more than one site, and the FDR was <0.05. We defined two criteria to select the significant DMSs: (1) a stringent criterion required both statistical and biological significance (Lacey et al., 2019), i.e., an absolute percent methylation difference (PMD) of >20% in the paired groups; (2) a less-stringent criterion only required statistical significance, i.e., it only depended on statistically significant FDR. We defined that the DMRs were overlapped if the width of the overlapped part exceeded 50% compared with the narrower region. We measured the overlapped conditions between KBD-associated DMG/DMR and OA-associated DMG/DMR using the Jaccard index:

$\text{Jaccard }\left(\text{A},\text{ B}\right)=\frac{\left|\text{A}\cap \text{B}\right|}{\left|\text{A}\cup \text{B}\right|}$

where |A| denotes the number of elements in set A.

Gene Set Enrichment Analysis

The DMGs are defined as the genes that were annotated to the differentially methylated regions with PMD larger than 20%, which had additional biological significance. With the significant DMGs, we then performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using enrichGO (a total of 10,271 GO terms) and enrichKEGG (a total of 229 pathways) functions implemented in R package ClusterProfiler (version 3.14.3) (Yu et al., 2012). We used the 19,553 genes covered by the microarray as the background genes, and set the minimal and maximal sizes of genes annotated by Ontology term for testing as 10 and 500, respectively. The significant GO terms or KEGGs were selected with the q-values <0.05.

Results

Clinical Characteristics of the Study Patients

The study involved five KBD patients, five OA patients, and five healthy controls, measured by genome-wide DNA methylation profiling. The average ages of KBD patients, OA patients, and normal individuals were 57.4 ± 7.13 years old, 64.6 ± 5.18 years old, and 51 ± 7.31 years old, respectively. The detailed clinical characteristics of the studied people are summarized in Table 1.

TABLE 1

Table 1. Basic characteristics of study subjects for genome-wide DNA methylation profiling in this study.

DNA Methylation Levels Is Highly Correlated With Age in Osteoarthritis but Not in Kashin–Beck Disease

We first examined the relationship between OA/KBD and ages. To do so, we performed a standard linear regression model on the methylation levels as the response variable, and the disease status as the predictor variable, and controlled ages as the covariate to examine whether the number of identified DMSs will reduce or not. As expected, for KBD vs. Controls, when age was excluded in the regression model, 12,666 DMSs were detected; however, if age was included as a covariate, the number of DMSs were reduced to 4,325 (i.e., roughly 65.9% decreased). In contrast, for OA vs. Controls, when age was not included as covariate, 3,499 sites were detected, while only one site was detected if age was included as covariate (i.e., roughly 100% decreased). No DMSs detected in OA vs. Controls was presumably due to OA is an age-related disorder (Anderson and Loeser, 2010). In addition, we performed the Pearson correlation between age and the genome-wide DNA methylation sites. The median absolute correlation was 0.34 (mean value was 0.34) for KBD, while the median absolute correlation was 0.67 (mean value was 0.66) for OA (Figure 1A and Supplementary Figure 2). Note that the age between OA and Control showed significant difference (i.e., 64.6 ± 5.18 vs. 51 ± 7.31; two sample t-test p-value = 1.1 × 10⁻²). Consequently, if performing the differential methylation analysis did not account for the age, most of the DMSs were detected in OA vs. Controls due to the age discrepancy between OA vs. Controls.

FIGURE 1

Figure 1. The differences of Kashin–Beck disease (KBD) and osteoarthritis (OA) in differentially methylated regions (DMRs) and differentially methylated genes (DMGs). (A) Violin plots show the Pearson correlations between age and the significant DMSs without accounting for age as covariate. OA shows the higher correlations while KBD displays the lower correlations (either positive or negative). (B) Venn diagram shows the overlap of the DMRs between KBD vs. Controls, OA vs. Controls with large absolute percent methylation difference (PMD) (i.e., larger than 20%). (C) Venn diagram shows the overlap of the corresponding DMRs annotated to DMGs with large PMD (i.e., larger than 20%) between KBD vs. Controls, OA vs. Controls. (D) The bubble plot to show CpG sites harbored in the DMRs of KBD vs. Controls. (E) The bubble plot to show CpG sites harbored in the OA vs. Controls. PMD, absolute percent methylation difference.

In addition, we extracted the top 10 methylation principal components based on M-values to examine the batch differences among KBD, OA, and Control samples. A heatmap showed that the batch differences among KBD, OA, and Control samples were not significant (Supplementary Figure 3).

Kashin–Beck Disease Associated Differentially Methylated Regions Has Higher Proportion Than Osteoarthritis-Associated DMRs in Gene Body Regions

Next, to investigate the differences of OA and KBD in DMRs, we performed three groups of DMR analyses, i.e., KBD vs. Controls, OA vs. Controls, and KBD vs. OA. Consequently, 3,009 DMRs (19,408 CpG sites harbored) were identified for KBD vs. Controls with a threshold FDR < 0.05, including 1,012 hypermethylated regions and 1,997 hypomethylated regions; for OA vs. Controls, 1,063 DMRs (6,467 CpG sites harbored) were identified, including 536 hypermethylated regions and 527 hypomethylated regions; and for KBD vs. OA, 399 DMRs (2,939 CpG sites harbored) were identified, including 290 hypermethylated regions and 109 hypomethylated regions. Most of the CpG sites were located at the body region and intergenic region (IGR). Such results were consistent with previous studies (van den Dungen et al., 2017). For example, there were 3,958 total CpG sites in open sea of the body region and 3,527 CpG sites in open sea of the IGR region (Figures 1D,E and Supplementary Figure 4). Methylation that occurred in the gene body regions could affect gene expression by stimulating transcription elongation or alternative splicing (Maor et al., 2015), and generally positively correlated with gene expression (Yang et al., 2014); therefore, it was crucial for gene regulation and genomic stability (Jones, 2012). The overlapping of DMRs between KBD and OA was 719 (Jaccard index = 0.103; Fisher's exact test: p − value < 2.2 × 10⁻¹⁶). The KBD- and OA-specific regions were 2,290 and 344, respectively (Supplementary Figure 5A). We observed that the proportion of hypomethylated regions of KBD vs. Controls (66.4%) located in the gene body regions was higher than the OA vs. Controls (50.4%), presumably implying that the epigenetic factors may repress more gene expression in KBD than OA. Besides, for KBD vs. OA, 49.5% hypomethylated regions were located in gene body regions, while 33.3% hypermethylated regions were located in gene body regions.

In addition, we further detected DMRs using the stringent criterion that required both statistical and biological significance. With this criterion, we identified 131 DMRs (548 CpG sites harbored) in the KBD vs. Controls, including 62 hypermethylated regions and 69 hypomethylated regions (Supplementary Table 1); 58 DMRs (293 CpG sites harbored) in the OA vs. Controls, including 37 hypermethylated regions and 21 hypomethylated regions (Supplementary Table 2), and four DMRs (24 CpG sites harbored) in the KBD vs. OA group, including three hypermethylated regions and one hypomethylated region (Supplementary Table 3). The overlapping of DMRs between KBD and OA was 44 (Jaccard index = 0.303; Fisher's exact test: p − value < 2.2 × 10⁻¹⁶). The KBD- and OA-specific DMRs were 87 and 14, respectively (Figure 1B). We observed that the proportion of KBD-associated DMRs (67.9%) located in the gene body regions is higher than that of OA-associated DMRs (32.8%), presumably implying that the epigenetic factors played a more important role in KBD than OA. We did not further examine the proportion of DMRs of KBD vs. OA since we only identified four DMRs.

Kashin–Beck Disease and Osteoarthritis-Associated Differentially Methylated Genes Were Enriched in Different Functional Pathways

The third analysis was performed to examine the crucial role of regulation mechanism of methylation regions associated with expressed genes. To do so, we inferred the detected DMRs to the genes covered by the array. With no stringent criterion, the associated DMRs of KBD vs. Controls (including 813 hypermethylated and 1,549 hypomethylated), OA vs. Controls (including 437 hypermethylated and 448 hypomethylated) were 2,362 genes and 885 genes, and 386 genes (KBD vs. OA, including 285 hypermethylated genes and 101 hypomethylated genes), respectively (Supplementary Figure 5B). Compared with OA vs. Control, more CpG sites were differentially methylated with KBD vs. Control. As a consequence, more genes were involved. Our results were largely consistent with previous studies (Duan et al., 2010; Wang et al., 2017). This was because OA was a classic age-related disorder (Anderson and Loeser, 2010; Loeser, 2017), while epigenetic factors largely contributed to the development of KBD (Lü et al., 2011). With the stringent criterion, the associated DMRs of KBD vs. Controls (including 51 hypermethylated and 59 hypomethylated), OA vs. Controls (including 14 hypermethylated and 31 hypomethylated), and KBD vs. OA were 110 genes, 45 genes, and three genes, respectively (Figure 1C). A total of 75 unique DMGs were associated with KBD, while 10 unique DMGs were associated with OA, 35 genes were overlapped between KBD and OA (Jaccard index = 0.291; Fisher's exact test: p − value < 2.2 × 10⁻¹⁶) and one gene overlapping between KBD samples and KBD vs OA samples (Jaccard index = 0.009; Fisher's exact test: p − value = 2.2 × 10⁻²; Figure 2C). All of the common DMGs shared the same methylation direction in KBD vs. Control and OA vs. Control, implying the similar pathogenic mechanism of KBD and OA. Among those common shared 35 DMGs, Homebox (HOX) family genes play a key role in regulating the chondrogenesis and, therefore, was required for chondrocyte proliferation (Yu et al., 2007); GATA3 was a crucial gene regulating the chondrocyte differentiation contributed to both KBD (Wang et al., 2017) and OA (Goldring, 2012). For KBD-specific DMGs, HDAC4 is a crucial gene that regulates growth plate chondrocyte differentiation (Chen et al., 2016), while one of the basic pathological feature of KBD is a focal necrosis of chondrocytes in the hypertrophic zone of growth plate cartilage (Ren et al., 2007). Another example is SMOC2, which plays an important role in bone mineralization, cell–matrix interactions, collagen binding, and bone remodeling (Brekken and Sage, 2000). The DMRs harbored in HDAC4 and SMOC2 are shown in Figures 2A,B. Besides, the KBD-specific DMGs, the OA-specific DMGs, and SULF1 encodes the enzyme Sulf1, which helps in establishing GAG sulfation in the endoplasmic reticulum and also regulates Wnt signaling through desulfation of cell surface HSPGs (Dhoot et al., 2001), while OA is associated with changes in GAG expression levels, sulfation patterns (Sauerland et al., 2003), and Wnt signaling pathway (Wang et al., 2019). A previous study also showed that SULF1 expressed significantly higher in the OA cartilage compared with normal human articular cartilage (Otsuki et al., 2008). All 14 putative marker genes with strong evidence for KBD and OA are shown in Table 2.

FIGURE 2

Figure 2. The distribution of methylation levels in DMGs and enriched significant pathways of both KBD and OA. (A) DMR annotated to DMG HDAC, which is located in chromosome 2 (position: 240,014,758–240,028,276). (B) DMR annotated to DMG SMOC2, which is located in chromosome 6 (position: 168,925,150–169,002,604). The x-axis represents the coordinate of chromosome and y-axis represents PMD after smoothing. The green line represents the PMD between KBD vs. Controls and the orange line represents the PMD between OA vs. Controls. The blue bars represent the CpG sites. Note that the green line and orange line would be negative at the location where the CpG sites are sparse. (C) Top 10 significantly enriched gene ontology (GO) terms of the DMGs for KBD vs. Controls, and OA vs. Controls. The size of dots represents the q-values of the GO, and the color of the blocks represents the ratio of genes enriched in the GO terms. PMD, the absolute percent methylation difference.

TABLE 2

Table 2. Kashin–Beck disease (KBD) and osteoarthritis (OA) associated differentially methylated genes (DMGs).

Finally, we further performed gene set enrichment analysis (GSEA) of DMGs (details in the Materials and Methods section). The KBD-associated DMGs were enriched in a total of 51 unique GO terms and one KEGG pathway (Supplementary Table 4), while the OA-associated DMGs were enriched in a total of 126 unique GO terms and one KEGG pathway (Supplementary Table 5). The KBD- and OA-shared DMGs were mainly enriched in GO terms or KEGG pathways that are directly related to the skeletal system and limb-associated pathways among the GO biological processes including embryonic organ development, skeletal system morphogenesis, limb morphogenesis, chondrocyte development, etc. (Supplementary Tables 4, 5, Figure 2C, and Supplementary Figure 6). Besides, we found that KBD-specific DMGs were enriched in wound response and healing pathways, and coagulation associated pathways (Supplementary Table 4). For example, regulation of response to wounding (GO:1903034) is affected by trace elements (Abou Shady et al., 1991; Chen et al., 2017). These results were consistent with the previous differential gene expression analysis finding that several KBD-associated genes were also associated with wound healing and repair (Duan et al., 2010; Wang et al., 2012). However, OA-specific DMGs were enriched in ion-associated GO molecular functions such as gated channel activity and calcium-activated potassium channel activity, for example, potassium channel activity (GO:0005267) and calcium-activated cation channel activity (GO:0005227). Besides, our results showed that DMGs associated with OA were enriched in several ion channel pathways, which was also consistent with the previous studies that ion channels were potential biomarkers for OA (Lewis and Barrett-Jolley, 2015) (Supplementary Table 5).

Discussion

In this paper, we performed extensive differentially methylated region (DMR) analysis to reveal the pathogenetic similarities and differences between KBD and OA. We presented the novel differentially methylated signals associated with articular cartilage degeneration in KBD and OA in terms of DMRs, DMGs, and functional pathways. Previous studies mainly focused on differentially methylated site analysis, while this work was preliminarily performing DMR analysis.

Age is usually strongly correlated with changes in DNA methylation. We found that 66% significant KBD-associated DMSs were strongly correlated with age and 100% significant OA-associated DMSs were strongly correlated with age. Further, we found that there were 4,325 DMSs for KBD that were correlated with disease status, while only one DMS for OA was correlated with disease status after removing age effects. These results confirmed that the DNA methylation pattern in OA samples was altered by age but not in KBD.

With controlling the age and sex as covariates, we identified 2.25-fold more KBD-associated DMRs and 1.5-fold more KBD-associated DMGs than the OA-associated DMRs and DMGs. Such results confirmed that epigenetic factors may contribute to KBD more than OA, while age contributed to OA more than KBD. Compared with the previous study, which identified 1,212, and 656 DMSs for KBD vs. Controls, and OA vs. Controls, respectively (Wang et al., 2017), we found that only 212 out of 1,212 DMSs (17.5%) were harbored in the DMRs we identified, and the 103 out of 656 DMSs (15.7%) were harbored in the DMRs we identified. This discrepancy was presumably due to the previous study that did not account for the confounding factors, such as age and sex. In addition, 74 out of 110 DMGs (67.2%) annotated for KBD vs. Controls and 31 out of 45 DMGs (68.9%) annotated for OA vs. Controls were also identified in previous study. Overall, more than 30% DMGs we identified were novel DMGs. Further examination of those non-overlapped DMGs would yield additional insights.

In this paper, we found several potential DMGs to be crucial for KBD. For example, the methylation region overlapped with gene HDAC4 was significantly differentially methylated between KBD and healthy subjects. HDAC4 encoded HDAC4 protein, which encode histone deacetylase 4 protein. Previous studies have reported that HDAC4 not only downregulates a series of cartilage-related gene expressions in human chondrocytes such as Runx2, MMP-1, MMP-3, MMP-13, ADAMTS4, and ADAMTS5 but also partially blocks the catabolic events in chondrocytes stimulated by 1L-1β (Cao et al., 2014). 1L-1β induced cartilage degradation by inducing the expression of MMP-1, MMP-13 (Lee et al., 2012), and aggrecanase-1/-2 (Sylvester et al., 2012) in chondrocytes and synovial fibroblasts. Further, HADC4 is a direct target of microRNA miR-365, which could repress chondrocyte differentiation in growth plate chondrocytes by targeting the 3′UTR of HDAC4 mRNA in osteoblasts (Yang et al., 2016; Xu et al., 2017). DMR analysis and previous studies suggested that the role of HDAC4 in KBD was crucial as it might serve as a potential drug target for clinical therapy.

Another example is SMOC2 that encodes SMOC2 protein, which belongs to the SPARC family (selected protein acidic and rich in cysteine/osteonectin/BM-40). SPARC is then matricellular proteins highly expressed in bone matrix that can modulate both matrix metalloprotease (MMP) and growth factor activates (Alford and Hankenson, 2006). Therefore, SMOC2 plays an important role in bone mineralization, cell–matrix interactions, collagen binding, and bone remodeling (Brekken and Sage, 2000). The SPARC family consists of several proteins, including SMOCs (secreted modular calcium-binding proteins), which are extracellular glycoproteins of the SPARC-related modular calcium-binding protein family (Gao et al., 2019). SMOC2 was observed that it negatively affects chondrogenic differentiation of precursor cells (Peeters et al., 2017). In addition, a positive correlation was observed between SMOC2 and 1L-1β. However, it may still lack a full understanding of the molecular interactions between SMOC2 and KBD, as well as how it affected KBD progression. Therefore, further understanding of SMOC2 is needed to better understand how SMOC2 affects KBD clinical outcome.

There are also several limitations and potential extensions for this study. First, articular cartilage specimens were collected from five KBD patients, five OA patients, and five healthy subjects for DMR analysis due to the limited number of KBD patients who need surgical treatment. Second, in addition to age and sex, DNA methylation patterns are also influenced by other potential factors such as body mass index (BMI). However, due to the small sample size, we did not consider other confounding factors in our model, which would cause a substantial loss of degree of freedom. Therefore, some other factors would remain probable confounding factors for the DMR analysis. Third, we only performed differentially methylated region analysis for KBD and OA, while integrative analysis of DNA methylation data and gene expression data from the same individuals has become an increasingly common approach in DNA methylation analysis (Hotta et al., 2018). Because the existing gene expression data (Wang et al., 2009; Duan et al., 2010) and DNA methylation data were not matched with the same individuals, we did not attempt to extend our analysis to gene expression data. Nevertheless, measuring the DNA methylation and gene expression data from the same individuals and performing integrating analysis is an important future direction. Finally, our analysis did not account for cellular heterogeneity present in microarray data. Indeed, articular cartilage is a physiologically non-self-renewing avascular tissue, consisting of only chondrocytes (Jiang and Tuan, 2015). Chondrocytes could be further classified to several articular cartilage chondrocyte subtypes (St-Jacques et al., 1999; Saito et al., 2010; Prein et al., 2016). Therefore, the cellular composition of cartilage samples is heterogeneous and varying. Failed accounting for the cellular heterogeneity would affect the accuracy of the methylation–disease association estimation.

In conclusion, we have presented a comprehensive DMR analysis for KBD vs. health controls, OA vs. health controls, and identified a group of DMRs and DMGs for KBD and OA. To the best of our knowledge, this is also the first DMR analysis for KBD and OA, and provides new insights about the similarities and differences between KBD and OA at the methylation level. We hope that our results could provide novel clues for better understanding the pathogenesis and diagnosis for KBD cartilage disease.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

This study was approved by the Institutional Review Board (IRB) of Xi'an Jiaotong University. All subjects gave their informed written consent by signing a document that had been carefully reviewed by the IRB of Xi'an Jiaotong University.

Author Contributions

SS conceived the idea, provided funding support, and designed the experiments. YF performed the real data analyses. FZ and YW provided the real data set. YF, DG, YW, JZ, FZ, LW, and XG revised the manuscript. SS and YF wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (Grant No. 61902319) and the Natural Science Foundation of Shaanxi Province (Grant No. 2019JQ127).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.636291/full#supplementary-material

References

Abou Shady, E. A., Farrag, H., el-Damarawy, N., Mohammed, F., Kamel, A., and Massoud, A. (1991). In vitro effects of trace elements on blood clotting and platelet function. B-Zinc and magnesium. J. Egypt. Public Health Assoc. 66:49.

PubMed Abstract | Google Scholar

Alford, A. I., and Hankenson, K. D. (2006). Matricellular proteins: extracellular modulators of bone development, remodeling, and regeneration. Bone 38, 749–757. doi: 10.1016/j.bone.2005.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, A. S., and Loeser, R. F. (2010). Why is osteoarthritis an age-related disease? Best practice and research. Clin. Rheumatol. 24, 15–26. doi: 10.1016/j.berh.2009.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Aryee, M. J., Jaffe, A. E., Corrada-Bravo, H., Ladd-Acosta, C., Feinberg, A. P., Hansen, K. D., et al. (2014). Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363–1369. doi: 10.1093/bioinformatics/btu049

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc Ser B 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Brekken, R. A., and Sage, E. H. (2000). SPARC, a matricellular protein: at the crossroads of cell–matrix. Matrix Biol. 19, 569–580. doi: 10.1016/S0945-053X(00)00105-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Butcher, L. M., and Beck, S. (2015). Probe Lasso: a novel method to rope in differentially methylated regions with 450K DNA methylation data. Methods 72, 21–28. doi: 10.1016/j.ymeth.2014.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, K., Wei, L., Zhang, Z., Guo, L., Zhang, C., Li, Y., et al. (2014). Decreased histone deacetylase 4 is associated with human osteoarthritis cartilage degeneration by releasing histone deacetylase 4 inhibition of runt-related transcription factor-2 and increasing osteoarthritis-related genes: a novel mechanism of human osteoarthritis cartilage degeneration. Arthritis Res. Ther. 16, 1–13. doi: 10.1186/s13075-014-0491-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Wei, X., Wang, S., Jiao, Q., Zhang, Y., Du, G., et al. (2016). Compression regulates gene expression of chondrocytes through HDAC4 nuclear relocation via PP2A-dependent HDAC4 dephosphorylation. Biochim. Biophy. Acta 1863, 1633–1642. doi: 10.1016/j.bbamcr.2016.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Tellez, G., Escobar, J., and Vazquez-Anon, M. (2017). Impact of trace minerals on wound healing of footpad dermatitis in broilers. Sci. Rep. 7, 1–9. doi: 10.1038/s41598-017-02026-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Cobb, J., Dierich, A., Huss-Garcia, Y., and Duboule, D. (2006). A mouse model for human short-stature syndromes identifies Shox2 as an upstream regulator of Runx2 during long-bone development. Proc. Natl. Acad. Sci. U.S.A. 103, 4511–4515. doi: 10.1073/pnas.0510544103

PubMed Abstract | CrossRef Full Text | Google Scholar

Commission, N. H. (2016). China Health and Family Planning Statistical Yearbook 2016. Beijing: Peking Union Medical College.

Compere, S. J., and Palmiter, R. D. (1981). DNA methylation controls the inducibility of the mouse metallothionein-I gene in lymphoid cells. Cell 25, 233–240. doi: 10.1016/0092-8674(81)90248-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhoot, G. K., Gustafsson, M. K., Ai, X., Sun, W., Standiford, D. M., and Emerson, C. P. (2001). Regulation of Wnt signaling and embryo patterning by an extracellular sulfatase. Science 293, 1663–1666. doi: 10.1126/science.293.5535.1663

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, P., Zhang, X., Huang, C.-C., Jafari, N., Kibbe, W. A., Hou, L., et al. (2010). Comparison of beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics 11:587. doi: 10.1186/1471-2105-11-587

PubMed Abstract | CrossRef Full Text

Duan, C., Guo, X., Zhang, X. D., Yu, H. J., Yan, H., Gao, Y., et al. (2010). Comparative analysis of gene expression profiles between primary knee osteoarthritis and an osteoarthritis endemic to Northwestern China, Kashin-Beck disease. Arthritis Rheum. 62, 771–780. doi: 10.1002/art.27282

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Q., Mok, H.-P., and Zhuang, J. (2019). Secreted modular calcium-binding proteins in pathophysiological processes and embryonic development. Chin. Med. J. 132:2476. doi: 10.1097/CM9.0000000000000472

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldring, M. B. (2012). Chondrogenesis, chondrocyte differentiation, and articular cartilage metabolism in health and osteoarthritis. Ther. Adv. Musculoskelet. Dis. 4, 269–285. doi: 10.1177/1759720X12448454

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Ma, W.-J., Zhang, F., Ren, F.-L., Qu, C.-J., and Lammi, M. (2014). Recent advances in the research of an endemic osteochondropathy in China: Kashin-Beck disease. Osteoarthr. Cartilage 22, 1774–1783. doi: 10.1016/j.joca.2014.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Hotta, K., Kitamoto, A., Kitamoto, T., Ogawa, Y., Honda, Y., Kessoku, T., et al. (2018). Identification of differentially methylated region (DMR) networks associated with progression of nonalcoholic fatty liver disease. Sci. Rep. 8, 1–11. doi: 10.1038/s41598-018-31886-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaenisch, R., and Bird, A. (2003). Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat. Genet. 33, 245–254. doi: 10.1038/ng1089

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., and Tuan, R. S. (2015). Origin and function of cartilage stem/progenitor cells in osteoarthritis. Nat. Rev. Rheumatol. 11:206. doi: 10.1038/nrrheum.2014.200

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, P. A. (2012). Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 13, 484–492. doi: 10.1038/nrg3230

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, M., Schulze, J., Hansen, U., Ashwodt, T., Keene, D. R., Brunken, W. J., et al. (2004). A novel marker of tissue junctions, collagen XXII. J. Biol. Chem. 279, 22514–22521. doi: 10.1074/jbc.M400536200

PubMed Abstract | CrossRef Full Text | Google Scholar

Lacey, M., Baribault, C., Ehrlich, K. C., and Ehrlich, M. (2019). Atherosclerosis-associated differentially methylated regions can reflect the disease phenotype and are often at enhancers. Atherosclerosis 280, 183–191. doi: 10.1016/j.atherosclerosis.2018.11.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D.-S., Roh, S. Y., Choi, H., and Park, J.-C. (2020). NFI-C is required for epiphyseal chondrocyte proliferation during postnatal cartilage development. Mol. Cells 43:739. doi: 10.14348/molcells.2020.2272

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y.-A., Choi, H. M., Lee, S.-H., Hong, S.-J., Yang, H.-I., Yoo, M. C., et al. (2012). Hypoxia differentially affects IL-1β-stimulated MMP-1 and MMP-13 expression of fibroblast-like synoviocytes in an HIF-1α-dependent manner. Rheumatology 51, 443–450. doi: 10.1093/rheumatology/ker327

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, R., and Barrett-Jolley, R. (2015). Changes in membrane receptors and ion channels as potential biomarkers for osteoarthritis. Front. Physiol. 6:357. doi: 10.3389/fphys.2015.00357

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Zhang, X., Cao, Y., Hu, S., Peng, F., Zhou, J., et al. (2017). Association between EN1 rs4144782 and susceptibility of knee osteoarthritis: a case-control study. Oncotarget 8:36650. doi: 10.18632/oncotarget.16842

PubMed Abstract | CrossRef Full Text | Google Scholar

Loeser, R. F. (2017). The role of aging in the development of osteoarthritis. Trans. Am. Clin. Climatol. Assoc. 128:44.

PubMed Abstract | Google Scholar

Lü, A., Guo, X., Aisha, M., Shi, X., Zhang, Y. Z., and Zhang, Y. (2011). Kashin–Beck disease and Sayiwak disease in China: prevalence and a comparison of the clinical manifestations, familial aggregation, and heritability. Bone 48, 347–353. doi: 10.1016/j.bone.2010.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Maor, G. L., Yearim, A., and Ast, G. (2015). The alternative role of DNA methylation in splicing regulation. Trends Genet. 31, 274–280. doi: 10.1016/j.tig.2015.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda-Duarte, A. (2018). DNA methylation in osteoarthritis: current status and therapeutic implications. Open Rheumatol. J. 12:37. doi: 10.2174/1874312901812010037

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, J. C., Przkora, R., and Cruz-Almeida, Y. (2018). Knee osteoarthritis: pathophysiology and current treatment modalities. J. Pain Res. 11:2189. doi: 10.2147/JPR.S154002

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno-Reyes, R., Suetens, C., Mathieu, F., Begaux, F., Zhu, D., Rivera, M. T., et al. (1998). Kashin–Beck osteoarthropathy in rural Tibet in relation to selenium and iodine status. New Engl. J. Med. 339, 1112–1120. doi: 10.1056/NEJM199810153391604

PubMed Abstract | CrossRef Full Text | Google Scholar

Otsuki, S., Taniguchi, N., Grogan, S. P., D'Lima, D., Kinoshita, M., and Lotz, M. (2008). Expression of novel extracellular sulfatases Sulf-1 and Sulf-2 in normal and osteoarthritic articular cartilage. Arthritis Res. Ther. 10:R61. doi: 10.1186/ar2432

PubMed Abstract | CrossRef Full Text | Google Scholar

Peeters, T., Lories, R., and Cailotto, F. (2017). SMOC2, a secreted calcium-binding protein affects chondrogenic differentiation of precursor cells and homeostasis of adult cartilage. Osteoarthr. Cartilage 25:S154. doi: 10.1016/j.joca.2017.02.256

CrossRef Full Text | Google Scholar

Peng, A., Yang, C., Rui, H., and Li, H. (1992). Study on the pathogenic factors of Kashin-Beck disease. J. Toxicol. Environ. Health 35, 79–90. doi: 10.1080/15287399209531597

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickering, G., Simpson, J., Kiss-Toth, E., and Wilkinson, M. (2018). KIF26B is necessary for osteogenic transdifferentiation and mineralisation in an in vitro model of heterotopic ossification. Osteoarthr. Cartilage, 26:S33. doi: 10.1016/j.joca.2018.02.082

CrossRef Full Text | Google Scholar

Piva, S. R., Susko, A. M., Khoja, S. S., Josbeno, D. A., Fitzgerald, G. K., and Toledo, F. G. (2015). Links between osteoarthritis and diabetes: implications for management from a physical activity perspective. Clin. Geriatr. Med. 31, 67–87. doi: 10.1016/j.cger.2014.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Prein, C., Warmbold, N., Farkas, Z., Schieker, M., Aszodi, A., and Clausen-Schaumann, H. (2016). Structural and mechanical properties of the proliferative zone of the developing murine growth plate cartilage assessed by atomic force microscopy. Matrix Biol. 50, 1–15. doi: 10.1016/j.matbio.2015.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinonez, S. C., and Innis, J. W. (2014). Human HOX gene disorders. Mol. Genet. Metab. 111, 4–15. doi: 10.1016/j.ymgme.2013.10.012

CrossRef Full Text | Google Scholar

Ren, F., Guo, X., Zhang, R., Wang, S. J., Zuo, H., Zhang, Z., et al. (2007). Effects of selenium and iodine deficiency on bone, cartilage growth plate and chondrocyte differentiation in two generations of rats. Osteoarthr. Cartilage, 15, 1171–1177. doi: 10.1016/j.joca.2007.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynard, L. N., Ratnayake, M., Santibanez-Koref, M., and Loughlin, J. (2016). Functional characterization of the osteoarthritis susceptibility mapping to CHST11—a bioinformatics and molecular study. PLoS ONE 11:e0159024. doi: 10.1371/journal.pone.0159024

PubMed Abstract | CrossRef Full Text | Google Scholar

Rice, S. J., Beier, F., Young, D. A., and Loughlin, J. (2020). Interplay between genetics and epigenetics in osteoarthritis. Nat. Rev. Rheumatol. 16, 268–281. doi: 10.1038/s41584-020-0407-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, T., Fukai, A., Mabuchi, A., Ikeda, T., Yano, F., Ohba, S., et al. (2010). Transcriptional regulation of endochondral ossification by HIF-2α during skeletal growth and osteoarthritis development. Nat. Med. 16, 678–686. doi: 10.1038/nm.2146

PubMed Abstract | CrossRef Full Text | Google Scholar

Sauerland, K., Plaas, A. H., Raiss, R. X., and Steinmeyer, J. (2003). The sulfation pattern of chondroitin sulfate from articular cartilage explants in response to mechanical loading. Biochim. Biophys. Acta 1638, 241–248. doi: 10.1016/S0925-4439(03)00089-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Saxon, L., Finch, C., and Bass, S. (1999). Sports participation, sports injuries and osteoarthritis. Sports Med. 28, 123–135. doi: 10.2165/00007256-199928020-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Schepman, K., Engelbert, R. H., Visser, M. M., Yu, C., and De Vos, R. (2011). Kashin Beck disease: more than just osteoarthrosis. Int. Orthop. 35, 767–776. doi: 10.1007/s00264-010-1043-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, P. N. P., Yadav, U. S., Azad, K., Goswami, P., Kinare, V., and Bandyopadhyay, A. (2018). NFIA and GATA3 are crucial regulators of embryonic articular cartilage differentiation. Development 145:dev.156554. doi: 10.1242/dev.156554

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J. Y., Pineault, K. M., Dones, J. M., Raines, R. T., and Wellik, D. M. (2020). Hox genes maintain critical roles in the adult skeleton. Proc. Natl. Acad. Sci. U.S.A. 117, 7296–7304. doi: 10.1073/pnas.1920860117

PubMed Abstract | CrossRef Full Text | Google Scholar

St-Jacques, B., Hammerschmidt, M., and McMahon, A. P. (1999). Indian hedgehog signaling regulates proliferation and differentiation of chondrocytes and is essential for bone formation. Genes Dev. 13, 2072–2086. doi: 10.1101/gad.13.16.2072

PubMed Abstract | CrossRef Full Text | Google Scholar

Stone, R. (2009). A medical mystery in middle China. Am. Assoc. Adv.Sci. 324, 1378–1381. doi: 10.1126/science.324_1378

CrossRef Full Text | Google Scholar

Sylvester, J., El Mabrouk, M., Ahmad, R., Chaudry, A., and Zafarullah, M. (2012). Interleukin-1 induction of aggrecanase gene expression in human articular chondrocytes is mediated by mitogen-activated protein kinases. Cell. Physiol. Biochem. 30, 563–574. doi: 10.1159/000341438

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, Y., Morris, T. J., Webster, A. P., Yang, Z., Beck, S., Feber, A., et al. (2017). ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics 33, 3982–3984. doi: 10.1093/bioinformatics/btx513

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Dungen, M. W., Murk, A. J., Kok, D. E., and Steegenga, W. T. (2017). Persistent organic pollutants alter DNA methylation during human adipocyte differentiation. Toxicol. In Vitro 40, 79–87. doi: 10.1016/j.tiv.2016.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanderman, K., Tremblay, M., Zhu, W., Shimojo, M., Mienaltowski, M., Coleman, S., et al. (2011). Brother of CDO (BOC) expression in equine articular cartilage. Osteoarthr. Cartilage 19, 435–438. doi: 10.1016/j.joca.2011.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ventham, N., Kennedy, N. A., Adams, A., Kalla, R., Heath, S., O'leary, K., et al. (2016). Integrative epigenome-wide analysis demonstrates that DNA methylation may mediate genetic risk in inflammatory bowel disease. Nat. Commun. 7, 1–14. doi: 10.1038/ncomms13507

PubMed Abstract | CrossRef Full Text | Google Scholar

Visman, I., Van den Oever, I., de Rooij, M., Roorda, L., Lems, W., Nurmohamed, M., et al. (2019). Cardiovascular disease is associated with activity limitations in osteoarthritis patients. Int. J. Clin. Rheumatol. 14:99.

Google Scholar

Wang, S., Guo, X., Wang, W., and Wang, S. (2012). Genome-wide study identifies the regulatory gene networks and signaling pathways from chondrocyte and peripheral blood monocyte of K ashin–B eck disease. Genes Cells 17, 619–632. doi: 10.1111/j.1365-2443.2012.01620.x

CrossRef Full Text | Google Scholar

Wang, W., Yu, Y., Hao, J., Wen, Y., Han, J., Hou, W., et al. (2017). Genome-wide DNA methylation profiling of articular cartilage reveals significant epigenetic alterations in Kashin-Beck disease and osteoarthritis. Osteoarthr. Cartilage 25, 2127–2133. doi: 10.1016/j.joca.2017.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W.-Z., Guo, X., Duan, C., Ma, W. J., Zhang, Y., Xu, P., et al. (2009). Comparative analysis of gene expression profiles between the normal human cartilage and the one with endemic osteoarthritis. Osteoarthr. Cartilage 17, 83–90. doi: 10.1016/j.joca.2008.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Fan, X., Xing, L., and Tian, F. (2019). Wnt signaling: a promising target for osteoarthritis therapy. Cell Commun. Signal. 17, 1–14. doi: 10.1186/s12964-019-0411-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., Gao, Y., Hu, N., Wu, L., and Chen, Q. (2017). miR-365 ameliorates dexamethasone-induced suppression of osteogenesis in MC3T3-E1 cells by targeting HDAC4. Int. J. Mol. Sci. 18:977. doi: 10.3390/ijms18050977

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Q., Forno, E., Celedón, J. C., and Chen, W. (2020). A region-based method for causal mediation analysis of DNA methylation data. bioRxiv, 2020.2011.2003.366989. Available online at: https://www.biorxiv.org/content/biorxiv/early/2020/11/04/2020.11.03.366989.full.pdf (accessed November 4, 2020).

Google Scholar

Yang, L., Zhang, J., Zhao, G., Wu, C., Ning, Y., Wang, X., et al. (2017). Gene expression profiles and molecular mechanism of cultured human chondrocytes' exposure to T-2 toxin and deoxynivalenol. Toxicon 140, 38–44. doi: 10.1016/j.toxicon.2017.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Guan, Y., Tian, S., Wang, Y., Sun, K., and Chen, Q. (2016). Mechanical and IL-1β responsive miR-365 contributes to osteoarthritis development by targeting histone deacetylase 4. Int. J. Mol. Sci. 17:436. doi: 10.3390/ijms17040436

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Han, H., De Carvalho, D. D., Lay, F. D., Jones, P. A., and Liang, G. (2014). Gene body methylation can alter gene expression and is a therapeutic target in cancer. Cancer Cell 26, 577–590. doi: 10.1016/j.ccr.2014.07.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, L., Liu, H., Yan, M., Yang, J., Long, F., Muneoka, K., et al. (2007). Shox2 is required for chondrocyte proliferation and maturation in proximal limb skeleton. Dev. Biol. 306, 549–559. doi: 10.1016/j.ydbio.2007.03.518

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwolanek, D., Veit, G., Eble, J. A., Gullberg, D., Ruggiero, F., Heino, J., et al. (2014). Collagen XXII binds to collagen-binding integrins via the novel motifs GLQGER and GFKGER. Biochem. J. 459, 217–227. doi: 10.1042/BJ20130642

PubMed Abstract | CrossRef Full Text | Google Scholar