Single-cell RNA seq analysis identifies the biomarkers and differentiation of chondrocyte in human osteoarthritis

doi:10.21203/rs.3.rs-25051/v1

Download PDF

Research article

Single-cell RNA seq analysis identifies the biomarkers and differentiation of chondrocyte in human osteoarthritis

https://doi.org/10.21203/rs.3.rs-25051/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Background： Single-cell RNA sequencing (scRNA-seq) was recently adopted for exploring molecular programmes and lineage progression patterns of pathogenesis of important diseases. In this study, we use scRNA-seq to identify potential markers for chondrocytes in osteoarthritis (OA) and explore the function of different types of chondrocytes in OA.

Methods：Here we aimed to identifies the biomarkers and differentiation of chondrocyte by Single-cell RNA seq analysis. GO and KEGG analysis were used to prove the function of candidate marker genes in chondrocytes. Protein–protein interaction (PPI) network was constructed to found the hub genes in 3 types of chondrocyte respectively. We also used qRT-PCR to detect the expression level of the candidate marker genes in different types of chondrocyte.

Results： In this study, we characterized the single-cell expression profiling of 480 chondrocyte samples and found hypertrophic chondrocyte (HTC), homeostatic chondrocyte (HomC) and fibrocartilage chondrocyte (FC) respectively. The results of GO and KEGG analysis to the candidate marker genes of 3 types of chondrocyte. showed the candidate marker genes made specific function in these chondrocytes to regulate the development of OAs respectively. We further revealed the differential expression of top 10 marker genes of 3 types of chondrocyte in different type of chondrocytes. The marker genes of HTC and FC were mainly expressed in their respective cell. The marker genes of HomC did not have obviously differential expression among different types of chondrocyte. Last, we proved the key genes in each cell subset. CD44, JUN and FN1 were proved tightly related to the proliferation and differentiation of chondrocytes in OAs and could be regarded as biomarkers to estimate the development of OA.

Conclusion： Our results provide new insights into exploring the roles of different types of chondrocyte in OA. The biomarkers of chondrocyte were also valuable for estimating OA progression.

Epigenetics & Genomics

osteoarthritis

chondrocyte

single-cell RNA-seq

cell differentiation

biomarker

Osteoarthritis (OA) is the most common type of arthritis. It is a disability, delayed and degenerative disease characterized by disappearance of inflammation of articular cartilage and synovium, resulting in stiffness, swelling, pain and loss of mobility [1, 2]. An important pathological feature of OA is the degradation of articular cartilage. Cartilage is composed of chondrocytes, stroma and fibers [3]. Chondrocytes are the only cell type present in mature cartilage and are responsible for repairing damaged cartilage tissue [4]. Dysfunction of chondrocyte will result in degradation of the matrix, cell death and total loss of cartilage integrity [5]. Some reports found chondrocytes in aging or diseased cartilage become senescent with associated phenotypic changes contributing to development or progression of OA [6]. Apoptotic chondrocyte death also occurs more frequently in OA compared to normal cartilage and in OA lesional compared to OA non-lesional cartilage [7]. Owing to the important roles of chondrocytes in the pathogenesis of osteoarthritis, some researchers have tried to determine the development of osteoarthritis based on the differentiation and growth of chondrocytes [8, 9]. A few biochemical markers have been found and used for estimating growth stage of chondrocytes in OA. However, biochemical markers that can effectively predict cell subpopulations of chondrocytes remain largely unknown.

The formation of chondrocyte begins with the aggregation of mesenchymal cells [10]. Mesenchymal cells undergo a series of changes such as aggregation, migration, proliferation and fusion, and eventually differentiate into chondrocytes [11]. Chondrocytes undergo in four stages: resting stage, proliferation stage, prehypertrophy stage, and hypertrophy stage, and finally calcify into bone [12]. There exist some changes in gene and protein expression, cell morphology, cell surface markers and metabolic activity in the process of chondrocyte differentiation. Of which, the expression of some genes and proteins can not only regulate chondrogenic differentiation but also determine the stage of chondrogenic differentiation. Sox5 (L-Sox5), Sox6 and Sox9 have been used as biochemical markers of proliferative chondrocyte and cooperatively activated expression of the chondrocyte differentiation marker COL2A1 [13]. COL10A1 and ALP are markers of hypertrophic chondrocyte, when chondrocytes are transferred from proliferation stage to hypertrophy stage, the expression levels of these two genes are significantly increased [14]. Runx2 is mainly expressed in the perichondrium and prehypertrophic chondrocytes, which regulate the transformation of chondrocytes into hypertrophic chondrocytes [15]. These biomarkers play important roles in defining types of different chondrocytes. However, because of the limited number of available markers to identify the cell subtype of chondrocytes in OA, the exact process of chondrocyte differentiation in osteoarthritis is not clear.

