Excess FGFR3 signaling in achondroplasia disrupts turnover of resting zone chondrocytes via CREB signaling

TL;DR

This study shows that overactive FGFR3 signaling in achondroplasia disrupts cartilage cell function and suggests CREB as a potential treatment target.

Contribution

The study identifies CREB signaling as a novel mediator of FGFR3 effects in achondroplasia.

Findings

Excess FGFR3 signaling disrupts resting zone chondrocyte turnover and stem cell-like behavior.

CREB signaling contributes to dwarfism by impairing resting zone chondrocyte properties.

CREB inhibitor 666-15 restores growth plate pathology and bone length in a mouse model.

Abstract

Achondroplasia, associated with gain-of-function mutations in FGFR3, causes growth plate cartilage dysfunction, resulting in short-limb dwarfism. However, its precise molecular and cellular mechanisms remain unclear. To address this, we aimed to generate knock-in mice (Fgfr3Ach) harboring the achondroplasia mutation (p.Gly380Arg). In addition to previously reported abnormalities, we observe an expansion of the resting zone. EdU labeling and lineage tracing analyses indicate that disruption of turnover and impairment of stem cell-like behavior of resting zone chondrocytes results in accumulation of cells in the resting zone. Single-cell RNA-seq and immunohistochemical analysis identify a cell cluster that corresponds to the expanded resting zone. Pathway analysis and functional experiments reveal that CREB disrupts stem cell-like properties in resting zone chondrocytes and contributes to dwarfism. Administration of CREB inhibitor 666-15 restores growth plate pathology and bone length. These findings demonstrate that excess FGFR3 signaling disrupts resting zone chondrocyte properties and suggest potential therapeutic targets for achondroplasia.

Achondroplasia, associated with FGFR3 gain-of-function mutations, is a chondrodysplasia characterized by short-limb dwarfism. Here, Horike et al. discover that excess FGFR3 signaling is mediated by CREB and disrupts resting zone chondrocyte properties in a knock-in mouse model harboring the achondroplasia mutation.

Figures

