Dmrt2 promotes transition of endochondral bone formation by linking Sox9 and Runx2

Endochondral bone formation is fundamental for skeletal development. During this process, chondrocytes undergo multiple steps of differentiation and coordinated transition from a proliferating to a hypertrophic stage, which is critical to advance skeletal development. Here, we identified the transcription factor Dmrt2 (double-sex and mab-3 related transcription factor 2) as a Sox9-inducible gene that promotes chondrocyte hypertrophy in pre-hypertrophic chondrocytes. Epigenetic analysis further demonstrated that Sox9 regulates Dmrt2 expression through an active enhancer located 18 kb upstream of the Dmrt2 gene and that this enhancer’s chromatin status is progressively activated through chondrocyte differentiation. Dmrt2-knockout mice exhibited a dwarf phenotype with delayed initiation of chondrocyte hypertrophy. Dmrt2 augmented hypertrophic chondrocyte gene expression including Ihh through physical and functional interaction with Runx2. Furthermore, Dmrt2 deficiency reduced Runx2-dependent Ihh expression. Our findings suggest that Dmrt2 is critical for sequential chondrocyte differentiation during endochondral bone formation and coordinates the transcriptional network between Sox9 and Runx2.


Introduction
Endochondral bone formation is the fundamental process of skeletal development in vertebrates 1  Endochondral bone formation occurs by sequential steps of chondrocyte differentiation 2 . Chondrocytes rst arise from mesenchymal cells derived from cranial neural crest cells, sclerotomes, and lateral plate mesoderm. These cells undergo mesenchymal cell condensation and differentiate into early-stage chondrocytes, including round and proliferating chondrocytes that produce abundant chondrocytespeci c extracellular matrix proteins such as collagen type II alpha 1 chain (Col2a1) and Aggrecan (Acan) 2 . Proliferating chondrocytes then stop proliferating and enlarge their cell size to become pre-hypertrophic and hypertrophic chondrocytes 3 , which are characterized by the expression of Indian hedgehog (Ihh) and collagen type X alpha 1 chain (Col10a1), respectively 4,5 . These late-stage chondrocytes then undergo terminal differentiation and produce matrix metalloproteinase 13 (MMP13), which allows vascular invasion into cartilage 6 . Finally, terminal chondrocytes become apoptotic and are replaced by bone 2 . The sequential processes of early, late, and terminal chondrocyte differentiation construct a well-arranged columnar layer of chondrocytes called growth plate chondrocytes.
The function of stage-speci c chondrocytes is strictly controlled by critical transcription factors [7][8][9] . In particular, Sry-related HMG-Box gene 9 (Sox9) plays indispensable roles in chondrocyte development and endochondral bone formation 10 , and runt related transcription factor 2 (Runx2) and Runx3 play essential roles in chondrocyte hypertrophy by directly regulating Ihh expression 9 . Chondrocyte-speci c Sox9knockout mice display severe defects in skeletal development and abnormal craniofacial development 10,11 and double-knockout mice of Runx2 and Runx3 show the complete absence of hypertrophic chondrocytes 9 . Importantly, Sox9 regulates early-chondrocyte gene expressions, Col2a1 and Acan, in collaboration with Sox5 and Sox6 12 . Recent genome-wide analyses further demonstrated the genomewide cooperation of Sox5/6/9 through super-enhancers of chondrocyte genes [12][13][14] . Because Sox9 expression initiates and promotes Sox5 and Sox6 expression 10 , Sox5 and Sox6 are not essential for initiation of chondrogenesis but required for promotion of Sox9-regulated chondrogenesis.
Although it has long been accepted that Sox9 regulates early chondrogenesis, recent studies established that Sox9 is important for chondrocyte hypertrophy. Sox9, in association with myocyte enhancer factor 2C (Mef2c) and AP-1 family members, directly activates Col10a1 expression to promote chondrocyte hypertrophy 15,16 . Moreover, Sox9 is expressed in upper hypertrophic chondrocytes and maintains Runx2 expression 15 . These reports strongly indicate that Sox9 target genes are also involved in chondrocyte hypertrophy and unknown molecules mediate the transition from proliferating to hypertrophic chondrocytes and conduct transcriptional machinery of the Sox9-Runx2 axis. However, in contrast to the wealth of knowledge regarding the Sox9 target genes in early chondrogenesis, the roles of Sox9 target genes in chondrocyte hypertrophy remain poorly understood. Thus, uncovering the target genes of Sox9 and their functional roles in chondrocyte hypertrophy would deepen our understanding of endochondral bone formation.
In this study, we discovered that Dmrt2 coordinates endochondral bone formation as a molecule downstream of Sox9. Sox9 and its partners Sox5 and Sox6 increased Drmt2 gene expression along with chondrocyte differentiation. Drmt2-de cient mice showed the dwarf phenotype and delayed endochondral bone formation. Mechanistically, Dmrt2 promoted Ihh expression and chondrocyte hypertrophy through physical and functional collaboration with Runx2. Epigenetic analysis further revealed that Sox9 directly bound to H3K27ac positive enhancer located 18 kb upstream of the Dmrt2 transcriptional start site (TSS) and that the chromatin status of this enhancer became progressively more active through chondrocyte differentiation. Thus, our ndings provide novel insights into the transcription network controlling endochondral bone formation.