In this study, we derived and characterized the cell features and marker genes among different cell subtype of chondrocytes using scRNA-seq profiling from chondrocytes isolated from 3 patients with OA. We conducted computational analysis to define the cell subtype of chondrocytes and proved the function and new biomarkers of cell subtype of chondrocytes in OA. Our results could shed light on differentiation process of chondrocytes in OA.

Data downloading of single-cell transcriptome profiling of chondrocyte cell samples

The single-cell transcriptome data of 480 chondrocytes were obtained from GSE104782 via the Gene Expression Omnibus (GEO) database (https://www.ncbi. nlm.nih.gov/geo/). The chondrocytes from 3 osteoarthritis patients. We filtered out poor-quality cells at first. The transcriptome data were merged into one matrix and the normalization process were conducted by using the R package of limma.

Processing of single-cell RNA-seq data

We utilized R packages of limma, Seurat, dplyr and magrittr to filtered out the poor data of matrix to transcriptome and normalized the data. We first read the data of matrix and merged the repeated genes. We used R package of Seurat to generate the object and filtered out cells with poor quality[16]. The percentage of the gene numbers, cell counts and mitochondria sequencing count were calculated. In the data, the genes with less than only 3 cells detected could be deleted. Moreover, the cells whose data including less than 50 genes could also be excluded. The proportion of mitochondria was restricted to less than 5%. We also extracted genes that have a higher coefficient of variation between cells and constructed a characteristic variance diagram.

After the quality control of data for the 3 osteoarthritis patients and found the gene symbols, we performed a PCA analysis and the top principle component for the cells was obtained [17]. We estimated the top 20 principle component by P value and the significant components were selected for t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm. For visualization, we performed the t-SNE analysis based on that selected PC dimensions with the R package of Seurat [18]. For classifying all filtered cells, we set the clustering parameter resolution to 0.5 for the function FindClusters in Seurat. The standard of expression feature of marker genes was cutting off cell samples and marker genes with logFCfilter = 0.5 and adjPvalFilter = 0.05[19]. The R package of ggplot2 was used to found the candidate marker genes in different cluster. The clusters and cell categories were also annotated by candidate marker genes on the R package of SingleR. Monocle pseudotime trajectory of cells were also constructed by the R package of monocle [20].

Culturing the ATDC5 cells and L-929 cells

ATDC5 cells were used as a vitro model for chondrocytes. The ATDC5 cells (1 × 10⁴ cells, purchased from ShZeYe, Shanghai, China) were inoculated in a cell vial containing DMEM/F12 (1:1) and 5% FBS. When the cells were more than half of the wall, added 1% of ITS (Insulin, Transferrin, Selenium). Changed the cell culture medium every 2d. The chondrocytes of 7d,14d and 21d were selected to identify the growth stage. The chondrocytes were divided into proliferative chondrocyte (proC), hypertrophic chondrocyte (HTC) and resting chondrocyte (resC) by the expression of marker genes for different phase. All cells were cultured in a constant temperature incubator at 37 ° C with 5% CO₂.

L-929 cells (MICL-002, purchased from Shanghai yanjin biological, Shanghai, China) belonging to mouse fibroblast were separated from mouse. We regarded it as FCs. All cells were cultured in MEM medium (Gibco, USA) containing 10% fetal bovine serum (FBS; Gibco), streptomycin (100 mg/ml; Gibco) and penicillin (100 units/ml; Gibco). All cells were cultured in a constant temperature incubator at 37 ° C with 5% CO₂.

RNA extraction and reverse-transcription quantitative PCR (RT-qPCR)

Total RNA from tissues was extracted and was reverse-transcribed into cDNA using RNeasy Mini kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The expression of mRNAs was determined using SYBR Premix Ex Taq TM II (Takara Bio Inc., Japan). The primers using GAPDH as internal parameter. The relative expression level was calculated using the 2^−ΔΔCt method. The experiment was repeated 3 times. The primers were listed in Supplementary Table 1

Protein–protein interaction (PPI) network analysis

PPI information of candidate marker genes in 3 clusters was acquired from Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/). Then, cytoscape software was used for construction of PPI network. The top 20 candidate marker genes with high degree of connectivity were selected as the hub genes and visualized in histogram.

Statistical analysis

Statistical analysis was performed using R packages. Bar graphs show average expression levels ± SEMs. P values were calculated using a two-tailed Student’s t-test. All statistical calculations were processed using SPSS 22.0 statistical software. The evaluation of RT-qPCR data was performed using one-way ANOVA with Tukey’s post hoc test. P < 0.05 was considered significant.

Single-cell profiling of human OA cartilage chondrocytes

