Uncovering pathways regulating chondrogenic differentiation of CHH fibroblasts.

Mutations in the non-coding snoRNA component of mitochondrial RNA processing endoribonuclease (RMRP) are the cause of cartilage-hair hypoplasia (CHH). CHH is a rare form of metaphyseal chondrodysplasia characterized by disproportionate short stature and abnormal growth plate development. The process of chondrogenic differentiation within growth plates of long bones is vital for longitudinal bone growth. However, molecular mechanisms behind impaired skeletal development in CHH patients remain unclear. We employed a transdifferentiation model (FDC) combined with whole transcriptome analysis to investigate the chondrogenic transdifferentiation capacity of CHH fibroblasts and to examine pathway regulation in CHH cells during chondrogenic differentiation. We established that the FDC transdifferentiation model is a relevant in vitro model of chondrogenic differentiation, with an emphasis on the terminal differentiation phase, which is crucial for longitudinal bone growth. We demonstrated that CHH fibroblasts are capable of transdifferentiating into chondrocyte-like cells, and show a reduced commitment to terminal differentiation. We also found a number of key factors of BMP, FGF, and IGF-1 signalling axes to be significantly upregulated in CHH cells during the chondrogenic transdifferentiation. Our results support postulated conclusions that RMRP has pleiotropic functions and profoundly affects multiple aspects of cell fate and signalling. Our findings shed light on the consequences of pathological CHH mutations in snoRNA RMRP during chondrogenic differentiation and the relevance and roles of non-coding RNAs in genetic diseases in general.

Introduction

