The integrative analysis of DNA methylation and mRNA expression profiles confirmed the role of selenocompound metabolism pathway in Kashin-Beck disease

ABSTRACT

Kashin-Beck disease (KBD) is an endemic chronic osteochondropathy. The etiology of KBD remains unknown. In this study, we conducted an integrative analysis of genome-wide DNA methylation and mRNA expression profiles between KBD and normal controls to identify novel candidate genes and pathways for KBD. Articular cartilage samples from 17 grade III KBD patients and 17 healthy controls were used in this study. DNA methylation profiling of knee cartilage and mRNA expression profile data were obtained from our previous studies. InCroMAP was performed to integrative analysis of genome-wide DNA methylation profiles and mRNA expression profiles. Gene ontology (GO) enrichment analysis was conducted by online DAVID 6.7. The quantitative real-time polymerase chain reaction (qPCR), Western blot, immunohistochemistry (IHC), and lentiviral vector transfection were used to validate one of the identified pathways. We identified 298 common genes (such as COL4A1, HOXA13, TNFAIP6 and TGFBI), 36 GO terms (including collagen function, skeletal system development, growth factor), and 32 KEGG pathways associated with KBD (including Selenocompound metabolism pathway, PI3K-Akt signaling pathway, and TGF-beta signaling pathway). Our results suggest the dysfunction of many genes and pathways implicated in the pathogenesis of KBD, most importantly, both the integrative analysis and in vitro study in KBD cartilage highlighted the importance of selenocompound metabolism pathway in the pathogenesis of KBD for the first time.

KEYWORDS:

• Kashin-Beck disease
• DNA methylation
• mRNA expression profile

Introduction

Kashin-Beck disease (KBD) is an endemic chronic osteochondropathy, distributed in certain areas of Eastern Siberia of Russia, North Korea, and China [1]. According to China Health Statistics Yearbook, 37.9 million people living at KBD prevalent counties are at the risk of KBD, and 0.53 million patients suffered from KBD in China in 2017, including 8,540 children under the age of 16. KBD mainly damages cartilage, especially in growth plate and joints [2]. It usually occurs during childhood between 3 and 12 years and becomes symptomatic in adulthood. The typical manifestations of KBD include serious joint pain, joint deformities, irreversible aggravated osteoarthritis (OA), and growth retardation [3].

Various etiological hypotheses have been proposed for KBD, mainly including environmental selenium deficiency and cereal contamination by myco-toxin-producing fungi (mainly for T-2 toxin) [1]. In the 1970s, Chinese researchers discovered that a low environmental content of selenium resulted in a poor nutritional selenium uptake in the KBD prevalent areas [1]. Selenium is an essential element for human health and relates to many diseases, such as cancer, cardiovascular, diabetes, especially Keshan Disease (KD) and KBD [1]. KBD is distributed in selenium-deficient areas in a diagonal belt from northeast to southwest China in the KBD endemic areas, the level of selenium in both soil and human biological samples is in lower amounts than non-KBD areas. The hair selenium level of school children in KBD prevalent areas is significantly lower than that of children living in non-KBD prevalent areas [4]. Consequently, selenium deficiency is an endemic risk factor for KBD, but the mechanism of environmental selenium deficiency implicated in the pathogenesis of KBD remains elusive. Besides, T-2 toxin has been identified in endemic food samples in KBD prevalent areas. The chicks fed with food contaminated by T-2 toxin have shown the similar cartilage damages with KBD patients [5]. Therefore, T-2 toxin has been used as toxic treatment in human chondrocyte line (C28/I2) for years to simulate the similar mechanism in KBD chondrocytes [X. 6].

Benefiting from the rapid development of high throughput analyzing technologies, several genomic studies of KBD samples identified multiple genes associated with KBD. For instance, genetic association studies identified HLA-DRB1 [S. 7] and ITPR2 [F. 8] are susceptibility genes for KBD. Gene expression profile analysis identified 79 differently expressed genes between KBD and normal cartilage [W. Z. 9]. DNA methylation is the most extensively studied epigenetic modulation [Y. 10]. Recent studies demonstrated that abnormal DNA methylation was strongly correlated with cartilage damage in OA. Roach et al. [11] confirmed that epigenetic changes contributed to the altered synthesis of cartilage degrading enzymes in the chondrocytes of late-stage OA patients. Jeffries et al. [12] observed abnormal methylation of OA susceptibility genes associated with both OA subchondral bone and cartilage. However, the role and mechanism of DNA methylation in the development of KBD are still nebulous.

It is acknowledged that traditional single omic (for instance, genome, or transcriptome) mostly could not unveil the elusive molecular mechanisms of complex diseases. Therefore, integrative multi-omics datasets are becoming valuable. It could not only lead to lower false discovery rates, but also provide more biological insights than single omic analysis [13, M. 14]. Moreover, this method could reveal the biological complexity and overcome the weakness of single omics, the lack of comprehensiveness, and the information gap regarding interacting across multiple omic layers [15]. By using tans-omics integrative strategy, many studies have uncovered the implicate dysfunction of multiple molecular levels in complex diseases. Several tools for users to integrate the different omics are widely used now, such as InCroMAP [16,17], a powerful, high-level cross-platform microarray dataset analysis tool for generic or pathway-based analysis, visualization of heterogeneous, and cross-platform datasets.

In this study, we performed integrating analysis of DNA methylation profiles and mRNA expression profiling of KBD and normal cartilage to identify novel candidate genes and pathways for KBD. Our results suggest the dysfunction of many genes and pathways are implicated in the pathogenesis of KBD. Most importantly, both the integrative analysis and in vitro study in KBD cartilage highlighted the importance of selenocompound metabolism pathway in the pathogenesis of KBD for the first time.