To explore the genes closely related to the development of chondrocytes, we acquired 480 chondrocyte samples from 3 patients with osteoarthritis (OA1, OA2 and OA3). Each OA included 160 chondrocytes. The program Seurat were used to analyzed data. The distribution pattern chart showed the gene numbers and the sequencing count of each cell (Fig. 1a). Most cells of the 3OAs included 1000–5000 genes. The sequencing count of cells focus on the range of 950000–980000. We did not find the mitochondrial sequencing count from the data. By analyzed the correlation between nFeature and nCount, the significantly positive correlation between nFeature and nCount were found (Fig. 1b). We also identified the gene symbols with significant difference across cells (Fig. 1c), of which10 significantly differentially expressed genes (COL1A1, ZNF585B, SLC5A12, COL1A2, IBSP, SPARCL1, TMSB4X, SPP1, CCL2, COL14A1) across the cell samples have been found.

Furthermore, a Principal component analysis (PCA) was performed based on the 1500 most variable genes. The correlated genes in each component were screened. The expression profile of 30 top relative genes and the correlation analysis of top relative genes in each component were shown in Supplementary Fig. 1. The distribution of cells under the two principal components PC1 and PC2 showed an obvious dispersion, indicating that our principal component analysis could be used to mapped the cells into two dimensions (Fig. 1d). The other components were calculated with an estimated P value, and we selected the significant components for subsequent analysis (Fig. 1e). Subsequently, cells were clustered based on a graph-based clustering approach, and were visualized in dimensional space using t-distributed Stochastic Neighbor Embedding (tSNE). By tSNE algorithm we classified the cells into 6 clusters (Fig. 1f). However, subgroups of the cells could not be distinguished. The 6 clusters were all annotated as chondrocyte. Therefore, a deeper analysis was need to found the subgroups of cells.

Cell trajectory analysis and annotations

In order to find the subgroups of cells, we first revealed the tendency curve of different clusters by cell trajectory analysis. Before conducting the cell trajectory analysis, we performed differential analysis using the limma package and identified a total of 1010 candidate marker genes with | log fold change (FC) | > 0.5 and adj Pval < 0.05 (Supplementary Table 2). The top 10 differential genes between the six clusters in the heatmap plot were illustrated in Fig. 2a. Moreover, we constructed the trajectory analysis to reveal the tendency curve of the six clusters (Fig. 2b). The results showed the cells of cluster2, cluster3 and cluster5 may had significant difference between each other and may belonged to different subgroup of chondrocyte. In addition, we empirically defined the subgroup of chondrocyte in each cluster by candidate marker genes. By annotation, we found the cells of cluster2 belonged to hypertrophic chondrocytes (HTCs). The cells of cluster3 belonged to homeostatic chondrocytes (HomCs) which were defined by Ji and his partners [21]. The cells of cluster5 belonged to fibrocartilage chondrocytes (FCs). We failed to annotate the cells of cluster0, cluster1 and cluster4. By analyzed these results, we suggested the trajectory analysis combined with the annotation to candidate marker genes of tSNE was helpful to defined the types of cells between different clusters. For the top 10 candidate marker genes of HTC (cluster2), HomC (cluster3) and FC (cluster5), we illustrated the significantly differential expressions of these genes in six clusters via bubble plot (Fig. 2c-e). The plots showed most of these genes were specifically overexpressed in the cells of each cluster.

Function and pathway enrichment analysis

By cell trajectory analysis and annotation of candidate marker genes in 3 cluster (cluster2, cluster3 and cluster5), we defined HTCs, HomCs and FCs. To testify the validity of the cell annotation to 3 cluster and detect the function of candidate marker genes in the cells of 3 cluster, we carried out the GO and KEGG analysis to the candidate marker genes in the 3 cluster respectively. For the genes belonging to HTCs, most functions including oxidative stress, epithelial cell proliferation and growth factor activity were closely related to HTCs (Fig. 3a). Similar to GO analysis, the KEGG analysis to genes in HTCs also had some important pathways (ECM-receptor interaction, IL-17, TNF) (Fig. 3b). The genes of these pathways all play important roles in HTCs. The results of GO and KEGG analysis to genes in HomCs showed some functions such as unfolded protein, temperature stimulus and heat shock protein could modulate cellular homeostasis and human circadian clock rhythm (Fig. 3c-d). These functions and pathways illustrated HomCs may regulated cell cycle, development, RNA metabolism and biosynthesis. FCs have close relationship with collagen formation. The function and pathway enrichment analysis to genes in cluster 5 showed some notably functions with collagen (collagen metabolic, collagen trimer and collagen binding) (Fig. 3e-f). These results further evidence the correctness of the annotation to cells in 3 clusters. The results of GO and KEGG analysis also indicate these candidate marker genes may play significant roles during the process of proliferation and differentiation for chondrocytes.