Cartilage-hair hypoplasia (CHH; OMIM #250250) is an autosomal recessive skeletal metaphyseal chondrodysplasia characterized by a disproportionate short stature, abnormal growth plate development, fine and sparse hair, immunodeficiencies, haematological abnormalities, and an increased risk of developing certain types of cancer [[1], [2], [3], [4], [5]]. CHH is caused by mutations in RMRP, a non-coding snoRNA component of mitochondrial RNA processing endoribonuclease (RNase MRP) [2,3,5]. To date, more than one hundred CHH-causing mutations in RMRP have been reported [2]. Mutations localized within the transcribed sequence, especially in highly conserved regions, affect RMRP stability, RNase MRP complex assembly, functions, and subcellular localization [[6], [7], [8]]. In addition, mutations in the promoter region of the RMRP gene, i.a. insertions or duplications in the region between the TATA box and the transcription initiation site, interfere with the transcription of the RMRP gene [2,6]. RNase MRP is ubiquitously present in mitochondria and nucleoli of eukaryotic cells, cleaving several RNA substrates and playing an important role in various cellular processes, such as mitochondrial DNA replication [9], pre-rRNA processing and ribosome maturation [10,11], cell cycle regulation [12,13] as well as cell's innate immune response [14]. RMRP is also a source of short regulatory RNAs that silence the expression of several genes relevant for skeletal development (PTCH2, SOX4), and levels of RMRP itself in a self-regulatory negative feedback mechanism [15,16]. However, not all cellular functions of RMRP are known.

The hallmark of CHH is short-limbed dwarfism. Longitudinal bone growth is regulated by endochondral ossification in the metaphyseal growth plates of developing bones. During this process, chondroprogenitor cells residing in the resting zone of the growth plate give rise to proliferative chondrocytes. These are characterized by the expression of SRY-Box transcription Factor 9 (SOX9), and cartilage extracellular matrix proteins such as type II collagen (COL2A1), and aggrecan (ACAN) [17,18]. Proliferative chondrocytes gradually differentiate into pre-hypertrophic chondrocytes and differentiate further into mineralizing hypertrophic chondrocytes that express the runt-related transcription factor 2 transcription factor (RUNX2), type X collagen (COL10A1), matrix-degrading enzymes such as matrix metalloproteinase 13 (MMP-13), and osteopontin (SPP1). These terminally differentiated chondrocytes then either become apoptotic or more rarely, transdifferentiate into osteoblasts [19]. Finally, the remaining cartilage extracellular matrix serves as a scaffold for endochondral bone formation [[19], [20], [21], [22], [23]]. This cellular process is tightly regulated by a variety of transcription factors, molecular signalling pathways, cytokines, and hormones [19]. However, to date, the molecular mechanisms leading to impaired skeletal development in CHH patients remain largely unclear. We previously reported that RMRP expression is regulated in an in vitro model for endochondral ossification and presented the functional link between Rmrp and impaired growth of long bone typical for CHH patients [24]. Recently, a zebrafish model of CHH was generated and confirmed that rmrp mutations disrupt chondrogenesis and the process of endochondral bone ossification in vivo [25]. The study of chondrogenic differentiation in primary CHH cells is challenging. Stem cells and primary chondrocyte cultures of CHH patients are extremely rare, require invasive sampling, and are difficult to maintain in culture and expand. A FDC (fibroblast-derived chondrocyte) model offers an alternative [26]. In this transdifferentiation model, dermal fibroblasts are transdifferentiated into chondrocyte-like cells when plated on aggrecan and stimulated with insulin-like growth factor 1 (IGF-1) and transforming growth factor β3 (TFG-β3). During FDC transdifferentiation cells lose their characteristic fibroblast morphology, form dense aggregates, and express cartilage collagens and proteoglycans [26].

Here, we used primary human CHH and control fibroblasts in a FDC transdifferentiation model and conducted transcriptomic analysis by RNA sequencing to determine the chondrogenic capacity and differential pathway regulation of CHH cells.

Materials and methods

Sample collection and dermal fibroblast culture

Dermal fibroblasts from control donors (n = 4) and CHH patients (n = 4) were isolated from skin biopsies obtained during standard diagnostic procedures at Freiburg University Hospital, Germany. Ethical approval for use of dermal fibroblasts was received from the medical ethical Institutional Review Board of Freiburg University Hospital. Informed consent was acquired from all subjects. Genotyping of CHH patients confirmed mutations within the coding region of their RMRP gene (Table 1). Dermal fibroblasts were cultured in DMEM (Dulbecco's Modified Eagle Medium)/F-12 (Invitrogen, Waltham, Massachusetts, United States) supplemented with 10% fetal calf serum (FCS) (Sigma-Aldrich, St. Louis, Missouri, United States), 1% non-essential amino acids (NEAA) (Thermo Fisher Scientific, Waltham, Massachusetts, United States), and 1% Penicillin/Streptomycin (Thermo Fisher Scientific) at 37 °C and 5% CO2.

Table 1

Table of top canonical pathways, upstream regulators, molecular and cellular functions, diseases and disorders identified by IPA for DE genes from RNA sequencing analysis for control and CHH cells undergoing FDC transdifferentiation. Pathways analysis confirmed the chondrogenic potential of control and CHH dermal fibroblasts. Significantly deregulated genes (FC ≥ 2 or FC ≤ 0.5 and FDR ≤0.05) were used as input.

Classification		Control			CHH
		Day 1/ Day 0	Day 3/ Day 0	Day 3/ Day 1	Day 1/ Day 0	Day 3/ Day 0	Day 3/ Day 1
		p-value	p-value	p-value	p-value	p-value	p-value
Canonical pathways	Role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis	3.79E-11	4.90E-12	.	1.93E-10	.	.
	Human embryonic stem cell pluripotency	2.96E-09	5.17E-11	.	.	.	.
	Osteoarthritis pathway	6.81E-09	5.92E-11	.	.	1.45E-10	.
	Role of macrophages, fibroblasts and endothelial cells in rheumatoid arthritis	.	3.98E-09	.	.	.	.
	Hepatic fibrosis/ Hepatic stellate cell activation	.	3.27E-12	9.62E-08	4.26E-11	1.12E-11	6.18E-08
	Hepatic fibrosis signalling pathway	1.81E-08	.	4.19E-06	.	.	7.07E-05
Upstream regulators	TGFβ1	1.72E-29	3.78E-28	1.00E-20	2.79E-47	1.56E-56	.
	TNFα	2.20E-29	3.70E-33	5.02E-21	4.13E-43	1.11E-55	7.13E-10
	Dexamethasone	9.60E-27	6.04E-29	1.11E-15	2.60E-33	.	4.37E-10
Molecular and cellular functions	Cellular movement	1.58E-07 - 5.05E-30	4.81E-09 - 7.34E-40	2.25E-04 - 5.26E-20	4.81E-11 - 4.45E-41	5.45E-13 - 1.48E-45	2.85E-03 - 1.57E-09
	Cellular development	1.53E-07 - 3.69E-20	5.92E-09 - 3.52E-31	2.28E-04 - 4.78E-10	.	8.60E-13 - 2.80E-32	2.44E-03 - 5.60E-09
	Cell death and survival	1.54E-07 - 1.36E-16	4.39E-09 - 1.06E-18	1.55E-04 - 7.35E-10	6.62E-11 - 2.86E-28	3.25E-13 - 2.44E-42	2.65E-03 - 1.66E-06
	Cellular function and maintenance	5.62E-08 - 3.86E-16	.	.	.	.	.
	Cellular growth and proliferation	1.53E-07 - 1.11E-15	7.74E-09 - 1.04E-19	2.11E-04 - 4.78E-10		8.60E-13 - 2.80E-32	2.44E-03 - 5.60E-09
	Cell morphology	.	7.00E-09 - 2.18E-17	1.97E-04 - 9.24E-12	.	.	2.65E-03 - 1.66E-08
	Cell cycle	.	.	.	2.37E-13 - 5.86E-26	8.56E-13 - 1.14E-35	.
Diseases and disorders	Cancer	1.35E-07 - 1.20E-44	7.38E-09 - 1.16E-32	2.29E-04 - 3.91E-22	6.65E-11 - 6.87E-78	8.56E-13 - 6.82E-68	2.89E-03 - 9.20E-13
	Organismal injury and abnormalities	1.54E-07 - 1.20E-44	7.38E-09 - 1.16E-32	2.29E-04 - 3.91E-22	6.65E-11 - 6.87E-78	8.56E-13 - 6.82E-68	2.89E-03 - 9.20E-13
	Gastrointestinal disease	1.35E-07 - 1.18E-33	7.38E-09 - 1.15E-27	1.99E-04 - 1.32E-16	4.63E-11 - 5.07E-57	3.42E-13 - 1.40E-49	.
	Dermatological diseases	9.70E-08 - 4.44E-31	3.09E-09 - 2.92E-28	.	.	.	.
	Endocrine system disorders	9.98E-08 - 1.48E-27	4.52E-09 - 4.21E-28	1.24E-04 - 1.58E-16	5.88E-11 - 2.88E-35	1.43E-13 - 5.29E-46	2.06E-03 - 6.86E-08

Chondrogenic transdifferentiation of human dermal fibroblasts

For fibroblast-derived chondrocytes (FDC) transdifferentiation, dermal fibroblasts were plated at high density (200 000 cells/well) into 24-well plates coated with 5 μg aggrecan (Sigma-Aldrich) and cultured for three days in FDC transdifferentiation medium consisting of DMEM/F12 + Glutamax (Invitrogen), 10% FCS (Sigma-Aldrich), 1% antibiotic/antimycotic (Invitrogen), 1% NEAA (Thermo Fisher Scientific), 1% ITS-A (Insulin-Transferrin-Selenium–Sodium Pyruvate; Life Technologies, Carlsbad, California, United States), 50 μg/ml ascorbic acid 2-phosphate (Sigma-Aldrich) and 1 ng/ml human recombinant TGF-β3 (Life Technologies) as previously described [26]. The medium was refreshed every other day.

RNA isolation and quality control

Total RNA was isolated at three different time points during FDC transdifferentiation (Day 0, 1, and 3; in total 24 RNA samples) using the mirVana miRNA Isolation Kit (Thermo Fisher Scientific) according to the manufacturer's protocol. RNA quality was assessed with an Agilent's 2100 Bioanalyzer and the RNA 6000 Nano Kit (Agilent, Santa Clara, California, United States). Only samples with an RNA integrity (RIN) number greater than 9 were included for further analysis.

RNA sequencing and bioinformatic analysis

RNA sequencing and bioinformatic analysis were performed by VIB Nucleomics Core (Leuven, Belgium). RNA-seq was performed using the NextSeq platform, TruSeq Library Prep Kit, NextSeq Reagent Kit (High Output-75-Cycles sequencing kit) producing single end/fr-first strand fragments. Bioconductor packages were used for the accuracy statistics (ShortRead 1.20.0), adapters trimming (cutadapt 1.2.1), removing reads shorted than 35bp, ambiguous reads, polyA-reads, low-quality reads and artefact reads (FastX 0.013; ShortRead 1.16.3). The processed reads were then aligned to the reference genome of Homo sapiens (GRCh3773) using Tophat v2.0.8b, samtools 0.1.19–44428cd, and samtools 0.1.19–44428cd. Reads that overlap with genes were counted and used to derive expression levels. Counting per gene was performed using htseq count 0.6.1p1; genes with less than one count per million were removed. Within-sample and between-sample normalization were performed using the EDASeq package from Bioconductor [27,28]. Finally, the results were expressed in Fragments Per Kilobase of gene sequence and per Million fragments of library size (FPKM). Statistical comparative analysis was performed using edgeR 3.4.0 package of Bioconductor, by fitting a negative binomial generalized linear model (GLM) [29]. The resulting p-values were corrected for multiple testing with Benjamini-Hochberg to control the false discovery rate (FDR) [30]. RNA sequencing data were uploaded to ArrayExpress under accession E-MTAB-10996.

Pathway analysis

Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., Hilden, Germany) was used for functional analysis of differentially expressed genes. Significantly deregulated genes, upregulated when FC ≥ 2.0 and FDR ≤0.05, downregulated when FC ≤ 0.5 and FDR ≤0.05, were used as input. A separate analysis was performed for each condition at the different time points, for both up and downregulated genes. IPA was performed as described previously [31].

RT-qPCR

Total RNA was reverse transcribed using random hexamers and quantitative polymerase chain reaction (qPCR) was performed using Takyon™ No Rox SYBR Master Mix blue dTTP (Eurogentec, Seraing, Belgium), 300 nM forward and reverse primers and protocol: 50 °C for 2 min, denaturation at 95 °C for 10 min, followed by 50 cycles of amplification (15 s 95 °C and 1 min 60 °C) followed by a melting curve (Bio-Rad CFX96 Real-Time PCR Detection System). Data were analyzed using Bio-Rad CFX Manager Software version 1.1 and the standard curve method (Bio-Rad, Hercules, California, United States). Relative quantification of target gene expression was normalized to a reference gene expression (primer sequences are listed in Supplementary Material, Table S1). Statistical significance was assessed by two-way ANOVA, Bonferroni post-tests, *p < 0.05, **p < 0.01, ***p < 0.001.

ELISA

To determine FGF7, BMP6, and IGF-1 protein concentrations in culture media, the hFGF7, hBMP6, and hIGF-1 Quantikine ELISA kits (R&D, DKG00, DY507, DG100B, Minneapolis, Minnesota, United States) were used according to the manufacturer's protocol. Plates were measured on a ThermoScientific Multiskan FC plate reader (ThermoScientific Multiskan FC, Waltham, MA, USA).

Results

Transcriptomic signature of transdifferentiating CHH cells

To study the process of chondrogenic differentiation in CHH cells we employed the FDC model and performed transcriptome analysis of CHH and control cells (Fig. 1). To characterize the overall transcriptomic changes provoked during FDC transdifferentiation we first determined the number of differentially expressed (DE) genes in control and CHH cells (Fig. 2). In control cells, the highest number of DE genes was observed between Day 1 and Day 0, and more than a third of these genes remained DE until Day 3 (Fig. 2B and C). A noticeably lower number of genes was DE between Day 3 and Day 0 and between Day 3 and Day 1. Overall, this showed that in control fibroblasts the FDC transdifferentiation protocol evoked the strongest transcriptional response within 24 h upon its initiation and this was followed by less strong changes in later time points. In CHH cells the timing of the largest gene expression change in response to initiation of the FDC transdifferentiation seemed to be delayed towards Day 3 (Fig. 2D and E). Next, we compared CHH and control cells to determine the impact of CHH-related mutations in the RMRP gene on the transcriptome during the chondrogenic transdifferentiation. We evaluated RMRP expression levels. RMRP expression was significantly lower in CHH cells compared to controls at each measured time point, and RMRP levels did not change due to the FDC transdifferentiation in controls nor CHH cells (Supplementary Material, Fig. S1). When we compared CHH and control cells we found DE of many genes already at baseline Day 0 (431 upregulated and 197 downregulated genes) and the number of DE genes increased even further during the course of FDC transdifferentiation (Fig. 2A). A significant portion of genes was consistently up- or downregulated in CHH cells at all time points (Fig. 2F and G). Altogether, we deduced that CHH cells showed a slower engagement into the FDC process.

Fig. 1

Chondrogenic transdifferentiation of human dermal fibroblasts - study design. Dermal fibroblasts of four CHH patients and four controls were used for the transdifferentiation experiment. Genotyping confirmed mutations within the coding region of the RMRP gene of CHH patients. Primary dermal fibroblasts were isolated from skin biopsies, seeded at high density on aggrecan coated plates, and cultured in transdifferentiation medium supplemented with ITS, ascorbic acid, and TGF-β3. Total RNA was isolated at Day 0, Day 1, and Day 3. Transcriptome analysis was performed using the NextSeq platform and Bioconductor packages. RNA sequencing expression data were then validated on transcriptional (RT-qPCR) and protein (ELISA) levels.

Fig. 2

The number of significantly up and downregulated genes for all three comparisons (control and CHH cells over time and CHH vs control cells) and pattern of gene regulation during the FDC transdifferentiation. Table summarizing the number of significantly differentially expressed (DE) genes for each comparison (A). Venn diagrams of significantly DE genes in control (B, C) and CHH cells (D, E) over time, and in CHH vs control cells comparison (F, G). Each circle represents the number of significantly DE genes for different comparisons and their overlaps. Changes in expression were calculated as fold change (FC) to control Day 0 and genes were considered significantly upregulated when FC ≥ 2.0 and FDR ≤0.05, while downregulated when FC ≤ 0.5 and FDR ≤0.05.

To acquire an insight into regulatory networks and pathway regulation during FDC transdifferentiation , DE genes were analyzed in Ingenuity Pathway Analysis (IPA) separately for all three comparisons (control, CHH, and CHH vs control). The top canonical pathways, upstream regulators, diseases and disorders, molecular and cellular functions are summarized in Table 1 (control and CHH cells over time) and Table 2, Table 3 (CHH vs control). In control as well as CHH cells (Table 1) two chondrocyte-related pathways, “The role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis” and “Osteoarthritis pathway” were among the top regulated canonical pathways. TNFα, TGF-β1 as well as dexamethasone (a glucocorticoid required for chondrogenic differentiation of human mesenchymal stem cells (hMSCs) in vitro [32]), were identified as the top upstream regulators of these pathways. “Cellular morphology”, “Development”, “Growth and proliferation”, and “Cell death and survival” were among the top enriched molecular and cellular functions in the FDC model. Analysis of DE genes in the CHH vs control comparison (Table 2) revealed “Cell cycle control and chromosomal replication”, “Cell cycle: G2/M DNA damage checkpoint control”, “Role of BRCA1 in DNA damage control” and “Nucleotide excision repair pathway” to be among the 4 top regulated canonical pathways at Day 1 and Day 3. Consistent with this notion, the top regulated molecular and cellular functions were those related to cell cycle, DNA replication, recombination and repair, cell death, and survival. Interestingly, the number of deregulated molecules associated with diseases such as cancer and gastrointestinal disorders was considerably higher in CHH cells (Table 3). This reflects clinical symptoms typically associated with CHH pathology, such as a higher risk of developing cancer and gastrointestinal problems [2].

Table 2

Table of top canonical pathways identified by IPA in CHH vs control cells comparison. Significantly deregulated genes (FC ≥ 2 or FC ≤ 0.5 and FDR ≤0.05) were used as input.

Classification		CHH vs control
		Day 0	Day 1	Day 3
		p-value	p-value	p-value
Canonical pathways	Cell cycle control of chromosomal replication	.	3.74E-11	1.57E-13
	Role of BRCA1 in DNA damage response	.	1.26E-09	7.48E-10
	Cell Cycle: G2/M DNA damage checkpoint regulation	.	1.12E-07	.
	Hereditary breast cancer signalling	.	2.83E-07	.
	Mitotic roles of polo-like kinase	.	1.45E-06	.
	Nucleotide excision repair pathway	.	.	3.11E-08

Table 3

Table of top molecular and cellular functions and diseases and disorders identified by IPA in CHH vs control comparison. Significantly deregulated genes (FC ≥ 2 or FC ≤ 0.5 and FDR ≤0.05) were used as input.

Classification		CHH vs control
		Day 0		Day 1		Day 3
		p-value	Nb of molecules	p-value	Nb of molecules	p-value	Nb of molecules
Molecular and cellular functions	Cellular movement	1.18E-04 - 3.19E-14	177	.	.	2.35E-11 - 1.75E-32	455
	Cellular development	1.11E-04 - 2.17E-08	200	3.50E-07 - 5.86E-19	456	3.67E-11 - 4.24E-31	578
	Cell death and survival	1.16E-04 - 3.64E-10	205	2.69E-07 - 1.45E-23	506	3.59E-11 - 1.07E-42	606
	Cellular growth and proliferation	1.11E-04 - 2.17E-08	182	.		3.67E-11 - 4.24E-31	572
	Cell cycle	.	.	3.67E-07 - 1.69E-29	347	3.79E-11 - 5.61E-43	390
	DNA replication, recombination and repair	.	.	3.57E-07 - 3.73E-25	289	.	.
Diseases and disorders	Cancer	1.46E-04 - 4.57E-25	523	3.15E-07 - 7.03E-61	1308	3.98E-11 - 1.04E-59	1392
	Organismal injury and abnormalities	1.46E-04 - 4.57E-25	537	3.58E-07 - 7.03E-61	1314	3.98E-11 - 1.04E-59	1405
	Gastrointestinal disease	1.46E-04 - 3.90E-22	483	3.15E-07 - 4.31E-48	1180	2.03E-11 - 8.65E-41	1243
	Dermatological diseases and conditions	1.15E-04 - 8.18E-18	360	.	.	.	.
	Endocrine system disorders	.	.	3.08E-07 - 3.12E-39	1102	1.09E-11 - 1.13E-44	1152

Overall, IPA analysis confirmed the chondrogenic potential of transdifferentiated CHH dermal fibroblasts. At the same time, overall transcriptomic changes provoked by FDC transdifferentiation suggested that engagement of CHH cells into the transdifferentiation process is delayed when compared to control fibroblasts.

Transdifferentiated CHH chondrocytes show a reduced commitment to terminal differentiation

To examine the acquired phenotype of transdifferentiated CHH chondrocytes we evaluated the expression of several chondrocyte phenotype marker genes. Starting with chondrogenic markers, the expression of SOX9, the master regulator of cartilage development and maturation [17,18], was significantly higher in CHH cells at Day 3 of FDC transdifferentiation (Fig. 3A). Next, we evaluated the expression of aggrecan and proteoglycan 4, two critical extracellular matrix (ECM) components of growth plate cartilage important for skeletal development [33,34]. We found the expression levels of proteoglycan aggrecan (ACAN) to be significantly higher in CHH cells at Day 0 and Day 3 (Fig. 3B). Proteoglycan 4 (PRG4) was induced significantly in both groups at Day 1 and Day 3 of FDC transdifferentiation, however, its expression was significantly higher in CHH cells at Day 3 (Fig. 3C). This mirrored the expression patterns of SOX9. In neither of the groups we were able to detect expression of chondrogenic COL2A1. Nevertheless, we detected pericellular COL6A3 (Fig. 3D). While COL6A3 expression increased significantly in controls already at Day 1, in CHH cells this was delayed to Day 3. This resulted in significantly lower levels of COL6A3 at Day 1 in CHH cells compared to controls. The expression level of another connecting and anchoring collagen, COL14A1, (Fig. 3E) was significantly higher in CHH cells at all measured time points.

Fig. 3

Expression levels of transcription factors and extracellular matrix components in transdifferentiating CHH and control cells confirmed successful FDC transdifferentiation and indicated hindered terminal differentiation of CHH cells. The expression levels of chondrogenic markers SOX9 (A), ACAN (B), PRG4 (C), connecting and anchoring collagens COL6A3 (D), and COL14A1 (E), and hypertrophic markers COL10A1 (F), RUNX2 (G) and GLI1 (H) were evaluated. RNA-seq data (FPKM values) are plotted as mean + standard deviation and represented as FC to control Day 0. Significance is indicated as: * CHH vs control cells, + CHH or control cells Day 1 vs Day 0; # CHH or control cells Day 3 vs Day 0; ^ CHH or control cells Day 3 vs Day 1; * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001 (same for +, #, ^).

Cartilage ECM collagen and marker of chondrocyte hypertrophy, COL10A1 [35], was the strongest upregulated collagen in control and CHH cells (Fig. 3F; control av. 270.2x, CHH av. 391.3x). We also found expression levels of the transcription factor RUNX2 to be significantly increased in both groups at Day 1 and Day 3 (Fig. 3G). Nevertheless, RUNX2 levels were significantly lower in CHH cells at Day 3. Expression of another chondrocyte maturation and osteoarthritis-associated transcription factor GLI1 [36] was induced significantly only in control cells while showing no difference in CHH cells (Fig. 3H). This resulted in significantly lower levels of GLI1 in CHH cells at Day 1 and Day 3. Chondrocyte terminal differentiation is typically accompanied by increased expression of proteolytic enzymes [37]. Therefore we examined the expression of major ECM degradation enzymes. We detected the expression of all seven MMPs previously shown to be expressed by chondrocytes (MMP- 1, 2, 3, 8, 9, 13, 14) [38]. Among these, the most strongly induced was MMP-13 and at Day 3 this induction was significant in both CHH and control cells (Fig. 4A; control 116x, CHH 14x). Nevertheless, the MMP-13 expression level was significantly lower in CHH cells at Day 3 when compared to control cells. MMP-14 expression was also induced significantly in both groups (Fig. 4B), however, its expression levels were significantly higher in CHH cells at Day 3. In relation to MMP activity regulators, we observed significantly higher expression levels of tissue inhibitor of metalloproteinases TIMP-1 in CHH cells at Day 3 (Fig. 4C), and significantly higher levels of TIMP-2 at Day 1 and Day 3 (Fig. 4D). Taken altogether, these results suggest that CHH chondrocytes have difficulties to fully engage in the transition towards their terminal differentiation in FDC culture conditions. Expression patterns of chondrocyte phenotype marker genes in CHH cells (compared to controls) are summarized in the heat map (Fig. 5A). Gene expression levels of selected genes (RMRP, SOX9, ACAN, COL10A1) were validated by RT-qPCR analysis (Supplementary Material, Figs. S2A–D). Significantly DE collagens, proteoglycans, and ECM-degrading enzymes are listed in Supplementary Material, Table S2). Chondrocyte fate is on the gene expression level regulated by several transcription factors (TFs) [39]. To better understand the observed phenotype of transdifferentiated CHH cells we examined expression levels of TFs regulating chondrocyte differentiation [40]. Only a few TFs were differentially expressed in CHH cells compared to controls at Day 0 (Table 4). As expected, the number of DE TFs increased after the initiation of FDC transdifferentiation. Compared to controls, CHH cells showed significantly increased expression levels (at Day 1 and/or Day 3) of TFs promoting chondrocyte proliferation and chondrogenic differentiation (SOX9, SOX5, SMAD3, SMAD6, FOS) [17,[41], [42], [43], [44]]; as well as TFs driving hypertrophic differentiation (NFATC1, ZBTB20, SNAI1, ATF3, FOSL2, CEBPB, and CEBPD) and terminal chondrocyte maturation (FOXO3, FOXO4) [40,45]. Expression of TFs STAT5A and STAT5B, which have been linked to FGFR3 [46], and GH-IGF-1 signalling [47] was significantly upregulated in CHH cells at Day 1 and Day 3. We also found several TFs to be significantly lower in their expression in CHH cells (at Day 1 and/or Day 3). Among these were TFs promoting hypertrophy, RUNX2 and GLI1, as well as HEY, and SP7 which drive terminal chondrocyte maturation [40]. Overall, expression levels of TFs regulating chondrocyte differentiation support observed hampered terminal differentiation of transdifferentiated CHH cells.

