The cell viability and morphological characteristics of chondrocytes after T-2 treatment
Cell viability was firstly measured by the MTT assay to determine the treatment concentration in the preliminary experiment. Detailed information of results can be seen in Supplementary Fig. 1. According to the result, cell viability was found to be decreased gradually with the growing T-2 toxin concentration. Finally, the treatment condition of 5 ng/ml for 24 h and 72 h (with cell viability of 89.21 and 46.50%, respectively) was employed for the following experiment.
Additionally, the most control chondrocytes appeared an intact structure, with clear nuclear membrane and abundant mitochondria. But compared with the control chondrocytes, the T-2 toxin exposed-chondrocytes showed less ribosome and their mitochondria were smaller and denser with part of cristae dissolved. In addition, the engulfed necrotic cell debris can be seen in the cytoplasm. (Supplementary Fig. 2).
In the control group, the normal cell morphology could be clearly observed. However, cell necrosis, pyknosis and lysis of nuclei as well as cytoplasm with light staining could be observed in the 24 h group and 72 h group by hematoxylin and eosin (HE) staining (Supplementary Fig. 3).
The global DNA methylation level of chondrocytes
The 5-mc contents of chondrocytes treated with T-2 toxin were detected by ELISA, which could indicate the overall DNA methylation levels of chondrocytes. Compared with the control group, 5-mc content of chondrocytes in the 24 h group was dramatically increased (P < 0.05). Additionally, the 5-mc content of chondrocytes in the 24 h group was notably higher than that in the 72 h group (P < 0.05) (Supplementary Fig. 4).
The differentially expressed genes and enrichment analysis
In comparison to the results of the control group, 189 genes had altered expression levels in the chondrocytes of 24 h treatment group, of which 35 genes were up-regulated and 154 genes were down-regulated. By making a comparison between the results of the 72 h group and those of the 24 h group, 1671 genes were detected to be differentially expressed, including 1045 up-regulated genes and 626 down-regulated genes (Supplementary Table 2). The heat maps and volcano plots were exhibited in Fig. 1.
Then, the DEGs were subject to GO and KEGG enrichment analysis. For 24 h group vs control group, GO terms of DEGs primarily involved response to toxic substance (GO:0009636, P = 1.86 × 10− 10), positive regulation of cell death (GO:0010942, P = 3.09 × 10− 9) and cellular response to external stimulus (GO:0071496, P = 7.24 × 10− 8). KEGG enrichment pathways of DEGs included MAPK signaling pathway (hsa04010, P = 1.62 × 10− 7), PI3K-Akt signaling pathway (hsa04151, P = 6.76 × 10− 6) and TGF-beta signaling pathway (ko04350, P = 0.0166) (Fig. 2).
Moreover, the enriched GO terms of DEGs between the 72 h group and 24 h group mainly included extracellular matrix (GO:0031012, P = 1.66 × 10− 19), positive regulation of cell death (GO:0010942, P = 1.70 × 10− 12) and collagen trimer (GO:0005581, P = 4.90 × 10− 11). Additionally, KEGG pathway of DEGs among the two groups largely involved TNF signaling pathway (hsa04668, P = 9.77 × 10− 8), MAPK signaling pathway (hsa04010, P = 1.20 × 10− 7) and PI3K-Akt signaling pathway (hsa04151, P = 9.33 × 10− 7) (Fig. 2).
Additionally, we also compared the 72 h group and control group. To the results, 1443 genes were observed to have altered expression levels in the 72 h group compared with control group, mainly including 584 down-regulated genes and 859 up-regulated genes (seen in Supplementary Table 2 and Supplementary Fig. 5). The DEGs were enriched in various GO terms including extracellular matrix (GO:0031012, P = 3.80 × 10− 17) and cellular response to growth factor stimulus (GO:0071363, P = 6.46 × 10− 8). The KEGG pathway of DEGs between the two groups mostly included TGF-beta signaling pathway (hsa04350, P = 6.76 × 10− 4), ECM-receptor interaction (ko04512, P = 2.40 × 10− 3) and TNF signaling pathway (hsa04668, P = 1.78 × 10− 4) (Supplementary Fig. 5).
The differentially methylated genes and enrichment analysis
By comparing the results of the 24 h group with those of the control group, a total of 955 loci were detected to be differentially methylated, including 634 hypermethylated loci and 321 hypomethylated loci. All the sites were corresponding to 316 hypermethylated genes and 274 hypomethylated genes. Additionally, by a comparison between the results of the 72 h group and those of the 24 h group, there were 975 loci detected to be differentially methylated, including 723 hypermethylated loci and 252 hypomethylated loci, corresponding to 460 hypermethylated genes and 177 hypomethylated genes (Supplementary Table 2).
GO and KEGG functional enrichment analyses were conducted to the DMGs of 24 h vs control. The top 20 enriched GO terms were displayed in Fig. 3, which mainly included cell projection morphogenesis (GO:0048858, P = 2.14 × 10− 9), cAMP-mediated signaling (GO:0019933, P = 1.32 × 10− 5), actin cytoskeleton organization (GO:0030036, P = 3.09 × 10− 6). The top 20 KEGG pathways of DMGs in the 24 h group were also presented in Fig. 3, which mainly involved in cAMP signaling pathway (hsa04024, P = 0.0022), MAPK signaling pathway (hsa04010, P = 0.0056) and TGF-beta signaling pathway (hsa04350, P = 0.0112) (Fig. 3).
Furthermore, the top enriched GO terms of DMGs between the 72 h group and 24 h group mainly included response to growth factor (GO:0070848, P = 1.29 × 10− 7), regulation of transforming growth factor beta activation (GO:1901388, P = 6.76 × 10− 6) and extracellular matrix organization (GO:0030198, P = 5.13 × 10− 5). KEGG pathways of DMGs mainly included MAPK signaling pathway (ko04010, P = 3.80 × 10− 5), parathyroid hormone synthesis, secretion and action (hsa04928, P = 0.0007) and cAMP signaling pathway (hsa04024, P = 0.0239) (Fig. 3).
Additionally, for 72 h vs control, 2591 loci had significantly differential methylation levels, involving 1978 hypermethylated loci and 613 hypomethylated loci, corresponding to 954 hypermethylated genes and 466 hypomethylated genes (Supplementary Table 2). The enriched GO terms of DMGs between the 72 h group and control group mainly involved response to growth factor (GO:0070848, P = 5.62 × 10− 10), positive regulation of cell death (GO:0010942, P = 7.08 × 10− 9) and negative regulation of cellular component organization (GO:0051129, P = 2.88 × 10− 8). KEGG pathway of DMGs between the two groups included PI3K-Akt signaling pathway (hsa04151, P = 1.20 × 10− 4), Apoptosis (ko04210, P = 2.57 × 10− 4) and Ras signaling pathway (hsa04014, P = 5.37 × 10− 4) (Supplementary Fig. 5).
Integrative analysis of DNA methylation and mRNA expression
We compared the results of the DNA methylation and the gene expression via integrating the DMGs and the DEGs. Precisely, 4 DMEGs were identified in the 24 h group compared with the control group. Hypermethylated and downregulated genes mostly involved CXCL3 (log2FC = − 2.32, Δβ = 0.104) and LSAMP (log2FC = − 1.75, Δβ = 0.122); hypomethylated and upregulated gene included CCL2 (log2FC = 1.34, Δβ = − 0.138); hypomethylated and downregulated gene included SLC16A6 (log2FC = − 1.08, Δβ = − 0.101) (Supplementary Table 3).
Additionally, there were 45 DMEGs between 72 h group and 24 h group. Precisely, 9 genes were hypermethylated and downregulated such as HLA-DRB1 (log2FC = − 1.28, Δβ = 0.101), CCL2 (log2FC = − 3.74, Δβ = 0.112); 8 genes were hypomethylated and upregulated such as HDAC9 (log2FC = 2.98, Δβ = − 0.102) and MYC (log2FC = 2.25, Δβ = − 0.103); 6 genes were hypomethylated and downregulated such as COL4A2 (log2FC = − 1.40, Δβ = − 0.102) and COL5A2 (log2FC = − 2.66, Δβ = − 0.120); 24 genes were hypermethylated and upregulated such as PDE4B (log2FC = 1.43, Δβ = 0.135), MITF (log2FC = 1.85, Δβ = 0.101). The detailed information of the DMEGs could be obtained in the (Supplementary Table 3).
For 72 h vs control, 81 genes were showed to have both differentially methylated and expressed levels. Indeed, 27 genes were hypermethylated and downregulated such as CCL2 (log2FC = − 2.40, Δβ = 0.119) and IL15 (log2FC = − 1.76, Δβ = 0.102); 10 genes were hypomethylated and upregulated for instance, HDAC9 (log2FC = 2.17, Δβ = − 0.101) and PTPRO (log2FC = 1.46, Δβ = − 0.114); 34 genes were hypermethylated and upregulated including MITF (log2FC = 1.61, Δβ = 0.116) and SORCS3 (log2FC = 2.75, Δβ = 0.122); 13 genes were hypomethylated and downregulated such as CCL2 (log2FC = − 2.40, Δβ = − 0.106) and YPEL2 (log2FC = − 3.30, Δβ = − 0.132) (Supplementary Table 3).
The expression level validation of the DMEGs by qPCR
To further validate the results of DNA methylation and expression, Sequenom MassARRAY and qPCR were carried out for the DMEGs including CCL2, CXCL3, SLC16A6, HDAC9, PDE4B and HLA-DRB1.
The mRNA levels of the selected genes were tested using the qPCR. Additionally, the results of qPCR were largely consistent with the results of RNA-sequencing. The results showed that the mRNA level of CCL2 was significantly higher in the 24 h group than that in the control group (P = 0.007). However, there were no statistical differences in the expression levels of CXCL3 and SLC16A6 between the control group and the 24 h group. Moreover, compared with the 24 h group, the mRNA levels of CCL2 and HLA-DRB1 were significantly decreased in the 72 h group (P CCL2 = 0.001; P HLA-DRB1 = 0.001). Whereas, the expression levels of HDAC9 and PDE4B were statistically higher in the 72 h group than that in the 24 h group (P HDAC9 = 0.001; P PDE4B = 0.001) (Supplementary Fig. 6).
The DNA methylation level validation of the DMEGs by Sequenom MassARRAY
The DNA methylation levels of genes including CCL2, CXCL3, SLC16A6, HDAC9, PDE4B and HLA-DRB1. were further validated by the Sequenom MassARRAY system. The results showed that the DNA methylation level of CXCL3 was remarkably higher in the 24 h group than that in the control group (P < 0.05). In addition, compared with the 24 h group, DNA methylation levels of HLA-DRB1 and PDE4B were significantly increased in the 72 h group (P HLA-DRB1 < 0.05; P PDE4B < 0.05). Although there were no statistical differences in the DNA methylation levels of other genes in the T-2 treated chondrocytes, the trends of DNA methylation detected by the MassARRAY were consistent with the results of the 850 BeadChip. For example, the DNA methylation levels of CCL2 and SLC16A6 were decreased in the 24 h group compared with the control group. The DNA methylation levels of the genes (CCL2, HDAC9) were higher in the 72 h group than those in the 24 h group (Supplementary Fig. 7).