b). After the secondary ossification center was formed, the height of growth plate cartilage was not significantly different between *Fgfr3*^*Ach*^ mice and control mice at P18. However, close examination revealed that the height of resting zone was slightly increased, and the height of proliferative zone was slightly decreased in *Fgfr3*^*Ach*^ mice compared with those of control mice at P18 (Supplementary Fig. 1a,c). The height of growth plate cartilage was significantly reduced in *Fgfr3*^*Ach*^ mice compared to that in control mice at P21, P24, and P28 (Fig. [1f–h](#Fig1) and Supplementary

Detection of proliferating chondrocytes and slowly dividing chondrocytes in *Fgfr3*^*Ach*^ mice. Histological sections from growth plate cartilage in proximal tibia were subjected to detection of EdU.a Short-term EdU labeling assay to detect proliferating chondrocytes. The mice were administered EdU at P20, sacrificed at P21, and subjected to histological analyses. b Long-term EdU labeling assay to detect slowly dividing cells. The mice were administered EdU on consecutive days from P4 to P7, sacrificed at P21, and subjected to histological analysis. Blue color is DAPI. Scale bar, 50 μm. Quant

Lineage tracing analysis of chondrocytes in *Fgfr3*^*Ach*^ mice. Sections of growth plate cartilage were subjected to confocal laser-scanning microscopy analysis.a Mice were administered tamoxifen (100 mg/kg body weight) at P11 and sacrificed at P21. *Left*, histological images of the distal femur. Scale bars, 50 μm. *Right*, number of stacked cells in each column was counted. Each circle represents one column of cells. Error bars denote means ± s.d. *P = 0.0426 by two-sided unpaired Student’s *t*-test. Data were representative of two pairs of *Fgfr3*^*Ach*^ mice and control mice. b *Top*, the

Identification of the CREB pathway as downstream of FGF signaling in chondrocytes.a *Top*, Semi-serial sections of growth plate cartilage from the proximal tibia of *Fgfr3*^*Ach*^ mice were immunostained. Panels in the fourth row indicate magnification of boxed regions in the respective panels in the third row. Scale bar, 50 μm. Data for FGFR3, phosphorylated FRS2 (P-FRS2), CBP, and phosphorylated ERK (P-ERK) were representative of two control mice and two *Fgfr3*^*Ach*^ mice. Data for phosphorylated CREB (P-CREB) were representative of three control mice and three *Fgfr3*^*Ach*^ mice. *Bottom

a), and deposited the data in the Gene Expression Omnibus database ([GSE275171](GSE275171)). We analyzed a total of 6438 cells, with each cell having an average of 131,911 reads. After processing the data, we performed principal component analysis (PCA), data clustering, dimensionality reduction using Uniform Manifold Approximation and Projection (UMAP), and a two-dimensional projection in Seurat. The resolution parameter was set to 0.1, which revealed eight cell clusters (Supplementary Fig. 6b). Feature plot functions, using a chondrocyte marker (*Col2a1*) and hematopoietic markers (*Lyz2*, *

Differentially active signaling pathways downstream of *Fgfr3*^*Ach*^ between the proliferative/hypertrophic zones and resting zone.*Left*, as previously reported, *Fgfr3*^*Ach*^ signaling activates the STAT and ERK pathways that regulate the proliferation and hypertrophic differentiation of cells in the proliferative and hypertrophic zones, respectively. *Right*, our results showed that *Fgfr3*^*Ach*^ signaling activates CREB in resting zone cells. This signaling pathway is responsible for the impairment of stem cell-like properties of resting zone chondrocytes in *Fgfr3*^*Ach*^ mice.

Full text

Introduction

Achondroplasia (Ach) is the most frequent form of chondrodysplasia and is characterized by disproportionate short-limb dwarfism, occurring at a frequency of 1 in 15,000–25,000 births^1,2^. Individuals with Ach suffer from foramen magnum and spinal canal stenosis, which often cause neurological symptoms and debilitating conditions. Bone growth occurs through a process called endochondral bone formation, where chondrocytes undergo sequential differentiation into three spatially arranged layers in growth plate cartilage: resting, proliferative, and hypertrophic zones^3^. Chondrocytes in the resting zone are irregularly scattered in a bed of cartilage matrix, whereas chondrocytes in the proliferative and hypertrophic zones are arranged in columns parallel to the longitudinal axis of bone. Resting zone chondrocytes are believed to have stem cell-like properties in that they undergo asymmetrical cell divisions: one daughter cell retains its identity as a resting zone chondrocyte; the other differentiates into a proliferative zone chondrocyte that migrates toward the primary ossification center to form the proliferative zone^4^. This zone is characterized by active cell replication. When a proliferative zone chondrocyte divides, the resulting daughter cells align along the longitudinal axis of the bone. As a result, chondrocyte clones are arranged in columns parallel to this axis. This spatial orientation directs growth in a specific direction, leading to the elongated shape observed in many endochondral bones. Hypertrophic zone chondrocytes are generated through the terminal differentiation of proliferative zone chondrocytes that are located farthest from the epiphysis. These cells cease dividing and become enlarged, substantially contributing to the overall growth process^5^.

Mutations in the gene encoding fibroblast growth factor receptor 3 (FGFR3) have been identified in patients with Ach^6,7^. FGFR3 functions as a transmembrane receptor tyrosine kinase. Observations that mice deficient for FGFR3 show bone overgrowth^8,9^ revealed it as a negative regulator of endochondral bone formation in growth plate cartilage, with Ach-causing mutations defined as gain-of-function mutations. Over 97% of Ach cases result from an autosomal dominant missense mutation (p.Gly380Arg) in the transmembrane domain of FGFR3^7,10,11^. Patients with Ach who do not have the p.Gly380Arg mutation harbor other less common FGFR3 mutations, such as p.Ser217Cys, p.Ser279Cys, p.Ser344Cys, and p.Gly375Cys^1^.

Research using genetically modified mice have explored the molecular mechanisms responsible for the disturbance of growth plate cartilage in Ach. The p.Gly380Arg mutation is an ortholog of the p.Gly374Arg mutation in mice, as mouse Fgfr3 is shorter than human FGFR3 by six amino acids at the N-terminal region. Transgenic mice that overexpress Fgfr3 with a p.Gly374Arg mutation under the chondrocyte-specific type II collagen gene (Col2a1) promoter/enhancer sequences^12^ and knock-in mice with the p.Gly374Arg mutation introduced in the Fgfr3 locus^13^ exhibit short-limb dwarfism and decreased height of proliferative and hypertrophic zones in the growth plate. In other studies, knock-in mice with gain-of-function mutations such as p.Lys644Glu^14^, p.Tyr367Cys^15,16^, and p.Gly369Cys^17^ in the Fgfr3 locus were also generated. These genetically modified mice, together with in vitro experiments, revealed that gain-of-function mutations in Fgfr3 inhibit the proliferation of proliferative zone chondrocytes and the hypertrophy of hypertrophic zone chondrocytes by activating the MAPK/ERK^18^ and STAT signaling pathways^19^. Based on the finding that C-type natriuretic peptide (CNP) signaling pathways cross-talk with and downregulate MAPK/ERK signaling in Ach^20^, a stabilized form of CNP called vosoritide was developed, which was shown to rescue the Ach phenotype in mice^21^ as well as in humans^22,23^ and is now clinically available for the treatment of Ach. However, vosoritide administration does not completely restore bone length^22,23^. In addition to the incomplete rescue brought about by current drugs, the lack of a systemic approach in the study of Ach suggests that important pathophysiological mechanisms influencing growth plate pathology have not yet been discovered. While most studies have focused on the proliferative and hypertrophic zones, abnormalities in the resting zone remain understudied. Furthermore, other cell signaling cascades beyond MAPK/ERK and STAT pathways have not been comprehensively investigated.

In the present study, we generated knock-in mice bearing the p.Gly374Arg mutation in Fgfr3 (Fgfr3^Ach^). These mice exhibited an Ach phenotype, including dwarfism. The growth plate cartilage of Fgfr3^Ach^ mice showed an expanded resting zone that hosted a unique pathological cell cluster at postnatal day 21 (P21) when prominent growth disturbance occurred. A systemic approach, including lineage tracing and single-cell RNA sequencing (scRNA-seq) analysis, revealed the molecular mechanisms at play in the pathological cell cluster within the expanded resting zone. These mechanisms may represent novel therapeutic targets against Ach.

Results

Generation of Fgfr3Ach mice

We first introduced the p.Gly374Arg mutation into one allele of the Fgfr3 in mouse ES cells using the CRISPR/Cas9 system. After in vitro fertilization using sperm obtained from these chimeric mice, we generated nine knock-in mice. All nine died before 6 weeks of age, suggesting that the p.Gly374Arg mutation is lethal in mice, making it difficult to establish conventional knock-in mouse lines. Thus, we employed a conditional knock-in approach. Insertion of the neo cassette into intron 9 (formerly designated as intron 10) inhibits expression of the Fgfr3 allele^14,17^. Thus, we introduced a p.Gly374Arg mutation into exon 9 and inserted a neo cassette flanked by FRT sequences into intron 9 of the Fgfr3 allele (Fgfr3^Ach-neo^) in mouse ES cells using the CRISPR/Cas9 system (Fig. 1a). After generating mice from ES cells, we mated them with CAG-Flp transgenic mice. The neo cassette flanked by FRTs at both ends was excised from the oocytes, leaving a p.Gly374Arg mutation and an FRT sequence in the Fgfr3 locus of the progeny mice (Fgfr3^Ach^). We used Flp/FRT rather than the Cre/loxP system because we sought to perform lineage tracing analysis using the CreER/loxP system. Fgfr3^Ach^ mice developed dwarfism (Fig. 1b). The incisors of Fgfr3^Ach^ mice did not align normally because of changes in the skull. We shortened the teeth of the mice and softened the pelleted mouse chow by wetting it. Body weight was decreased as they aged in both males and females (Fig. 1c). Further, 85% of mice died by P36 (Fig. 1d), confirming the results from conventional knock-in mice described above. Skeletal elements, including the femur, tibia, and ulna, were significantly shorter in Fgfr3^Ach^ mice than in those of control mice at P21 (Fig. 1e). The degree of skeletal phenotypes of Fgfr3^Ach^ mice appears more severe than that of human patients with achondroplasia. These phenotypes are consistent with the previously reported phenotype of mice bearing the same mutation^13^.Fig. 1. Generation and phenotype of Fgfr3^Ach^ mice.a Introduction of the p.Gly374Arg mutation and a neocassette into the mouse Fgfr3 locus. Asterisk indicates a mutation that causes p.Gly374Arg; arrowheads denote FRP sequences; an arrow indicates the guide RNA sequence. b Gross appearance of Fgfr3^Ach^ mice at postnatal day 28 (P28). Scale bar, 5 mm. c Weight of Fgfr3^Ach^ mice after birth. N = 16 mice. d Survival rate of Fgfr3^Ach^ mice after birth. N = 43 mice. e. Length of skeletal elements of Fgfr3^Ach^ mice at P21. N = 6 mice. Lengths of each skeletal element (femur, tibia, and ulna) in Fgfr3^Ach^ mice were compared with sex-matched controls. f Histology of growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice at P18, P21, and P28. Hematoxylin-eosin staining and safranin O-fast green-iron hematoxylin staining. Data were representative of two, six, and six control mice and two, six, and six Fgfr3^Ach^ mice at P18, P21, and P28, respectively. g Semi-serial histological sections were immunostained with antibodies that recognize the proteins indicated on the left and stained with safranin O-fast green-iron hematoxylin at P21. N = 2 mice. h Top, magnification of histological images of growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice at P21. Dotted lines indicate the resting, proliferative, and hypertrophic zones. SOC, secondary ossification center; RZ, resting zone; PZ, proliferative zone; HZ, hypertrophic zone; POC, primary ossification center. Bottom, heights of the resting, proliferative, and hypertrophic zones (RZ, PZ, and HZ, respectively) were measured. The heights of the entire growth plate were also measured. N = 4 mice (2 males and 2 females). g, h Data were representative of four control mice and four Fgfr3^Ach^ mice. Data are presented as the means ± s.d. P values by two-sided unpaired Student’s t-test are indicated. Scale bars, 5 mm in (b) and 100 μm in (f–h). Source data are provided as a Source data file.

Resting zone expansion in the growth plate of Fgfr3Ach mice

Histological analysis indicated that primordial cartilage appeared normal in Fgfr3^Ach^ mice at 15.5 days postcoital (E15.5) (Supplementary Fig. 1). Although the formation of the secondary ossification center was slightly delayed at P16 (Supplementary Fig. 1), it was clearly formed as indicated by expression of COL10 at the periphery of the secondary ossification center in Fgfr3^Ach^ mice at P18 (Supplementary Figs. 1–3).

The height of proliferative zone was slightly reduced in Fgfr3^Ach^ mice compared with those of control mice at P8 (Supplementary Fig. 1b). After the secondary ossification center was formed, the height of growth plate cartilage was not significantly different between Fgfr3^Ach^ mice and control mice at P18. However, close examination revealed that the height of resting zone was slightly increased, and the height of proliferative zone was slightly decreased in Fgfr3^Ach^ mice compared with those of control mice at P18 (Supplementary Fig. 1a,c). The height of growth plate cartilage was significantly reduced in Fgfr3^Ach^ mice compared to that in control mice at P21, P24, and P28 (Fig. 1f–h and Supplementary Fig. 2). The proliferative and hypertrophic zone heights decreased at P21, P24, and P28 (Fig. 1h and Supplementary Fig. 2), consistent with previously reported Ach model mice^12,13,16,21^, suggesting the inhibition of both chondrocyte proliferation and terminal hypertrophic differentiation. Importantly, we found a unique abnormality in the resting chondrocyte zone, which has not been clearly mentioned in previously reported papers. The resting zone was expanded in Fgfr3^Ach^ mice compared to that in control mice at P18 and P21 (Fig. 1g, h and Supplementary Figs. 1 and 2). The expansion was slightly at P18 and became remarkable at P21. The expanded resting zone in Fgfr3^Ach^ mice contained limited amounts of proteoglycan, indicated by weak safranin O staining, and type II collagen (COL2), as evident by weak immunohistochemical staining (Fig. 1g, h). Because the secondary ossification center was clearly formed in Fgfr3^Ach^ mice at P18 (Supplementary Figs. 1 and 3), the expanded resting zone is not simply residual epiphyseal cartilage.

Abnormal behavior of Fgfr3Ach resting zone chondrocytes

To analyze chondrocyte proliferation, we performed 5-ethynyl-2’-deoxyuridine (EdU) labeling experiments. To detect actively dividing cells, we performed a short-term EdU chase assay by administering EdU at P20 and sacrificing mice at P21. The number of labeled chondrocytes in proliferative zone in Fgfr3^Ach^ mice significantly decreased (Fig. 2a), consistent with previous reports on Ach model mice^12,13,16,21^. The number of labeled chondrocytes in resting zones was very low both in Fgfr3^Ach^ mice and control mice, indicating that proliferation rates of resting zone chondrocytes were very low in both mice.Fig. 2. Detection of proliferating chondrocytes and slowly dividing chondrocytes in Fgfr3^Ach^ mice. Histological sections from growth plate cartilage in proximal tibia were subjected to detection of EdU.a Short-term EdU labeling assay to detect proliferating chondrocytes. The mice were administered EdU at P20, sacrificed at P21, and subjected to histological analyses. b Long-term EdU labeling assay to detect slowly dividing cells. The mice were administered EdU on consecutive days from P4 to P7, sacrificed at P21, and subjected to histological analysis. Blue color is DAPI. Scale bar, 50 μm. Quantifications of % EdU-positive cells in the resting, proliferative, and hypertrophic zones (RZ, PZ, and HZ, respectively) are indicated in the graphs at the bottom. N = 3 female mice. Data are presented as the means ± s.d. P values by two-sided unpaired Student’s t-test are indicated. Source data are provided as a Source data file.

Next, we analyzed slowly dividing cells in the growth plate cartilage by performing a long-term EdU labeling-chasing assay. We administered EdU to the mice on consecutive days from P4 to P7 and sacrificed the mice at P21. Only a few cells that lined the resting zone were specifically labeled with EdU in control mice (Fig. 2b). This result indicates that among cells that incorporated EdU at P4-7, the majority continuously underwent cell division, while a few underwent slow and limited cell division and retained EdU until P21. In contrast, a number of cells in the expanded resting zone were labeled in Fgfr3^Ach^ mice (Fig. 2b), suggesting that many of the cells that incorporated EdU at P4-7 underwent slow and limited cell division and remained in the resting zone until P21. To exclude the possibility that more chondrocytes were labeled in Fgfr3^Ach^ mice than in control mice during P4-7, we euthanized mice at P8. Chondrocytes were similarly labeled in control and Fgfr3^Ach^ mice at P8 (Supplementary Fig. 4). These results suggest that the resting zone consists of more slowly dividing chondrocytes in Fgfr3^Ach^ mice.

To further analyze behavior of resting zone chondrocytes, we performed lineage tracing analysis. We crossed mice to generate Fgfr3^Ach^ mice bearing chondrocyte-specific Col11a2-CreERT transgene^24^ and multicolor R26R-Confetti reporter gene^25^. Mice were administered tamoxifen at P11 and sacrificed at P21. Stacked cells in each column expressed a single color in the growth plate cartilage of control mice, suggesting that chondrocytes in the resting zone give rise to clonal progeny that migrate toward the primary ossification center (Fig. 3a). This result is consistent with that of a previous report^4^, suggesting the success of our experiments. In contrast, only a few clonal cells were stacked in Fgfr3^Ach^ mice (Fig. 3a). When tamoxifen was administered at P20 and mice were sacrificed at P24, the angles of the clonal cell stacks relative to the longitudinal axis of the growth plate cartilage were more variable in Fgfr3^Ach^ mice than in control mice (Fig. 3b). This suggests that clonal progeny cells moved in random directions and remained within the resting zone, and did not move into the proliferative zone in the Fgfr3^Ach^ mice. This resulted in the accumulation of cells in the resting zone, contributing to the expansion of the resting zone and a decreased population of proliferative zone chondrocytes.Fig. 3. Lineage tracing analysis of chondrocytes in Fgfr3^Ach^ mice. Sections of growth plate cartilage were subjected to confocal laser-scanning microscopy analysis.a Mice were administered tamoxifen (100 mg/kg body weight) at P11 and sacrificed at P21. Left, histological images of the distal femur. Scale bars, 50 μm. Right, number of stacked cells in each column was counted. Each circle represents one column of cells. Error bars denote means ± s.d. *P = 0.0426 by two-sided unpaired Student’s t-test. Data were representative of two pairs of Fgfr3^Ach^ mice and control mice. b Top, the mice were administered tamoxifen (100 mg/kg body weight) at P20 and sacrificed at P24. Histological images of proximal tibia. Bottom, mice were administered tamoxifen at P18, sacrificed at P23 (blue) or P20, and sacrificed at P24 (orange). The angles of clonal cell stacks relative to the longitudinal axis of the growth plate cartilage in the proximal tibia were measured. Scale bars, 50 μm. N = 2 mice. c Scheme explaining the behaviors of resting zone chondrocytes in the growth plate of control and Fgfr3^Ach^ mice. Left, in control mice, resting zone chondrocytes behave like stem cells, with one daughter cell continuously undergoing self-renewal (blue color cells), while the other differentiates into proliferative zone chondrocytes that move toward the primary ossification center (brown color cells). Right, in Fgfr3^Ach^ mice, the majority of resting zone chondrocytes undergo slow cell divisions with progenies moving toward random directions and accumulating in the resting zone. Color density of cells corresponds to intensity of EdU labeling after a long-term EdU chase (Fig. 2b). d Histological images of growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice at P28. Safranin O-fast green-iron hematoxylin staining. Scale bars, 50 μm. Data were representative of two control mice and two Fgfr3^Ach^ mice. e The histological sections were immunostained with antibodies against CD73 at P28. Scale bar, 50 μm. Data were representative of two control mice and two Fgfr3^Ach^ mice. Source data are provided as a Source data file.

Based on the results from EdU labeling and lineage tracing experiments, we propose the following model whose scheme is shown in Fig. 3c. Resting zone chondrocytes behave like stem cells in control mice but not in Fgfr3^Ach^ mice. In control mice, after cell division, one of the daughter cells from resting zone chondrocytes continuously underwent self-renewal, while the other differentiated into proliferative zone chondrocytes that migrated toward the primary ossification center to form columns (Fig. 3c, left). As a result, the number of resting zone chondrocytes does not change. In contrast, in Fgfr3^Ach^ mice, a number of resting zone chondrocytes underwent slow cell divisions, with progeny moving in random directions and accumulating in the resting zone. Consequently, although resting zone chondrocytes underwent cell division slowly, continuous accumulation of progeny within the resting zone resulted in the expansion of the resting zone in Fgfr3^Ach^ mice by P21 (Fig. 3c, right). In addition, reduced migration of progenies toward proliferative zone would contribute to the reduced population of proliferative zone chondrocytes and hypertrophic zone chondrocytes.

Notably, the height of the resting chondrocyte zone decreased in Fgfr3^Ach^ mice between P21 and P28 (Fig. 1f and 3d). This finding suggested that the cell population in the expanded resting zone was decreased after P21. Tunel assay indicated no significantly increased staining in resting zone chondrocytes in Fgfr3^Ach^ mice compared to that in control mice (Supplementary Fig. 5), excluding the possibility that significant apoptosis occurred in Fgfr3^Ach^ resting zone. Heights of proliferative and hypertrophic chondrocyte zones were still decreased in Fgfr3^Ach^ mice at P28, excluding a possibility that resting zone chondrocytes differentiated into proliferative zone chondrocytes after P21. Other possibilities include that slow cell division of resting zone chondrocytes could be one of the reasons why the expanded resting zone was reduced in Fgfr3^Ach^ mice. Overall, Fgfr3^Ach^ resting zone chondrocytes lost stem cell-like properties, being consistent with the observation that immunohistochemical staining for CD73, a chondrocyte stem cell-like marker^4^, indicated that resting zone chondrocytes did express CD73 in control mice at P28, whereas chondrocytes in the expanded resting zone in Fgfr3^Ach^ mice did not (Fig. 3e).

scRNA-seq analysis of the Fgfr3Ach growth plate cartilage

To investigate how excess Fgfr3 signaling causes resting zone chondrocyte abnormalities in Fgfr3^Ach^ mice, we dissected the growth plate cartilage from one pair of Fgfr3^Ach^ and control mice at P19 and two pairs at P22, subjected these to single-cell RNA sequencing (scRNA-seq) analysis (Supplementary Fig. 6a), and deposited the data in the Gene Expression Omnibus database (GSE275171). We analyzed a total of 6438 cells, with each cell having an average of 131,911 reads. After processing the data, we performed principal component analysis (PCA), data clustering, dimensionality reduction using Uniform Manifold Approximation and Projection (UMAP), and a two-dimensional projection in Seurat. The resolution parameter was set to 0.1, which revealed eight cell clusters (Supplementary Fig. 6b). Feature plot functions, using a chondrocyte marker (Col2a1) and hematopoietic markers (Lyz2, Hba-a1, Vpreb3, Mpo, and Emcn) revealed that clusters #0, #1, and #3 consisted of chondrocytes, while the other clusters were composed of blood cells present in the bone marrow (Supplementary Fig. 6c). We then selected cells in clusters #0, #1, and #3 using the subset function and subjected these to another clustering analysis. The dataset comprised 4413 cells, with an average of 123,355 reads per cell. The resolution parameter was set to 0.3, revealing five cell clusters (Supplementary Fig. 7a). Feature plots, using chondrocyte markers confirmed that all cells were chondrocytes (Supplementary Fig. 7b). Regarding the differentiation stages of chondrocytes, the feature plot function revealed that clusters #2 and #3 were enriched in resting zone chondrocytes, as indicated by the preferential expression of Clu^26^ and Pthlh^27^, marker genes for resting zone chondrocytes (Supplementary Fig. 7c). In contrast, cluster #0 was notably enriched in prehypertrophic and hypertrophic chondrocytes, as evident by the preferential expression of Sp7, Pth1r, Ihh, and Col10a1. Cluster #4 showed abundance of dividing cells, as reflected by the high G2M and S scores (high expression of cell cycle-related genes), suggesting that cluster #4 was enriched in proliferative zone chondrocytes. Some cells in cluster #1 expressed Sp7 and Pth1r, while other cells in cluster #1 did not, suggesting that cluster #1 is likely associated with prehypertrophic zone chondrocytes and proliferative zone chondrocytes.

When we compared the number of cells in each cluster between Fgfr3^Ach^ and control mouse samples, cluster #0 (prehypertrophic and hypertrophic chondrocytes) was enriched in control mouse cells (Fig. 4a), consistent with the histological finding that the height of the hypertrophic chondrocyte zone was reduced in Fgfr3^Ach^ mice (Fig. 1h). Conversely, cluster #3 was enriched in Fgfr3^Ach^ mouse cells (Fig. 4a). When the data were reduced to each mouse sample, these enrichment patterns were consistently observed across the three pairs of Fgfr3^Ach^ and control mice (Supplementary Fig. 8a, b).Fig. 4. Characterization of growth plate chondrocytes in the proximal tibia of Fgfr3^Ach^ mice and control mice via scRNA-seq.a Uniform Manifold Approximation and Projection (UMAP) plots (Supplementary Fig 7a) were separated by genotype. UMAP in control and Fgfr3^Ach^ mouse chondrocytes are shown. Five clusters are visualized. b Feature plots of Spon1, the top-ranking DEG of cluster #3 (Supplementary Fig. 9a). Control mouse cells and Fgfr3^Ach^ mouse cells are separately shown in the left and right column, respectively. c Violin plot for the Spon1 for each cluster. Control mouse cells and Fgfr3^Ach^ mouse cells are separately shown. d Left, Sections of the growth plate cartilage of the proximal tibia in control or Fgfr3^Ach^ mice at P21 were immunostained for SPONDIN1, which is encoded by the Spon1 gene. Scale bar, 50 μm. Data were representative of three control mice and three Fgfr3^Ach^ mice. Right, the % SPONDIN1-positive area in the resting zones was calculated. Error bars denote means ± s.d. **P = 0.0017 by two-sided unpaired Student’s t-test. N = 3 mice. Source data are provided as a Source data file.

Analysis of Fgfr3Ach resting zone chondrocytes via scRNA-seq

Results from histological and scRNA-seq analyses of the growth plate suggested that cluster #3 corresponds to cells in the expanded resting zone in Fgfr3^Ach^ mice. To test this hypothesis, we first identified the marker genes of cluster #3. Differentially expressed genes (DEGs) between cluster #3 and the remaining clusters were identified using the FindMarkers function in Seurat (Supplementary Fig. 9a and Supplementary Data 1). The feature plot function confirmed that DEGs were preferentially expressed in cluster #3 (Fig. 4b and Supplementary Fig. 9b). Among DEGs, the top-ranked gene was Spon1, which encodes SPONDIN1 (also known as F-spondin). The violin plot and feature plot indicate that cluster #3 preferentially expressed Spon1 (Fig. 4b, c). Immunohistochemical analysis using an anti-SPONDIN1 antibody indicated that it was specifically expressed in the expanded resting zone in Fgfr3^Ach^ mice, with weak expression detected in the resting zone of control mice (Fig. 4d). In addition, feature plot analysis revealed that cluster #3 preferentially expressed known markers for resting zone chondrocytes such as Clu1, Apoe, and Sfrp5 (Supplementary Fig. 10a). Immunohistochemical analysis using an anti-SFRP5 antibody indicated that it was specifically expressed in the expanded resting zone in Fgfr3^Ach^ mice and in the resting zone of control mice (Supplementary Fig. 10b). These results collectively confirm the hypothesis that cluster #3 corresponds to cells in the expanded resting zone of chondrocytes in Fgfr3^Ach^ mice.

In addition to cluster #3, cluster #2 was also enriched in Fgfr3^Ach^ mouse cells (Fig. 4a and Supplementary Fig. 8a, b). DEGs between cluster #2 and the remaining clusters included Col1a2 and Omd (Supplementary Fig. 11a). The Omd gene encodes OSTEOMODULIN (also known as OSTEOADHERIN), which localize in osteoid surface^28^. The feature plot function confirmed that DEGs were preferentially expressed in cluster #2 (Supplementary Fig. 11b). Immunohistochemical analysis using an anti-COL1 or anti-OSTEOMODULIN antibodies indicated that they were specifically expressed in the boundary between resting zone and secondary ossification center in both control and Fgfr3^Ach^ mice (Supplementary Fig. 11c). Cells in cluster #2 also expressed Col10a1 (Supplementary Fig. 11a), and cells at the boundary between the resting zone and the secondary ossification center were hypertrophic and expressed COL10 (Supplementary Fig. 3). These results suggest that cluster #2 corresponds to cells in the boundary between the resting zone and secondary ossification center. We selected clusters #2 and #3 cells and subjected them to trajectory inference and RNA velocity analysis. The result indicated that cluster #3 became cluster #2 (Supplementary Fig. 11d). These results collectively suggest that the primary lesion resides in cluster #3 rather than cluster #2.

CREB signaling in the resting chondrocytes in Fgfr3Ach mice

We analyzed whether FGFR3 signaling is activated in expanded resting zone chondrocytes in Fgfr3^Ach^ mice. scRNA-seq analysis indicated that Fgfr3 was highly expressed in the cluster #3 cells in Fgfr3^Ach^ mice than those in control mice (Supplementary Fig. 12). Immunohistochemistry confirmed expression of FGFR3 in the expanded resting zone chondrocytes in Fgfr3^Ach^ mice (Fig. 5a). FRS2, a direct substrate of fibroblast growth factor receptors (FGFR), was phosphorylated in chondrocytes within the expanded resting zone of Fgfr3^Ach^ mice (Fig. 5a). It has been reported that, activated FGFR3 signaling increases expression STAT1 and STAT5 in chondrocytes^29^. scRNA-seq analysis revealed that Stat1 was slightly and Stat5a and Stat5b were highly expressed in cluster #3 in Fgfr3^Ach^ mice (Supplementary Fig. 13). Immunohistochemistry confirmed high expression of STAT5a and/or STAT5b in the expanded resting zone chondrocytes in Fgfr3^Ach^ mice (Supplementary Fig. 14). These results indicate that FGFR3 is highly expressed and FGFR3 signal is activated in the expanded resting zone chondrocytes in Fgfr3^Ach^ mice.Fig. 5. Identification of the CREB pathway as downstream of FGF signaling in chondrocytes.a Top, Semi-serial sections of growth plate cartilage from the proximal tibia of Fgfr3^Ach^ mice were immunostained. Panels in the fourth row indicate magnification of boxed regions in the respective panels in the third row. Scale bar, 50 μm. Data for FGFR3, phosphorylated FRS2 (P-FRS2), CBP, and phosphorylated ERK (P-ERK) were representative of two control mice and two Fgfr3^Ach^ mice. Data for phosphorylated CREB (P-CREB) were representative of three control mice and three Fgfr3^Ach^ mice. Bottom, The percentage of the number of P-CREB–positive cells in the resting zone was caliculated. ***P = 0.0002 by two-sided unpaired Student’s t-test. (N = 3 mice). b Ingenuity pathway analysis on DEGs for cluster #3 (Supplementary Fig. 4a) detected enrichment of “CREB Signaling in Neurons.” c Primary chondrocytes from wild-type mice were treated with or without 10 μM forskolin, 10 ng/ml bFGF, or 30 nM BGJ398 (FGFR inhibitor). Forskolin, a cAMP analog, was used as a positive control for CREB activation. Cell lysates were subjected to immunoblot analysis using antibodies against the proteins indicated on the left side. Molecular weights (kDa) as determined by the protein ladder are shown on the right. Data were representative of two independent experiments. d Rat chondrosarcoma (RCS) cells were transduced with a CRE-luciferase reporter vector and treated with the growth factors indicated at the bottom. Each circle indicates one sample. e RCS cells were transduced with the CRE-luciferase reporter vector and treated with bFGF at the indicated concentrations. N = 5. f ATDC5 cells were treated with or without 10 μM forskolin or 666-15 (CREB inhibitor). Spon1 mRNA expression was analyzed using real-time RT-PCR. Each circle indicates one sample. g Primary chondrocytes were treated with vehicle, 10 µM forskolin, or 1 µM 666-15 (CREB inhibitor) for 3 days and subjected to western blot analysis. Left: blots representative of three independent experiments are shown. Molecular weights (kDa) as determined by the protein ladder are shown on the right. Right: quantification of SPONDIN1 to the β-ACTIN ratio. N = 3. The data were representative of three independent experiments. Error bars denote means ± s.d. (d, f, g) P values by one-way ANOVA with Tukey’s multiple comparisons test are indicated. Source data are provided as a Source data file.

Ihh expressed in prehypertrophic chondrocytes and PTHrP expressed in the resting zone chondrocytes, regulate the recruitment of resting zone cells into the proliferative zone^30,31^. scRNA-seq analysis and immunohistochemical analysis indicated that expression levels of Ihh and PTHrP were similar between Fgfr3^Ach^ mice and control mice (Supplementary Fig. 12). These results collectively suggest that the abnormalities of the Fgfr3^Ach^ resting zone are the direct effect of mutant FGFR3 signaling in these cells rather than the indirect effect through the Ihh / PTHrP pathway.

To explore the molecular mechanisms responsible for abnormal resting zone chondrocytes in Fgfr3^Ach^ mice, we performed ingenuity pathway analysis (IPA) on DEGs between cluster #3 and the residual clusters (Supplementary Fig. 9a and Supplementary Data 1) and found enrichment of the CREB pathway (“CREB Signaling in Neurons”) (Fig. 5b). Single-cell regulatory network inference and clustering (SCENIC)^32^ analysis also detected Creb1 regulons enriched in Fgfr3^Ach^ mice (Supplementary Fig. 15a) and cluster #3 (Supplementary Fig. 15b). The Creb1 target genes enriched in cluster #3 compared with cluster #1 included Sdc4, Trak1 and Rras (Supplementary Fig. 15c). Immunohistochemistry using an antibody that recognizes phospho-CREB, an activated form of CREB, revealed phospho-CREB expression in the resting zone chondrocytes and hypertrophic zone chondrocytes of control mice (Fig. 5a), as reported previously^33^. In Fgfr3^Ach^ mice, cells in the expanded resting zone also expressed phospho-CREB (Fig. 5a). Phospho-CREB promotes the recruitment of coactivator CREB-binding protein (CBP) to stimulate the transcription of CREB-dependent genes. CBP was also expressed in the expanded resting zone chondrocytes of Fgfr3^Ach^ mice (Fig. 5a), suggesting an interaction between phospho-CREB and CBP for CREB activation in these cells.

To analyze whether excessive FGFR3 signaling activates the CREB pathway, we performed in vitro experiments using mouse primary chondrocytes. The addition of bFGF increased CREB phosphorylation (Fig. 5c), which was abolished by further addition of FGFR inhibitor BGJ398. In rat chondrosarcoma (RCS) cells, bFGF, but not IGF1, IGF2, EGF, or HGF, significantly activated the reporter vector under the control of cAMP response element (CRE)-luciferase (Fig. 5d). We found that bFGF2 activated the CRE-luciferase reporter vector in RCS cells in a dose-dependent manner (Fig. 5e). These results collectively suggest that excess FGFR3 signaling upregulates CREB activity in the expanded resting zone chondrocytes of Fgfr3^Ach^ mice. Regarding Spon1, a marker for the expanded resting zone chondrocytes of Fgfr3^Ach^ mice (Fig. 4), the addition of forskolin increased Spon1 mRNA expression in chondrogenic ATDC5 cells, which was abolished by further addition of CREB inhibitor 666-15 (Fig. 5f). Forskolin is a cAMP analog and activates CREB. The addition of bFGF increased expression of SPONDIN1 protein in primary chondrocytes, which was abolished by further addition of 666-15 (Fig. 5g). These results suggest that excess FGFR3 signal increases Spon1 expression through CREB in the resting zone chondrocytes of Fgfr3^Ach^ mice.

Immunohistochemical analysis indicated that ERK was phosphorylated in hypertrophic chondrocytes, but not in chondrocytes within the expanded resting zone of Fgfr3^Ach^ mice (Fig. 5a). This finding suggested that CREB phosphorylation in the expanded resting zone was independent of the ERK pathway in Fgfr3^Ach^ mice.

To test whether excess CREB activity mediates an abnormal phenotype in Fgfr3^Ach^ mice, we administered the CREB inhibitor 666-15 at 10 mg*/kg from P7 to P27 and sacrificed them at P28. Administration of 666-15 partially and significantly restored the weight (Fig. 6a) and bone growth (Fig. 6b) of Fgfr3^Ach^ mice. Histological analysis revealed that the administration of 666-15 partially and significantly rescued abnormalities in structures of growth plate cartilage in Fgfr3^Ach^ mice at P28 (Fig. 6c). Administration of 666-15 increased areas that express COL2 or COL10 in Fgfr3^Ach^ mice at P28 (Fig. 6d). Heights of the resting zone was decreased and heights of proliferative and hypertrophic zones were increased. The height of the entire growth plate was increased (Fig. 6e). Immunohistochemical analysis indicated that the administration of 666-15 to Fgfr3^Ach^ mice reduced the expression of phosphor-CREB, SPONDIN1, and STAT5 in resting zone chondrocytes in Fgfr3^Ach^ mice (Fig. 6f). Immunohistochemical staining also indicated that 666-15 administration enhanced CD73 expression in resting zone chondrocytes (Fig. 6g), suggesting the recovery of stem cell-like properties. On the other hand, administration of 666-15 significantly changed neither weight, femur length, nor expression of CD73 in the resting zone in control mice (Supplementary Fig. 16a–c), suggesting that the 666-15 at a dose of 10 mg/*kg is not effective in a physiological condition. Collectively, these results suggest that activated CREB is at least partially responsible for the abnormalities in chondrocytes within the resting zone and dwarfism in Fgfr3^Ach^ mice.Fig. 6. Systemic administration of CREB inhibitor 666-15 partially rescued the phenotype of Fgfr3^Ach^ mice. Fgfr3^Ach^ mice were treated with 666-15 or vehicle via intraperitoneal injection from P7 to P27 and sacrificed and analyzed at P28.a Body weights of Fgfr3^Ach^ mice treated with 666-15 (10 mg/kg body weight) or vehicle. N = 10 Fgfr3^Ach^ mice (4 males and 6 females) with vehicle and 6 Fgfr3^Ach^ mice (3 males and 3 females) with 666-15 were recruited. Weights of untreated control mice (N = 11 (6 males and 5 females)) are also indicted. N = 7 Fgfr3^Ach^ mice with vehicle (3 males and 4 females), all 6 Fgfr3^Ach^ mice with 666-15, and all 11 control mice survived at P28. Data are presented as the mean ± s.d. P values by One-way ANOVA with Tukey’s multiple comparisons test at P28 are indicated. b Length of femur of Fgfr3^Ach^ mice treated with 666-15 or vehicle at P28. Each circle represents one mouse. Data are presented as the mean ± s.d. P values by One-way ANOVA with Tukey’s multiple comparisons test are indicated. c Histology of growth plate cartilage in knee joints of Fgfr3^Ach^ mice treated with 666-15 or vehicle. Hematoxylin-eosin or safranin O-fast green-iron hematoxylin staining. Data were representative of 11 control mice, seven Fgfr3^Ach^ mice treated with vehicle, and six Fgfr3^Ach^ mice with 666-15. d Sections of the growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice treated with 666-15 or vehicle were stained with hematoxylin-eosin or immunostained for COL2 or COL10. Data were representative of two Fgfr3^Ach^ mice treated with vehicle and two Fgfr3^Ach^ mice with 666-15. e Heights of the resting, proliferative, and hypertrophic zones (RZ, PZ, and HZ, respectively) in the proximal tibia of Fgfr3^Ach^ mice treated with 666-15 or vehicle. The heights of the entire growth plate were also measured. Each circle represents one mouse. N = 7 Fgfr3^Ach^ mice with vehicle and 6 Fgfr3^Ach^ mice with 666-15. As for control mouse sample, six mice were randomly selected from 11 mice. Data are presented as mean ± s.d. P values by one-way ANOVA followed by Tukey’s multiple comparisons test are indicated. f Left, Sections of the growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice treated with 666-15 or vehicle were immunostained. Data for phosphorylated CREB (P-CREB) and STAT5a/b were representative of two Fgfr3^Ach^ mice treated with vehicle and two Fgfr3^Ach^ mice with 666-15. Data for SPONDIN1 were representative of three Fgfr3^Ach^ mice treated with vehicle and three Fgfr3^Ach^ mice with 666-15. Right, the % SPONDIN1-positive area in the resting zones was calculated. Error bars denote means ± s.d. ***P = 0.0007 by two-sided unpaired Student’s t-test. N = 3 mice. g Sections of growth plate cartilage in the proximal tibia of Fgfr3^Ach^ mice treated with 666-15 or vehicle were immunostained for CD73. Data were representative of two Fgfr3^Ach^ mice treated with vehicle and two Fgfr3^Ach^ mice with 666-15. Scale bar, 100 μm. Each circle represents one mouse. Source data are provided as a Source data file.

We finally analyzed CREB activity in growth cartilage at earlier stages by immunohistochemistry. Expression of SPONDIN1, phospho-CREB, and STAT5a/b were not so much different between Fgfr3^Ach^ mice and control mice at P8 (Supplementary Fig. 17b). The number of resting zone chondrocytes expressing SPONDIN1, phospho-CREB, and STAT5a/b was slightly increased in Fgfr3^Ach^ mice compared to those in control mice at P18 (Supplementary Fig. 17b). Based on these results, we speculated that increased CREB activity and impaired behavior of resting zone chondrocytes might begin at earlier stages and become prominent at P21. The reason why CREB is activated after the formation of the resting zone (after formation of the secondary ossification center) remains to be investigated.

Discussion

Decreased numbers of proliferative and hypertrophic zone chondrocytes in Ach growth plates are thought to result from the decreased proliferation and hypertrophic differentiation of cells within each zone^12,13,16,21^. FGFR3 signaling inhibits proliferation by overexpression and activation of STATs in proliferative zone chondrocytes^1,19,29^, while suppressing hypertrophic differentiation through MAPK/ERK activation in hypertrophic chondrocytes^1,18^. Current treatments for Ach, including vosoritide, are intended to target these signaling molecules in chondrocytes, exclusively in the proliferative and hypertrophic zones^1,20^. In this study, we found expansion of the resting zone in Fgfr3^Ach^ mice harboring the Ach mutation. Similar expansion of the resting zone can be recognized in some of the figures of the previous papers that reported Ach model mice^13^, although the authors did not mention it in the text. Our comprehensive examination of Fgfr3^Ach^ mice led to the discovery of an aberrant stem cell-like population within the expanded resting zone. We found that excess FGFR3 signaling disrupts the stem cell-like properties of resting zone chondrocytes. In addition, excess FGFR3 signaling in resting zone chondrocytes activates the CREB pathway, which is distinct from the previously identified STAT or MAPK/ERK pathways in proliferative and hypertrophic zones (Fig. 7). One reason vosoritide, which downregulates MAPK/ERK signaling, does not correct bone length completely^22,23^ may be that it does not target CREB signaling in resting zone chondrocytes.Fig. 7. Differentially active signaling pathways downstream of Fgfr3^Ach^ between the proliferative/hypertrophic zones and resting zone.Left, as previously reported, Fgfr3^Ach^ signaling activates the STAT and ERK pathways that regulate the proliferation and hypertrophic differentiation of cells in the proliferative and hypertrophic zones, respectively. Right, our results showed that Fgfr3^Ach^ signaling activates CREB in resting zone cells. This signaling pathway is responsible for the impairment of stem cell-like properties of resting zone chondrocytes in Fgfr3^Ach^ mice.

Thus, CREB may represent a new therapeutic target for Ach. Administration of CREB inhibitor 666-15 reverted growth plate abnormalities and elongated bone in Fgfr3^Ach^ mice. Systemic inhibition of CREB by 666-15 is well-tolerated in mice^34^, although CREB plays an important role in multiple organs throughout the body and can cause side effects when administered systemically. One approach to prevent these side effects is through local drug delivery to the cartilage. Cystine-dense peptides (CDPs)^35^ and octaarginine^36^ preferentially accumulate in cartilage and can be used as carriers for the selective delivery of drugs to cartilage.

Vosoritide targets MAPK/ERK pathway in proliferative zone chondrocytes and hypertrophic zone chondrocytes, while 666-15 targets CREB pathway in resting zone chondrocytes. The difference in target pathways and the difference in target cells raise a hypothesis that simultaneous administration of vosoritide and 666-15 would cause additive effects on elongation of bone of Fgfr3^Ach^ mice. This hypothesis has yet to be tested in future studies.

We found that CREB activates the transcription of Spon1 in the resting zone chondrocytes of Fgfr3^Ach^ mice. This is consistent with a previous finding that addition of FGF2 increased the expression of Spon1 in articular chondrocytes^37^. Although we did not analyze whether SPONDIN1 mediates CREB-induced abnormalities in resting zone chondrocytes and dwarfism phenotype in Fgfr3^Ach^ mice, it is noteworthy that inactivation of SPONDIN1 causes elongation of limb bones in organ culture^38^, suggesting that SPONDIN1 can be a drug target for Ach.

The stem cell-like properties of resting zone chondrocytes are maintained by Hedgehog^4,30^, Wnt^39^, and IGF1^26^. However, the upstream and downstream regulatory mechanisms have not yet been clarified. In this study, we found that excess FGFR3 signaling due to a gain-of-function mutation prevents resting zone chondrocytes from behaving like stem cells by limiting the number of cell divisions and by inhibiting the provision of proliferative chondrocytes that migrate toward the primary ossification center. Additionally, we found that the disruption of stem cell-like properties is caused by CREB activation. These mechanisms are at least partly responsible for the dwarfism in Ach. The current findings will contribute to a deeper understanding of the molecular mechanisms that maintain the stem cell-like properties of resting zone chondrocytes and will facilitate drug discovery for Ach.

Methods

Animal model: creation of Fgfr3Ach conditional knock-in mice

First, we introduced a p.Gly374Arg (c.1120G > A) mutation in Fgfr3 on one allele of mouse C57BL/6 ES cells (EGR-G101) using the CRISPR/Cas9 system. The target sequence of the guide RNA was 5′-tacgcaggcgtcctcagcta-3′. We also introduced the c.1107C > A, c.1110C > A, and c.1113C > G mutations to prevent the cleavage of KI alleles and vectors by sgRNA/Cas9, establishing a ScaI site without changing the amino acid sequence. We then PCR-cloned the 5′ homology arm sequence containing c.1120G > A and the 5′ homology arm sequence from the genomic DNA of mutant ES cells. We inserted these homology arms into pNT1.1, which contains a neocassette, to create the targeting vector (Fig. 1a, top). We transfected EGR-G101 cells with the targeting vector, plasmids encoding Cas9, and guide RNAs designed for intron 9. The target sequences of the guide RNAs were 5′-gatgcctccttatctccata-3′ and 5′-ggctccacaccttgcctcat-3′. Mutant ES cell clones were established and used to generate chimeric mice, which were crossed with C57BL/6 mice to produce knock-in mice (Fgfr3^Ach-neo^). Then, the Fgfr3^Ach-neo^ knock-in mice were mated with B6-Tg(CAG-FLPe)36 mice, which express flippase (Flp) under the control of CAG promoter (provided by the RIKEN BRC through the National Bio-Resource Project of MEXT, Japan), to delete the neo cassettes flanked by FRT sequences and generate progeny mice bearing a mutant allele (Fgfr3^Ach^) (Fig. 1a, bottom). Fgfr3^Ach^ mice had malalignment of incisors, and we shortened the teeth of the mice and softened the pelleted mouse chow by wetting to allow for food intake.

Mice were housed under a controlled 14:10 h light–dark cycle, with the light phase from 07:00 to 21:00 and the dark phase thereafter. The ambient temperature was maintained at 23 ± 1.5 °C with a relative humidity of 45 ± 15%. The mice were euthanized in a carbon dioxide (CO_2_) chamber for experiments.

Mouse strains for lineage tracing analysis

Col11a2-CreERT mice^24^ express CreERT under the control of a chondrocyte-specific Col11a2 promoter/enhancer. Cre-mediated DNA recombination is dependent on tamoxifen. Rosa26R-Confetti^25^ (referred to as R26R-Confetti or Confetti, and obtained from the Jackson Laboratory) is a reporter mouse strain that contains the brainbow 2.1 construct. Upon Cre-mediated DNA recombination, each allele randomly expressed one of four different fluorescent proteins (nuclear green, cytoplasmic red, cytoplasmic yellow, and membrane-bound cyan) in a stochastic manner, enabling clonal identification.

For clonal genetic tracing, mice were crossed to generate the Col11a2-CreERT; Rosa26R-Confetti; CAG-Flp; Fgfr3^Ach-neo^ strain. Mice were administered a single intraperitoneal injection of tamoxifen at a dose of 100 μg per g body weight, to achieve recombination of varying efficiency for the spatial separation of individual clones.

Short-term and long-term EdU chasing assays

EdU (Invitrogen) dissolved in phosphate buffered saline (PBS) was administered to mice at indicated postnatal days. Mice were sacrificed at P21. A Click-iT Imaging Kit with Alexa Fluor 488-azide (C10337; Invitrogen, Carlsbad, CA, USA) was used to detect EdU in cryosections.

Histological analysis

Samples were fixed with 4% paraformaldehyde and embedded in paraffin. Semi-serial sections were stained with hematoxylin-eosin or safranin O-fast green and hematoxylin.

The heights of growth plate cartilage zones were measured at 20 positions across the growth plates in the proximal tibia using BZ-X800 analyzer software (KEYENCE, Osaka, Japan). Individual growth plate zones were defined following previously described morphological criteria^40^. The resting zone was defined as the region distal to the secondary ossification center containing round, single chondrocytes. The proliferative zone adjacent to the resting zone was characterized by flattened cells packed into multicellular clusters to form columns of chondrocytes perpendicular to the growth axis. The hypertrophic zone comprised all enlarged chondrocytes distal to the proliferative zone and proximal to the primary center of ossification. Measurements were averaged to yield a single value for each zone for each animal.

For immunohistochemistry, the antigens were unmasked by treatment with hyaluronidase and EDTA. Supplementary Table 1 lists the antibodies used in this study. Anti-type I collagen, anti-type II collagen, anti-type X collagen, anti-f-spondin, anti-Erk-p, anti-CREB-p, anti-FRS2-p, and anti-CD73 antibodies were detected using a CSA II Biotin-free Tyramide Signal Amplification System Kit (Agilent Technologies, Santa Clara, CA, USA) with 3,3′-diaminobenzidine (DAB) as the chromogen. Antigen–antibody binding was detected using an ImmPACT AMEC Red Peroxidase Substrate Kit (SK-4285; Vector Laboratories, Newark, CA, USA) for CD73 and CBP.

For TUNEL assay, apoptotic cells were detected using the In situ Cell Death detection kit, TMR red (fluorescein, Roche Diagnostics).

Single-cell preparation for scRNA-seq analysis

Sample collection, cellular barcoding, and flow cytometry were performed as previously reported^41^. Briefly, the mice were sacrificed, and tissues around the growth plate cartilage in the proximal tibia were dissected. Samples were cut into small pieces (1–2 mm) and incubated in digestion medium (preparation medium with 0.2 mg/ml Liberase^TM^ [Roche, Basel, Switzerland] and 2 kU/ml DNase Ⅰ [Merck, Rahway, NJ, USA]) for 4 h. The samples were then filtered using a 70 μm cell strainer (BD Biosciences, Franklin Lakes, NJ, USA).

For cell surface labeling, cell surface proteins were biotinylated as previously described^42^. Cells were resuspended in 1 ml of PBS supplemented with 1% fetal bovine serum (FBS) and 1 ng EZ-Link Sulfo-NHS-Biotin (Thermo Fisher Scientific, Waltham, MA, USA) for 10 min at 4 °C and washed. Next, the cells were stained with 0.6 μg/ml of TotalSeq (A0951-A0955, and A0437 [BioLegend, San Diego, CA, USA]) for 20 min at 4 °C, washed, and resuspended in PBS supplemented with 10% FBS and 1 μM Sytox Blue Dead Cell Stain (Invitrogen, Waltham, MA, USA) for 5 min at room temperature.

PE-positive and Sytox blue-negative cells were sorted using a FACS Aria II flow cytometer (BD Bioscience) with BD FACS Diva 9.0.1 (BD Biosciences) and suspended in PBS supplemented with 20% FBS. Gating strategy is shown in Supplementary Fig. 18.

Library preparation, sequencing, and FASTQ file preprocessing for scRNA-seq analysis

Library preparation and the first processing of paired-end FASTQ files were performed according to the workflow outlined in the previously reported TAS-seq method^43^. Single-cell suspensions were subjected to cDNA synthesis using the BD Rhapsody Express Single-Cell Analysis System (BD Biosciences) and a BD Rhapsody Targeted & Abseq Reagent Kit (BD Biosciences). After reverse transcription, the resulting BD Rhapsody beads were treated with exonuclease I at 37 °C for 60 min at 1200 rpm on an Eppendorf Thermomixer C with Thermotop. The beads were then chilled on ice, whereafter the supernatant was removed. The beads were washed, resuspended, and stored at 4 °C. During the washing step, bead-containing DNA LoBind tubes were replaced twice.

cDNAs were amplified via the TAS-seq method at Immunogeneteqs Inc. (Noda City, Chiba, Japan) and sequenced with an Illumina Novaseq 6000 sequencer (Illumina, San Diego, CA, USA) using a Novaseq 6000 S4 Reagent Kit v1.0 or v1.5.

To assign cDNA reads to each transcript, bowtie2-indexes built from reference RNA sequences (cDNA and ncRNA fasta files from the Ensembl database; GRCm38.p6 Ensembl release 102^44^) were used. The inflection point in the knee plot (total read count versus rank of the read count) was detected using the DropletUtils package^45^ in R 3.6.3 (https://cran.r-project.org/) from the resulting single-cell gene expression matrix files. Cells whose total read count exceeded the inflection point were considered valid. Demultiplexing of single cells via expression of TotalSeq streptavidin/anti-biotin was performed as described previously^46^.

ScRNA-seq data processing and analysis

The sequencing depth was approximately 132,000 reads per cell. Data were imported into the Seurat R package (version 4.0.2)^47^. Quality control was performed on all datasets, retaining cells with 2000–7000 detected genes and <6% mitochondrial transcripts. We applied log-normalization to all datasets, using a “scale factor” of 1000,000 molecules for each cell, identifying the top 5000 highly variable genes as reported previously^43^. We calculated the cell cycle phase scores from the integrated data using the CellCycleScoring function and scaled the data by regressing them to mitigate cell cycle heterogeneity. We performed PCA on the scaled expression values using the first 30 principal components to build an SNN graph for data clustering, dimension reduction with UMAP, and two-dimensional data projection.

DEGs between clusters were identified using the FindMarkers function in Seurat. DEGs were subjected to IPA (QIAGEN Inc., version 81348237)^48^.

SCENIC analysis

To infer transcription factor (TF) regulatory networks, we applied SCENIC^30^ using pySCENIC (package version 0.12.1) in Python. The count matrix from the preprocessed scRNAseq data (as described above) was used as input, together with auxiliary datasets, including the cisTarget database (mc9nr), the motif-to-TF annotation file (v9), and the TF list (allTFs_mm). Only genes present in the cisTarget database (approximately 18,000) were considered for the analysis. First, GRNBoost2 was employed to infer co-expression modules using the “grnboost2” function, from which candidate regulons were derived. Next, cisTarget analysis was performed using the “df2regulons” function to identify regulons whose regulator binding motifs were significantly enriched across their target genes. This resulted in the identification of 723 regulons, with the Creb1 regulon ranked 99^th^ by normalized enrichment score (NES) in the cisTarget analysis. Finally, AUCell scores were computed using the “aucell” function to quantify regulon activity for each cell barcode.

For Supplementary Fig.15a, AUCell scores for the Creb1 regulon were averaged across biological replicates (d19g.cont.3, d22d.cont.1, and d22e.cont.6 for Control; d22d.ach.2, d22e.ach.5, and d19g.ach.4 for Fgfr3^Ach^). For Supplementary Fig. 15b, AUCell scores for the Creb1 regulo n were averaged across both biological replicates and clusters (clusters #0-4 for d19g.cont.3, d22d.cont.1, d22e.cont.6, d22d.ach.2, d22e.ach.5, and d19g.ach.4).

RNA velocity analysis

We performed the RNA velocity analysis as described previously^41^. TAS-Seq data cDNA reads were mapped to the reference genome (GRCm38 release 101) using HISAT2-2.2.1^49^ and the following parameters: -q -p 6 –rna-strandness F –very-sensitive –seed 656565 –reorder –omit-sec-seq –mm. For the HISAT2 index build, a corresponding Ensembl GTF file was filtered to retain protein-coding RNA, long non-coding RNA, and T cell chain/immunoglobulin chain annotations according to the 10X Genomics’s method (https://support.10xgenomics.com/single-cell-gene expression/software/pipelines/latest/advanced/references#mkgtf). Then, the cell barcode information of each read was added to the HISAT2-mapped BAM files, and associated gene annotations were assigned using featureCounts v2.0.2^50^ and the following parameters: -T 2 -Q 0 -s 1 -t gene -g gene_name –primary -M -O –largestOverlap –fraction -R BAM. In the featureCounts analysis, a “gene” annotation was used to capture unspliced RNA information for the RNA velocity analysis, and primary annotations were maintained. The resulting BAM file was split using valid cell barcodes and nim 1.0.6 and hts-nim v0.2.23. The split files were processed into loom files using velocyto run (version 0.17.17) with the -c and -U options, and the loom files were concatenated using the loompy. combine function (version 3.0.6)^51^.

Then, we performed RNA velocity analysis using scVelo^52^. We read the loom files to an AnnData object. After estimating the RNA velocity, we inferred the trajectory using PAGA^53^, which was extended by the velocity-inferred directionality. We also performed dynamical modeling and velocity analysis.

Culture of mouse primary chondrocytes, rat chondrosarcoma (RCS) cells, and ATDC5 cells

Primary chondrocytes were prepared from C57BL/6 mice at 18.5 postcoitus as described previously^54^. Briefly, epiphyseal cartilage was dissected from the knee, elbow, shoulder joints, and femoral heads of mice and digested with 3 mg/ml collagenase D (Roche) in Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F12; Invitrogen) containing 5% FBS and 1% penicillin-streptomycin (Life Technologies, Carlsbad, CA, USA) at 37 °C overnight. Approximately 5 × 10^5^ primary chondrocytes were obtained from each mouse and cryopreserved in LaboBanker (Kurabo Industries, Ltd., Osaka, Japan). Before the experiments, the cells were thawed, plated, and cultured in DMEM/F12 supplemented with 5% FBS and 1% penicillin-streptomycin for less than 10 days. Chondrocytes were seeded (0.5 × 10^5^ cells/well) into 12-well tissue culture plates (Corning, Corning, NY, USA). Chondrocytes were cultured in DMEM/F12 containing 5% FBS and 1% penicillin-streptomycin (Invitrogen) at 5% CO_2_ in humidified air. Before the experiments, cells were pretreated overnight with a serum-free starvation medium (serum-free). For experiments, cells were treated with 10 ng/ml bFGF (PeproTech, Cranbury, NJ, USA), 10 μM forskolin (F3917, Sigma-Aldrich, St. Louis, MO, USA), or 30 nM NVP-BGJ398 (ChemScene LLC) in dimethyl sulfoxide (DMSO) for 30 min for immunoblot analysis. Primary chondrocytes were pre-incubated with or without 30 nM NVP-BGJ398 for 1 h before stimulation with 10 ng/ml bFGF. Primary chondrocyte cells were treated with 10 µM forskolin, or 30 nM NVP-BGJ398, or 1 µM 666-15 (S8846, Selleck Chem) for 72 h were lysed and analyzed by western blotting.

Rat chondrosarcoma (RCS) cells^55^ were provided by Dr. James H. Kimura (Section of Biochemistry, Bone and Joint Center, Henry Ford Hospital, Detroit, Michigan). RCS cells were cultured in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin at 37 °C at 5% CO2. RCS cells were used for luciferase assay.

ATDC5 cells were obtained from Riken BRC and cultured in the maintenance medium consisting of a 1:1 mixture of DME and Ham’s F-12 medium (Nacalai Tesque, Cat. 08460-95) containing 5% FBS. ATDC5 cells were treated with 10 µM forskolin or 1 µM 666-15 (S8846, Selleck Chem) for 48 h, and subjected to mRNA expression analysis by real-time quantitative RT-PCR.

Immunoblot analysis

Mouse primary chondrocytes were lysed in RIPA buffer (10 mM Tris-HCl [pH 7.5], 150 mM NaCl, 0.1% SDS, 0.1% sodium deoxycholate, 1 mM EDTA, 1% NP-40, complete protease inhibitors from Roche, and phosphatase inhibitor cocktails 1 and 2 from Sigma-Aldrich). The growth plate cartilage was dissected from the proximal tibiae of mice at P21, frozen in liquid nitrogen, crushed using a Multi-Beads Shocker (Yasui Kikai, Osaka, Japan), and lysed in RIPA lysis buffer. These samples were separated on a 4–12% gradient sodium dodecyl sulfate-polyacrylamide gel electrophoresis, and the proteins were electroblotted onto a polyvinylidene difluoride membrane (Invitrogen). The membranes were immunostained with rabbit anti-pCREB (Cell Signaling Technology, Danvers, MA, USA; #9198 s, 1:1000), rabbit anti-CREB (Cell Signaling Technology #9197 s, 1:1000), rabbit anti-pERK (Cell Signaling Technology # 4370 s, 1:1000), rabbit anti-ERK (Cell Signaling Technology # 9102 s, 1:1000), rabbit anti-histone (Cell Signaling Technology #9715 s), and rabbit anti-β-actin (Cell Signaling Technology #4967 s, 1:1000) antibodies. Goat anti-rabbit IgG-HRP (Santa Cruz, Santa Cruz, CA, USA; 1:5000) was used as the secondary antibody. An ECL system and LAS4000 (GE Healthcare, Chicago, IL, USA) were used for chemiluminescent immunodetection. Uncropped blots are provided in the Source data file. Antibodies are listed in Supplementary Table 1.

Luciferase assay

The following plasmids were obtained from commercial sources: pTAL, pTAL-CRE, and pM from Clontech (Palo Alto, CA, USA); pGL4.10 and pRL-TK from Promega (Madison, WI, USA). RCS cells^55^ in a 12-well plate were cotransfected with the pTAL-CRE vector (0.5 μg/well) with an internal reporter, pRL-TK (0.03 μg). After 48 h of transfection, the cells were treated with or without bFGF (1, 3, 10, 30, or 100 ng/ml; WAKO, cat. no. 064-05381) and FGF-1 (10 ng/ml; PeproTech, cat. no. 100-11), IGF2 (10 ng/ml; PeproTech, cat. no. 100-12), EGF (10 ng/ml; PeproTech, cat. no. AF-100-15), and HGF (10 ng/ml; PeproTech, cat. no. 100-39H) for 6 h, then harvested to measure luciferase activity using a Dual-luciferase Reporter Assay System (cat. #E1910, Promega Corp.). The specific promoter activities of CRE genes were expressed as fold change compared to the reporter activity of the empty vector. Luciferase activity was measured and normalized to Renilla luciferase activity.

mRNA expression analysis

Total RNA was extracted using the RNeasy kit (Qiagen, Hilden, Germany). For quantitative reverse transcription PCR (RT-PCR), total RNA was reverse-transcribed into first-strand cDNA using ReverTra Ace (Toyobo, Osaka, Japan) and an oligo(dT)20 primer. PCR amplification was performed using the KAPA SYBR FAST qPCR Master Mix ABI Prism kit (KAPA Biosystems, Wilmington, DE, USA) and StepOnePlus Real-Time PCR System (Thermo Fisher Scientific). Sequences of the PCR primers used are listed in Supplementary Table 2. The RNA expression levels of target genes were normalized to those of β-actin, and the results indicate relative expression.

Administration of CREB inhibitor 666-15 to mice

Ten milligrams of 666-15 (Selleck; S8846) were dissolved in 100 ml DMSO. The prepared solution was preserved at −20 °C. Control and Fgfr3^Ach^ mice were intraperitoneally injected with 666-15 at a dose of 10 mg*/*kg for 5 days each week. After three weeks of treatment, Fgfr3^Ach^ mice were sacrificed at 28-days-old and subjected to X-ray imaging (Faxitron X-ray DX-50). The hind limbs were dissected, and the femur was isolated and subjected to X-ray imaging. Femur lengths were measured using the FIJI ImageJ software (National Institutes of Health). We measured the distance between the lateral metaphyseal end of the femoral head and the bottom of the concave between the femoral condyles and designated it as the femoral length, as reported previously^56^.

Statistical analysis

Data are presented as the mean ± standard deviation. Two-sided unpaired Student’s t-test was used for two-sample comparisons, and one-way analysis of variance with post-hoc Tukey’s Honest Significant Difference test was used for multiple comparisons. GraphPad Prism 8 (GraphPad Software, La Jolla, CA, USA) was used for the statistical analysis. Assumptions of normality and equal variance were not formally tested; however, data distributions appeared approximately normal based on visual inspection of histograms, and group variances were considered comparable. All statistical tests were two-sided unless otherwise noted. Each data point represents an independent biological sample unless otherwise specified. No statistical method was used to predetermine sample size. The number of animals to be used was estimated from the previously reported experiment^56^ and our pilot study. The experiments were not randomized. The Investigators were not blinded to allocation during experiments and outcome assessment.

Ethics statement

All methods were performed in accordance with relevant guidelines and regulations. Experiments using recombinant DNA were approved by the Recombinant DNA Experiments Safety Committee of Osaka University (Approval number: #04794). All animal experiments were approved by the Institutional Animal Committee of Osaka University (Approval number: #Doi- 03-044-027, #Biken-AP-H30-01, and #Biken-AP-R03-01).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information

Supplementary Information Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Reporting Summary Transparent Peer Review file

Source data

Source Data