Determining the expression level of candidate marker genes among resCs, proCs, HTCs and FCs

The GO and KEGG analysis proved the annotation of clusters was accurate. In order to confirmed the candidate marker genes indeed had aberrant expression in different subtypes of chondrocyte, we detected the expression level of top 10 candidate marker genes in 3 clusters among resCs, proCs, HTCs and FCs. We first used ATDC5 cells to separate resCs, proCs and HTCs. We separated the 3 subtypes of chondrocyte by comparing the differential expression of marker genes (resC: SOX9 and Fn1, proC: COL2A1 and Runx2, HTC: COL10A1 and COL1A1) among 7d, 14d and 21d in ATDC5 cells (Supplementary Fig. 2). For FCs, we used L-929 cells which were proximal cartilage fibroblasts of mouse as FCs.

After finishing the preparation of subtypes of chondrocyte, we detected the expression level of top 10 candidate marker genes in 3 clusters among resCs, proCs, HTCs and L-929 cells. For HTCs, the results showed the 10 genes all had higher expression level in recCs, proC, HTCs compared with L-929 cells. The genes of TNFAIP6, GAS1, SLC7A2 and TNC were markedly upregulated in HTCs (Fig. 4a-b). The results of top 10 genes in HomCs showed the 10 genes all had higher expression level in recCs, proC, HTCs. However, compared with top 10 genes in HTCs the expression level of top 10 genes in HomCs did not have obviously differential expression among resCs, proCs, HTCs (Fig. 4c-d). When examined the expression level of top 10 genes in FCs, we found the expression of most genes up-regulated in L-929 except for CFD, SFRP2 and SPARCL1. The gene expression of COL1A1 and COL1A2 also showed up-regulated in HTC compared with resC and proC (Fig. 4e-f). These results implied the most candidate marker genes may had special expression in their own subgroup.

Identification of hub genes from candidate marker genes in HTCs, HomCs and FCs

To further prove the key genes in each cell subset. Protein–protein interaction (PPI) network of candidate marker genes in 3 clusters was constructed by Cytoscape software respectively, based on STRING database (Supplementary Fig. 3). The top 20 candidate marker genes with high degree of connectivity were selected as the hub genes of each cell subset (Fig. 5a-c). These top genes might play critical roles in chondrocyte progression.

CD44, JUN and FN1 were the hub genes with highest degree of connectivity in each cell subset. In order to test the relationship between candidate marker genes and cell subset. The expression of the 3 genes in the 4 chondrocyte subsets was detected. The expression of CD44 was gradually increased with proliferation and differentiation of chondrocytes. The expression of JUN was contrary to CD44 among resC, proC and HTC. A low expression level of JUN and CD44 in L-929 was observed. The expression level of FN1 was markedly upregulated in L-929 (Fig. 5d). These results combined with previous analysis confirmed some candidate marker genes can reveal the progression of chondrocytes. CD44, JUN and FN1 may play important roles in different subtypes of chondrocyte.

Osteoarthritis (OA) is the most common disease of joints and currently lacks treatment options that modify structural pathology [22, 23]. It is usually defined as the result of an imbalance in the degradation and synthesis of the normal coupling of chondrocytes, cell matrix, and subchondral bone under the influence of mechanical and biological factors [23, 24]. It is now generally accepted that the chondrocytes are the target of these abnormal biomechanical factors, and that biochemical and genetic factors also contribute to alterations in the normal functional activities of these cells [25]. Clinically, some researchers can judge the progress of osteoarthritis patients by the degradation of cartilage tissue proteoglycan and collagen fibers [26, 27]. Several genes have also been found to play important roles in the development of chondrogenesis in osteoarthritis and can serve as markers for the development of osteoarthritis. Zhu and his partners found lncRNA FAS-AS1 could promote the degradation of extracellular matrix of cartilage in osteoarthritis by increased the expressions of MMP1 and MMP13 and decreased the expression of COL2A1 in chondrocytes [27]. These genes could also be used for predicting the internal state of chondrocytes in OA pathogenesis. However, due to the complexity of chondrocyte differentiation in osteoarthritis and lack of specific biomarker to characterize chondrocyte differentiation in osteoarthritis, the differentiation and cell subtype of chondrocytes in osteoarthritis is not known. Here, we investigated OA chondrocytes at different stages by using scRNA-seq analysis. We found 3 populations of chondrocytes (HTCs, HomCs and FCs), new markers of chondrocyte populations and signaling pathways involved in OA pathogenesis based on scRNA-seq analysis. We also found the hub genes which may play important roles in the 3 populations respectively. There results may contribute to the early diagnosis and therapeutic treatment of human OA.