Results
Dmrt2 is a target gene of Sox5/6/9 in primary chondrocytes To uncover the transcription factors involved in chondrocyte differentiation, we attempted to identify genes induced by Sox9 and its transcriptional cofactors, Sox5 and Sox6 (Sox5/6/9), in primary chondrocytes by performing RNA-seq based-cloning approach (Figure 1a). RNA-seq analysis with three biological replicates identi ed 1295 downregulated and 1209 upregulated genes using thresholds of false discovery rate (FDR) < 0.05 and fold change > 2 ( Figure 1b). Relative to the levels in the control, primary chondrocytes overexpressing Sox5/6/9 exhibited the upregulation of many known Sox9 target genes, including Col2a1, Col11a1, Acan, and Matn1 (Figure 1c). Enrichment analysis of GO molecular function revealed 66 DNA-binding transcription factors upregulated by Sox5/6/9 (Figure 1d and   Supplementary Table 1). Among these Sox5/6/9 target transcription factors, we were intrigued by the transcription factor Dmrt2. Although severe skeletal defects associated with DMRT2 mutations have been reported in humans and mice [17][18][19] , little is known about the functional roles of Dmrt2 in endochondral ossi cation. RT-qPCR con rmed that Sox5/6/9 upregulated Col2a1 and Dmrt2 expression in newborn rib chondrocytes, limb bud mesenchyme, and C3H10T1/2 cells (Figure 1e-g). Additionally, among the seven Dmrt family transcription factors, only Dmrt2 showed strong expression in the rib cartilage of newborn mice (Supplementary Figure 1a), as well as Sox5/6/9-dependent induction in differentiating C3H10T1/2 cells (Supplementary Figure 1b). Collectively, these ndings suggest that the collaboration of Sox5, Sox6, and Sox9 induces the expression of Dmrt2 in chondrocytes.

Dmrt2 is expressed in chondrocytes during endochondral bone formation in vivo
To examine whether Dmrt2 is expressed in chondrocytes during endochondral bone formation, we rst determined tissue distribution of Dmrt2 by perfuming RT-qPCR. We observed relatively high expression of Dmrt2, as well as the Sox9 target gene Col2a1, in rib cartilage ( Figure 2a). Additionally, Dmrt2 expression increased during in vitro chondrocyte differentiation of ATDC5 cells in the presence of insulin-transferrinselenium (ITS), monitored by Alcian blue staining and Col2a1 expression ( Figure 2b). Furthermore, immuno uorescence analysis of E15.5 mouse tibia sections revealed that Dmrt2 was highly expressed in pre-hypertrophic chondrocytes, but weakly expressed in proliferating or hypertrophic chondrocytes ( Figure   2c). Taken together, these ndings suggest the possibility that Dmrt2 plays stage-speci c roles in prehypertrophic chondrocytes during endochondral bone formation.

