Abstract
Aims
To assess the alterations in cell-specific DNA methylation associated with chondroitin sulphate response using peripheral blood collected from Kashin-Beck disease (KBD) patients before initiation of chondroitin sulphate treatment.
Methods
Peripheral blood samples were collected from KBD patients at baseline of chondroitin sulphate treatment. Methylation profiles were generated using reduced representation bisulphite sequencing (RRBS) from peripheral blood. Differentially methylated regions (DMRs) were identified using MethylKit, while DMR-related genes were defined as those annotated to the gene body or 2.2-kilobase upstream regions of DMRs. Selected DMR-related genes were further validated by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) to assess expression levels. Tensor composition analysis was performed to identify cell-specific differential DNA methylation from bulk tissue.
Results
This study revealed 21,060 hypermethylated and 44,472 hypomethylated DMRs, and 13,194 hypermethylated and 22,448 hypomethylated CpG islands for differential global methylation for chondroitin sulphate treatment response. A total of 12,666 DMR-related genes containing DMRs were identified in their promoter regions, such as CHL1 (false discovery rate (FDR) = 2.11 × 10-11), RIC8A (FDR = 7.05 × 10-4), and SOX12 (FDR = 1.43 × 10-3). Additionally, RIC8A and CHL1 were hypermethylated in responders, while SOX12 was hypomethylated in responders, all showing decreased gene expression. The patterns of cell-specific differential global methylation associated with chondroitin sulphate response were observed. Specifically, we found that DMRs located in TESPA1 and ATP11A exhibited differential DNA methylation between responders and non-responders in granulocytes, monocytes, and B cells.
Conclusion
Our study identified cell-specific changes in DNA methylation associated with chondroitin sulphate response in KBD patients.
Cite this article: Bone Joint Res 2024;13(5):237–246.
Article focus
-
To investigate the alterations in cell-specific DNA methylation associated with chondroitin sulphate response for endemic arthritis.
Key messages
-
Cell-specific DNA methylation alterations are associated with chondroitin sulphate response in Kashin-Beck disease (KBD) patients.
-
The study revealed hypermethylated and hypomethylated differentially methylated regions (DMRs) and CpG islands for chondroitin sulphate treatment response.
-
Specific genes (e.g. CHL1 and RIC8A) showed differential methylation and gene expression in chondroitin sulphate responders.
Strengths and limitations
-
This is the first ever investigation to report cell-specific differential DNA methylation pattern for chondroitin sulphate response in endemic arthritis.
-
There were only ten samples received, which reduced representation bisulphite sequencing, thus restricting our ability to detect smaller alterations in DNA methylation at genome-wide significance.
Introduction
Kashin-Beck disease (KBD) is a chronic osteochondral disease caused by environmental risk factors, which leads to irreversible pathological and clinical signs.1 The primary damage targets of KBD are epiphyseal cartilage and articular cartilage, resulting in endochondral ossification disorders and secondary osteoarthropathy.2 Patients with advanced stages of the disease may experience joint deformity, short stature, permanent disability, and reduced quality of life.3 Up to now, treatments for KBD have been limited to palliative measures such as steroids,4 selenium supplementation,5 physiotherapy,6 and various joint debridement techniques,7 along with medications for symptomatic relief including chondroitin sulphate,8 non-steroidal anti-inflammatory drugs,9 sodium selenite,9 and intra-articular hyaluronic acid.10 It should be noted that there are substantial individual differences in the therapeutic response of these drugs among KBD patients.
Chondroitin sulphate is a natural sticky polysaccharide derived from animal cartilage, which possesses anti-inflammation effects11 and enhances immunity.12,13 Cartilage injury in adult KBD patients could be related to the modification of chondroitin sulphate structure.14 Luo et al15 indicated that alterations of enzymes involved in cartilage chondroitin sulphate metabolism on proteoglycan play an essential role in KBD pathogenesis among adolescent children. Chondroitin sulphate is often used to alleviate pain and promote cartilage regeneration, thereby improving joint functions.16 A previous study showed that chondroitin sulphate increases type II collagen and proteoglycan synthesis in human articular chondrocytes.17 Additionally, it can reduce the production of some proinflammatory mediators and proteases, improving the anabolic/catabolic balance of the extracellular cartilage matrix.17 Clinical trials have demonstrated that chondroitin sulphate is effective in reducing Western Ontario and McMaster Universities osteoarthritis index (WOMAC)18 total score and stiffness score over a six-month duration in KBD patients when compared with placebo controls.8 A cluster-randomized study also suggested that chondroitin sulphate may play a protective role in preserving articular cartilage.19 These findings suggest that maintaining adequate levels of chondroitin sulphate may help to prevent or reduce the severity of joint damage associated with KBD.
DNA methylation is catalyzed by the methyltransferase enzymes that affect gene regulation without altering the DNA sequence.20,21 Approximately 80% of all CpGs across the genome are methylated, making DNA methylation a highly stable biomarker compared with gene expression.22 Aberrant methylation has been reported in several diseases, indicating its potential as a diagnostic and prognostic tool.23 For example, Nair et al24 indicated that DNA methylation plays a mediating role between disease genetic variants and patient prognosis in rheumatic arthritis. Furthermore, the predictive capability of DNA methylation in response to therapy has been established, underscoring its significance in personalized medicine approaches.25 Specifically, a study investigating DNA methylation patterns in relation to etanercept responses identified two methylation sites that served as predictive biomarkers in rheumatic arthritis.26 So far, several studies have constructed DNA methylation profiles of KBD, osteoarthritis, and normal control articular cartilage, and have identified chondroitin sulphate-related epigenetic differences.27,28 Peripheral blood is composed of diverse cellular populations, each exhibiting unique methylation signatures.29 The epigenome-wide association study of peripheral blood samples at bulk level may obscure the influence of cell-specific differential DNA methylation. However, no investigations have explored the effect of cell-specific DNA methylation patterns on the response of chondroitin sulphate in KBD patients.
In this study, we conducted reduced representation bisulphite sequencing (RRBS) on peripheral blood samples obtained from KBD patients before initiation of chondroitin sulphate treatment. Paired differentially methylated region (DMR) analyses were performed to demonstrate the distinct epigenetic alternations between chondroitin sulphate responders and non-responders. The DMRs were then annotated with their associated differentially methylated genes (DMGs), followed by the enrichment analysis. Additionally, tensor composition analysis was performed to assess the characteristics of cell-specific DNA methylation associated with chondroitin sulphate response.
Methods
Patients and samples
A total of 120 KBD patients were enrolled in this study, and each subject underwent a thorough clinical and radiological examination. Diagnosis of KBD was based on the clinical diagnosis criteria of China (WS/T 207-2010). The KBD patients received treatment with chondroitin sulfate (Hengxin Pharmaceutical, China; national drug approval number H32025748) at a dosage of three tablets thrice daily for a duration of one month. Clinical data and joint dysfunction index scores were measured at baseline before treatment initiation and again at one month after treatment initiation. The information included self-reported ethnicities, lifestyle characteristics, health statuses, and family and medical histories, as well as other relevant information to assess patient demographics and disease severity. Subjects with missing blood and questionnaire data, with a joint dysfunction index score of less than 6 before treatment, and diagnosed with genetic bone and cartilage diseases, rheumatoid arthritis, or other skeletal disorders were excluded from the study. Peripheral blood samples were collected at baseline before chondroitin sulphate treatment initiation.
Clinical characteristics of study participants
Detailed characteristics of the study participants are summarized in Table I. This study involved peripheral blood samples from ten KBD patients, which were measured by RRBS. The mean age of chondroitin sulphate responders was 65.0 years (standard deviation (SD) 2.76), while that of non-responders was 65.6 years (SD 6.50). All KBD participants self-reported as Chinese Han ethnicity.
Table I.
Characteristics of the study subjects used in reduced representation bisulphite sequencing.
| Characteristic | Chondroitin sulphate responder | Chondroitin sulphate non-responder |
|---|---|---|
| Number of participants | 5 | 5 |
| Male, n (%) | 3 (60) | 3 (60) |
| Mean age at baseline, yrs (SD) | 65.0 (2.76) | 65.6 (6.50) |
| Mean baseline SJDI score (SD) | 8.2 (0.40) | 7.8 (1.16) |
| Mean 3-mth follow-up SJDI score after CS treatment (SD) | 5.4 (0.80) | 8.8 (0.75) |
-
CS, chondroitin sulphate; SD, standard deviation; SJDI, sum of joint dysfunction index.
Definition of chondroitin sulphate response
The joint dysfunction index scores for KBD were measured using a document from the National and Family Planning Commission of China titled “Assessment for Therapeutic Efficacy on Kashin-Beck disease (WS/T 79-2011)” as the primary measure (http://www.nhc.gov.cn/ewebeditor/uploadfile/2014/10/20141022174419733.pdf). The joint dysfunction index consists of five items: Q1 quantifies the degree of arthralgia during periods of nocturnal rest; Q2 evaluates arthralgia experienced during ambulation; Q3 measures the severity of morning stiffness; Q4 measures the maximum distance ambulated; and Q5 assesses activities of lower limb. Each item is allocated a score ranging from 0 to 2, with a total score range of 0 to 10 (Supplementary Table i).30 These scores provide an objective assessment of disease severity in KBD patients. At the pre-treatment baseline visit and at the follow-up visit (one month after baseline), the sum of joint dysfunction index (SJDI) scores was measured. The response to chondroitin sulphate was determined based on the difference between pre-treatment SJDI scores and post-treatment SJDI scores: the value ≥ 2 indicated responder, while the value ≤ 0 indicated non-responder.
Library construction and sequencing
DNA was extracted from peripheral blood specimens obtained from five response and five non-response patients using QIAamp Fast DNA Tissue Kit (Qiagen, Germany) following the manufacturer’s instructions. The quantity of extracted DNA was determined by agarose gel electrophoresis and A260/280 ratios measured by a spectrophotometer. The fragmented DNA samples were treated with MspI (New England Biolabs (NEB), USA) and then subjected to bisulphite conversion for RRBS. The Accel-NGS Methyl-Seq DNA Library Kit (Swift Biosciences, Integrated DNA Technologies, USA) was used to attach adapters to the single-stranded DNA fragments. This involved a highly efficient Adaptase step that performed end repair, tailing of 3' ends, and ligation of the first truncated adapter complement to 3' ends simultaneously. To incorporate truncated adapter 1, a primer extension reaction was carried out during the extension step. The ligation step added the second truncated adapter to the bottom strand only. To enhance yield and facilitate the integration of full-length adapters, the indexing PCR step was implemented. Bead-based solid-phase reversible immobilization (SPRI) clean-ups were used to eliminate oligonucleotides and small fragments, as well as to modify enzymatic buffer composition. Ultimately, pair-end 2 × 150 base pair sequencing was conducted utilizing an Illumina HiSeq 4000 platform (Illumina, USA) housed in LC Sciences (USA).
Differentially methylated analysis
Initially, Cutadapt (Marcel Martin, Sweden) and in-house Perl scripts (Lianchuan, China) were employed to remove reads containing adapter contamination, low-quality bases, and undetermined bases.31 Subsequently, the sequence quality was assessed using FastQC (Babraham Bioinformatics, UK). Reads that passed the quality control were mapped to reference genome hg38 using Bismark.32 Following alignment, reads were de-duplicated using SAMtool.33 For each cytosine site (or guanine corresponding to a cytosine on the opposite strand) in the reference genome sequence, we determined the DNA methylation level by calculating the ratio of reads supporting C (methylated) to that of total reads (methylated and unmethylated) using MethPipe (Andrew D. Smith, USA).34 DMRs were calculated using default parameters in R package MethylKit (R Foundation for Statistical Computing, Austria), which included 1,000 bp slide windows with 500 bp overlap and a p-value < 0.05 for statistical significance.35
Total RNA extraction, reverse transcription, and quantitative reverse transcriptase PCR
The quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) was carried out as described previously.2 Briefly, total RNA samples from peripheral blood (ten responders and ten non-responders) were extracted using the TRIzol reagent (Invitrogen Thermo Fisher Scientific, USA). Next, 500 ng of the extracted RNA was reverse transcribed into complementary DNA according to the manufacturer’s protocol of Maxima H Minus cDNA Synthesis Master Mix with dsDNase (Thermo Fisher Scientific).36 Then, 25 μl of reaction volume was used for the qRT-PCR on CFX96 Real-Time PCR System (Bio-Rad Laboratories, USA) using SYBR Green qPCR Master Mix (Bimake.com, USA) to detect expression levels of DMR-related genes.37 All reactions were carried out at an annealing temperature of 60°C; cycling conditions were: 95°C-30 s (one cycle), 95°C-5 s, followed by 60°C-30 s (40 cycles). The primers used for cell adhesion molecule L1 like (CHL1), RIC8 guanine nucleotide exchange factor A (RIC8A), SRY-box transcription factor 12 (SOX12), acid phosphatase 1 (ACP1), and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were purchased from Takara Bio (Japan) (Supplementary Table ii). All primers span intron-exon boundaries, and used at a final concentration of 500 nM and validated using a standard curve of five serial dilutions so that all primer efficiencies were between 90% and 105%. The 2-ΔΔCt method was employed to determine relative messenger RNA (mRNA) levels and standardize them for each transcript.38 The GAPDH gene was used to normalize the data. All the qRT-PCR experiments were performed as triplicate replications.
Statistical analysis
GraphPad Prism 7.03 (GraphPad Software, USA) was used for statistical analysis. qRT-PCR data were presented as means with SDs, and differences in expression levels were assessed using the two-tailed independent-samples t-test. Statistical significance was defined as p < 0.05.
Functional enrichment analysis of DMR-related genes
DMR-related genes were defined as those annotated to the DMRs, including gene body or 2.2-kilobase upstream regions of gene body, which had additional biological significance. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used for annotating both DMR-related genes and background genes. To analyze the enrichment of GO terms and KEGG pathways in genes related to DMRs, we used the GOseq R package (Matthew Young, UK)39 and KOBAS software (Dechao Bu, China),40 respectively. The enriched pathways with a corrected p-value < 0.05 were considered statistically significant.
Global cell proportion estimation
The input methylation data were normalized based on the chromosome position of CpGs in Infinium HumanMethylation450 BeadChip array (Illumina, USA).41 Briefly, the chromosome position of DMR was mapped to CpG, and the methylation level values of different DMRs at each CpG were then averaged to obtain their corresponding methylation values. Furthermore, global cell type proportions for B cells, CD8+ and CD4+ T cells, monocytes, natural killer (NK) cells, and granulocytes were estimated separately within the methylation level matrix using centDHSbloodDMC.m () in EpiDISH R package.42
Cell type-specific analyses
We used tensor composition analysis to estimate differential DNA methylation between chondroitin sulphate treatment responses and cell type-specific DNA methylation. Tensor composition analysis is a method that estimates cell type-specific associations between DNA methylation and disease phenotypes using bulk DNA methylation data.43 The input included estimated global cell proportions and bulk DNA methylation data related to chondroitin sulphate treatment response. To estimate differential DNA methylation in tensor composition analysis, we implemented a two-step pipeline: a joint model and a marginal conditional model. The joint model is used to detect differential DNA methylation within any cell type at a CpG, while the marginal conditional model is employed to detect differential DNA methylation within a specific cell type adjusted for the other cell types at a CpG. The joint model can be analogized to an analysis of variance (ANOVA) test, and provides evidence for differential DNA methylation in at least one cell type. The p-values of joint model were adjusted for multiple tests using the Benjamini-Hochberg method.43 All CpGs with p-values < 0.05 were tested for cell-specific differential DNA methylation.
Results
Genome-wide DNA methylation profiles
The specific distribution for CpG, CHH, and CHG sites across the whole genome are shown in Supplementary Figure a. Mean methylation levels and patterns were almost identical between chondroitin sulphate responders and non-responders for CpG, CHH, and CHG contexts in promoter, exon, intron, and downstream regions (Supplementary Figures b to d). The promoter and downstream regions exhibited higher mean methylation levels than gene bodies such as internal exons and introns.
DMR results and DMR-related genes between responders and non-responders
To identify methylated regions and genes potentially involved in chondroitin sulphate response, we compared the DNA methylation profiles of chondroitin sulphate responder versus non-responder (Figure 1a). A total of 65,532 DMRs (false discovery rates (FDRs) < 0.05) were identified and annotated with genomic features; however, over 90% of these DMRs were found in promoter regions and intron regions. Compared with non-responder, 21,060 DMRs were hypermethylated, and the remaining 44,472 DMRs were hypomethylated (Figure 1b). In functional elements such as promoter, exon, intron, and intergenic regions (Table II), the number of hypomethylated regions was higher than that of hypermethylated areas. DMRs in promoters were mainly observed in low CpG-containing promoter (LCP) (50.81%), followed by high CpG-containing promoter (HCP) (25.49%) and intermediate CpG-containing promoter (ICP) (23.70%). Furthermore, more DMRs exhibited hypomethylation in LCP, HCP, and ICP regions. We identified a total of 12,666 DMR-related genes that contained DMRs within their promoter regions, such as CHL1 (FDR = 2.11 × 10-11), RIC8A (FDR = 7.05 × 10-4), SOX12 (FDR = 1.43 × 10-3), and ACP1 (FDR = 1.73 × 10-5).
Fig. 1
Overview of DNA methylation profile. a) A heatmap was used to assess the DNA methylation level of the top 50 methylated sites. Each methylated site is represented by a single row of coloured boxes, while each sample is represented by a single column. b) Volcano plots for differentially methylated regions (DMRs). The colour scheme uses red to denote high expression and blue to indicate low expression. c) Confirmation of the differential expression of DMR-related genes. Messenger RNA (mRNA) fold change for cell adhesion molecule L1 like (CHL1), RIC8 guanine nucleotide exchange factor A (RIC8A), SRY-box transcription factor 12 (SOX12), and acid phosphatase 1 (ACP1) in peripheral blood, measured by quantitative reverse transcriptase polymerase chain reaction with triplicate replication. Symbols represent each individual donor. Bars show the mean and standard deviation. *p < 0.05, ***p < 0.001. GAPDH, glyceraldehyde 3-phosphate dehydrogenase.
Table II.
Statistics of differentially methylated regions in promoters for chondroitin sulphate responder versus non-responder.
| Total promoters | Hypermethylated, n (n = 87,058) |
Hypomethylated, n (n = 117,025) |
|---|---|---|
| Proximal | ||
| HCP | 8,567 | 10,378 |
| ICP | 7,996 | 10,115 |
| LCP | 17,072 | 21,428 |
| Intermediate | ||
| HCP | 8,572 | 10,046 |
| ICP | 7,200 | 9,654 |
| LCP | 15,499 | 20,825 |
| Distal | ||
| HCP | 6,206 | 8,255 |
| ICP | 5,021 | 8,373 |
| LCP | 10,925 | 17,951 |
-
HCP, high CpG-containing promoter; ICP, intermediate CpG-containing promoter; LCP, low CpG-containing promoter.
To investigate the relationship between DMR-related genes and chondroitin sulphate response, we conducted an analysis of four DMR-related genes based on previous studies involving chondrogenesis and/or drug treatment. Research has indicated that DNA methylation on functional regions can regulate gene expression by modifying chromatin structure or transcription efficiency.44 Thus, we examined the transcriptional expressions of CHL1, RIC8A, SOX12, and ACP1, of which DMR-related genes were expressed differently between chondroitin sulphate responders and non-responders (Figure 1c). Specifically, the expression levels of CHL1, RIC8A, and SOX12 were significantly higher in non-responders than those in responders (p < 0.05, two-tailed independent-samples t-test).
Functional implication of DMR-related genes
In order to investigate the potential biological significance of differential methylation between responders and non-responders to chondroitin sulphate therapy, we conducted pathway enrichment analyses utilizing annotated DMR-related genes. GO enrichment analysis of DMG identified 571 significant GO terms for biological process, such as phosphorylation (FDR = 1.94 × 10-20), protein phosphorylation (FDR = 2.46 × 10-18), and regulation of transcription, DNA-templated (FDR = 3.84 × 10-18). For molecular function, 193 significant GO terms were identified, such as ATP binding (FDR = 5.15 × 10-30), nucleotide binding (FDR = 6.87 × 10-26), and metal ion binding (FDR = 1.23 × 10-20). For cellular component, 168 significant GO terms were identified, such as cytoplasm (FDR = 4.88 × 10-32), nucleus (FDR = 8.64 × 10-30), and cytosol (FDR = 1.84 × 10-27) (Figure 2a). Furthermore, our investigation of KEGG annotated pathways revealed the top five hits as glutamatergic synapse (FDR = 6.00 × 10-4), type II diabetes mellitus (FDR = 8.83 × 10-4), focal adhesion (FDR = 1.13 × 10-3), circadian entrainment (FDR = 1.50 × 10-3), and axon guidance (FDR = 1.51 × 10-3) (Figure 2b).
Fig. 2
Function analysis of differentially methylated region-related genes. a) Gene Ontology (GO) analysis of differentially methylated region (DMR)-related genes. b) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DMR-related genes. ATP, adenosine triphosphate.
Differentially methylated region within CpG island between responders and non-responders
Our analysis of methylation sites within CpG island (CGI) were collectively analyzed using EMBOSS (Rice P, UK).45 In total, 35,642 significant CGIs with a FDR < 0.05 were identified, including 13,194 hypermethylated CGIs and 22,448 hypomethylated CGIs. In particular, our investigation revealed notable associations between specific CGIs and various genes such as WDR43 (hypermethylated, FDR = 3.92 × 10-208) and CLIC4 (hypomethylated, FDR = 3.71 × 10-199).
Global methylation changes are cell-specific between responders and non-responders
Joint test revealed two regions that differed across all cell types, chr12:93963963-93967395 (p = 0.026) and chr13:113424898-113425105 (p = 0.006). Marginal conditional test identified that chr12:93963963-93967395 in monocytes and B cells was associated with differential DNA methylation in chondroitin sulphate responders compared with non-responders at baseline. This region is located within the transcription sites of KIAA0748 (TESPA1) on chromosome 12. Moreover, chr13:113424898-113425105 in granulocytes, monocytes, and B cells was associated with differential DNA methylation in chondroitin sulphate responders compared with non-responders at baseline, which was located within the transcription sites of ATP11A on chromosome 13 (Table III).
Table III.
Results of the tensor composition analysis of differentially methylated regions.
| chr: stratpos-endpos (hg38) | Gene annotation | Cell type | Δ methylation | Joint model p-value* | Marginal conditional model p-value |
|---|---|---|---|---|---|
| chr12:93963963-93967395 | KIAA0748 | Granulocyte | Hyper | 0.026 | 0.282 |
| Monocyte | 0.034 | ||||
| B cell | 0.034 | ||||
| chr13:113424898-113425105 | ATP11A | Granulocyte | Hypo | 0.006 | 0.005 |
| Monocyte | 0.006 | ||||
| B cell | 0.006 |
-
*
Joint model in tensor composition analysis tested for evidence of differential DNA methylation within any cell type at each differentially methylated region.
Discussion
In this study, we conducted the first genome-wide RRBS of peripheral blood to reveal the epigenetic characteristics between chondroitin sulphate responders and non-responders in KBD patients. We presented novel differentially methylated signals associated with chondroitin sulfate treatment response among KBD patients in terms of 756 DMRs, 48 CGIs, and multiple functional pathways. We also performed tensor composition analysis to clarify the cell type-specific DNA methylation of chondroitin sulphate treatment response in KBD. Our findings provide insight into the epigenetic mechanisms underlying chondroitin sulphate treatment response and the specific DNA methylation patterns among cell types in KBD patients.
Pre-treatment DNA methylation patterns may serve as predictive markers of treatment response. DNA methylation that occurs in the gene body regions can interfere with the ability of transcription factors to interact with the DNA, thereby preventing access by the transcription machinery and altering gene expression.46 The cytosine-guanine (C-G) dinucleotides (CpG) are commonly found in clusters in gene promoter regions, where they are known as CpG islands; these islands tend to be hypomethylated when gene transcription is active.23 Conversely, DNA hypermethylation is generally associated with a decrease in gene expression.23
Our analysis revealed multiple novel genes of differential DNA methylation associated with chondroitin sulphate that were not present in the Comparative Toxicogenomics Database,47 such as CHL1, RIC8A, SOX12, and ACP1. CHL1 is a bone morphogenetic protein antagonist,48 and could serve as an important biomarker for both diagnosis and treatment outcome of complex disease.49RIC8A is a highly conserved cytosolic protein that plays a regulatory role in asymmetric cell divisions.50 Knocking down RIC8A has been shown to result in anomalous radial migratory behaviour, characterized by reduced cell spreading and focal adhesion formation.51SOX12 plays a crucial role in establishing non-canonical wingless/integrated (WNT) signal crosstalk, which positively regulates SOX4 and SOX11 in growth plate chondrocytes.52ACP1 plays a pivotal role in the biosynthesis of cytosolic low molecular weight protein tyrosine phosphatase,53 which serves as a key regulator in regulation of growth factors and cellular adhesion. A previous study confirmed an association between the ACP1 gene and treatment-related osteonecrosis in children.54
Our GO enrichment analysis of DMR-related genes revealed functional categories mainly related to regulation of transcription, RNA polymerase, and protein binding. Additionally, these genes were enriched in five KEGG main classes potentially associated with environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems, such as focal adhesion, long-term potentiation, glutamatergic synapse, and insulin signalling pathway. A study by Han et al55 investigated differences in N-glycoproteins present in knee cartilage from individuals with normal control, OA, and KBD, and found that key N-glycoproteins were closely related to focal adhesion pathway for OA and KBD. However, there is limited evidence regarding the treatment response of chondroitin sulphate in KBD patients. Given that epigenetic factors largely promote the development of KBD,56 studying chondroitin sulphate treatment-related pathways could have considerable translational value and impact on developing therapeutic strategies for KBD.
Tensor composition analysis indicated that chondroitin sulphate treatment was associated with cell-specific DNA methylation changes between non-responders and responders. The observed distinctions at baseline suggest that gene expression influenced by DNA methylation levels located in chr12:93963963-93967395 and chr13:113424898-113425105 may differ between chondroitin sulphate response groups. Our findings contribute to the understanding of the underlying mechanisms of KBD, particularly in relation to the response to chondroitin sulphate treatment. By identifying cell-specific changes in DNA methylation associated with chondroitin sulphate response, we uncovered potential biomarkers that could be used to monitor treatment response. This could lead to improved diagnostic accuracy and more personalized treatment strategies for KBD patients. Cell-specific DNA methylation patterns will be invaluable for epigenome-wide association studies and KBD epigenome studies to help improve biological interpretation, for prioritizing candidates that require functional validation, and to help elucidate causal pathways to KBD.
Furthermore, our findings on the patterns of cell-specific differential global methylation associated with chondroitin sulphate response, including one located in ATP11A that was significantly associated with differential DNA methylation between responders and non-responders in granulocytes, monocytes, and B cells, provide novel insights into the cellular mechanisms of KBD. Given that a CGI within the ATP11A gene is susceptible to methylation, it is plausible that the downregulation of ATP11A expression in B cell progenitors is a consequence of DNA methylation.57 It has been shown that an increased expression level of ATP11A can confer protection to malignant lymphoblastic leukaemia cells against a variety of novel small molecule signal transduction inhibitors.58 Therefore, determining the gene expression level of ATP11A may have prognostic value,58 and potentially serve as a novel therapeutic target for KBD. Our findings could be helpful in the design of targeted therapies that aim to modulate these specific cellular responses, thereby enhancing the efficacy of chondroitin sulphate treatment in KBD patients.
To the best of our knowledge, this is the first investigation to report cell-specific differential DNA methylation for chondroitin sulphate response in KBD patients. However, our study also has several limitations. First, there were only ten samples receiving RRBS, which restricted our ability to detect smaller alterations in DNA methylation at genome-wide significance. Second, we used GAPDH as the internal reference gene in PCR, which may not be the most suitable internal reference gene. The stability of multiple reference genes in KBD samples should be compared in future studies. Third, estimates of differential DNA methylation derived from tensor composition analysis were not based on cell-specific DNA methylation measurements obtained from our study participants. The ethnic differences may influence DNA methylation patterns and subsequently affect cell type prediction accuracy. It is imperative to obtain DNA methylation data from sorted cells for more accurate results. Additionally, there exists a greater capacity to detect differential DNA methylation in more abundant cell types compared with less abundant ones.
In summary, we presented a comprehensive analysis of DMR in response to chondroitin sulphate in KBD patients, and estimated alterations in DNA methylation using deconvolute cell-specific DNA methylation methods. We identified a group of differentially methylated genes related to chondroitin sulphate response, and identified several cell-specific DMRs. Our findings suggest that there is a strong association between cell-specific differential DNA methylation and chondroitin sulphate response. Our study provides novel insights into understanding the characteristics of target cells affected by chondroitin sulfate therapy and identifying potential biomarkers for evaluating its efficacy in KBD patients.
References
1. Wang X , Ning Y , Zhang P , et al. Comparison of the major cell populations among osteoarthritis, Kashin-Beck disease and healthy chondrocytes by single-cell RNA-seq analysis . Cell Death Dis . 2021 ; 12 ( 6 ): 551 . Crossref PubMed Google Scholar
2. Cheng B , Liang C , Yang X , et al. Genetic association scan of 32 osteoarthritis susceptibility genes identified TP63 associated with an endemic osteoarthritis, Kashin-Beck disease . Bone . 2021 ; 150 : 115997 . Crossref PubMed Google Scholar
3. Ning Y , Hu M , Chen S , et al. Investigation of selenium nutritional status and dietary pattern among children in Kashin-Beck disease endemic areas in Shaanxi Province, China using duplicate portion sampling method . Environ Int . 2022 ; 164 : 107255 . Crossref PubMed Google Scholar
4. Liu W , Liu G , Pei F , et al. Kashin-Beck disease in Sichuan, China: report of a pilot open therapeutic trial . J Clin Rheumatol . 2012 ; 18 ( 1 ): 8 – 14 . Crossref PubMed Google Scholar
5. Xie D , Liao Y , Yue J , et al. Effects of five types of selenium supplementation for treatment of Kashin-Beck disease in children: a systematic review and network meta-analysis . BMJ Open . 2018 ; 8 ( 3 ): e017883 . Crossref PubMed Google Scholar
6. Mathieu F , Suetens C , Begaux F , De Maertelaer V , Hinsenkamp M . Effects of physical therapy on patients with Kashin-Beck disease in Tibet . Int Orthop . 2001 ; 25 ( 3 ): 191 – 193 . Crossref PubMed Google Scholar
7. Jin Z-K , Xu C-X , Dong X-H , et al. Long-term outcomes of arthroscopic debridement of the knee in adults with Kashin-Beck disease: an 18-year follow-up . J Int Med Res . 2021 ; 49 ( 10 ): 3000605211050781 . Crossref PubMed Google Scholar
8. Yue J , Yang M , Yi S , et al. Chondroitin sulfate and/or glucosamine hydrochloride for Kashin-Beck disease: a cluster-randomized, placebo-controlled study . Osteoarthritis Cartilage . 2012 ; 20 ( 7 ): 622 – 629 . Crossref PubMed Google Scholar
9. Jirong Y , Huiyun P , Zhongzhe Y , et al. Sodium selenite for treatment of Kashin-Beck disease in children: a systematic review of randomised controlled trials . Osteoarthritis Cartilage . 2012 ; 20 ( 7 ): 605 – 613 . Crossref PubMed Google Scholar
10. Yu FF , Xia CT , Fang H , Han J , Younus MI , Guo X . Evaluation of the therapeutic effect of treatment with intra-articular hyaluronic acid in knees for Kashin-Beck disease: a meta-analysis . Osteoarthritis Cartilage . 2014 ; 22 ( 6 ): 718 – 725 . Crossref PubMed Google Scholar
11. Melgar-Lesmes P , Garcia-Polite F , Del-Rey-Puech P , et al. Treatment with chondroitin sulfate to modulate inflammation and atherogenesis in obesity . Atherosclerosis . 2016 ; 245 : 82 – 87 . Crossref PubMed Google Scholar
12. Hatano S , Watanabe H . Regulation of Macrophage and Dendritic Cell Function by Chondroitin Sulfate in Innate to Antigen-Specific Adaptive Immunity . Front Immunol . 2020 ; 11 : 232 . Crossref PubMed Google Scholar
13. Qi SS , Shao ML , Sun Z , et al. Chondroitin Sulfate Alleviates Diabetic Osteoporosis and Repairs Bone Microstructure via Anti-Oxidation, Anti-Inflammation, and Regulating Bone Metabolism . Front Endocrinol (Lausanne) . 2021 ; 12 : 759843 . Crossref PubMed Google Scholar
14. Han J , Li D , Qu C , et al. Altered expression of chondroitin sulfate structure modifying sulfotransferases in the articular cartilage from adult osteoarthritis and Kashin-Beck disease . Osteoarthritis Cartilage . 2017 ; 25 ( 8 ): 1372 – 1375 . Crossref PubMed Google Scholar
15. Luo M , Chen J , Li S , et al. Changes in the metabolism of chondroitin sulfate glycosaminoglycans in articular cartilage from patients with Kashin-Beck disease . Osteoarthritis Cartilage . 2014 ; 22 ( 7 ): 986 – 995 . Crossref PubMed Google Scholar
16. Li L , Li Y , Feng D , et al. Preparation of Low Molecular Weight Chondroitin Sulfates, Screening of a High Anti-Complement Capacity of Low Molecular Weight Chondroitin Sulfate and Its Biological Activity Studies in Attenuating Osteoarthritis . Int J Mol Sci . 2016 ; 17 ( 10 ): 1685 . Crossref PubMed Google Scholar
17. Henrotin Y , Marty M , Mobasheri A . What is the current status of chondroitin sulfate and glucosamine for the treatment of knee osteoarthritis? Maturitas . 2014 ; 78 ( 3 ): 184 – 187 . Crossref PubMed Google Scholar
18. Bellamy N , Buchanan WW , Goldsmith CH , Campbell J , Stitt LW . Validation study of WOMAC: a health status instrument for measuring clinically important patient relevant outcomes to antirheumatic drug therapy in patients with osteoarthritis of the hip or knee . J Rheumatol . 1988 ; 15 ( 12 ): 1833 – 1840 . PubMed Google Scholar
19. Zhang Y , Dong W , Liu H , Cicuttini F , de Courten M , Yang J . Effects of chondroitin sulfate and glucosamine in adult patients with Kaschin-Beck disease . Clin Rheumatol . 2010 ; 29 ( 4 ): 357 – 362 . Crossref PubMed Google Scholar
20. Nam AS , Chaligne R , Landau DA . Integrating genetic and non-genetic determinants of cancer evolution by single-cell multi-omics . Nat Rev Genet . 2021 ; 22 ( 1 ): 3 – 18 . Crossref PubMed Google Scholar
21. He CP , Jiang XC , Chen C , et al. The function of lncRNAs in the pathogenesis of osteoarthritis . Bone Joint Res . 2021 ; 10 ( 2 ): 122 – 133 . Crossref PubMed Google Scholar
22. Bird A . DNA methylation patterns and epigenetic memory . Genes Dev . 2002 ; 16 ( 1 ): 6 – 21 . Crossref PubMed Google Scholar
23. van der Maarel SM . Epigenetic mechanisms in health and disease . Ann Rheum Dis . 2008 ; 67 Suppl 3 : iii97 – 100 . Crossref PubMed Google Scholar
24. Nair N , Wilson AG , Barton A . DNA methylation as a marker of response in rheumatoid arthritis . Pharmacogenomics . 2017 ; 18 ( 14 ): 1323 – 1332 . Crossref PubMed Google Scholar
25. Ouchi K , Takahashi S , Yamada Y , et al. DNA methylation status as a biomarker of anti-epidermal growth factor receptor treatment for metastatic colorectal cancer . Cancer Sci . 2015 ; 106 ( 12 ): 1722 – 1729 . Crossref PubMed Google Scholar
26. Plant D , Webster A , Nair N , et al. Differential Methylation as a Biomarker of Response to Etanercept in Patients With Rheumatoid Arthritis . Arthritis Rheumatol . 2016 ; 68 ( 6 ): 1353 – 1360 . Crossref PubMed Google Scholar
27. Wang W , Yu Y , Hao J , et al. Genome-wide DNA methylation profiling of articular cartilage reveals significant epigenetic alterations in Kashin-Beck disease and osteoarthritis . Osteoarthritis Cartilage . 2017 ; 25 ( 12 ): 2127 – 2133 . Crossref PubMed Google Scholar
28. Wen Y , Li P , Hao J , et al. Integrating genome-wide DNA methylation and mRNA expression profiles identified different molecular features between Kashin-Beck disease and primary osteoarthritis . Arthritis Res Ther . 2018 ; 20 ( 1 ): 41 . Crossref PubMed Google Scholar
29. Houseman EA , Accomando WP , Koestler DC , et al. DNA methylation arrays as surrogate measures of cell mixture distribution . BMC Bioinformatics . 2012 ; 13 : 86 . Crossref PubMed Google Scholar
30. Yu F-F , Xia C-T , Fang H , Wang D-M , Guo X . Reliability and validation of the joint dysfunction index as a new assessment instrument for therapeutic efficacy for Kashin-Beck disease . Clin Rheumatol . 2016 ; 35 ( 11 ): 2815 – 2821 . Crossref PubMed Google Scholar
31. Martin M . Cutadapt removes adapter sequences from high-throughput sequencing reads . EMBnet J . 2011 ; 17 ( 1 ): 10 . Crossref Google Scholar
32. Krueger F , Andrews SR . Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications . Bioinformatics . 2011 ; 27 ( 11 ): 1571 – 1572 . Crossref PubMed Google Scholar
33. Li H , Handsaker B , Wysoker A , et al. The Sequence Alignment/Map format and SAMtools . Bioinformatics . 2009 ; 25 ( 16 ): 2078 – 2079 . Crossref PubMed Google Scholar
34. Song Q , Decato B , Hong EE , et al. A reference methylome database and analysis pipeline to facilitate integrative and comparative epigenomics . PLoS One . 2013 ; 8 ( 12 ): e81148 . Crossref PubMed Google Scholar
35. Akalin A , Kormaksson M , Li S , et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles . Genome Biol . 2012 ; 13 ( 10 ): R87 . Crossref PubMed Google Scholar
36. Liu Z-M , Shen P-C , Lu C-C , Chou S-H , Tien Y-C . Suramin enhances chondrogenic properties by regulating the p67phox/PI3K/AKT/SOX9 signalling pathway . Bone Joint Res . 2022 ; 11 ( 10 ): 723 – 738 . Crossref PubMed Google Scholar
37. Du X , Jing Z , Wen G , et al. Meniscus cell lysate induces mitochondrial dysfunction of fibroblast-like synoviocytes via upregulating ANT3 in osteoarthritis: ANT3 is a potentially important novel mediator and drug target in osteoarthritis . Bone Joint Res . 2023 ; 12 : 274 – 284 . Crossref PubMed Google Scholar
38. Schmittgen TD , Livak KJ . Analyzing real-time PCR data by the comparative C(T) method . Nat Protoc . 2008 ; 3 ( 6 ): 1101 – 1108 . Crossref PubMed Google Scholar
39. Young MD , Wakefield MJ , Smyth GK , Oshlack A . Gene ontology analysis for RNA-seq: accounting for selection bias . Genome Biol . 2010 ; 11 ( 2 ): R14 . Crossref PubMed Google Scholar
40. Mao X , Cai T , Olyarchuk JG , Wei L . Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary . Bioinformatics . 2005 ; 21 ( 19 ): 3787 – 3793 . Crossref PubMed Google Scholar
41. Teschendorff AE , Zhu T , Breeze CE , Beck S . EPISCORE: cell type deconvolution of bulk tissue DNA methylomes from single-cell RNA-Seq data . Genome Biol . 2020 ; 21 ( 1 ): 221 . Crossref PubMed Google Scholar
42. Zheng SC , Breeze CE , Beck S , Teschendorff AE . Identification of differentially methylated cell types in epigenome-wide association studies . Nat Methods . 2018 ; 15 ( 12 ): 1059 – 1066 . Crossref PubMed Google Scholar
43. Rahmani E , Schweiger R , Rhead B , et al. Cell-type-specific resolution epigenetics without the need for cell sorting or single-cell biology . Nat Commun . 2019 ; 10 ( 1 ): 3417 . Crossref PubMed Google Scholar
44. Moore LD , Le T , Fan G . DNA methylation and its basic function . Neuropsychopharmacology . 2013 ; 38 ( 1 ): 23 – 38 . Crossref PubMed Google Scholar
45. Rice P , Longden I , Bleasby A . EMBOSS: the European Molecular Biology Open Software Suite . Trends Genet . 2000 ; 16 ( 6 ): 276 – 277 . Crossref PubMed Google Scholar
46. Lev Maor G , Yearim A , Ast G . The alternative role of DNA methylation in splicing regulation . Trends Genet . 2015 ; 31 ( 5 ): 274 – 280 . Crossref PubMed Google Scholar
47. No authors listed . The Comparative Toxicogenomics Database (CTD) . 2024 . https://ctdbase.org/ ( date last accessed 12 April 2024 ). Crossref PubMed Google Scholar
48. Sakuta H , Suzuki R , Takahashi H , et al. Ventroptin: a BMP-4 antagonist expressed in a double-gradient pattern in the retina . Science . 2001 ; 293 ( 5527 ): 111 – 115 . Crossref PubMed Google Scholar
49. Yang CR , Ning L , Zhou FH , et al. Downregulation of Adhesion Molecule CHL1 in B Cells but Not T Cells of Patients with Major Depression and in the Brain of Mice with Chronic Stress . Neurotox Res . 2020 ; 38 ( 4 ): 914 – 928 . Crossref PubMed Google Scholar
50. Boularan C , Hwang I-Y , Kamenyeva O , et al. B Lymphocyte-Specific Loss of Ric-8A Results in a Gα Protein Deficit and Severe Humoral Immunodeficiency . J Immunol . 2015 ; 195 ( 5 ): 2090 – 2102 . Crossref PubMed Google Scholar
51. Fuentealba J , Toro-Tapia G , Arriagada C , et al. Ric-8A, a guanine nucleotide exchange factor for heterotrimeric G proteins, is critical for cranial neural crest cell migration . Dev Biol . 2013 ; 378 ( 2 ): 74 – 82 . Crossref PubMed Google Scholar
52. Kato K , Bhattaram P , Penzo-Méndez A , Gadi A , Lefebvre V . SOXC Transcription Factors Induce Cartilage Growth Plate Formation in Mouse Embryos by Promoting Noncanonical WNT Signaling . J Bone Miner Res . 2015 ; 30 ( 9 ): 1560 – 1571 . Crossref PubMed Google Scholar
53. Teruel M , Martin J-E , González-Juanatey C , et al. Association of acid phosphatase locus 1*C allele with the risk of cardiovascular events in rheumatoid arthritis patients . Arthritis Res Ther . 2011 ; 13 ( 4 ): R116 . Crossref PubMed Google Scholar
54. Gagné V , Aubry-Morin A , Plesa M , et al. Genes identified through genome-wide association studies of osteonecrosis in childhood acute lymphoblastic leukemia patients . Pharmacogenomics . 2019 ; 20 ( 17 ): 1189 – 1197 . Crossref PubMed Google Scholar
55. Han J , Deng H , Lyu Y , et al. Identification of N-Glycoproteins of Knee Cartilage from Adult Osteoarthritis and Kashin-Beck Disease Based on Quantitative Glycoproteomics, Compared with Normal Control Cartilage . Cells . 2022 ; 11 ( 16 ): 2513 . Crossref PubMed Google Scholar
56. Lü AL , Guo X , Aisha MMT , Shi XW , Zhang YZ , Zhang YY . Kashin-Beck disease and Sayiwak disease in China: prevalence and a comparison of the clinical manifestations, familial aggregation, and heritability . Bone . 2011 ; 48 ( 2 ): 347 – 353 . Crossref PubMed Google Scholar
57. Zhao S , Geybels MS , Leonardson A , et al. Epigenome-Wide Tumor DNA Methylation Profiling Identifies Novel Prognostic Biomarkers of Metastatic-Lethal Progression in Men Diagnosed with Clinically Localized Prostate Cancer . Clin Cancer Res . 2017 ; 23 ( 1 ): 311 – 319 . Crossref PubMed Google Scholar
58. Zhang B , Groffen J , Heisterkamp N . Resistance to farnesyltransferase inhibitors in Bcr/Abl-positive lymphoblastic leukemia by increased expression of a novel ABC transporter homolog ATP11a . Blood . 2005 ; 106 ( 4 ): 1355 – 1361 . Crossref PubMed Google Scholar
Author contributions
B. Cheng: Data curation, Methodology, Writing – original draft, Funding acquisition.
C. Wu: Data curation, Methodology, Writing – review & editing.
W. Wei: Data curation.
H. Niu: Data curation.
Y. Wen: Data curation.
C. Li: Methodology.
P. Chen: Methodology.
H. Chang: Methodology.
Z. Yang: Data curation, Methodology, Writing – review & editing.
F. Zhang: Funding acquisition, Writing – review & editing.
Funding statement
The authors disclose receipt of the following financial or material support for the research, authorship, and/or publication of this article: support from the National Key Research and Development Program of China (2022YFC2503100), National Natural Science Foundation of China (82304265), and Fundamental Research Funds for the Central Universities (xzy012023092), as reported by B. Cheng and F. Zhang, and support from the Natural Science Basic Research Program of Shaanxi Province (2021JCW-08), as reported by F. Zhang.
ICMJE COI statement
B. Cheng reports grants from Fundamental Research Funds for the Central Universities (xzy012023092) and the National Natural Science Foundation of China (82304265), related to this study.
F. Zhang reports grants from the National Key Research and Development Program of China (grant number 2022YFC2503100) and Natural Science Basic Research Program of Shaanxi Province (2021JCW-08), related to this study.
Data sharing
The data that support the findings for this study are available to other researchers from the corresponding author upon reasonable request.
Acknowledgements
The authors thank Zhen Zhang for the English language editing and review assistance.
Ethical review statement
This study was approved by the Human Ethics Committees of Xi'an Jiaotong University. Written informed consent was obtained from all study participants.
Open access funding
The authors report that they received open access funding for their manuscript from Fundamental Research Funds for the Central Universities (xzy012023092).
Supplementary material
Figures displaying distribution trend of chromosome-wide methylation sequencing data and mean methylation level in the gene body. Tables listing the joint dysfunction index criterion and the primers used in quantitative reverse transcriptase polymerase chain reaction detection.
© 2024 Cheng et al. This is an open-access article distributed under the terms of the Creative Commons Attribution Non-Commercial No Derivatives (CC BY-NC-ND 4.0) licence, which permits the copying and redistribution of the work only, and provided the original author and source are credited. See https://creativecommons.org/licenses/by-nc-nd/4.0/