In this study, we found 3 cell populations (HTCs, HomCs and FCs) in chondrocytes of 3 OAs. Hypertrophic chondrocytes are divided into prohypertrophy and hypertrophy [28]. In recent years, the research on hypertrophic chondrocytes has become a hot topic in the research of osteoarthritis. Many hypertrophic chondrocytes were found in the degenerative articular cartilage [29]. These chondrocytes have different types of collagen which have different roles in the development of osteoarthritis. It is important to study the specific role of hypertrophic chondrocytes in osteoarthritis and the genes that regulate chondrocytes. In our study, we also found HTCs were existed in chondrocytes of OAs. The GO analysis of HTCs concentrated on oxidative stress, epithelial cell proliferation and growth factor. Some researchers found SIN-1 could induce oxidative stress leads to necrosis in hypertrophic chondrocytes [30, 31]. For cell proliferation and growth factor, some growth factor such as COL10A1 and ALP were high expression in hypertrophic chondrocytes and promote the proliferation of HTCs [32, 33]. Therefore, the candidate marker genes of HTCs may have close relationships with proliferation and differentiation of HTCs in OAs. The expression analysis of top 10 candidate marker genes in different types of chondrocyte also showed TNFAIP6, GAS1, SLC7A2 and TNC were high expression in HTCs. TNFAIP6, GAS1 and TNC have been found associated with the biological function of chondrocyte or cartilage by other researchers [34–36]. The change in expression level of these genes in HTCs of OAs can be used to judge the development status of OA.

Compared with HTCs, there are few reports on the role of HomCs and FCs in osteoarthritis. FCs are new cell types that have been discovered to differentiate from chondrocytes in recent years. FCs, that was mainly found in the late stages of OA and possessed a high ratio of genes related to unfavorable OA outcomes and a capacity for vascularization [21, 37]. The GO and KEGG analysis showed FCs have close relationship with collagen formation. The macromolecular framework of articular cartilage consists of collagens (predominantly type II collagen), proteoglycans, and non-collagenous proteins [3, 38]. Of which, collagen is the matrix that provides the cartilage with its tensile stiffness and strength. It is possible for researchers to find new treatments of OA by regulating the collagen formation and proliferation of fibrocartilage chondrocytes which could secret collagen. We also found 6 candidate marker genes (S100A4, COL1A2, TMSB4X, THY1, COL14A1 and IFI27) were high expression in FCs. These genes may be used as biomarkers for distinguish FCs from other chondrocytes. HomCs can help regulate the homeostatic balance of chondrocyte metabolic activity in response to environmental signals. This type of chondrocyte may play an important role in maintaining rhythmic behaviors and daily cycles of metabolism during OA cartilage degeneration [39, 40]. The functions and pathways analysis of HomCs in our results showed regulated cell cycle, development, RNA metabolism and biosynthesis. We also found all top 10 candidate marker genes of HomCs were high expression among among resCs, proCs, HTCs. These results imply these candidate marker genes of HomCs might regulate the circadian clock system of chondrocyte in OAs.

Three hub genes were also found among these candidate genes in HTCs, HomCs and FCs respectively. CD44 is the membrane receptor for hyaluronan. CD44 is the primary receptor of retention and anchoring of proteoglycan aggregates to the chondrocyte plasma membrane [41, 42]. This receptor also paly important roles in cartilage such as maintains cartilage homeostasis and affects both chondrocyte survival pathways as well as apoptotic pathways [43, 44]. In this study we regarded CD44 as the hub gene of HTCs and found this receptor was high expression in HTCs. These results show CD44 is a biomarker of HTCs in OA and this receptor exerts more influence on HTCs than other types of chondrocyte in OA. JUN is also known as AP-1. The abnormal expression of JUN in chondrocyte can affect the differentiation of chondrocyte in growth cartilage. Motomura and his partners found specific inhibition of c-Fos/AP-1 could regulate differentiation of chondrocyte in OA which effectively prevented cartilage destruction and osteophyte formation [45]. Some researchers also found the suppression of c-Fos/AP-1 expression could regulate activity of OA chondrocyte under pathological conditions [46]. In our study, we found JUN which was the hub gene of HomCs was high expression among different types of chondrocyte except for FCs. The function of HomCs mainly focus on maintain rhythmic behaviors and daily cycles of metabolism during OA cartilage degeneration. We speculate JUN may inhibits the differentiation of chondrocyte by affecting cell cycle of chondrocyte. Fibronectin1 is a glycoprotein present in a soluble dimeric form in plasma, and in a dimeric or multimeric form at the cell surface and in extracellular matrix. Some researchers found FN1 could increase chondrocyte cell death [47]. The high expression of FNI in fibrocartilage chondrocytes in our study implied this gene could be used for testing the growth of fibrocartilage.