Materials and methods

Ethical approval

This investigation was approved by the Human Ethics Committee of Xi’an Jiaotong University (Project Number: 2015–237). The study was conducted in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. A written informed consent was obtained from all the patients or relatives of donors.

Sample collection

Human articular cartilage samples from 17 grade III KBD and 17 normal controls were used in this study (Table 1). Both KBD and normal subjects were of Chinese Han lineage and matched with age. KBD was diagnosed according to the diagnosis criteria of KBD in China (WS/T 207–2010). Clinical and radiography examinations were taken carefully for each subject. Other demographic characteristics of human subjects such as lifestyle habits, family history, health status, and ethnicities were collected from nurse-administered questionnaires. The subjects with genetic bone and cartilage diseases, rheumatoid arthritis (RA), and other skeletal disorders were excluded from this study.

Table 1. Basic characteristics of study samples.

Method	DNA methylation profile		Gene expression profile		Western Blot		IHC
Samples	Control	KBD	Control	KBD	Control	KBD	Control	KBD
Number	5	5	4	4	5	5	3	3
Age(Mean±SD)	51.0 ± 7.3	57.4 ± 7.1	48.5 ± 8.3	49.5 ± 5.6	66.2 ± 13.6	57.4 ± 8.0	54.0 ± 16.1	53.0 ± 3.0
Male/Female	3/2	3/2	3/1	3/1	1/4	1/4	1/2	2/1

Genome-wide DNA methylation profiling

The cartilage specimens of 5 KBD and control pairs were harvested from knee joints [W. 18]. QIAamp DNA Blood Mini Kit (QIAGEN, Germany) was used to extract DNA from the cartilage specimens according to the manufacturer’s recommendations. Samples were dissected and frozen in liquid nitrogen rapidly for genome-wide DNA methylation profiling. The quality of extracted DNA was evaluated by agarose gel electrophoresis.

Illumina Infinium HumanMethylation450 BeadChip was used for genome-wide DNA methylation profiling. A total of 500 ng of DNA were used for bisulfite conversion according to the standard protocol for the EZ DNA Methylation Kit (Zymo Research, USA). Then, the samples were incubated overnight at 37°C followed by the addition of whole genome amplification. After that, we use HumanMethylation450 array for hybridization, the array was washed and stained following the protocol guide of Infinium HD Assay Methylation. Then IScan SQ scanner (iScan System, Illumina, USA) was used for scanning the raw image intensities. The obtained raw data of image intensity were processed using GenomeStudio software (Illumina, USA).

Differential methylation analysis

A β value, varying from 0 (completely unmethylated) to 1 (completely methylated), was defined as the expression of the average percentage of methylated cytosine at a given CpG site. Differentially methylated CpG sites were identified using the empirical Bayes moderated t-test in limma [19], which was described by Illumina Methylation Analyzer package. A false discovery rate (FDR) adjusted P-value for each CpG site was calculated by using Benjamini–Hochberg method. The Pearson correlation coefficients were calculated to evaluate the correlations of cartilage specimens.

For quality control, the quality of extracted DNA was evaluated by agarose gel electrophoresis. The DNA specimens passing the quality control were used for DNA methylation profiling. The definition of significant CpG met the following two conditions: (1) FDR adjusted P-value (P_adj) ≤0.05; (2) |β difference (∆β)| ≥0.2. The CpG sites with P_adj>0.05, or P_adj ≤ 0.05 with missing values in more than 90% of cartilage specimens accompanied by P > 0.05 were eliminated.

Gene expression profiling of KBD and normal cartilage

Our previous gene expression profiling of KBD and normal cartilage was used here [W. Z. 9]. In brief, the articular cartilages of Four KBD patients and four normal controls were collected and matched according to age and sex. Four KBD patients from Yongshou County and Linyou County of Xi’an city in China underwent free total knee replacement surgery, and their articular cartilage tissues were collected. Four matched healthy control knee cartilage samples were collected from the donors within 8 hours of death caused by traffic accidents, the status of healthy samples was diagnosed by histological examination with H&E staining, for excluding KBD, genetic bone, cartilage diseases, OA and RA. Detailed procedures for the experiment were described in our previous study [W. Z. 9].

Integrated analysis of DNA methylation and mRNA expression profiles