Fig. 4

Expression levels of matrix metalloproteinases and MMP-related factors suggested lower ECM-degradation potential of transdifferentiated CHH cells. The expression levels of MMP-13 (A), MMP-14 (B) as well as tissue inhibitors of metalloproteinases TIMP-1 (C) and TIMP-2 (D) were evaluated. RNA-seq data (FPKM values) are plotted as mean + standard deviation and represented as FC to control Day 0. Significance is indicated as: * CHH vs control cells, + CHH or control cells Day 1 vs Day 0; # CHH or control cells Day 3 vs Day 0; ^ CHH or control cells Day 3 vs Day 1; * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001 (same for +, #, ^).

Fig. 5

Heat map summarizing differential expression of chondrocyte phenotype marker genes and signalling molecules in transdifferentiating CHH cells compared to controls. Significance (CHH vs control cells) is indicated as * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001.

Table 4

Summary of transcription factors regulating chondrocyte fate and differentiation found significantly DE in CHH cells compared to controls during FDC transdifferentiation. Changes in expression were calculated as fold change (FC) to control Day 0. The FDR values are indicated as follows: the significance of CHH patients against controls for each day of differentiation: * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001.

Gene	CHH vs control cells									Function
	Day 0			Day 1			Day 3
	FC	FDR	Direction	FC	FDR	Direction	FC	FDR	Direction