In conclusion, our single-cell analysis identified 3 types of chondrocyte with their candidate marker genes in OA. Our analysis identified the function and discriminative markers of HTCs, HomCs and FCs in OA. The expression level of hub genes and the growth of 3 types of chondrocyte have been proved with close relationships with the development of OA. These biomarkers provide important diagnostic and therapeutic strategies for patients with OA.

scRNA-seq: Single-cell RNA sequencing; OA: osteoarthritis; HTC: hypertrophic chondrocyte (HTC); HomC: homeostatic chondrocyte; FC: fibrocartilage chondrocyte; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: Protein–protein interaction; GEO: Gene Expression Omnibus; t-SNE: t-Distributed Stochastic Neighbor Embedding; proC: proliferative chondrocyte; HTC: hypertrophic chondrocyte; resC: resting chondrocyte; RT-qPCR: RNA extraction and reverse-transcription quantitative PCR; STRING: Search Tool for the Retrieval of Interacting Genes; PCA: Principal component analysis

Acknowledgments

Not applicable.

Authors’ contributions

Xiaolu Zhang, Nianlai Huang, and Rongfu Huang designed the algorithms. Liangming Wang and Qingfeng Ke provided bio-medical analysis. Liquan Cai and Xiaolu Zhang cultured cells and conducted the experiments of RT-qPCR. Xiaolu Zhang and Shiqiang Wu wrote the manuscript. All the authors have read and approved the manuscript.

Funding

This work was supported by grants from the Natural Science Foundation of Fujian Province (2019J01050826) to Xiaolu Zhang.

Availability of data and materials

The single-cell transcriptome data of 480 chondrocytes were obtained from GSE104782 via the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104782)

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

Not applicable.

Author details

¹Department of Orthopaedic, The Second Affiliated Hospital of Fujian Medical University, Quanzhou/362000, Fujian, China.²Department of Clinical Medicine, Fujian Medical University, Fuzhou/350000, Fujian, China.³Department of Clinical Laboratory, The Second Affiliated Hospital of Fujian Medical University, Quanzhou/362000, Fujian, China.