The DNA methylation and mRNA expression profiling datasets were input into InCroMAP for integrated analysis. Before running analysis, the filtering thresholds of InCroMAP were set as p-value <0.05 for DNA methylation profile data and fold-change (FC) ≥1.5 for mRNA expression profile data. After input, the uniform data required by InCroMAP, integration of multiple data types was performed by mapping each probe to every individual gene. This process is straightforward for mRNA datasets, while DNA methylation datasets are region based on genes by defining the upstream and downstream of each gene’s transcription start site [17]. On the basis of data-pairing, the common genes altered both in DNA methylation profiles and mRNA expression profiles were identified. InCroMAP supports gene sets from the KEGG PATHWAY database. It could also visualize all the pathways by using KEGG translator. Therefore, the mRNA expression or protein modification’s expression value was highlighted by node color or small box. After that, KEGG pathway database provided by InCroMAP was used for pathway enrichment analysis in this study [17]. At last, the identified common genes were submitted to the DAVID 6.7 functional annotation analysis tool (http://david.abcc.ncifcrf.gov) for Gene Ontology (GO) enrichment analysis.

Primary chondrocyte isolation and culture

In order to confirm the involvement of specific genes determined by the integrated analysis of DNA methylation and mRNA expression profiles, we used an in vitro cell model of KBD. The human chondrocyte line (C28/I2), obtained from Mary B. Goldring [20], treated with T-2 toxin[X. 6] for detecting the mRNA expression of significant genes. Then, another independent human articular chondrocytes samples from eight grade III KBD patients and 8 normal controls of femoral neck fracture (FNF) were used for Western blot (5 pairs) and immunohistochemistry (IHC) (3 pairs) to identify the expression of SEPSECS protein in KBD articular cartilage (Table 1).

KBD and control cartilage samples were cut into 5–10 mm³ slices, treated with 0.25% trypsin (Xi’an GuoAn Biological Technology Co.) for 30 minutes at room temperature (RT) and immersed into 0.2% type II collagenase solution (Gibco, 17,101–015) at 37°C for 10 hours. After that, the isolated chondrocytes were cultured at 37°C in 5% CO₂ in DMEM/F12 supplemented with 10% fetal calf serum (Sijiqing, Zhejiang Tianhang Biotechnology Company, Hangzhou, China) and 1% penicillin/streptomycin.

3-(4,5-dimethyl-2-thiazolyl)-2,5-diphenyl-2-H-tetrazolium bromide, thiazolyl blue tetrazolium bromide (MTT) assay

C28/I2 cell line was cultured at 37°C in 5% CO₂ in DMEM/F12 supplemented with 10% fetal calf serum and 1% penicillin/streptomycin. C28/I2 and T-2 toxin were used to make the KBD cell models as described previously [X. 6]. We used MTT assay to determine the cytotoxicity in C28/I2 human chondrocytes after treated with T-2 toxin. C28/I2 human chondrocytes were seeded in 96-well plates and treated with T-2 toxin at eight different concentrations (1, 2, 5, 8, 10, 12, 15, 20 ng/ml) for 12 h, 24 h, 48 h and 72 h separately. Each concentration has five repeats. Then, 20ul of MTT (5 ng/ml) was added in each wells for 4 hours. After that, the medium was discarded and 200ul dimethyl sulfoxide (DMSO) was added to dissolve the formazan. The results were measured by an automatic enzyme-linked immunosorbent assay reader at 490 nm [X. 6].

Quantitative real-time polymerase chain reaction (qPCR)

Total RNA from C28/I2 and T-2 toxin (treated C28/I2 cells 2 ng/ml, 5 ng/ml, 8 ng/ml for 24 h) isolated according to TRIzol protocol and reverse-transcribed to cDNA by One step TaKaRa PrimeScript^TM RT reagent Kit. GAPDH was chosen as reference gene [X. 6]. mRNA expression of identified genes (CTH, SEPSECS and PAPSS2) from Selenocompound metabolism pathway was measured by qPCR in triplicate analysis (Table 2). The cycle of threshold (CT) of each sample was averaged and normalized to that of GAPDH. The triplicate results were analyzed by comparative CT equation (2^−ΔΔCT).

Table 2. The primers used for qPCR.

Gene	Forward	Reverse
SEPSECS	5ʹ-GGCTGTGGCAGTGTGCTTAGAA-3′	5ʹ-TCCTAAGGCTGGCGTCACATATC-3′
CTH	5ʹ-ATCCACAGCATGAGTTGGTGAAG-3′	5ʹ-CCAAGCTCTCGGCCAGAGTAAATA-3′
PAPSS2	5ʹ-AATGTGGCGATGTCGTGGA-3′	5ʹ-TGCCTTGCTAGTCTGATAGTGTGAA-3′
MMP13	5ʹ-TCCCAGGAATTGGTGATAAAGTAGA-3′	5ʹ-GCATGACGCGAACAATACGG-3′
SOX9	5ʹ-GGAGATGAAATCTGTTCTGGGAATG-3′	5ʹ-TTGAAGGTTAACTGCTGGTGTTCTG-3′
ADAMTS5	5ʹ-TCAGCTTTGGCAACCTTCAAATC-3′	5ʹ-AGGTCCCAAACGGCTCAGTTC-3′
GAPDH	5ʹ-GCACCGTCAAGGCTGAGAAC-3′	5ʹ-TGGTGAAGACGCCAGTGGA-3′

Western blot

Radio immunoprecipitation assay reagent (Xi’an Hat Biotechnology Co., Ltd) was used for isolating the protein from KBD and control chondrocyte. The concentration of the protein was measured by BCA protein assay kit (Xi’an Hat Biotechnology Co., Ltd). Then, proteins were equally added to the SDS-polyacrylamide gel electrophoresis (SDS-PAGE) and ran separately before transferred to polyvinylidene difluoride (PVDF) membranes. The PVDF membranes were blocked in 5% nonfat milk solution and incubated with the SEPSECS primary antibody (1:1000, ab56033, abcam) at 4°C overnight and the secondary antibody (1:5000, Xi’an Zhuangzhi Biotechnology Co., Ltd, Cat: EK020) for 1 hour at RT. Gene Gnome XRQ (Gene Co., Ltd) was used to analyze the protein bands.

Immunohistochemistry (IHC)

The independent three grade III KBD patients and three control articular cartilage samples were used here. First, samples were embedded in paraffin for IHC. The procedure includes cutting sections, dewax, rehydrate, and antigen retrieval using citrate buffer (pH 6.0) overnight at 37°C and 12.5% trypsin (Xi’an GuoAn Biological Technology Co.) at 37°C for 20–30 minutes. Zhong Shan Gold Bridge Rabbit SP reagent was used in the following procession (Beijing Zhong Shan Gold Bridge Biological Technology Co., 18112A11). Three percent H₂O₂ was used to prevent endogenous peroxidase reactions for 10 minutes at RT (white bottle), and serum was used to block the samples for 20 minutes at RT to avoid nonspecific background (blue bottle). Then, primary antibody of SEPSECS (1:60, ab56033, abcam) treated the samples at 4°C overnight. The slices were rewarmed at 37°C the next day for 1 h at RT before treating with secondary antibody (yellow bottle) and horseradish peroxidase (red bottle) for 20 minutes at RT, respectively. After that, DAB solution treated the sections until the appropriate color could be seen before dehydrating, mounting the sections, and dyeing the nucleus.

Lentiviral vectors transfection

In order to further detect the role of SEPSECS in the mechanism of KBD cell model. The recombinant lentivirus of SEPSECS acquired from Genechem (Shanghai, China) was used here. C28/I2 human chondrocytes were first treated with 8 ng/ml T-2 toxin for 24 h, after that, the SEPSECS lentivirus and negative lentivirus were transfected into C28/I2 for 96 h. Then, the mRNA expression of SEPSECS and the selected genes (MMP13, SOX9, and ADAMTS5) were measured by qPCR (Table 2).

Statistical analysis

A one-way analysis of variance (ANOVA) and independent samples t-test were used for comparisons of the results. Statistical analyses were performed by SPSS 19.0. The threshold was defined as P < 0.05.

Results

DNA methylation profile

In total 473,607 methylation sites analyzed across the genome by using Illumina Infinium HumanMethylation450 BeadChip and 10,161 methylated CpG sites were identified significantly changed (P < 0.05) between KBD and normal samples, including 4,027 hypermethylation CpG cites (corresponding to 1,047 genes) and 6,134 hypomethylation CpG cites (corresponding to 1,747 genes). Table 3 shows the top 20 significant hypermethylation and hypomethylation CpG sites, among which, cg07687119 (P-value = 3.11 × 10⁻⁹, ∆β = 0.345) with the target genes (HOXC4, HOXC6, and HOXC5) is the most significant hypermethylation locus for KBD. Meanwhile, cg03223172 (P-value = 1.28 × 10⁻⁹, ∆β = −0.294) is the most significant hypomethylation locus for KBD; however, no target gene was identified with this locus. Multiple genes were identified significant differentially methylated between KBD and normal controls, such as SEPSECS (cg20177285, P-value = 3.47 × 10⁻⁵, ∆β = −0.4166), TNFAIP6 (cg01189638, P-value = 1.33 × 10⁻³, ∆β = −0.2732). Some genes are related more than one CpG site. For instance, COL4A1 is related to cg00050792 (P-value = 2.07 × 10⁻³, ∆β = 0.2958), cg08234256 (P-value = 9.19 × 10⁻³, ∆β = 0.2148) and cg03564009 (P-value = 1.09 × 10⁻³, ∆β = 0.2948).

Table 3. The top 20 significant hypermethylation and hypomethylation CpG sites.

Type	CpG cite	P-value	Padj-value	∆β	Gene
Hypermethylation	cg07687119	3.11 × 10^−9	3.69 × 10^−4	0.345	HOXC4;HOXC6;HOXC5
	cg13763232	6.29 × 10^−9	5.96 × 10^−4	0.2436	SLC6A6
	cg18503031	3.09 × 10^−8	2.44 × 10^−3	0.2958	MITF
	cg05401764	5.89 × 10^−8	2.94 × 10^−3	0.5756	SHOX2
	cg01214458	9.79 × 10^−8	3.26 × 10^−3	0.2028	MN1
	cg01410314	1.03 × 10^−7	3.26 × 10^−3	0.2048	C9orf41
	cg01575836	1.23 × 10^−7	3.26 × 10^−3	0.4202
	cg24937727	1.58 × 10^−7	3.26 × 10^−3	0.3358	RGL3
	cg07464206	2.15 × 10^−7	3.60 × 10^−3	0.3282
	cg06359632	2.33 × 10^−7	3.60 × 10^−3	0.442
	cg15964309	2.34 × 10^−7	3.60 × 10^−3	0.3868
	cg10928103	2.73 × 10^−7	3.60 × 10^−3	0.3786
	cg17848054	3.09 × 10^−7	3.60 × 10^−3	0.241	RGL3
	cg00862290	4.08 × 10^−7	3.86 × 10^−3	0.4134	KCNMB3
	cg20266894	4.08 × 10^−7	3.86 × 10^−3	0.401
	cg22799141	5.05 × 10^−7	4.05 × 10^−3	0.3512
	cg04403850	5.13 × 10^−7	4.05 × 10^−3	0.3178
	cg02532853	5.68 × 10^−7	4.13 × 10^−3	0.2698	NCRNA00114
	cg03255182	5.86 × 10^−7	4.13 × 10^−3	0.301	HOXC4;HOXC6;HOXC5
	cg04875697	6.00 × 10^−7	4.13 × 10^−3	0.2044	RXRA
Hypomythylation	cg03223172	1.28 × 10^−9	3.04 × 10^−4	−0.294
	cg16549994	6.88 × 10^−8	2.94 × 10^−3	−0.2026
	cg19291576	6.90 × 10^−8	2.94 × 10^−3	−0.3832
	cg27159719	7.45 × 10^−8	2.94 × 10^−3	−0.2292	TMEM71
	cg21400851	1.40 × 10^−7	3.26 × 10^−3	−0.308	KCTD2
	cg22786333	1.52 × 10^−7	3.26 × 10^−3	−0.2962	HEATR2
	cg24193140	1.84 × 10^−7	3.53 × 10^−3	−0.3038	KCTD16
	cg15058852	1.99 × 10^−7	3.53 × 10^−3	−0.3516	GNS
	cg00347729	2.62 × 10^−7	3.60 × 10^−3	−0.292	MMP10
	cg14378231	2.69 × 10^−7	3.60 × 10^−3	−0.4056
	cg27050407	2.84 × 10^−7	3.60 × 10^−3	−0.611	FAM65B
	cg20746459	3.00 × 10^−7	3.60 × 10^−3	−0.4686	C1orf125;NPHS2
	cg09501687	3.05 × 10^−7	3.60 × 10^−3	−0.3636	C5orf62
	cg01420491	3.12 × 10^−7	3.60 × 10^−3	−0.2736	HS1BP3
	cg07954691	3.23 × 10^−7	3.65 × 10^−3	−0.514
	cg09060914	4.14 × 10^−7	3.86 × 10^−3	−0.4484	SELP
	cg12302621	4.52 × 10^−7	3.91 × 10^−3	−0.2294	SLC28A1
	cg06133897	4.58 × 10^−7	3.91 × 10^−3	−0.212
	cg22659564	4.66 × 10^−7	3.91 × 10^−3	−0.24	FAT3
	cg21332304	5.23 × 10^−7	4.06 × 10^−3	−0.3174	NPEPPS

Note: ∆β > 0,hypermethlation; ∆β < 0, hypomethlation.

Integrative analysis of DNA methylation profile and mRNA expression

To reveal the complicated mechanism of KBD, we conducted a genome-wide integrative analysis of DNA methylation profiles and mRNA expression profiles of KBD and normal cartilage. After considering Gene Bank ID, P_methyl <0.05 (DNA methylated CpG sites) and fold-change (FC) ≥1.5 (mRNA expression files), 7,526 methylated CpG cites and 2,400 mRNA expression profiles between KBD and normal controls were subjected to the InCroMAP. Two hundred ninety-eight commonly affected genes in both DNA methylation profiles and mRNA expression profiles between KBD and normal controls were identified (Additional File: Table S1). These genes include the collagen family members, such as COL1A1 (NM_000088, P_methyl = 4.42 × 10⁻³, FC = 4.92), COL1A2 (NM_000089, P_methyl = 2.26 × 10⁻³, FC = 2.05), COL4A1 (NM_001845, P_methyl = 4.12 × 10⁻³, FC = 2.16), COL4A2 (NM_001846, P_methyl = 8.05 × 10⁻³, FC = 2.42), COL5A2 (NM_000393, P_methyl = 2.09 × 10⁻³, FC = 5.66), COL9A3 (NM_001853, P_methyl = 0.02, FC = 2.81) and COL14A1 (NM_021110, P_methyl = 4.47 × 10⁻³, FC = 1.61), and the homeobox family, such as HOXA13 (NM_000522, P_methyl = 0.013, FC = 1.84), HOXB2 (NM_002145, P_methyl = 1.03 × 10⁻³, FC = 5.08), HOXC9 (NM_006897, P_methyl = 9.15 × 10⁻³, FC = 2.22), and HOXD10 (NM_002148, P_methyl = 0.028, FC = 1.59), altered both in DNA methylation and mRNA expression in KBD cartilage compared to normal controls. Moreover, some other genes, such as TNFAIP6 (NM_007115, P_methyl = 1.34 × 10⁻³, FC = 5.94) and TGFBI (NM_000358, P_methyl = 3.21 × 10⁻³, FC = 7.86), which have been found as significant genes in KBD patients previously [21, S. 22] also identified significantly in this study.

Enrichment analysis

GO enrichment analysis of these 298 genes identified 36 significantly overrepresented GO terms (FDR<0.05), including those functionally involved in collagen, skeletal system development, growth factor, and apoptosis (Figure 1, Table S2).

Figure 1. The barplot of GO enrichment analysis results of common genes altered in both DNA methylation profiles and mRNA expression profiles. BP stands for biological process. CC stands for cellular component. MF stands for molecular function.

Integrative pathway enrichment analysis of DNA methylation profiles and mRNA expression profiles demonstrated 32 pathways significantly associated with KBD (Table 4), such as Selenocompound metabolism pathway (KEGG path: hsa00450, P = 0.0123, Figure 2), PI3K-Akt signaling pathway (KEGG path: hsa04151, P = 3.40 × 10⁻³, Figure 3), TGF-beta signaling pathway (KEGG path: hsa04350, P = 4.30 × 10⁻³), and Regulation of actin cytoskeleton (KEGG path: hsa04810, P = 3.95 × 10⁻³).

Table 4. Integrated pathway enrichment analysis results.

KEGG ID	Pathway name	P-value	NO. of KEGG Genes	NO. of Significant Genes in KBD
hsa05200	Pathways in cancer	6.14 × 10^−4	526	91
hsa04510	Focal adhesion	7.03 × 10^−4	199	56
hsa04151	PI3K-Akt signaling pathway	3.40 × 10^−3	354	82
hsa04810	Regulation of actin cytoskeleton	3.95 × 10^−3	213	55
hsa04350	TGF-beta signaling pathway	4.30 × 10^−3	84	25
hsa04512	ECM-receptor interaction	5.43 × 10^−3	82	25
hsa04144	Endocytosis	6.50 × 10^−3	244	50
hsa00010	Glycolysis/Gluconeogenesis	6.98 × 10^−3	68	17
hsa05205	Proteoglycans in cancer	8.00 × 10^−3	203	56
hsa01200	Carbon metabolism	8.35 × 10^−3	63	26
hsa05169	Epstein-Barr virus infection	8.64 × 10^−3	201	53
hsa04141	Protein processing in endoplasmic reticulum	1.01 × 10^−2	165	47
hsa05222	Small cell lung cancer	1.21 × 10^−2	93	26
hsa00450	Selenocompound metabolism	1.23 × 10^−2	17	7
hsa04390	Hippo signaling pathway	1.27 × 10^−2	154	38
hsa04530	Tight junction	1.41 × 10^−2	170	27
hsa05215	Prostate cancer	1.41 × 10^−2	97	27
hsa05202	Transcriptional misregulation in cancer	1.43 × 10^−2	186	50
hsa04730	Long-term depression	1.55 × 10^−2	60	18
hsa04670	Leukocyte transendothelial migration	1.62 × 10^−2	112	28
hsa04710	Circadian rhythm	1.63 × 10^−2	31	11
hsa00620	Pyruvate metabolism	2.33 × 10^−2	39	12
hsa00380	Tryptophan metabolism	2.40 × 10^−2	40	14
hsa00562	Inositol phosphate metabolism	2.40 × 10^−2	74	15
hsa05213	Endometrial cancer	2.91 × 10^−2	58	16
hsa05218	Melanoma	2.91 × 10^−2	72	16
hsa05014	Amyotrophic lateral sclerosis (ALS)	3.39 × 10^−2	51	17
hsa05219	Bladder cancer	3.66 × 10^−2	41	14
hsa05100	Bacterial invasion of epithelial cells	3.72 × 10^−2	74	19
hsa00020	Citrate cycle (TCA cycle)	3.81 × 10^−2	30	10
hsa05221	Acute myeloid leukemia	4.59 × 10^−2	66	17
hsa00532	Glycosaminoglycan biosynthesis – chondroitin sulfate/dermatan sulfate	4.88 × 10^−2	20	8

Figure 2. InCroMAP analysis results of Selenocompound metabolism pathway. The genes highlighted in yellow denote the genes analyzed by this study within the Selenocompound metabolism pathway.

Figure 3. InCroMAP analysis results of PI3K-ATK signaling pathway. The genes highlighted in yellow denote the genes analyzed by this study within the PI3K-ATK signaling pathway.

In vitro confirmation of selected genes

To ensure there are sufficient cells to perform mRNA expression analysis and simulate the similar mechanism of KBD chondrocytes, we used MTT assay to follow cell survival rate to response to different concentrations of T-2 toxin. Our results showed that C28/I2 chondrocyte survival rate was decreased with the increasing concentration of T-2 toxin. The median inhibitory concentration (IC50) of T-2 toxin treated on C28/I2 for each time point was 14.90 ng/ml, 11.98 ng/ml, 7.05 ng/ml, and 6.52 ng/ml for 12 h, 24 h, 48 h, and 72 h, respectively. The least significant difference (LSD) method detected significant differences between the four different time groups and control (F = 4.237, P = 0.006). 24 h, 48 h, 72 h showed significant differences when compared with control (P = 0.028, 0.001, 0.001, respectively). Concentration analysis identified that T-2 toxin-induced significant differences among the eight concentration groups (F = 3.631, P = 0.004). In this study, we chose the concentration of T-2 toxin as 2 ng/ml, 5 ng/ml, and 8 ng/ml for 24 h as our best protocol to ensure adequate cell response and sufficient cell survival for further analysis.

The mRNA expression of selected genes in KBD cell model

A set of 3 distinct genes from the significantly affected Selenocompound metabolism pathway was selected and gene expression levels quantified in the KBD in vitro model that C28/I2 cell line treated with T-2 toxin (2 ng/ml, 5 ng/ml, and 8 ng/ml) for 24 h. We observed that the mRNA expression levels of three genes (SEPSECS, CTH, and PAPSS2) from the Selenocompound metabolism pathway were all significantly increased in C28/I2 chondrocytes treated with T-2 comparing with control (P_SEPSECS<1.0 × 10⁻⁴, P_CTH = 4.0 × 10⁻³, P_PAPSS2 = 7.2 × 10⁻³) (Figure 4).

Figure 4. The mRNA expression of SEPSECS (a), CTH (b), PAPSS2 (c) in C28/I2 in KBD cell models. The x-axis indicates different T-2 toxin intervention concentrations. The y-axis represents the mRNA expression levels of the genes compared to normal control. **** means P < 0.0001, and ** means P < 0.01.

After treated by 8 ng/ml T-2 toxin for 24 h in C28/I2 human chondrocytes, the mRNA expression of SEPSECS was significantly increased (P < 0.001, Figure 5(a)). Then, after knocking out SEPSECS by lentivirus transfection, the mRNA expression of SEPSECS was significantly decreased (P < 0.001, Figure 5(a)), which means the effectiveness of knocking out SEPSECS in C28/I2 human chondrocytes treated with T-2 toxin. The mRNA expression of MMP13 was decreased notably (P < 0.001, Figure 5(b)), while the SOX9 was overexpressed after knocking out SEPSECS in C28/I2 human chondrocytes (P < 0.01, Figure 5(c)), but the mRNA expression of ADMATS5 had no significant difference between the lentivirus group and T-2 toxin group (Figure 5(d)).

Figure 5. The mRNA expression of SEPSECS (a), MMP13 (b), SOX9 (c) and ADAMTS5 (d) in C28/I2 in KBD cell models. The y-axis represents the mRNA expression levels of the genes compared to normal control. *** means P < 0.001, ** means P < 0.01, * means P < 0.05.

The expression of SEPSECS in KBD cartilage

The expression of SEPSECS protein by IHC was observed in the superficial, middle, and deep layer of KBD and control articular cartilage. The results showed that the expression of SEPSECS protein was significantly higher in superficial layer of KBD samples compared with control samples (P = 0.021), while the middle and deep layer showed no significance (Figure 6). In addition, western blot results showed that the expression of SEPSECS protein (54 kDa) was 2.19 times higher in KBD chondrocytes than that in normal controls (P = 7.10 × 10⁻³) (Figure 7).

Figure 6. The expression of SEPSECS protein in the KBD and control cartilage. Immunostaining for SEPSECS in the superficial, middle and deep zone of KBD cartilage (a). The percentage of positive of SEPSECS in the superficial (b), middle (c), deep layer (d) of KBD and control cartilage. *P = 0.021 in the superficial zone, n = 3. Scale bar stands for 50 μm.

Figure 7. The expression of SEPSECS protein in the KBD and control chondrocytes. The expression of SEPSECS protein in KBD chondrocytes and normal chondrocytes (a). **P = 7.10 × 10⁻³ (b).

Discussion

In this study, a trans-omic analysis of DNA methylation and mRNA expression profiling between KBD and normal samples was performed by using the cross-platform microarray dataset analysis tool InCroMAP. Two hundred ninety-eight common genes, 36 GO terms, and 32 KEGG pathways were identified associated with the pathogenesis of KBD.

Among all the 298 common genes, we found several genes of the collagen type family that have significantly altered, such as COL1A1, COL1A2, COL4A1, COL4A2, COL5A2, COL9A3, and COL14A1. The Collagen family is the most abundant proteins in the extracellular matrix (ECM) and is essential for bone adhesion and formation. Different collagen types based on molecular defects of collagen genes are related to many disorders, such as chondrodysplasias, osteogenesis imperfecta, epidermolysis bullosa, OA, and osteoporosis [23]. COL4A1 promotes osteoblastic cells survival and ECM formation and is related to osteoporosis in KBD patients [24]. Additionally, Wang et al. [S. 25] demonstrated COL5A2 is a significant gene in KBD chondrocytic but not in peripheral blood monocytic networks. Type I collagen (COL1) is participated in the repair of articular cartilage [26]. Therefore, we suggested that the candidate genes of COL1A1 and COL1A2 may be involved in the cartilage repair of KBD. Since few studies have showed collagen family genes (COL1A1, COL1A2, COL4A2, COL9A3, COL14A1) identified in this study related to KBD, our findings indicated the important roles of collagen family in the molecular mechanism of KBD articular cartilage.

Another interesting series of genes we found related to KBD is homeobox family (HOXA13, HOXB2, HOXC9, and HOXD10). Homeobox genes are a group of transcription genes that function as the formation of many body structures during early embryonic development. HOXA13 encodes a DNA-binding transcription factor and may regulate gene expression, morphogenesis, and differentiation. It has been reported that mutated HOXA13 associated with the limb and genitourinary deformity development and can cause hand-foot-genital syndrome [27]. Additionally, HOXD10 is expressed in the development of limb buds and is involved in differentiation and limb development. These limb deformities related genes show the similar syndrome in KBD patients, indicating their possible relation to KBD. Therefore, the roles of these homeobox family genes involved in the development of KBD should be attended to.

Moreover, we also confirmed other significant genes that have been demonstrated associated with KBD, such as TNFAIP6 and TGFBI. TNFAIP6, known as tumor necrosis factor alpha-induced protein 6, is a member of the hyaluronan-binding protein family and plays an important role in ECM formation, inflammation and cell proliferation [S. 25]. Wang et al. [S. 22] identified TNFAIP6 as a significant gene associated with the biological function of chondrocyte and cartilage of KBD. In addition, TGFBI, transforming growth factor beta-induced, which encodes a collagen-related protein to bind the type I, II, and IV collagens, plays a critical role in osteoblast differentiation, tissue regeneration, bone remodeling, cell collagen interaction, even affects bone resorption, recovery, and cartilage development [S. 25, P. 28]. TGFBI has been identified as a significant gene from the chondrocytic canonical pathway involved in KBD [S. 25].

Further GO enrichment analysis of the 298 common genes observed numerous GO terms involved in the development of KBD, such as cell adhesion, collagen, skeletal system development, cell death, apoptosis, cytoskeleton, ECM. Some of them have been identified playing important roles in the cartilage damage of KBD, like collagen and ECM. Cartilage contains chondrocytes as well as their surrounding ECM, which is a complex network formed by water, collagen, proteoglycans, and other proteins. Collagens are the most abundant family of ECM proteins, making up two-thirds of the dry mass of adult articular cartilage [Y. 29]. Multiple collagen subtypes have been identified in articular cartilage, including the most major collagen (type II collagen) and some minor collagens (type IX, X, XI, VI, VII collagens). It has been known that the cartilage degeneration caused by abnormal catabolism in the ECM is the chief characteristic in the pathogenesis of KBD [X. 6]. Wang et al. demonstrated that Wnt/β-catenin signaling pathway involved in disrupting anabolism and catabolism in the ECM in KBD pathogenesis [X. 6]. Taken together, the identified GO terms may provide us important biological process, cellular component, and molecular function involved in the pathogenesis of KBD.

We identified 32 significantly modified KEGG pathways statistically associated with KBD, such as focal adhesion and PI3K-Akt signaling pathway. The focal adhesion complex is a signaling center regulating dynamic cellular behaviors, such as cell adhesion, spreading, migration, proliferation, and initiated [30]. Shin et al. identified inhibition of the focal adhesion complex could result in a significant reduction of chondrocyte dedifferentiation [30]. Another study performing whole-exome sequencing found focal adhesion is a major pathway for KBD pathogenesis. Interestingly, the focal adhesion complex is linked to the intracellular signaling pathway such as PI3K [30], which is a pathway that also identified significant in this study. PI3K is the primary signaling pathway regulating chondrocyte survival and apoptosis [31]. Du et al. demonstrated that SNP of Selenoprotein S SEPS1 − 105 G > A associated with the expression of PI3K-Akt signaling pathway played important roles in the development of KBD [31]. This evidence indicated that chondrocyte apoptosis was partly regulated by PI3K-Akt.

It has been known that low environmental selenium has a close relation to KBD. Selenium, whose anti-inflammatory, antioxidant, and antiviral properties are well known, shows biologically ambivalent characteristics in animals and exerts its biological effects mainly through enzymatically active selenoproteins [32]. Selenium is an integral part of selenoproteins and contains the redox-active amino acid selenocysteine in their catalytically active center [33]. A targeted deletion of tRNA gene for the selenocysteine in the osteo-chondroprogenitor cells causes many features similar to KBD, suggesting its potential importance for KBD [34]. In addition, Speckmann et al. [35] demonstrated that selenium influenced epigenetic modifications of DNA. In rat liver and colon, it was suggested that selenium deficiency resulted in less DNA methylation mainly because selenium provides the methylation donor for DNA methylation [36]. In contrast, other studies showed that high selenium exposure decreased total DNA methylation levels in the liver of rats, perhaps because of selenium inhibiting DNA methyltransferase expression [37,38]. Evidences showed that monomethylated selenium compounds are crucial for the tumor preventive effects of the selenium. The methylselenolate was identified to interact with H₂O₂ and participate in the ROS of tumor formation [39]. The methylselenol is the key intermediate contributing to selenide in the case of selenium-methylated selenocompounds [39,40].

For the first time, we have identified the selenocompound metabolism pathway associated with KBD via genetic study. It has been known selenium mainly exerts its biological effects through enzymatically active selenoproteins. The transport and storage of selenium are mainly depended on liver-derived circulating glycoprotein selenoprotein P [32], while selenoprotein S, a selenoprotein localized in the plasma membrane and endoplasmic reticulum (ER) stress, could control ER homeostasis and participate in adipogenesis [C. Y. 41]. Du et al. demonstrated that SNP of Selenoprotein S SEPS1 − 105 G > A played important roles in KBD pathogenesis [31]. Selenocompound metabolism pathway includes 17 genes and 27 compound of chemical substance, related to 4 diseases (homocystinuria, hypermethioninemia, spondyloepimetaphyseal dysplasia[SEMD], Methylcobalamin deficiency type G), among which, SEMD is the most similar disease to KBD. SEME is characterized by short, bowed lower limbs, mild brachydactyly, enlarged knee joints, and early-onset osteoarthropathy, which are nearly identical to the symptom of KBD, suggesting that KBD and SEMD share common pathways of disease that may include the selenocompound metabolism pathway. Additionally, seven genes in selenocompound metabolism pathway were identified significant association with KBD in this study, such as CTH, SEPSECS, PAPSS2. CTH, known as cystathionine gamma-lyase, is a novel osteoclast regulator. Itou et al. identified CTH may accelerate osteoclast differentiation and has a potent calcium resorption function, which may serve as a molecular switch regulating between osteoclast and osteoblast balance [42]. SEPSECS, whose official full name is Sep (O-phosphoserine) tRNA:Sec (selenocysteine) tRNA synthase, is involved in the synthesis of selenocysteine. It is identified as a DNA hypomethylation-related gene related to KBD and showed increased protein expression level in KBD cartilage tissue and chondrocytes in this study. Further functional studies were warranted to identify the interaction between SEPSECS and KBD. PAPSS2, whose official full name is 3ʹ-phosphoadenosine 5ʹ-phosphosulfate synthase 2, is also known as SK2, BCYM4, ATPSK2. Mutations of PAPSS2 can cause an autosomal recessive form of SEMD in humans [43]. Another study showed that TGF-β maintains biomechanical properties and regulates expression of PAPSS2 in mouse articular cartilage [44]. The three genes mentioned above are all up-regulated in the KBD cell models in our study. Additionally, the mRNA expression of MMP13 was significantly decreased, while SOX9 expression was increased notably after knocking out SEPSECS in C28/I2 human chondrocytes treated with T-2 toxin. MMP13 is known in the breakdown of extracellular matrix in OA development. Little et al. found that OA progression is prohibited in MMP13 knockout mice [45]. SOX9 is found in chondrocyte differentiation. The overexpression of SOX9 could alleviate the progression of human OA [46]. All these results indicated that the lower expression of SEPSECS may act as a protective role in the development of KBD. Moreover, the increased expression of SEPSECS protein, which is recognized as 54kDa band by the company, also showed 43 kDa notably in KBD chondrocytes. Same results could be found in human hepatoma cells [47]. The in vitro study in KBD cartilage highlight for the first time the importance of selenocompound metabolism pathway in the pathogenesis of KBD. All together, the identified KEGG pathways may have some potential interactions and give us novel insight in the development of KBD.

Although we identified multiple genes and pathways associated with the development of KBD, only one pathway was confirmed by further biological experiments. Therefore, more functional studies are warranted to detect the potential mechanism of these candidate genes as well as the identified pathways involved in the pathogenesis of KBD.

Conclusions

Up to now, the pathogenesis of KBD as well as the biological effects of selenium in vivo is still under debate. By trans-omics analysis of genome-wide DNA methylation profiling and mRNA expression profiling, we identified 298 common genes, 36 GO terms, and 32 KEGG pathways associated with KBD pathogenesis. Most importantly, we first observed the selenocompound metabolism pathway related to KBD via genetic study and verified it by validation tests. Our results have emphasized the dysfunction of many genes and pathways implicated in the pathogenesis of KBD and further confirmed the importance of selenium and selenocompound metabolism pathway in the development of KBD.

Acknowledgments

The study was supported by the National Natural Scientific Foundation of China [grant numbers 81472925, 81673112, 81703177]; Key projects of international cooperation among governments in scientific and technological innovation [grant number 2016YFE0119100]; and the Fundamental Research Funds for the Central Universities. We appreciate Zengtie Zhang in the Pathology Department for the guidance of processing the tissues.

Disclosure statement

The authors report no conflict of interest.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [81472925]; National Natural Science Foundation of China [81673112]; National Natural Science Foundation of China [81703177]; Key projects of international cooperation among governments in scientific and technological innovation [2016YFE0119100].