SOX9	.	.	.	.	.	.	5.8	*	UP	Chondrocyte proliferation, differentiation
SOX5	.	.	.	.	.	.	4.8	*	UP	Chondrocyte proliferation, differentiation
SMAD3	.	.	.	.	.	.	1.7	*	UP	Chondrocyte proliferation, differentiation
SMAD6	1.9	*	UP	1.8	*	UP	2.3	***	UP	Chondrocyte proliferation, differentiation
NFIX	.	.	.	2.1	**	UP	1.7	*	UP	Chondrocyte differentiation, endochondral ossification
RARB	.	.	.	.	.	.	0.1	*	DOWN	Chondrocyte proliferation and ECM production
RARG	.	.	.	1.7	*	UP				Chondrocyte proliferation and ECM production
RUNX2	0.4	**	DOWN	.	.	.	0.5	*	DOWN	Chondrocyte hypertrophic maturation
NFATC1	.	.	.	2.1	*	UP	2.5	**	UP	Chondrocyte hypertrophic maturation
STAT5A	.	.	.	2.2	**	UP	3.0	***	UP	Chondrocyte hypertrophic maturation, suppress proliferation
STAT5B	.	.	.	1.8	***	UP	1.5	*	UP	Chondrocyte hypertrophic maturation, suppress proliferation
GLI1	.	.	.	0.1	**	DOWN	0.2	**	DOWN	Chondrocyte proliferation, maturation
SP7	.	.	.	0.0	***	DOWN	0.2	*	DOWN	Terminal chondrocyte maturation
ZBTB20	1.8	*	UP	.	.	.	1.8	**	UP	Chondrocyte hypertrophic, terminal maturation
SNAI1	.	.	.	.	.	.	2.6	**	UP	Chondrocyte hypertrophic maturation
FOXO3	.	.	.	2.1	***	UP	2.1	***	UP	Terminal chondrocyte maturation
FOXO4	2.0	*	UP	2.1	**	UP	2.1	***	UP	Terminal chondrocyte maturation
HIF1A	.	.	.	3.4	*	UP	.	.	.	Growth plate chondrocyte survival and differentiation
HEY	.	.	.	0.2	*	DOWN	.	.	.	Chondrocyte hypertrophic maturation
ATF4	0.5	***	DOWN	.	.	.	.	.	.	Chondrocyte proliferation, hypertrophic maturation
ATF3	.	.	.	4.3	***	UP	5.4	***	UP	Chondrocyte proliferation, hypertrophic maturation
FOS	.	.	.	2.4	*	UP	2.9	**	UP	Chondrocyte hypertrophic maturation (suppress)
FOSL2	.	.	.	1.6	*	UP	.	.	.	Chondrocyte hypertrophic maturation
CEBPB	.	.	.	2.1	**	UP	2.9	***	UP	Chondrocyte hypertrophic maturation
CEBPD	.	.	.	4.1	***	UP	3.4	***	UP	Chondrocyte hypertrophic maturation

