Skip to main content
Download PDF

Thyroid hormone induces ossification and terminal maturation in a preserved OA cartilage biomimetic model

Arthritis Research & Therapy volume 26, Article number: 91 (2024) Cite this article

Abstract

Objective

To characterize aspects of triiodothyronine (T3) induced chondrocyte terminal maturation within the molecular osteoarthritis pathophysiology using the previously established T3 human ex vivo osteochondral explant model.

Designs

RNA-sequencing was performed on explant cartilage obtained from OA patients (n = 8), that was cultured ex vivo with or without T3 (10 ng/ml), and main findings were validated using RT-qPCR in an independent sample set (n = 22). Enrichment analysis was used for functional clustering and comparisons with available OA patient RNA-sequencing and GWAS datasets were used to establish relevance for OA pathophysiology by linking to OA patient genomic profiles.

Results

Besides the upregulation of known hypertrophic genes EPAS1 and ANKH, T3 treatment resulted in differential expression of 247 genes with main pathways linked to extracellular matrix and ossification. CCDC80, CDON, ANKH and ATOH8 were among the genes found to consistently mark early, ongoing and terminal maturational OA processes in patients. Furthermore, among the 37 OA risk genes that were significantly affected in cartilage by T3 were COL12A1, TNC, SPARC and PAPPA.

Conclusions

RNA-sequencing results show that metabolic activation and recuperation of growth plate morphology are induced by T3 in OA chondrocytes, indicating terminal maturation is accelerated. The molecular mechanisms involved in hypertrophy were linked to all stages of OA pathophysiology and will be used to validate disease models for drug testing.

Introduction

Osteoarthritis (OA) is a prevalent and chronic joint disease causing pain and disability, particularly among the elderly. Despite intense research efforts, no disease-modifying drug has been identified, hence patients have to face an invasive joint replacement surgery at end stage disease. As a result, OA has a major impact on the quality of life. And because of the ageing world population, the healthcare burden is expected to increase even further [1, 2].

To discover pathways underlying OA etiopathophysiological pathways, large-scale genome-wide association studies have been performed, identifying key genes robustly conferring risk to OA [3]. The role of these OA risk genes in the major causal routes to OA onset highlight aberrant articular cartilage maintenance processes in response to intrinsic and environmental challenges such as ageing. To study the role of OA risk genes in OA pathophysiology further, gene expression studies of cartilage and bone have been performed. These studies, on one hand, focused on assessing allelic expression imbalances (AEI) associated with OA risk alleles, which confirmed that OA risk alleles frequently act by changing the expression of position genes [4,5,6]. On the other hand, they focused on transcriptome-wide differential expression studies in both early (non-OA versus preserved-OA joint tissues) and ongoing (preserved versus lesioned OA joint tissues) OA. The non-OA versus preserved-OA studies highlighted that age-related changes in chondrocyte gene expression are marked by the differential expression of among others, COL1A1 and MMP13, likely reflecting hypertrophic chondrocytes prone to OA onset. The preserved versus lesioned OA studies showed that, with ongoing OA pathophysiology, chondrocytes recuperate a detrimental growth plate morphology and undergo terminal maturation [7,8,9]. During this maturational process of endochondral ossification, growth plate chondrocytes start to proliferate, become hypertrophic and eventually break down cartilage matrix and undergo apoptosis to allow the transition to bone. Relevant to OA, in this respect, is that active thyroid (T3) is known to signal terminal maturation in growth plate chondrocytes, while we have previously shown by in vitro and in vivo functional follow-up studies [10,11,12,13,14] that upregulation of thyroid signalling makes cartilage prone to an OA disease state upon environmental challenges such as mechanical loading.

In the absence of reliable human age-related biomimetic OA models of the interacting osteochondral compartment, the exact OA-related cellular processes downstream of thyroid signalling in cartilage are, however, poorly defined. In that light, we have recently shown that T3 induces hypertrophy in our ex vivo human biomimetic model of OA as marked by modest upregulation of conventional hypertrophy markers matrix metalloproteinase 13 (MMP13), Endothelial PAS domain-containing protein 1 (EPAS1), alkaline phosphatase (ALPL) and collagens COL1A1 and COL10A1 in cultured osteochondral explants [15]. Nonetheless, expression of these genes in preserved cartilage samples is highly variable and not always robustly changing with ongoing pathophysiology [16, 17]. In order to characterize robust changes in cellular processes upon T3 exposure we here exploited our human ex vivo biomimetic model of OA by using the same established method and subsequently performing RNA sequencing (RNA-seq) analysis. Moreover, subsequent integration of our results with previously published genetic and transcriptome wide datasets, allowed identification of relevant markers of the OA pathophysiological process [3, 5, 16, 18,19,20].

Materials and methods

Osteochondral explant culture

Material from 20 patients who underwent total knee replacement surgery as a consequence of OA, was included in this study. Material was collected and processed anonymously under waste material protocols from local hospitals and the RAAK study protocol (ethical permission number P08.239/P19.013). No traceable data was collected, therefore no informed consent was required, in agreement with the declarations of Helsinki and Taipei.

Osteochondral explants were harvested from the macroscopically preserved areas of the condyle, washed in PBS, and cultured in chondrogenic differentiation medium as described before [15]. After 3-4 days, explants were treated for 10 days with 10 ng/ml triiodothyronine (T3, Sigma-Aldrich), matching control samples were left untreated. Medium was refreshed every 3-4 days. After the culture period, supernatant was stored at -80 degrees Celsius and a 1mm section from the middle of the explant was stored in 4% formaldehyde. Of the remaining explant, cartilage was separated from bone, snap-frozen in liquid nitrogen and stored at -80 degrees Celsius. To extract RNA, cartilage was pulverized using the Retsch MM200 bead miller under liquid nitrogen cooling with Trizol reagent. RNA was extracted using chloroform and isolated using Qiagen RNeasy Mini Kit (Qiagen GmbH, Hilden, Germany).