Hunter DJ, Bierma-Zeinstra S, Osteoarthritis. Lancet. 2019;393(10182):1745–59.
Edwards JJ, Khanna M, Jordan KP, Jordan JL, Bedson J, Dziedzic KS. Quality indicators for the primary care of osteoarthritis: a systematic review. Ann Rheum Dis. 2015;74(3):490–8.
Buckwalter JA, Mankin HJ. Articular cartilage: tissue design and chondrocyte-matrix interactions. Instr Course Lect. 1998;47:477–86.
Monfort J, Garcia-Giralt N, Lopez-Armada MJ, Monllau JC, Bonilla A, Benito P, Blanco FJ. Decreased metalloproteinase production as a response to mechanical pressure in human cartilage: a mechanism for homeostatic regulation. Arthritis Res Ther. 2006;8(5):R149.
Stiller D. [Pathology of osteoarthritis. Histology, histochemistry and ultrastructure]. Zentralbl Allg Pathol. 1982;126(5–6):425–46.
Price JS, Waters JG, Darrah C, Pennington C, Edwards DR, Donell ST, Clark IM. The role of chondrocyte senescence in osteoarthritis. Aging Cell. 2002;1(1):57–65.
Kim HA, Lee YJ, Seong SC, Choe KW, Song YW. Apoptotic chondrocyte death in human osteoarthritis. J Rheumatol. 2000;27(2):455–62.
Guilak F, Nims RJ, Dicks A, Wu CL, Meulenbelt I. Osteoarthritis as a disease of the cartilage pericellular matrix. Matrix Biol. 2018;71–72:40–50.
Singh P, Marcu KB, Goldring MB, Otero M. Phenotypic instability of chondrocytes in osteoarthritis: on a path to hypertrophy. Ann N Y Acad Sci. 2019;1442(1):17–34.
Harada K, Oida S, Sasaki S, Enomoto S. Chondrocyte-like colony formation of mesenchymal cells by dentin extracts in agarose gel culture. J Dent Res. 1990;69(9):1555–9.
Cao Y, Lv C. [Effect of Wnt6 in Proliferation, Differentiation, and Migration of Bone Marrow Mesenchymal Stem Cells]. Zhongguo Xiu Fu Chong Jian Wai Ke Za Zhi. 2015;29(1):92–6.
Hall BK, Miyake T. The membranous skeleton: the role of cell condensations in vertebrate skeletogenesis. Anat Embryol (Berl). 1992;186(2):107–24.
Lefebvre V, Li P, de Crombrugghe B. A new long form of Sox5 (L-Sox5), Sox6 and Sox9 are coexpressed in chondrogenesis and cooperatively activate the type II collagen gene. EMBO J. 1998;17(19):5718–33.
Nishida T, Kubota S, Aoyama E, Takigawa M. Impaired glycolytic metabolism causes chondrocyte hypertrophy-like changes via promotion of phospho-Smad1/5/8 translocation into nucleus. Osteoarthritis Cartilage. 2013;21(5):700–9.
Enomoto H, Enomoto-Iwamoto M, Iwamoto M, Nomura S, Himeno M, Kitamura Y, Kishimoto T, Komori T. Cbfa1 is a positive regulatory factor in chondrocyte maturation. J Biol Chem. 2000;275(12):8695–702.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20.
Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, Tirosh I, Bialas AR, Kamitaki N, Martersteck EM, et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell. 2015;161(5):1202–14.
Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33(5):495–502.
Pont F, Tosolini M, Fournie JJ. Single-Cell Signature Explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets. Nucleic Acids Res. 2019;47(21):e133.
Hou R, Denisenko E, Forrest ARR. scMatch: a single-cell gene expression profile annotation tool using reference datasets. Bioinformatics. 2019;35(22):4688–95.
Ji Q, Zheng Y, Zhang G, Hu Y, Fan X, Hou Y, Wen L, Li L, Xu Y, Wang Y, et al. Single-cell RNA-seq analysis reveals the progression of human osteoarthritis. Ann Rheum Dis. 2019;78(1):100–10.
Eckstein F, Le Graverand MP. Plain radiography or magnetic resonance imaging (MRI): Which is better in assessing outcome in clinical trials of disease-modifying osteoarthritis drugs? Summary of a debate held at the World Congress of Osteoarthritis. 2014; Semin Arthritis Rheum 2015, 45(3):251–256.
Han J, Guo X, Tan W, Zhang F, Liu J, Wang W, Xu P, Lammi MJ. The expression of p-ATF2 involved in the chondeocytes apoptosis of an endemic osteoarthritis, Kashin-Beck disease. BMC Musculoskelet Disord. 2013;14:209.
Karsdal MA, Bay-Jensen AC, Lories RJ, Abramson S, Spector T, Pastoureau P, Christiansen C, Attur M, Henriksen K, Goldring SR, et al. The coupling of bone and cartilage turnover in osteoarthritis: opportunities for bone antiresorptives and anabolics as potential treatments? Ann Rheum Dis. 2014;73(2):336–48.
Goldring MB. The role of the chondrocyte in osteoarthritis. Arthritis Rheum. 2000;43(9):1916–26.
Bay-Jensen AC, Wichuk S, Byrjalsen I, Leeming DJ, Morency N, Christiansen C, Karsdal MA, Maksymowych WP. Circulating protein fragments of cartilage and connective tissue degradation are diagnostic and prognostic markers of rheumatoid arthritis and ankylosing spondylitis. PLoS One. 2013;8(1):e54504.
Zhu JK, He TD, Wei ZX, Wang YM. LncRNA FAS-AS1 promotes the degradation of extracellular matrix of cartilage in osteoarthritis. Eur Rev Med Pharmacol Sci. 2018;22(10):2966–72.
Zhang F, Xu L, Xu L, Xu Q, Li D, Yang Y, Karsenty G, Chen CD. JMJD3 promotes chondrocyte proliferation and hypertrophy during endochondral bone formation in mice. J Mol Cell Biol. 2015;7(1):23–34.
Akkiraju H, Nohe A. Role of Chondrocytes in Cartilage Formation, Progression of Osteoarthritis and Cartilage Regeneration. J Dev Biol. 2015;3(4):177–92.
He Y, Zhang Y, Zhang D, Zhang M, Wang M, Jiang Z, Otero M, Chen J. 3-morpholinosydnonimine (SIN-1)-induced oxidative stress leads to necrosis in hypertrophic chondrocytes in vitro. Biomed Pharmacother. 2018;106:1696–704.
He Y, Yao W, Zhang M, Zhang Y, Zhang D, Jiang Z, Ma T, Sun J, Shao M, Chen J. Changes in osteogenic gene expression in hypertrophic chondrocytes induced by SIN-1. Exp Ther Med. 2018;16(2):609–18.
Chen J, Chen F, Bian H, Wang Q, Zhang X, Sun L, Gu J, Lu Y, Zheng Q. Hypertrophic chondrocyte-specific Col10a1 controlling elements in Cre recombinase transgenic studies. Am J Transl Res. 2019;11(10):6672–9.
Diederichs S, Klampfleuthner FAM, Moradi B, Richter W. Chondral Differentiation of Induced Pluripotent Stem Cells Without Progression Into the Endochondral Pathway. Front Cell Dev Biol. 2019;7:270.
Wang S, Wang W, Zhao J, Zhang F, He S, Guo X. To screen the effective software for analysing gene interactions from Kashin-Beck disease genome profiling pathway and network, according to the tool of GeneMANIA. Int J Data Min Bioinform. 2015;12(1):100–14.
Lee KK, Leung AK, Tang MK, Cai DQ, Schneider C, Brancolini C, Chow PH. Functions of the growth arrest specific 1 gene in the development of the mouse embryo. Dev Biol. 2001;234(1):188–203.
Chen H, Yang J, Tan Z. Upregulation of microRNA-9-5p inhibits apoptosis of chondrocytes through downregulating Tnc in mice with osteoarthritis following tibial plateau fracture. J Cell Physiol. 2019;234(12):23326–36.
Huang L, Li M, Li H, Yang C, Cai X. Study of differential properties of fibrochondrocytes and hyaline chondrocytes in growing rabbits. Br J Oral Maxillofac Surg. 2015;53(2):187–93.
Pearle AD, Warren RF, Rodeo SA. Basic science of articular cartilage and osteoarthritis. Clin Sports Med. 2005;24(1):1–12.
Guo B, Yang N, Borysiewicz E, Dudek M, Williams JL, Li J, Maywood ES, Adamson A, Hastings MH, Bateman JF, et al. Catabolic cytokines disrupt the circadian clock and the expression of clock-controlled genes in cartilage via an NFsmall ka, CyrillicB-dependent pathway. Osteoarthritis Cartilage. 2015;23(11):1981–8.
Berenbaum F, Meng QJ. The brain-joint axis in osteoarthritis: nerves, circadian clocks and beyond. Nat Rev Rheumatol. 2016;12(9):508–16.
Knudson CB. Hyaluronan receptor-directed assembly of chondrocyte pericellular matrix. J Cell Biol. 1993;120(3):825–34.
Chow G, Knudson CB, Homandberg G, Knudson W. Increased expression of CD44 in bovine articular chondrocytes by catabolic cellular mediators. J Biol Chem. 1995;270(46):27734–41.
Chow G, Nietfeld JJ, Knudson CB, Knudson W. Antisense inhibition of chondrocyte CD44 expression leading to cartilage chondrolysis. Arthritis Rheum. 1998;41(8):1411–9.
Knudson CB, Knudson W. Hyaluronan and CD44: modulators of chondrocyte metabolism. Clin Orthop Relat Res. 2004;(427 Suppl):S152-162.
Motomura H, Seki S, Shiozawa S, Aikawa Y, Nogami M, Kimura T. A selective c-Fos/AP-1 inhibitor prevents cartilage destruction and subsequent osteophyte formation. Biochem Biophys Res Commun. 2018;497(2):756–61.
Khan NM, Haseeb A, Ansari MY, Haqqi TM. A wogonin-rich-fraction of Scutellaria baicalensis root extract exerts chondroprotective effects by suppressing IL-1beta-induced activation of AP-1 in human OA chondrocytes. Sci Rep. 2017;7:43789.
Williams A, Smith JR, Allaway D, Harris P, Liddell S, Mobasheri A. Carprofen inhibits the release of matrix metalloproteinases 1, 3, and 13 in the secretome of an explant model of articular cartilage stimulated with interleukin 1beta. Arthritis Res Ther. 2013;15(6):R223.