Based on expression patterns of major chondrocyte TFs, ECM collagens, proteoglycans as well as matrix degradation enzymes and their regulators, we demonstrated that the FDC transdifferentiation model is relevant in vitro model of chondrogenic differentiation, with an emphasis on the chondrocyte terminal differentiation phase. We also showed that transdifferentiated CHH cells present with a reduced commitment to terminal differentiation.

Expression of FGF7, BMP6, and IGF-1 is higher in CHH cells

Since we confirmed transdifferentiation of fibroblasts into chondrocyte-like cells and found differences in expression of chondrogenic differentiation marker genes between control and CHH cells, we next examined the expression of growth factors known to regulate chondrogenic differentiation within the growth plate [48]. The expression of several members of the fibroblast growth factor family, notably FGF2, FGF7, and FGF9 was differentially regulated in CHH cells during the FDC transdifferentiation (Fig. 6A–C). In control cells, the expression levels of all three factors dropped significantly at Day 1 and Day 3. In CHH cells, however, the expression remained stable (FGF2 and FGF7) or decreased as well (FGF9). Overall, the expression levels of all three factors were significantly higher in CHH cells during FDC transdifferentiation (Day 1 and Day 3). Another signalling pathway crucial for chondrocyte proliferation and differentiation during endochondral bone development is the bone morphogenic protein (BMP) signalling pathway. We found BMP2 expression to be strongly upregulated at Day 1 and Day 3 in both groups (Fig. 6D; control av. 82.3x, CHH av. 50.3x). Expression of BMP6 increased significantly in CHH cells at Day 1 and then dropped almost to its basal (Day 0) level at Day 3 (Fig. 6E). However, overall, BMP6 expression levels were significantly higher in CHH cells at both time points (Day 1 and Day 3). BMPR1A, BMPR1B, and BMPR2 were all expressed in control and CHH cells, however, none of the BMP receptors was regulated during FDC transdifferentiation (data not shown). BMP endothelial cell precursor-derived regulator (BMPER) expression levels were higher in CHH cells already at baseline and remained significantly higher even as BMPER expression dropped in both groups during the transdifferentiation process (Fig. 6F). For the growth hormone/insulin-like growth factor 1 (GH/IGF-1) axis, we found expression levels of growth hormone receptor (GHR) to be significantly higher in CHH cells compared to controls at baseline Day 0 and at Day 1 of FDC transdifferentiation (Fig. 7A). Expression of IGF-1 was induced in both, CHH and control cells, in response to FDC transdifferentiation, however, IGF-1 levels were consistently significantly higher in CHH cells at all time points (Day 0, 1, and 3) (Fig. 7B). Next, we looked at expression levels of IGF binding proteins (IGFBP), which regulate the availability of IGF-1 for binding to its receptor. Levels of IGFBP2 and IGFBP4 were significantly higher in CHH cells at Day 1 and Day 3 (Fig. 7C and D). However, expression levels of PAPPA and PAPPA2, enzymes cleaving the IGFBPs, and thus releasing IGF-1, were also significantly higher in CHH cells at Day 3 of FDC transdifferentiation (Fig. 7E and F). Expression patterns of these FGF, BMP, and IGF-1 signalling molecules in CHH cells are summarized in the heat map (Fig. 5B).