Epigenetic regulation of Dmrt2 in chondrocytes
To investigate the molecular mechanism by which Sox5/6/9 regulates Dmrt2 gene expression, we analyzed an epigenetic dataset of chondrocytes including ATAC-seq, which allows the genome-wide pro ling of the open chromatin region 20 , and a ChIP-seq dataset 21 . We performed combination analysis of ATAC-seq pro les of growth plate chondrocytes (GSE100585) and ChIP-seq pro les of rib chondrocytes for Sox9, H3K27ac, and IgG (GSE69108). As shown in Figure 3a Figure  3a). It should be noted that this region extensively overlaps with the peak of ChIP-seq for Sox9 and H3K27ac, an active enhancer mark of transcription ( Figure 3a). Consistent with the bioinformatics analyses, reporter assay indicated that Sox9 and Sox5/6/9 signi cantly upregulated the transcriptional activity on the enhancer region ( Figure 3b). These ndings suggest that Sox9 upregulates Dmrt2 expression through an 18 kb upstream Sox9-bound enhancer.
We next examined whether the chromatin status of this Sox9-bound enhancer changes during chondrocyte differentiation. To achieve this, we obtained ChIP-seq datasets of E12.5 limb buds for Sox9 (GSE73225) 22 and H3K27ac (GSE45456) 23 , and compared them with that of mature chondrocytes isolated from newborn ribs 21 . We found that Sox9 occupancy and active enhancer mark (H3K27ac) were very weak in E12.5 limb buds compared with the levels in rib chondrocytes (Figure 3c). ChIP-qPCR analysis further demonstrated that the enrichment of H3K27ac in differentiated ATDC5 cells became signi cantly higher than that of undifferentiated ATDC5 cells (Figure 3d). Taken together, these ndings suggest that Sox9 directly regulates Dmrt2 expression through epigenetic regulation of the active enhancer.

Dmrt2 is critical for endochondral bone formation in vivo
We next tested whether Dmrt2 regulates endochondral bone formation in vivo by examining Dmrt2de cient (Dmrt2 -/-) mice. Dmrt2 -/mice died soon after birth, as reported previously 18 . Newborn Dmrt2 -/mice showed a dwarf phenotype in contrast to wild-type (WT) and Dmrt2 heterozygous mice, as determined by skeletal preparations and microcomputed tomography (microCT) analysis (Supplementary  Figure 5b). Next, we quantitatively assessed in situ hybridization sections to measure the lengths of the Col2a1 expressing zone (resting and proliferating chondrocyte zone), Ihh expressing zone (pre-hypertrophic chondrocyte zone) and Col10a1 expressing zone (hypertrophic zone) in the tibiae of E15.0 mice (Figure 5c). We did not observe signi cant differences in the lengths of resting-proliferating zone between WT and Dmrt2 -/embryos ( Figure 5c). However, the Ihhpositive pre-hypertrophic and Col10a1-positive hypertrophic zones were signi cantly shorter in Dmrt2 -/embryos than in WT embryos ( Figure 5c). These ndings suggest that impaired endochondral bone formation is partly responsible for the skeletal abnormalities in Dmrt2 -/mice.