Supplementary Figure 1 The top 4 components from PCA procedure were shown and we exhibited the correlated genes in each component.

a: Cluster analysis across each component. The colors ranging from purple to golden yellow represent the expression levels of correlated genes from low to high; b: Correlation analysis of top relative genes in each component.

Supplementary Figure 2 The expression level of marker genes of resC, proC and HTC in ATDC5 cells during different stage.

SOX9 and Fn1 were the marker genes of resC, COL2A1 and Runx2 were the markers genes of proC, COL10A1 and COL1A1 were the marker genes of HTC.

Supplementary Figure 3 Protein–protein interaction (PPI) network of candidate marker genes in 3 clusters.

a-c: Protein–protein interaction (PPI) network of candidate marker genes in HTCs (a), HomCs (b) and FCs (c).

Download PDF

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Single-cell RNA seq analysis identifies the biomarkers and differentiation of chondrocyte in human osteoarthritis

Status:

Version 1

Abstract

Figures

Background

Methods

Data downloading of single-cell transcriptome profiling of chondrocyte cell samples

Processing of single-cell RNA-seq data

Culturing the ATDC5 cells and L-929 cells

RNA extraction and reverse-transcription quantitative PCR (RT-qPCR)

Protein–protein interaction (PPI) network analysis

Statistical analysis

Results

Single-cell profiling of human OA cartilage chondrocytes

Cell trajectory analysis and annotations

Function and pathway enrichment analysis

Determining the expression level of candidate marker genes among resCs, proCs, HTCs and FCs

Identification of hub genes from candidate marker genes in HTCs, HomCs and FCs

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Figures

Supplementary Files

Status:

Version 1