Fig. 6

Multiple regulators of chondrogenesis were differentially expressed in CHH cells during FDC transdifferentiation. Expression of FGF signalling factors FGF2 (A), FGF7 (B), and FGF9 (C) was significantly higher in CHH cells during the course of FDC. For BMP signalling, expression levels of BMP2 (D), BMP6 (E), and BMPER (F) are showing differential expression in CHH and control cells. RNA-seq data (FPKM values) are plotted as mean + standard deviation and represented as FC to control Day 0. Significance is indicated as: * CHH vs control cells, + CHH or control cells Day 1 vs Day 0; # CHH or control cells Day 3 vs Day 0; ^ CHH or control cells Day 3 vs Day 1; * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001 (same for +, #, ^).

Fig. 7

Several components of the growth hormone/insulin-like growth factor 1 (GH/IGF-1) axis are differentially regulated in CHH cells. Expression of GHR (A), IGF-1 (B), IGFBP2 (C), IGFBP4 (D), PAPPA (E), and PAPPA2 (F) was higher in CHH cells at Day 1 and Day 3 indicating deregulation of the GH/IGF-1 axis. RNA-seq data (FPKM values) are plotted as mean + standard deviation and represented as FC to control Day 0. Significance is indicated as: * CHH vs control cells, + CHH or control cells Day 1 vs Day 0; # CHH or control cells Day 3 vs Day 0; ^ CHH or control cells Day 3 vs Day 1; * FDR ≤0.05; ** FDR ≤0.01, *** FDR ≤0.001 (same for +, #, ^).

IPA upstream regulator analysis (URA) identified FGF2, 7 and 9, BMP2, and 6 as well as GHR and IGF-1 as significantly regulated (p-value of overlap) upstream molecules in CHH cells (compared to controls) at Day 3 of FDC transdifferentiation (Fig. 8). Comparing the expected effects of these upstream regulators on their downstream targets with the differential gene expression data URA predicted that upstream regulators FGF2, BMP2 and BMP6 are activated in CHH cells at Day 3 of FDC.

Fig. 8