Dmrt2 promotes late chondrogenesis by functional interaction with Runx2
Given the delay in chondrocyte hypertrophy observed upon the loss of Drmt2, we determined whether Dmrt2 promotes the initiation of this process. Importantly, Dmrt2 overexpression in primary chondrocytes signi cantly upregulated the expression of Ihh ( Figure 6a). We noticed that the promoter of mouse Ihh, a speci c marker of pre-hypertrophic chondrocytes, contains a consensus Dmrt2-binding element (BE, GnTACA) ( Figure 6b). We found that Flag-tagged Dmrt2 bound to this element within the Ihh gene promoter in primary chondrocytes as determined by ChIP assay (Figure 6c). A DNA-pulldown assay using a biotinylated Dmrt2-binding element also demonstrated that Dmrt2 directly bound to this element, and the Dmrt2 binding activity was decreased by a non-biotinylated probe in a dose-dependent manner ( Figure 6d). These ndings collectively suggest that Dmrt2 directly promotes Ihh expression during endochondral bone formation.
The transcription factor Runx2 directly regulates Ihh expression, consequently stimulating chondrocyte hypertrophy 9 . We found that the Runx2 binding element in the Ihh proximal promoter is located close to the Dmrt2 binding element ( Figure 6b) and con rmed the direct binding of Runx2 to the Ihh gene promoter by ChIP assay (Figure 6e). Thus, we hypothesized that Dmrt2 functionally collaborates with Runx2 to regulate Ihh expression. To prove this hypothesis, we ectopically expressed Dmrt2 and Runx2 in primary chondrocytes and examined Ihh expression. Notably, primary chondrocytes overexpressing both Dmrt2 and Runx2 displayed higher Ihh mRNA levels than cells overexpressing either Drmt2 or Runx2 alone, suggesting a synergistic interaction (Figure 6f). Dmrt2 and Runx2 also synergized to induce other Runx2 targets in primary chondrocytes, including Alpl (alkaline phosphatase), Col10a1, and Tcf7 (Supplementary Figure 3). Further, coimmunoprecipitation experiments indicated physical association of Dmrt2 with Runx2 ( Figure 6g). Venus-tagged Dmrt2 colocalized with DsRed-tagged Runx2 in the nucleus of C3H10T1/2 cells (Figure 6h), consistent with Dmrt2-Runx2 protein-protein interaction. Taken together, these ndings suggest that Dmrt2 interacts with and enhances Runx2 function at the target gene, Ihh.