RNA-seq

In total, 26 samples from 12 donors were sequenced. NEBNext Ultra II Directional RNA Library Prep Kit for Illumina was used for sample preparation. The sample preparation was performed according to the protocol "NEBNext Ultra II Directional RNA Library Prep Kit for Illumina" (NEB #E7760S/L). Briefly, rRNA was depleted from total RNA using the rRNA depletion kit (NEB#E6310). After fragmentation of the mRNA, a cDNA synthesis was performed. This was used for ligation with the sequencing adapters and PCR amplification of the resulting product. The quality and yield after sample preparation was measured with the Fragment Analyzer. Clustering and DNA sequencing using the NovaSeq6000 was performed according to manufacturer's protocols. NovaSeq control software NCS v1.7 was used. The experiments were performed at GenomeScan B.V., Leiden.

Differential gene expression analysis was performed using our in-house pipeline as follows: Resulting RNA-seq reads were aligned using GSNAP [21] against HRCh38 and abundances were estimated using HTseq count v0.11.1 [22]. Only uniquely mapping reads were used and the quality was checked using MultiQC v1.7 [23]. Principal component analysis (PCA) was performed to detect outliers. Outliers were further functionally checked in a pairwise manner by using normalized counts in edgeR and assessing the effect of removing each donor-pair on the outcome. Samples from four donors were removed based on these criteria, leaving 16 paired control and T3 samples, from 8 donors, for final analysis. Differential gene expression analysis was performed using DESeq2_1.30.0 in R version 4.0.2. Variance stabilizing transformation (VST) was applied and only genes that had at least 4 reads in 50 % of samples were included in the analysis. Differentially expressed genes were found using a pair-wise approach in a general linear model with donor and treatment as contrasts. False discovery rate (FDR) correction was applied and FDR < 0.05 was considered statistically significant. Enrichment analysis was performed using clusterProfiler_3.18.0 and the function enrichGO was applied.

Quantitative Real-Time PCR

RT-qPCR validation was performed on independent cartilage RNA samples (n = 22) from 8 donors. cDNA was synthesized using Transcriptor First strand cDNA synthesis kit (Roche, Basel, Switzerland). Quantative gene expression analysis was performed for genes shown in Supplemental Table 1, using Fast Start Sybr Green Master mix (Roche Applied Science) with the QuantStudio 6 Real-Time PCR system or using the Biomark 96.96 Dynamic Arrays Fluidigm RT-qPCR platform. Gene expression was normalized using housekeeping genes SDHA, GAPDH and RSP11 by calculating the -deltaCT. Data was checked for normality using Shapiro-Wilk test and evaluation of Q-Q plots. RT-qPCR results were compared using the generalized estimating equation (GEE) with robust variance estimators to account for donor effects in SPSS v26. A linear model was used when data was normally distributed, otherwise a gamma distribution with log-link was used. A P-value < 0.05 was considered statistically significant.

Histology

Histological evaluation was performed on paraffin embedded sections (5 um thickness) using hematoxylin and eosin straining (H&E) performed by pathology department and toluidine blue staining as reported previously [15]. In short, osteochondral explants were fixed in formaldehyde and decalcified using EDTA (12.5%). Following paraffin embedding, sectioning and rehydration steps, sections were incubated with toluidine blue (Sigma-Aldrich). Histological grading according to the 14-point Mankin score was performed [24, 25].

Results

Sample characteristics

In order to increase understanding of how T3 induces hypertrophy and/or terminal maturation and how these processes mark OA pathophysiology, we performed cartilage RNA-seq analysis of control and T3 treated human osteochondral explants that were harvested from total knee replacement surgery material. The dataset consisted of 16 paired cartilage samples from 8 donors, 6 males and 2 females and with an average age of 69.9 years (range 59-76). Histological evaluation of cartilage samples after 2 weeks culture with or without T3, confirmed that relatively preserved cartilage areas were included, with an average Mankin score of 4.2 (SD = 0.8), and there was no significant difference in Mankin score between control and T3 treated cartilage.

Previously identified T3-induced genes related to hypertrophy

Prior to genome wide expression analyses between controls and T3 exposed human articular cartilage, we extracted differential expression of MMP13, EPAS1, COL1A1, COL10A1 and ALPL genes, consistently associated with hypertrophy in the past, and compared them to previously published differential expression datasets of early [20] and ongoing [16] OA pathophysiology (Table 1).

Table 1 RNA-seq results for known hypertrophic genes
Full size table

Notable in Table 1 are COL1A1 and COL10A1 that were significantly downregulated upon T3 exposure (FC = -3.0, P = 2.9 x 10-03 and FC= -3.9, P = 9.22 x 10-4, respectively) yet highly upregulated between non-OA versus OA (FC 126 and 8.6, respectively) and not significantly changed with ongoing OA pathophysiology (preserved versus lesioned OA). Moreover, we showed that EPAS1 was significantly upregulated by T3 exposure (FC = 1.6 P = 9.6 x 10-03) yet significantly downregulated with ongoing OA pathophysiology. MMP13 was very lowly expressed in preserved cartilage explants and did not change upon T3 exposure. ALPL was upregulated but did not reach significance due to a wider range in expression levels. Together these data suggest that T3 exposure is not manifested by upregulation of CO10A1, known as an early hypertrophy marker, but rather, as marked by mineralizing marker EPAS1, triggers late hypertrophy towards early maturation and mineralization.

Differential expression analysis

In our RNA sequencing dataset, 10.060 genes exceeded the gene count threshold (VST > 4 counts) in at least 50% of the samples and could be included in the differentially expression (DE) analysis. Genome wide analysis with T3 exposure resulted in 247 FDR significant DET3 genes (Supplemental Table 2). Of these DET3 genes, 182 were upregulated and 65 were downregulated (Fig. 1).

Fig. 1
figure 1

Volcano plot of differentially expressed genes after T3 stimulation. Circles represent all expressed genes in OA explant cartilage with genes significantly differentially expressed after T3-induced hypertrophy indicated in red (up) or green (down) with an FDR > 0.05. X-axis represents fold change (FC) of gene expression in T3-treated compared to control cartilage

Full size image

The most significantly upregulated gene was CCDC80, coiled-coil domain containing 80, also known as URB (FC = 9.7, FDR = 1.3 x 10-43 ), predicted to bind to glycosaminoglycans and associated with cartilage extracellular matrix development and pre-hypertrophic chondrocytes [26]. The most significantly downregulated gene was ANGPTL7 (FC = -8.1, FDR = 2.4 x 10-25), also known as CDT6, a member of the angiopoietin-like family. Notable among the differential 247 genes were COL2A1 and ACAN as important extracellular matrix (ECM) components that were significantly upregulated after T3 treatment (FC = 6.4, FDR = 8.2x10-8 and FC = 1.5, FDR = 1.2 x 10-3, respectively), likely reflecting increased metabolic activity. As shown in Supplemental Table 3 and Fig. 2, differential expression of a selection of genes was confirmed by RT-qPCR in an independent set of 22 cartilage samples from 8 donors consisting of 4 males and 4 females with an average age of 63.3 years (range 50-81).

Fig. 2
figure 2

Gene expression in control and T3 stimulated explant cartilage, figure shows connected paired samples. VST data from RNA-seq analysis (paired samples from n = 8 donors). Statistical difference determined by DESeq analysis and with application of false discovery rate (FDR) correction. Validation data shown as -ΔCT values obtained by RT-qPCR analysis (paired samples from 8 donors). Statistical difference reflects linear generalized estimation equation (GEE). ** < 0.01;**** < 0.0001; n.s = not significant

Full size image

To identify genes that are directly downstream of T3, the presence of thyroid hormone receptor response elements in the 247 differentially expressed genes was retrieved using TFLink [27] and is shown in Supplemental Table 4. The results indicate that a substantial number of the DET3 genes contain THRA, THRB and RXRA response elements, including the most upregulated genes CCDC80 and MLXIPL. Among the transcription factors in our dataset, KLF9 was identified as not containing these thyroid response elements indicating other signalling pathways may be involved in T3 signalling in cartilage. To identify potential upstream transcription factors (TF) regulating T3 induced gene expression, TF enrichment analysis was performed using DAVID [28, 29] (Supplemental Table 5). Although ZIC3 and HEN1 were most significant, notable was enrichment for NFKAPPAB65, critical in thyroid-specific gene regulation [30].

Pathway enrichment analysis

To explore the biological pathways in which the n = 247 FDR significant DET3 genes act, we performed ClusterProfiler analyses. As shown in Supplemental Table 6, there were 119 enriched GO pathways of which the majority (n = 70) was in the Biological Processes (BP) subontology. The top 25 most significant genes in the BP subontology and their most commonly associated genes are shown in Fig. 3. The two most significant pathways with the highest gene ratios in the BP subontology (both Padjusted = 7.7 x 10-16) were ECM organization (GO:0030198) and extracellular structure organization (GO:0043062) and these pathways contained the same 34 genes (Fig. 3, Supplemental Table 6).

Fig. 3
figure 3

Genes present in the top 25 pathways enriched in the GO biological process subontology. Only genes present in at least 3 pathways are shown. Color reflects fold change T3 compared to control

Full size image

Notably, in these ECM pathways were also CCDC80, as most significantly upregulated gene, and besides several collagens, other components of cartilage ECM that were upregulated were SPARC (FC = 3.7, FDR = 1.1 x 10-15) and BGN (FC = 1.6, FDR = 1.9 x 10-2). In contrast, the most significantly downregulated gene, ANGPTL7, is involved in the expression and assembly of several ECM components such as fibronectin (FN1) myocilin (MYOC), COL1A1 and versican (VCAN) and VCAN expression was also downregulated by T3 (FC = -4.8, FDR = 0.049), suggesting impairment of proper matrix assembly.

Additionally, Fig. 3 showed that ossification (GO:0001503) had, with 22 genes, the next highest gene ratio in the BP subontology (Padjusted = 8 x 10-6), highlighting pathways related to organophosphorus (GO:0046683) and endochondral bone morphogenesis (GO:0060350) (see also Supplemental Table 6). Highest differential expression after stimulation with T3 for genes within these pathways was seen for CLEC3A (FC = 14.6, FDR = 2.9 x 10-29) and SPARC, genes that are known to be expressed in hypertrophic growth-plate cartilage. Other upregulated genes in this pathway known to be related to cartilage calcification were ANKH (FC = 1.9, FDR = 2.9 x 10-4), PHOSPHO1 (FC = 6.4, FDR = 2.0 x 10-14), AKT1 (FC = 2.1, FDR = 2.0 x 10-5) and RARG (FC = 1.7, FDR = 5.0 x 10-3). Downregulation within the ossification pathway was observed for TNFRSF11B (FC = -3.3, FDR = 6.2 x 10-4), encoding osteoprotegerin (OPG), the decoy receptor for TNFSF11/RANKL and TRAIL, involved in osteoclast activation. Of final note is RSPO2, R-spondin 2 (FC = -1.8, FDR = 1.7 x 10-2), a key regulator of Wnt signalling, facilitating the differentiation of proliferating chondrocytes into hypertrophic chondrocytes by enhancing Wnt/β-catenin signalling in endochondral ossification.

Together, pathway analyses suggested that chondrocyte proliferation and calcification are present in our T3 treated explants, indicating the transition to late-stage hypertrophy is accelerated by T3.

T3 induced genes relevant for OA pathophysiology

To identify which DET3 genes we identified show overlap with OA pathophysiology, we first performed a comparison with the DE gene dataset in knee joints by Karlsson et al. who compared early OA cartilage to cartilage from healthy individuals (DEOAH) [20]. We identified 40 genes with the same direction of effect in DET3 and DEOAH (Fig. 4, Supplemental Table 7). Notable are the upregulated CLEC3A and CCDC80 as well as risk genes COL2A1 and SPARC, suggesting the same activation of extracellular matrix remodelling is induced in early OA. Next we looked at the RAAK data set of paired preserved (P) and lesioned (OA) cartilage samples of OA patients (DEOAP) [16]. We identified 39 genes with same direction of effect in DET3 and DEOAP, with most notable upregulated genes being again CCDC80, as well as KLF2 (FC = 2.6, FDR = 1.6 x 10-2) and LOXL3 (FC = 2.0, FDR = 3.8 x 10-4) that play a role in cartilage development. Among the genes upregulated in both DET3 and DEOAP were also transcription factors Hairless (HR) and Krueppel-like factor (KLF)9 (FC = 9.0 and FC = 3.8, respectively) which are known to be downstream of T3. Suggesting thyroid hormone signalling plays a significant role in OA pathophysiology.

Fig. 4
figure 4

Illustrating the number of overlapping genes between T3 cartilage RNA-seq dataset, Karlsson OAH dataset [20], and RAAK OAP dataset [16]

Full size image

When comparing all three transcriptome wide datasets, there were 10 genes in total with the same direction of effect in DET3, DEOAH and DEOAP, as such marking consistent expression changes with ongoing OA pathophysiology. As already observed, CCDC80 is one of these robustly upregulated genes, suggesting CCDC80 is a strong quantitative marker of ongoing chondrocyte hypertrophy with OA pathophysiology. Furthermore, ANKH and CDON (FCT3 = 2.7, FDR = 3.6 x 10-4) were also found to be robustly upregulated during early, ongoing and terminally matured OA. Additionally, transcription factor ATOH8 (FCT3 = -1.8, FDR = 2.8 x 10-2) was downregulated all three datasets and was also previously shown to be downregulated after mechanical loading [31], indicating it could be a robust marker for cartilage stress. We also identified the risk gene PAPPA, which is upregulated in OAP, and by T3 (FC = 3.4 and 2.7, respectively), but conversely, is significantly downregulated in OA compared to healthy cartilage (FC = -2.1).

OA risk genes responding to T3 treatment

To investigate which OA risk genes are among the DET3 genes, we compiled the genes reported as OA risk genes in genome wide association studies (GWAS) [3, 18, 19], while incorporating allelic expression imbalances of genes, to infer OA risk effects [5]. We identified 37 DET3 genes as potential OA risk genes (see Supplemental Table 7), among which were many collagens as well as other cartilage matrix associated genes such as SPARC, COMP, CHST3 and CILP. For most of these OA risk genes no AEI effect has been shown, some of the few exceptions being COL11A1 (increased risk linked to an upregulation of type XI collagen; upregulated by T3 FC = 2.2 FDR = 7.7 x 10-3), COL12A1 (increased risk linked to a downregulation of type XII collagen; downregulated by T3 FC = 0.35 FDR = 2.0 x 10-3) and TNC (increased risk linked to a downregulation of tenascin-C; downregulated by T3 FC = 0,31 FDR = 6.0 x 10-3) [4]. Other OA risk genes that fit pathological events in OA are SMO, encoding the protein smoothened, upregulated by T3 (FC = 1.9, FDR = 1.9 x 10-2), known to mediate cartilage degeneration and PAPPA (FC = 2.7, FDR = 5.9 x 10-4), an important regulator of the IGF pathway. A selection of risk genes were included in the validation (Fig. 2, Supplemental Table 3). Fig. 2 shows results for risk genes SPARC, PAPPA, TNC and PLA2G2A. Upregulation of COL27A1, and SPARC was confirmed as well as downregulation of COMP, PLA2G2A and TNC. However, upregulation of PAPPA (P-value 0.14) and COMP (P-value 0.10) was not significant in the independent validation dataset.

Discussion

This study aimed to elucidate the cellular processes related to hypertrophy and/or terminal maturation of articular chondrocytes with OA pathophysiology by T3 exposure to human ex vivo explants of preserved OA cartilage. In doing so we identified 247 DET3 genes with CCDC80 as most consistently upregulated gene. Subsequent pathway enrichment analyses highlighted that DE genes upon T3 exposure were enriched within the pathways involved in extracellular matrix production and ossification, indicating acceleration of the chondrocyte transition to late-stage hypertrophy. Notable genes related to extracellular matrix, ossification and late hypertrophy included CCDC80, SPARC, CLEC3A, ANKH and PHOSPHO1. Upon comparing T3 responsive genes to those previously identified in early OA [20] and/or ongoing OA [16] pathophysiology we identified 10 genes, among which were CCDC80, ANKH, CDON and ATOH8, that consistently mark hypertrophy in various stages (early, ongoing, and terminal) OA pathophysiology (Fig. 4, Supplemental Table 7). We further identified relevant genes that are involved in the initiation and progression of OA hypertrophic disease stages, such as TNC, SMO and PAPPA. Together, we provide valuable insights into hypertrophic events in human OA articular cartilage.

Notably, CCDC80, the gene that was most significantly differentially expressed after T3 exposure, was also among the 10 genes that consistently mark the early and ongoing OA pathophysiological process. CCDC80 is known to bind extracellular matrix components and is suggested to play a role in skeletal development [32]. It is increased during cartilage development in mice, where it is expressed in the ECM and in chondrocytes of different maturation stages, particularly in hypertrophic chondrocytes [26]. However, not much research has been done on CCDC80 in human OA. Our results indicate that CCDC80 could be a very robust and sensitive marker for hypertrophy in OA. ATOH8 for that matter, was the only gene that was downregulated by T3 and also downregulated in the OAH and OAP datasets, and this gene is also significantly downregulated after exposure to mechanical stress [31]. Transcription factor ATOH8 (Atonal BHLH transcription factor 8), acts as a regulator of chondrocyte proliferation and differentiation in endochondral bones and is found in the region of pre-hypertrophic chondrocytes [33].

Thyroid hormone signalling pathways are complex and may involve multiple mechanisms both via nuclear and non-genomic signalling pathways [34]. Transcription factor analysis of the T3 perturbed ex vivo cartilage explants indicated that a large number of the DET3 genes responded via other elements than the nuclear thyroid hormone response elements (Supplemental Table 4). Notable in this respect was the regulation of genes via KLF9 which has been directly linked to thyroid hormone receptor signalling [35] as well as the enrichment of transcription factors ZIC3 and NFKAPPAB65. As such they represent non-genomic pathways which mediate thyroid hormone signalling and have partial or no similarities to thyroid hormone receptors [36]. A weakness of our study is that we have not experimentally validated the nuclear and/or non-signalling e.g. via immunohistochemical data of the thyroid receptors.

To understand how T3 signalling leads to OA, the recently published GWAS study by Boer et. al. [3] provides a valuable set of genes for researchers in the OA research field to probe for genes consistently conferring risk to OA. Importantly, it was shown by functional genomic studies [4,5,6] that these risk alleles frequently act by changing expression of positional genes in disease relevant tissues such as cartilage and bone. Upon comparing the here identified DE genes after T3 exposure with OA risk genes, in total 37 genes were overlapping. In a few of these genes such as TNC, encoding Tenacin and COL12A1 we could subsequently compare via genome wide AEI data sets [4] the direction of effect of the risk allele in relation to the gene expression effect upon T3 exposure. For TNC and COL12A1, we showed that T3 exposure resulted in a unbeneficial downregulation of these genes e.g. preventing attempts to cartilage repair [37, 38]. Among the risk genes that are also DE after T3 exposure is PAPPA, encoding pappalysin-1. PAPPA was significantly upregulated by T3 and is also upregulated during ongoing OA pathophysiology. Pappalysin-1 cleaves IGFBP-2, IGFBP-4 and IGFBP-5 thereby increasing the activity of insulin growth factor IGF [39].

A large proportion of the extracellular matrix related genes that were upregulated by T3, for example COL2A1 and ACAN, are considered normal, healthy constituents of cartilage. Our previous study also observed upregulation of COL2A1 and ACAN in osteochondral explants after T3 treatment [15], and similar findings were observed in primary chondrocyte 3D culture [40]. Upregulation of these genes highlights that T3 causes chondrocytes to leave their maturationally arrested state, become metabolically active, and increase production of matrix proteins, possibly implicating a loosening of epigenetic control that sets them on the path towards recuperation of a growth plate-like morphology. Notably, many of the cartilage extracellular matrix related genes that are upregulated by T3 such as collagens (COL2A1, COL3A1, COL11A2, COL15A1, COL27A1), are risk genes for OA as well as being upregulated in early OA cartilage compared to healthy cartilage, indicating their upregulation is, indeed, not generally beneficial but might be a response to injury that triggers an anabolic rescue response, that causes the chondrocytes to leave their maturationally arrested state, leading to hypertrophy and the recuperation of a growth plate phenotype. Although it must be noted that there was no macroscopic effect on the cartilage noted at this stage, as shown by the Mankin score.

The relationship between thyroid hormone and OA has been long established but the processes involved that lead to cartilage degeneration remained unclear. It has previously been shown that risk variants of the DIO2 gene, the catalytic activator of thyroid hormone, increase the risk of OA by regulating local activity of T3 in the cartilage [11]. Furthermore, hyperthyroidism, hypothyroidism and systemic sensitivity to thyroid hormone have also been associated with increased risk of developing OA [41, 42]. Our study indicates that CCDC80 may prove to be a useful diagnostic marker for OA that could directly link the thyroid hormone pathway and joint disease.

In the current study we have shown that T3 induced upregulation of the hypertrophy genes EPAS1 and ANKH, yet resulted in downregulation of the well-known hypertrophy markers COL10A1 and COL1A1. Since COL10A1 and COL1A1 are upregulated during early OA (i.e. non-OA versus preserved OA cartilage) [20], and EPAS1 and ANKH particularly in ongoing OA (i.e. preserved versus lesioned OA cartilage) [16], we suggest that COL10A1 and COL1A1 are markers for early hypertrophy in cartilage while the upregulation of EPAS1 and ANKH is more reflective of late stage hypertrophy induced by T3. This is further confirmed by the study of van der Kraan and van den Berg (2012) which showed COL10A1 expression was higher in less degenerated than in more degenerated areas of the cartilage [17]. Moreover, HIF-2α (encoded by the EPAS1 gene) is essential for endochondral ossification and acts by enhancing promoter activities of MMP13 and VEGFA [43] whereas ANKH actively contributes to tissue calcification and mineralization by enhancing inorganic phosphate (Pi) concentrations and modulating osteopontin (OPN) phosphorylation [44]. Together, this confirms that T3 particularly induces late stage hypertrophy and/or terminal maturation.

Conclusion

In conclusion, this study identified key players involved in cartilage degeneration induced by T3 and provides new insights into the role of T3 in OA pathophysiology. T3-induced changes in genes related to the extracellular matrix and ossification. Upregulation of extracellular matrix genes by T3 seems to reflect the activation of OA chondrocytes, which could initiate their recuperation to a growth plate-like morphology. Our analysis shows that there are several genes with high expression that are affected by T3 and potentially reflect late-stage hypertrophic OA disease progression, chondrocyte terminal maturation and cartilage mineralization. These results will improve the identification of chondrocyte terminal maturation in patients and the development of targeted therapeutics. Our established human OA disease model of chondrocyte hypertrophy and terminal maturation will be further exploited for preclinical studies on drug target discovery and efficacy testing.

Availability of data and materials

The main data this article are available in the article and in its online supplementary material or were derived from sources in the public domain as referenced in the article. Any further underlying data and code will be shared on reasonable request to the corresponding author.

Abbreviations

AEI:

Allelic Expression Imbalances

ALPL:

Alkaline Phosphatase gene

BP:

Biological Processes

CCDC80:

Coiled-coil domain containing 80

COL1A1:

Collagen I, alpha 1

COL10A1:

Collagen X, alpha 1

DE:

Differential Expression

ECM:

Extracellular Matrix

EPAS1:

Endothelial PAS domain-containing rotein 1

FC:

Fold Change

FDR:

False Discovery Rate

FN1:

Fibronectin

GEE:

Generalized Estimating Equations

HR:

Hairless

KLF9:

Krueppel-like factor 9

MMP13:

Matrix Metalloproteinase 13

MYOC:

Myocilin

NA:

Not Available

OA:

Osteoarthritis

PCA:

Principal Component Analysis

Pi:

Inorganic phosphate

RNA-seq:

RNA sequencing

RT-qPCR:

Reverse Transcription Quantitative Polymerase Chain Reaction

T3:

Triiodothyronine

VCAN:

Versican

VST:

Variance Stabilizing Transformation

References

  1. Vitaloni M, Botto-van Bemden A, Sciortino Contreras RM, Scotton D, Bibas M, Quintero M, et al. Global management of patients with knee osteoarthritis begins with quality of life assessment: a systematic review. BMC Musculoskeletal Disorders. 2019;20(1):493.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Turkiewicz A, Petersson IF, Björk J, Hawker G, Dahlberg LE, Lohmander LS, et al. Current and future impact of osteoarthritis on health care: a population-based study with projections to year 2032. Osteoarthritis and Cartilage. 2014;22(11):1826–32.

    Article  CAS  PubMed  Google Scholar 

  3. Boer CG, Hatzikotoulas K, Southam L, Stefánsdóttir L, Zhang Y, Coutinho de Almeida R, et al. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations. Cell. 2021;184(18):4784-818.e17.

  4. Coutinho de Almeida R, Tuerlings M, Ramos Y, Den Hollander W, Suchiman E, Lakenberg N, et al. Allelic expression imbalance in articular cartilage and subchondral bone refined genome-wide association signals in osteoarthritis. Rheumatology (Oxford). 2022(1462-0332 (Electronic)).

  5. den Hollander W, Pulyakhina I, Boer C, Bomer N, van der Breggen R, Arindrarto W, et al. Annotating Transcriptional Effects of Genetic Variants in Disease-Relevant Tissue: Transcriptome-Wide Allelic Imbalance in Osteoarthritic Cartilage. Arthritis Rheumatol. 2019;71(4):561–70.

    Article  Google Scholar 

  6. Shepherd C, Reese AE, Reynard LN, Loughlin J. Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP. Arthritis Res Ther. 2019;21(1):149.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Drissi H, Zuscik M, Rosier R, O’Keefe R. Transcriptional regulation of chondrocyte maturation: potential involvement of transcription factors in OA pathogenesis. Mol Aspects Med. 2005;26(3):169–79.

    Article  CAS  PubMed  Google Scholar 

  8. Coutinho de Almeida R, Mahfouz A, Mei H, Houtman E, den Hollander W, Soul J, et al. Identification and characterization of two consistent osteoarthritis subtypes by transcriptome and clinical data integration. Rheumatology (Oxford). 2021;60(3):1166-75.

  9. Vincent TL. Of mice and men: converging on a common molecular understanding of osteoarthritis. Lancet Rheumatol. 2020;2(10):e633–45.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Bos SD, Bovee JV, Duijnisveld BJ, Raine EV, van Dalen WJ, Ramos YF, et al. Increased type II deiodinase protein in OA-affected cartilage and allelic imbalance of OA risk polymorphism rs225014 at DIO2 in human OA joint tissues. Ann Rheum Dis. 2012;71(7):1254–8.

    Article  CAS  PubMed  Google Scholar 

  11. Bomer N, den Hollander W, Ramos YF, Bos SD, van der Breggen R, Lakenberg N, et al. Underlying molecular mechanisms of DIO2 susceptibility in symptomatic osteoarthritis. Ann Rheum Dis. 2015;74(8):1571–9.

    Article  CAS  PubMed  Google Scholar 

  12. Williams GR. Thyroid hormone actions in cartilage and bone. Eur Thyroid J. 2013;2(1):3–13.

    CAS  PubMed  Google Scholar 

  13. Nagase H, Nagasawa Y, Tachida Y, Sakakibara S, Okutsu J, Suematsu N, et al. Deiodinase 2 upregulation demonstrated in osteoarthritis patients cartilage causes cartilage destruction in tissue-specific transgenic rats. Osteoarthritis Cartilage. 2013;21(3):514–23.

    Article  CAS  PubMed  Google Scholar 

  14. Bomer N, Cornelis FM, Ramos YF, den Hollander W, Storms L, van der Breggen R, et al. The effect of forced exercise on knee joints in Dio2(-/-) mice: type II iodothyronine deiodinase-deficient mice are less prone to develop OA-like cartilage damage upon excessive mechanical stress. Ann Rheum Dis. 2016;75(3):571–7.

    Article  CAS  PubMed  Google Scholar 

  15. Houtman E, van Hoolwerff M, Lakenberg N, Suchiman HED, van der Linden-van der Zwaag E, Nelissen R, et al. Human Osteochondral Explants: Reliable Biomimetic Models to Investigate Disease Mechanisms and Develop Personalized Treatments for Osteoarthritis. Rheumatol Ther. 2021;8(1):499-515.

  16. Coutinho de Almeida R, Ramos YFM, Mahfouz A, den Hollander W, Lakenberg N, Houtman E, et al. RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage. Ann Rheum Dis. 2019;78(2):270-7.

  17. van der Kraan PM, van den Berg WB. Chondrocyte hypertrophy and osteoarthritis: role in initiation and progression of cartilage degeneration? Osteoarthritis Cartilage. 2012;20(3):223–32.

    Article  PubMed  Google Scholar 

  18. Styrkarsdottir U, Lund SH, Thorleifsson G, Zink F, Stefansson OA, Sigurdsson JK, et al. Meta-analysis of Icelandic and UK data sets identifies missense variants in SMO, IL11, COL11A1 and 13 more new loci associated with osteoarthritis. Nat Genet. 2018;50(12):1681–7.

    Article  CAS  PubMed  Google Scholar 

  19. Zengini E, Hatzikotoulas K, Tachmazidou I, Steinberg J, Hartwig FP, Southam L, et al. Genome-wide analyses using UK Biobank data provide insights into the genetic architecture of osteoarthritis. Nat Genet. 2018;50(4):549–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Karlsson C, Dehne T, Lindahl A, Brittberg M, Pruss A, Sittinger M, et al. Genome-wide expression profiling reveals new candidate genes associated with osteoarthritis. Osteoarthritis Cartilage. 2010;18(4):581–92.

    Article  CAS  PubMed  Google Scholar 

  21. Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005;21(9):1859–75.

    Article  CAS  PubMed  Google Scholar 

  22. Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  PubMed  Google Scholar 

  23. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32(19):3047–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Pauli C, Whiteside R Fau - Heras FL, Heras Fl Fau - Nesic D, Nesic D Fau - Koziol J, Koziol J Fau - Grogan SP, Grogan Sp Fau - Matyas J, et al. Comparison of cartilage histopathology assessment systems on human knee joints at all stages of osteoarthritis development. (1522-9653 (Electronic)).

  25. Mankin HJ. Biochemical and metabolic aspects of osteoarthritis. (0030-5898 (Print)).

  26. Wilson R, Norris EL, Brachvogel B, Angelucci C, Zivkovic S, Gordon L, et al. Changes in the chondrocyte and extracellular matrix proteome during post-natal mouse cartilage development. Mol Cell Proteomics. 2012;11(1):M111.014159.

  27. Liska O, Bohár B, Hidas A, Korcsmáros T, Papp B, Fazekas D, et al. TFLink: an integrated gateway to access transcription factor–target gene interactions for multiple species. Database. 2022;2022:baac083.

  28. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  PubMed  Google Scholar 

  29. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 2022;50(W1):W216-w21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Giuliani C, Bucci I, Napolitano G. The Role of the Transcription Factor Nuclear Factor-kappa B in Thyroid Autoimmunity and Cancer. Front Endocrinol (Lausanne). 2018;9:471.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Houtman E, Tuerlings M, Riechelman J, Suchiman HED, van der Wal RJP, Nelissen R, et al. Elucidating mechano-pathology of osteoarthritis: transcriptome-wide differences in mechanically stressed aged human cartilage explants. Arthritis Res Ther. 2021;23(1):215.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Liu Y, Monticone M, Tonachini L, Mastrogiacomo M, Marigo V, Cancedda R, et al. URB expression in human bone marrow stromal cells and during mouse development. Biochem Biophys Res Commun. 2004;322(2):497–507.

    Article  CAS  PubMed  Google Scholar 

  33. Schroeder N, Wuelling M, Hoffmann D, Brand-Saberi B, Vortkamp A. Atoh8 acts as a regulator of chondrocyte proliferation and differentiation in endochondral bones. PLoS One. 2019;14(8): e0218230.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Brent GA. Mechanisms of thyroid hormone action. J Clin Invest. 2012;122(9):3035–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Cvoro A, Devito L, Milton FA, Noli L, Zhang A, Filippi C, et al. A thyroid hormone receptor/KLF9 axis in human hepatocytes and pluripotent stem cells. Stem Cells. 2015;33(2):416–28.

    Article  CAS  PubMed  Google Scholar 

  36. Davis PJ, Goglia F, Leonard JL. Nongenomic actions of thyroid hormone. Nat Rev Endocrinol. 2016;12(2):111–21.

    Article  CAS  PubMed  Google Scholar 

  37. Matsui Y, Hasegawa M, Iino T, Imanaka-Yoshida K, Yoshida T, Sudo A. Tenascin-C Prevents Articular Cartilage Degeneration in Murine Osteoarthritis Models. Cartilage. 2018;9(1):80–8.

    Article  CAS  PubMed  Google Scholar 

  38. Unno H, Hasegawa M, Suzuki Y, Iino T, Imanaka-Yoshida K, Yoshida T, et al. Tenascin-C promotes the repair of cartilage defects in mice. J Orthop Sci. 2020;25(2):324–30.

    Article  PubMed  Google Scholar 

  39. Laursen LS, Overgaard MT, Søe R, Boldt HB, Sottrup-Jensen L, Giudice LC, et al. Pregnancy-associated plasma protein-A (PAPP-A) cleaves insulin-like growth factor binding protein (IGFBP)-5 independent of IGF: implications for the mechanism of IGFBP-4 proteolysis by PAPP-A. FEBS Lett. 2001;504(1–2):36–40.

    Article  CAS  PubMed  Google Scholar 

  40. Lee JK, Gegg CA, Hu JC, Reddi AH, Athanasiou KA. Thyroid hormones enhance the biomechanical functionality of scaffold-free neocartilage. Arthritis Res Ther. 2015;17(1):28.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Chen S, Sun X, Zhou G, Jin J, Li Z. Association between sensitivity to thyroid hormone indices and the risk of osteoarthritis: an NHANES study. Eur J Med Res. 2022;27(1):114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Kim BY, Kim SS, Park HK, Kim HS. Assessment of the relationship between knee ultrasound and clinical symptoms in patients with thyroid dysfunction. J Int Med Res. 2020;48(1):300060519897701.

    Article  PubMed  Google Scholar 

  43. Saito T, Fukai A, Mabuchi A, Ikeda T, Yano F, Ohba S, et al. Transcriptional regulation of endochondral ossification by HIF-2alpha during skeletal growth and osteoarthritis development. Nat Med. 2010;16(6):678–86.

    Article  CAS  PubMed  Google Scholar 

  44. Boraldi F, Lofaro FD, Quaglino D. Apoptosis in the Extraosseous Calcification Process. Cells. 2021;10(1):131.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Participation of patients and doctors at Alrijne Hospital Leiderdorp is gratefully acknowledged.

Funding

The FOACUS research collaboration project (project number LSHM19036) is financially supported by Health~Holland, Top Sector Life Sciences & Health, to stimulate public-private partnerships. The funding body had no influence on study design, reporting or the decision to submit the manuscript for publication.

Author information

Authors and Affiliations

  1. Department Biomedical Data Sciences, Section of Molecular Epidemiology, LUMC, Einthovenweg 20, Postzone S05-P, 2333 ZC, Leiden, The Netherlands

    N. M. Korthagen, E. Houtman, I. Boone, R. Coutinho de Almeida, R. G. H. H Nelissen, Y. F. M. Ramos & I. Meulenbelt

  2. Galapagos BV, Willem Einthovenstraat 13, Oegstgeest, 2342 BH, The Netherlands

    K. Sivasubramaniyan & M. A Tessari

  3. Alrijne hospital, Simon Smitweg 1, Leiderdorp, 2353 GA, The Netherlands

    R. Mahdad

Authors
  1. N. M. Korthagen
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. E. Houtman
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. I. Boone
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. R. Coutinho de Almeida
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. K. Sivasubramaniyan
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. R. Mahdad
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. R. G. H. H Nelissen
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Y. F. M. Ramos
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. M. A Tessari
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. I. Meulenbelt
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

NK, IM, RC, EH, YR and IB acquired and analyzed data. IM, YR and MT obtained funding. RM, RN provided patient material. All authors contributed to conceptual design of the study as well as critical revision of the article and have approved the final version.

Corresponding author

Correspondence to I. Meulenbelt.

Ethics declarations

Ethics approval and consent to participate

Left over surgical material was collected and processed anonymously under waste material protocols and approval from local hospitals. No traceable data was collected, therefore no informed consent was required, in agreement with the declarations of Helsinki and Taipei. Material was processed under the same protocols as the RAAK study (approved by the medical ethical committee of the Leiden University Medical Center P08.239/P19.013).

Consent for publication

Not applicable

Competing interests

The authors have declared no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Korthagen, N.M., Houtman, E., Boone, I. et al. Thyroid hormone induces ossification and terminal maturation in a preserved OA cartilage biomimetic model. Arthritis Res Ther 26, 91 (2024). https://doi.org/10.1186/s13075-024-03326-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13075-024-03326-5

Keywords

Arthritis Research & Therapy

ISSN: 1478-6362

Contact us