Target molecules of selected upstream regulators and their interactions depicted by IPA upstream regulator analysis. Networks were generated from the dataset of significantly deregulated genes (FC ≥ 2 or FC ≤ 0.5 and FDR ≤0.05) in CHH vs control cells comparison at Day 3 of FDC transdifferentiation. IPA upstream regulator analysis (URA) predicted the upstream molecules which might explain observed changes in gene expression. An overlap p-value was calculated using the right-tailed Fisher's Exact Test, based on the significant overlap between genes known from the literature (Ingenuity Knowledge Base) to be regulated by the upstream regulator and experimental gene expression data. P-value ≤ 0.01 is considered significant. Correlation between the expected effect of the upstream regulator on its downstream targets and differential gene expression data is indicated as activation z-score. Activation z-score ≥ 2 predicts activation, while activation z-score ≥ −2 predicts inhibition.

To validate gene expression data, we determined protein levels of FGF7, BMP6, and IGF-1 in the culture supernatant of these FDC cultures by ELISA (Fig. 9A–C). We measured significantly higher protein levels of all three growth factors in CHH cells supernatants at Day 3. This is in concert with the observed increase in gene expression levels of these factors measured by RNA sequencing and activation of FGF7 and BMP6 predicted by IPA URA. In conclusion, these results showed that several factors of FGF, BMP, and IGF-signalling pathways are regulated in CHH cells during FDC transdifferentiation.

Fig. 9

Levels of secreted factors FGF7, BMP6, and IGF-1 are following the changes observed on the gene expression level. Supernatants were collected from the same samples that were used for RNA sequencing and levels of secreted proteins FGF7 (A), BMP6 (B), and IGF-1 (C) in FDC culture supernatants were determined by ELISA. Statistical significance was assessed by two-way ANOVA, Bonferroni post-tests, *p < 0.05, **p < 0.01.

Discussion

The hallmark of CHH is short-limbed dwarfism. However, until now it is not clear how mutations in the non-coding snoRNA RMRP lead to the skeletal metaphyseal chondrodysplasia characteristic for CHH. Previously, we showed that Rmrp expression is regulated during different stages of chondrogenic differentiation in vitro, and spatiotemporally regulated in developing growth plates in vivo [24]. To study the process of chondrogenic differentiation and its regulation in primary CHH fibroblasts, we used a chondrogenic transdifferentiation model as described by French et al. [26], coupled to whole transcriptome analysis. At baseline, we noted a dominant global upregulation of gene expression in CHH cells compared to control cells. An opposite observation, a general downregulation of gene expression in CHH cells, was reported in a recently published CHH transcriptome and cell cycle study [13]. However, several relevant distinctions between the studies, such as genotype of CHH cells (3 different mutations in the transcribed region of the RMRP vs homozygous 70A > G RMRP) and passage number could explain these differences. During chondrogenic transdifferentiation, differences in gene expression between CHH and control cells became even more pronounced. Pathway analysis identified cellular processes such as cell cycle control, chromosomal replication and DNA damage checkpoints, cell death, and survival to be enriched in these genes. This is consistent with previously described functions of RMRP in cell cycle regulation, proliferation [2,[10], [11], [12], [13]], as well as telomere homeostasis [49,50], and it confirms that RMRP plays a role in a variety of cellular processes.