Discussion
Endochondral bone formation is achieved through successive chondrocyte differentiation steps that are strictly regulated in a spatiotemporal manner by various transcription factors. Previous biochemical and mouse genetic studies have shown that Sox9 regulates multiple steps of chondrocyte differentiation including early chondrogenesis and chondrocyte hypertrophy 15,16 . Although the mechanisms by which Sox9 and its target genes regulate early chondrogenesis have been well studied, Sox9 target genes and their functional roles in chondrocyte hypertrophy remain poorly understood. In the present study, we discovered that the transcription factor Dmrt2 is selectively expressed in pre-hypertrophic chondrocytes and is induced by Sox5/6/9. We also found that Dmrt2 contributes to endochondral bone formation by promoting Runx2 functions. Our results suggest that Dmrt2 coordinates successive chondrocyte differentiation processes during skeletal development.
We identi ed Dmrt2 as a transcription factor functioning downstream of Sox5/6/9 in primary chondrocytes. Several studies have reported the targets of Sox9 or Sox5/6/9 in chondrocytes, but these targets promote early chondrogenesis and negatively regulate chondrocyte hypertrophy 24 25 . Additionally, Saito et al. found that S100A1 and S100B1 are directly regulated by Sox5/6/9 and they suppress chondrocyte hypertrophy and maturation 24 . We also identi ed the transcription factor Sp6 as a Sox9 target and showed that the deletion of Sp6 in mice results in a dwarf phenotype and impaired limb development 26,27 . These results t the idea that Sox9 plays essential role in chondrogenesis. However, it is now established that Sox9 is required for chondrocyte hypertrophy 15,16 and our ndings suggest that target genes of Sox5/6/9 positively regulate chondrocyte hypertrophy. In addition to Dmrt2, FoxA2 and AP-1 family members, which were identi ed as Sox5/6/9-inducible genes (Supplemental Data 1), promote chondrocyte hypertrophy through functional interaction with Sox9 16,28 . Thus, other Sox5/6/9 target genes and their biological interaction with hypertrophic transcription factors warrant further investigation.
We found that Sox9 directly bound to the enhancer region of the Dmrt2 gene located 18 kb upstream of its TSS by epigenetic analysis using ATAC-seq and ChIP-seq datasets. This region showed both open chromatin and H3K27ac, a histone mark for active enhancers (Fig. 3a). Ohba  Although the ChIP-seq dataset and reporter assay demonstrated that Sox9 is responsible for the 18 kb upstream enhancer of the Dmrt2 gene, Sox9 alone failed to increase Dmrt2 expression in primary chondrocytes, limb bud mesenchyme, and C3H10T1/2 cells (Fig. 1e-g). We propose several possibilities to explain this. First, appropriate amounts of Sox5 and Sox6 are necessary to induce the transcription of Dmrt2. Previous studies have shown that Sox5 and Sox6 secure the DNA binding activity of Sox9 and cooperatively promote chondrocyte gene expression through super-enhancers at the genome-wide level 14,29 . Liu et al. also reported that Sox6 and Sox9 bind genomic regions in the vicinity of each other 14 . It would be interesting to examine whether Sox6 binds a genomic region located close to the Sox9 binding region 18 kb upstream of Dmrt2 TSS. Second, an unknown epigenetic mechanism controlled by Sox5 and Sox6 is necessary for induction of the Dmrt2 gene. Many histone modi cation enzymes including demethylases and acetyltransferases are necessary for activating transcription in chondrocytes. Previous biochemical and epigenetic studies showed that Sox9 functionally interacts with CBP/P300, which works as a histone acetyltransferase 21,30 . Our group has shown that the histone demethylase PHF2 associates with Sox9 to promote Sox9 target gene expression in chondrocytes 31 . Although we did not provide direct evidence showing that Sox5 and Sox6 control histone modi cation, it is likely that they upregulate Dmrt2 expression through epigenetic remodeling. More studies are necessary to uncover the precise molecular mechanism underlying the effects of Sox5 and Sox6 to activate Dmrt2 gene expression in chondrocytes.
Because Sox5/6/9 increase Dmrt2 mRNA and are strongly expressed in resting and proliferating chondrocytes, we rst predicted that Dmrt2 is widely expressed in resting and proliferating chondrocytes. Unexpectedly, immunohistochemical analysis revealed that Dmrt2 protein is selectively expressed in prehypertrophic chondrocytes (Fig. 2c). These ndings raise the possibility that other unknown transcription factors that are selectively expressed in pre-hypertrophic chondrocytes are required to accelerate Dmrt2 expression. Alternatively, post-transcriptional regulation of Dmrt2 controls pre-hypertrophic speci c protein expression. This complexity found to be exhibited for many genes expressed in chondrocytes of the growth plate. For instance, Ihh expression is known to be limited to pre-hypertrophic chondrocytes, but Runx2, which directly regulates Ihh, is widely expressed in both pre-hypertrophic and hypertrophic chondrocytes 9,32 . Additionally, PTHrP is only detectable in the periarticular region, even though the Gli family, signaling molecules of Ihh, are diffusely distributed in round and proliferating chondrocytes 5,33 . The molecular mechanisms that control the speci c expression of chondrocyte genes are not fully understood, so further studies are necessary to clarify these.
Dmrt2 contains a highly conserved DNA binding domain called the DM domain, but its sequence similarity with other Dmrt family members is low outside the DM domain 34 . The DM domain recognizes a consensus sequence and physically interacts with DNA in the minor groove 35 . Dmrt proteins are predicted to bind DNA as heterodimers or homodimers with other Dmrt family proteins 36 . Importantly, Dmrt proteins act as bifunctional transcriptional regulators to activate or repress transcription 34 . For instance, Dmrt1 inhibits Stra8 but activates Sohlh1 in germ cells, which prevents meiosis and promotes spermatogonial development 37 . Whether Dmrt1 activates or represses transcription appears to depend on motifs around Dmrt1 binding sites, which suggests that the function of Dmrt1 depends on coactivators and co-repressors 38 . These ndings raise the possibility that Dmrt2 also exerts a reciprocal function in chondrocytes in addition to the increase of Runx2 function. Notably, pre-hypertrophic chondrocytes are the transition stage from a proliferating to a hypertrophic state, at which the expression of proliferating chondrocyte-speci c genes such as Col2a1 and Aggrecan should be inhibited at the transcriptional level. It would be interesting to determine whether Dmrt2 controls chondrocyte differentiation by suppressing Sox9-dependent gene expression but promoting Runx2-dependent gene expression in pre-hypertrophic chondrocytes. The regulation of Sox9 function by Dmrt2 during endochondral bone formation awaits further investigation.
In conclusion, our work suggests a novel role for Dmrt2 during endochondral bone formation as a transcriptional coactivator of Runx2. These ndings increase our understanding of the molecular mechanisms of endochondral bone formation and provide new insights into the transcription factor network controlling skeletal development.