Looking at the gene expression of chondrocyte-specific markers such as SOX9, RUNX2, and other TFs, collagens, and MMPs we confirmed that control, as well as CHH fibroblasts, transdifferentiated into chondrocytes-like cells. This was further supported by pathway analysis, as several pathways related to cartilage and chondrocytes were enriched in both groups. In control and CHH cells, hypertrophy-related COL10A1 and MMP-13 were the strongest induced collagen and matrix-degrading enzyme. In our previous work, we also observed a strong induction of COL10A1 expression in the FDC model and found COL10A1 expression at Day 3 and Day 5 to be significantly lower in CHH cells compared to controls [24]. However, we then followed cells for 7 days and measured a limited number of genes by RT-qPCR. In the present study, we did not detect significantly lower levels of COL10A1 expression in CHH cells. However, it should be noted that it is due to a particularly high expression of COL10A1 measured in a single CHH donor that this did not reach significance. We did detect significantly lower expression levels of MMP-13 in CHH cells compared to controls. This could be explained by accompanying lower expression levels of transcription factors RUNX2 and GLI1. RUNX2 is the master regulator of chondrocyte maturation and hypertrophy [51,52], driving the expression of COL10A1 [53] and MMP-13 [54]. MMP-13 activity is regulated on multiple levels [55], one of which is represented by tissue inhibitors of metalloproteinases (TIMPs). We found increased expression of TIMP1 and TIMP2 in CHH cells, which indicates that not only was MMP-13 expressed less in CHH cells, but it may also be less active at the protein level. Mmp13-null mice are characterized by elongated growth plates and delayed terminal differentiation [56]. In humans, mutations in MMP13 cause Missouri variant of spondyloepimetaphyseal dysplasia (SEMDM; OMIM #602111), an autosomal dominant disorder which manifestations largely overlap with those observed in CHH, including metaphyseal changes, shortened limbs, and abnormal growth plate development [57]. Therefore we speculate that decreased expression and activity of MMP-13 in CHH chondrocytes contribute to impaired skeletal development of CHH patients.

Next, we looked more closely at the expression of TFs regulating chondrocyte differentiation. We observed that the expression of SOX9 and SOX5, driving chondrocyte proliferation, and chondrogenic differentiation [17,43] was significantly higher in CHH cells. This was accompanied by higher expression of negative regulators of chondrocyte hypertrophy and maturation such as SMAD3, SMAD6 and FOS [41,42,44] in CHH cells when compared to controls. SOX9 is a pivotal transcription factor for cartilage development and maturation [17,43]. It promotes chondrocyte specification and early differentiation [17,18], drives the expression of two major cartilage ECM proteins, COL2A1 and ACAN [17,58], and suppresses hypertrophy by blocking the expression of RUNX2 [59,60]. While we did not detect the expression of COL2A1, the higher expression of SOX9 in CHH cells is in line with the observed higher expression of ACAN and reduced expression of MMP-13. We also detected the expression of the chondrocyte pericellular matrix collagen COL6A3. While FDC transdifferentiation provoked induction of COL6A3 expression in both groups, the response was slightly delayed in CHH cells. Together, our results demonstrated that transdifferentiated control chondrocytes acquired a terminal differentiation-like phenotype. Employing the FDC model as an in vitro model for chondrogenic differentiation uncovered that CHH chondrocytes have difficulties to fully commit to terminal differentiation, and they are balancing between a chondrogenic and terminally differentiated phenotype.

Survival, proliferation, and differentiation of chondrocytes within the growth plate are tightly regulated by signalling pathways (IGF, FGF, BMP, IHH/PTHrP, Wnt, TGFβ) [48]. We found that several factors of BMP, FGF, and IGF-signalling pathways were affected in CHH cells during FDC transdifferentiation. FGF/FGFR signalling plays an essential role in bone development and growth. Activating mutations of FGF3R in humans cause achondroplasia, a form of short-limbed dwarfism [[61], [62], [63]]. FGF-FGFR3 signalling stimulates expression and activation of STAT proteins, including STAT5a and STAT5b, resulting in abnormal growth plate development and inhibited bone growth [46,64]. In the current study, we found higher expression levels of FGF ligands 2, 7, and 9 as well as STAT proteins STAT5A and STAT5B in transdifferentiating CHH cells. We speculate that more active FGF signalling contributes to metaphyseal chondrodysplasia typical for CHH patients. However, this needs to be specifically addressed in follow-up work.

The two most abundant BMPs in the growth plate are BMP-2 and BMP-6 [65,66]. These BMPs drive chondrocytes from a quiescent to a proliferative state, and finally towards their terminal differentiation [67]. Here we found an increase in expression of BMP2 in both control and CHH cells, which fits with the observed mature and differentiated phenotype of the FDC chondrocytes. In CHH cells, this was accompanied by a higher expression of BMP6 as well as SMAD3 and SMAD6. Both Smad6 and Smad3 have been shown to inhibit BMP-induced chondrocyte maturation and Col10a1 expression in mice [41,42], and they might attenuate BMP-induced terminal differentiation of CHH cells in the FDC model.

The GH/IGF-axis is a key endocrine mechanism regulating longitudinal bone growth [68]. Locally produced IGF-1 is particularly important as it promotes chondrocyte proliferation and hypertrophy [[68], [69], [70], [71]]. Several factors of the IGF signalling pathway were differentially expressed in CHH cells. First of all, expression of IGF-1 was massively induced in CHH cells and this translated into increased levels of secreted IGF-1 in CHH cells at Day 3 of FDC transdifferentiation. Other factors of the GH/IGF-axis, such as IGF-binding proteins IGFBP2 and IGFBP4 and proteinases PAPP-A and PAPP-A2 were less affected by FDC transdifferentiation. Nevertheless, they showed an overall higher expression in CHH cells when compared to controls. IGFBP2 and IGFBP4 are produced by differentiating and maturing chondrocytes and negatively regulate IGF-mediated proliferation and matrix synthesis [72,73]. They are regulated by proteinases PAPP-A and PAPP-A2, which cleave IGFBPs thereby releasing IGF-1 [72,74,75]. Our results show that several aspects of the IGF-1 axis are regulated in CHH cells and altogether it appears that IGF-1 signalling is more active in CHH chondrocytes. Cross-talks between BMP, FGF, and IGF signalling could also play a role here, however, relationships between these pathways are complex and context-dependent [[76], [77], [78], [79]].

Elements of all three pathways (BMP, FGF, and IGF signalling) were affected in CHH cells during FDC transdifferentiation, and we speculate that this deregulation might explain the differences in the commitment of CHH cells to terminal chondrogenic differentiation. However, due to its predominantly descriptive character, our study has limitations. Understanding how RMRP snoRNA regulates the identified pathways and elucidating the consequences of these deregulated signalling pathways calls for a deeper functional analysis that should address the biological relevance of the proposed mechanisms [24]. Despite this limitation, our study offers important novel insights into key deregulated aspects of chondrogenic differentiation in CHH cells and uncovers potential molecular mechanisms of defective growth plate development in CHH patients. For future studies, it would be beneficial to increase the number of study subjects to cover an even wider range of CHH-causing mutations. The sample size presented in this study is limited to four individuals per group, however, considering the rarity of CHH it is a reasonable number. In addition, it would be of interest to follow the FDC transdifferentiation in these cells for a longer time to investigate how the FDC differentiation of CHH cells develops at later stages.

Conclusions

In conclusion, this is the first report describing the global transcriptomic changes in CHH fibroblasts transdifferentiated into chondrocytic cells. Our data demonstrate that CHH cells present with a hampered terminal chondrogenic differentiation, accompanied by altered expression of key factors of BMP, FGF, and IGF-1 signalling axes. This supports conclusions that RMRP snoRNA has pleiotropic cellular functions and that CHH-related mutations in snoRNA RMRP profoundly affect multiple aspects of cell fate and signalling. Our findings highlight the relevance and roles of non-coding RNAs in human genetic diseases in general.

Sources of funding

This work was financially supported by a bi-lateral cooperation grant (DN 82–304) from the and the , by the Dutch Arthritis Foundation (, grants LLP14 and 17-2-401) and by a grant from Stichting de Weijerhorst (Bewegen Zonder Pijn). Mandy Peffers is funded by a Welcome Trust Clinical Intermediate Fellowship (107471/Z/15/Z).

CRediT authorship contribution statement

Alzbeta Chabronova: Investigation, Formal analysis, Writing – original draft, Writing – review & editing. Guus G.H. van den Akker: Investigation, Supervision, Formal analysis, Writing – original draft, Writing – review & editing. Mandy M.F. Meekels-Steinbusch: Conceptualization, Methodology, Resources, Investigation, Formal analysis, Writing – original draft, Writing – review & editing. Franziska Friedrich: Conceptualization, Methodology, Resources, Investigation, Formal analysis, Writing – original draft, Writing – review & editing. Andy Cremers: Resources. Don A.M. Surtel: Resources, Investigation. Mandy J. Peffers: Investigation, Formal analysis, Writing – original draft, Writing – review & editing. Lodewijk W. van Rhijn: Writing – review & editing. Ekkehart Lausch: Resources, Writing – review & editing. Bernhard Zabel: Conceptualization, Methodology, Supervision, Resources, Writing – review & editing, Funding acquisition, Project administration. Marjolein M.J. Caron: Conceptualization, Methodology, Resources, Investigation, Formal analysis, Writing – original draft, Writing – review & editing. Tim J.M. Welting: Conceptualization, Methodology, Resources, Supervision, Investigation, Formal analysis, Writing – original draft, Writing – review & editing, Funding acquisition, Project administration.

Declaration of competing interest

The authors declare no conflict of interest.