Cell culture and reagents
The mouse broblast-like cell line C3H10T1/2 and mouse teratocarcinoma cell line ATDC5 were purchased from the RIKEN Cell Bank (Ibaraki, Japan). These cells were cultured at 37°C in a humidi ed 5% CO 2  RNA-seq data were analyzed using iDEP (integrated Differential Expression and Pathway analysis) 40 . Brie y, read count data of three replicates for control and Sox5/6/9 were generated and uploaded to the iDEP website (http://bioinformatics.sdstate.edu/idep/). Differentially expressed genes (DEGs) were identi ed using a threshold of false discovery rate (FDR) < 0.05 and minimal fold change > 2. The raw data have been deposited in the NCBI Gene Expression Omnibus database (GEO GSE155118).

Reverse-transcription polymerase chain reaction (RT-PCR) and RT-qPCR
Total RNA was isolated using a Nucleo Spin RNA Plus kit. cDNA was synthesized using ReverTra Ace ® qPCR RT Master Mix (TOYOBO, Osaka, Japan). For RT-PCR analysis, cDNA was ampli ed using KOD FX (TOYOBO) and then the PCR products were electrophoresed in a 1.6% agarose gel with ethidium bromide. Primer pairs used for RT-PCR analysis are listed in Supplemental Table S2. For RT-qPCR analysis, cDNA was ampli ed with EagleTaq Universal Master Mix (ROX) using a StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Primers and TaqMan probes used for cDNA ampli cation are listed in Supplementary Table S2. The mRNA expression was normalized to β-actin expression levels.

Epigenetic datasets and analysis
ChIP-seq and ATAC-seq datasets were downloaded from the GEO database. We obtained GSE69109 ChIPseq pro les of newborn mouse rib chondrocytes for Sox9 and H3K27ac 21 . ChIP-seq pro les of E12.5 limb buds for Sox9 were obtained from GSE73225 23 and H3K27ac from GSE45456 22 . ATAC-seq datasets of growth plate chondrocytes were downloaded from GSE100585 20 .
FASTQ data of ChIP-seq and ATAC-seq were aligned to the mouse genome (mm10) using Bowtie and peak calling was performed using MACS software with the default settings (p-value cut-off = 1.00e-05). ChIP-seq and ATAC-seq data were visualized using the Integrative Genomics Viewer (IGV).

Reporter assay
Dmrt2 minimal promoter (-100 to +10) and Sox9 binding region located 18 kb upstream of the Dmrt2 TSS were introduced upstream of the luciferase gene. Reporter genes were co-transfected with the expression vectors and Renilla into HEK293 cells using the FuGENE6 reagent. After 48 h of transfection, the cells were lysed and luciferase activity was measured using speci c substrates in a luminometer All mice were maintained on the C57BL/6 background. Littermate embryos were used for histological analysis. All animal experiments were performed using protocols approved by the Animal Committee of Osaka University Graduate School of Dentistry.

Skeletal preparation
The skin of the mice was removed and xed in 95% ethanol overnight. Cartilage tissues were stained with 1.5% Alcian blue followed by staining of bone tissues with 0.02% Alizarin red S. Skeletal samples were photographed under a stereoscopic microscope.
In situ hybridization The protocol for in situ hybridization has been described in a previous report 41  buffer. The magnetic beads were resuspended with SDS sample buffer, boiled for 5 min, and then subjected to western blotting.

Statistical analysis
Randomization and blinding were not performed in the animal studies. Sample sizes were estimated based on previous studies of endochondral bone formation 31,41 . Data were statistically analyzed by Student's t-test for comparison between two groups. For more than two groups, we used one-way analysis of variance (ANOVA) or two-way ANOVA followed by the Tukey test. At least ve mice (n = 5-6) were used for the phenotypic analysis. P-values of less than 0.05 were considered statistically signi cant. infected with control (Cont), Sox9 (Sox9), or Sox5+Sox6+Sox9 (Sox5/6/9) adenoviruses. Col2a1 and Dmrt2 mRNA expression levels were analyzed by RT-qPCR. The RNA level is indicated as the fold increase compared with that of the control. Data are shown as the mean ± s.d. (n = 3). **p < 0.01; one-way analysis of variance (ANOVA) followed by the Tukey test.