Replacement of the cartilaginous callus with woven bone is associated with increased cell proliferation and immune cell infiltration

We performed tibial diaphyseal osteotomy and stabilized the fracture using an intramedullary nail, this is an established model of secondary fracture healing that proceeds through the endochondral ossification process.1,2,3 As we and others have reported for this model,23,24,25 at the end of the inflammatory phase and the early repair phase, i.e., at Day 7 (d7) postfracture, the fracture gap was filled with and surrounded by a fibrocartilaginous (soft) callus (Fig. S1). The soft callus increased in size as the repair phase advanced to d10 (Fig. S1).23 As healing further proceeded to d14, the densely packed chondrocytes that filled the soft callus progressed into a hypertrophic state (Figs. 1a and S2a, b) and secreted the hypertrophy markers type 10 collagen (Col X) (Fig. 1b) and matrix metallopeptidase 13 (MMP13) (Fig. 1c). Additionally, a subset of soft-callus chondrocytes became terminally hypertrophic, as shown by the secretion of type I collagen (Col I) (Figs. 1d and S2c, d), which is an essential protein for soft-callus mineralization. At the 3 timepoints (i.e., d7 to d14), the callus areas distal to the fracture gap were filled with newly formed bone that was interwoven with bone marrow (Figs. S1 and S2b, c, f). As expected, immune cells were enriched in the bone marrow that populated the woven-bone regions, while fewer immune cells were observed in the soft-callus area (Fig. S2e, g). As healing further progressed from d14 to d21, endochondral ossification mediated the replacement of the cartilaginous callus with newly formed woven bone, which was accompanied by the re-establishment of bone marrow elements, including immune cell populations (Figs. 1e, f and S3a–f). The quantitative measurement of gene expression using reverse transcriptase combined with quantitative polymerase chain reaction (RT‒qPCR) corroborated our histological analysis and demonstrated very low expression of chondrocyte and hypertrophic chondrocyte markers at d21 (Fig. 1g); these results indicated the nearly complete resorption of the cartilaginous callus.

Fig. 1

Cells with higher proliferation capacity populate the callus as the cartilaginous callus that bridges the fracture gap is replaced by woven bone. a–d Staining of d14 calli. a Safranin O/fast green staining. The cartilaginous extracellular matrix of the soft callus was stained reddish orange. The white arrows point to the fracture line. The dashed black line outlines the callus. The scale bar = 500 μm. IF staining of Col X (purple) (b) and MMP13 (green) (c). Both proteins are secreted by hypertrophic chondrocytes. d IF costaining of Col I (red) and Col II (green). Areas where both proteins are expressed (yellow) surround mineralized chondrocytes. IF images shown in b–d were captured in the region indicated by the yellow box in a. e As in a except that d21 callus samples were stained. The absence of a reddish orange stain indicates resorption of the cartilaginous callus. f As in d except that d21 callus samples were stained. The image was captured in the region indicated by the yellow box in e. The trabecular structure of woven bone bridging the fracture gap is obvious. In all the IF images, DAPI stains nuclei (blue), and the scale bar = 200 μm. All the images are representative of five callus tissues. g RT‒qPCR quantification of the indicated transcripts in the RNA that was purified from d14 and d21 callus tissues. Both Col1a2 and Aggrecan (Acan) are chondrocyte markers, while Col10a1 is a hypertrophic chondrocyte marker. The expression level of each transcript was normalized to that of β-actin mRNA, and the normalized level on d14 was defined as 100. n = 3. The bar graphs present the average ± SEM. ***P < 0.001 using unpaired Student’s t test. h A volcano plot of RNA-seq data displaying the fold change in gene expression values on d21 relative to d14 (Table S1). Significantly differentially expressed genes (FDR < 0.05) are highlighted in red and blue to indicate upregulated and downregulated genes, respectively. RNA-seq data were generated from three biological replicates at each timepoint. i GO analysis of significantly differentially expressed genes shown in h. Representative examples of functional categories enriched by upregulated (red) or downregulated (blue) genes are shown, and the adjusted P value (Padj) for each pathway is given (see also Figs. S5 and S6a, b and Tables S2 and S3 for comprehensive lists of significantly enriched functional categories)

To gain insight into the biological processes involved in the calcification and resorption of the cartilaginous callus as well as the formation of woven bone, we performed RNA-seq using total RNA that was isolated from callus tissues on d14 and d21. We identified 9 330 differentially expressed genes (Padj < 0.05; DESeq226) at d21 relative to d14 (Figs. 1h and S4 and Table S1). As expected, the expression of genes involved in cartilage development and ossification was downregulated at d21 compared to d14 (Gene Ontology (GO) analysis;27,28 Figs. 1i and S5 and S6a, b and Tables S2 and S3). On the other hand, d21 showed the upregulation of genes in several pathways that are involved in immune cell responses (Figs. 1i and S5 and S6a, b and Tables S2 and S3), which was consistent with the re-establishment of the bone marrow population and the replacement of the dense fibrocartilaginous callus, which has relatively low immune cell infiltration, with the marrow-infiltrated woven bone (Figs. S2 and S3). RT‒qPCR analysis of the relative expression of the immune cell marker CD45 between Days 7 and 21 demonstrated the lowest CD45 expression levels on Days 10 and 14 (Fig. S6c), which are the 2 timepoints that are distant from the initial immune response and are characterized by the largest soft callus areas. On the other hand, the expression level of CD45 was highest on d21 among the 4 analyzed timepoints, and it was increased by ∼7 fold on d21 relative to d14 (Fig. S6c), which was consistent with the RNA-seq data (Table S1) and GO analysis (Figs. 1i and S5 and S6a, b and Tables S2 and S3). Importantly, GO analysis also showed that cell division, DNA replication, and chromosome segregation were among the most upregulated pathways on d21 relative to d14 (Figs. 1i and S5 and S6a, b and Tables S2 and S3). Accordingly, as endochondral ossification proceeds, rapidly proliferating cells, including immune cells, populate the callus to replace cells with lower proliferation capacities. These observations are important to our studies, as rapidly proliferating cells are reported to exhibit global 3′ UTR shortening8,10 (see below).

The transcriptome of the healing callus exhibits widespread utilization of alternative PAs

When the last exon of a gene contains multiple PAs (Fig. 2a), utilization of the distal PA (dPA) during mRNA maturation results in the incorporation of the alternative UTR (aUTR) sequences, generating an APA isoform with a longer 3′ UTR (long APA isoform, or “lAPA”) (Fig. 2a). On the other hand, utilization of the proximal PA (pPA) removes the aUTR region and limits the size of the 3′ UTR to the constitutive UTR (cUTR), generating a short APA isoform (sAPA) (Fig. 2a). To study the extent of APA in the healing callus, we first performed transcriptomic analysis of the APA isoforms in the callus on d14 using the DaPars-based InPAS algorithm.17,18 Quantification of the relative expression of APA isoforms is based on calculating the Percent of Distal Utilization Index (PDUI) (Fig. 2a), which indicates the fraction of each analyzed transcript that is expressed as a lAPA isoform (Fig. 2a). The results indicated that ∼14% of the analyzed genes exclusively utilized the dPA (PDUI = 1; Group 1 in Fig. 2b; Table S4), while ∼2% exclusively utilized the pPA (PDUI = 0, Group 7 in Fig. 2b; Table S4). Importantly, ∼84% of the analyzed genes utilized both the dPA and pPA to variable extents (Groups 2-6 in Fig. 2b; Table S4). Hereafter, this group will be called APA genes/events. The same analysis was performed using APAlyzer, which utilizes a different approach to quantify APA isoform expression, and this analysis yielded comparable data, as APA events were identified in ∼78% of the genes, which showed ∼81% overlap with the APA events identified by InPAS (Fig. 2c and Table S5). Importantly, in ∼50% of the APA genes, the minor APA isoform constituted ≥20% of the total gene reads (Fig. 2b). GO analysis29 of the APA genes indicated significant enrichment in a wide range of biological pathways that play essential roles in fracture healing (Figs. 2d and S7 and Table S6). In fact, a large number of transcripts that are known to play roles in ossification, angiogenesis, osteoblast differentiation, endochondral ossification, and extracellular matrix deposition generated a PDUI of 40%–50%, suggesting substantial utilization of both the pPA and dPA (Table S7). These genes included the bone matrix components collagen, type I, alpha 1 (Col1a1) and alpha 2 (Col1a2); the angiogenic factor vascular endothelial growth factor A (Vegfa); the osteoclast marker tartrate resistant acid phosphatase 5 (Acp5; also known as TRAP); and the bone/cartilage development inducer bone morphogenic protein 5 (Bmp5) (Table S7). These results indicate that thousands of genes are expressed as multiple APA isoforms in the callus and suggest critical roles of APA in bone regeneration.

Fig. 2

APA is commonly found in transcripts that are expressed in the callus on day 14 postfracture. a Schematic of the short APA (sAPA) and long APA (lAPA) isoforms of the same gene. Utilization of the proximal polyadenylation site (pPA) confines the size of the 3′ UTR to the constitutive UTR (cUTR), generating the sAPA isoform. On the other hand, utilization of the distal PA (dPA) results in the inclusion of the alternative UTR (aUTR) of the 3′ UTR, generating the lAPA isoform. The InPAS-based calculation of PDUI values is provided. b Pie chart presentation of the distribution of the indicated 7 PDUI groups in the transcriptome of the d14 callus. c Presentation of the overlap between the APA events in the transcriptome of the d14 callus that were identified by InPAS and APAlyzer. Events identified by both programs are shown in red; events identified by InPAS alone are shown in blue. d GO analysis of APA events in d14 calli. Significantly enriched functional categories are shown, and the false discovery rate (FDR) of enrichment for each pathway is given (comprehensive lists of significantly enriched functional categories are provided in Fig. S7 and Table S6)

Utilization of the pPA generates shortened 3′ UTRs in the Col1a1 and Col1a2 transcripts

Type I Collagen, which is the main component of the organic bone matrix, is a heterotrimer comprising two α1 chains (encoded by Col1a1) and one α2 chain (encoded by Col1a2).30 According to our RNA-seq data analyses, Col1a1 and Col1a2 are the most abundant transcripts in the d14 callus, and both are expressed as 2 main APA isoforms (Fig. 3a and Tables S4, S5, S7). To confirm these results, we mapped the 3′ end of both the Col1a1 and Col1a2 mRNAs by performing rapid amplification of cDNA 3′ ends (3′ RACE) (Fig. S8a, b) using RNA purified from d14 calli. Indeed, we identified two APA isoforms of the Col1a1 and Col1a2 mRNAs (Fig. 3a–c) and mapped their 3′ ends to CPSs that were consistent with the bioinformatics predictions (Fig. S9a, b). Interestingly, we identified the same Col1a1 and Col1a2 APA isoforms when we performed 3′ RACE on RNA that was isolated from undifferentiated or differentiated murine MC-3T3-E1 osteoblastic cells (Fig. 3d–f).

Fig. 3

Long and short APA isoforms of the Col1a1 and Col1a2 mRNAs are expressed in the fracture callus as well as osteogenic and chondrogenic cell lines. a Schematic of the lAPA and sAPA isoforms of Col1a1 (left) and Col1a2 (right). The determined size of the 3′ UTR of each isoform is shown. b Representative gel image of Col1a1 3′ RACE-PCR products showing the 2 APA isoforms. 3′ RACE was performed using RNA purified from d14 calli (3′ RACE procedure is shown in Fig. S8a, b). The 3′ end of each APA isoform was mapped by excising and sequencing the corresponding PCR band (Fig. S9a). c As in b, except 3′ RACE-PCR was performed on Col1a2 mRNA. Sequencing of the PCR products confirmed the lAPA and sAPA isoforms that were identified by InPAS (Fig. S9b). d, e As in b and c, respectively, except that RNA was purified from the indicated cell lines. ATDC5 cells were differentiated under conditions that promoted hypertrophic differentiation and mineralization (Fig. S10). f RT‒qPCR analysis of the expression of the osteoblast and mineralization marker bone gamma carboxyglutamate protein (Bglap) in undifferentiated and differentiated MC3T3 cells. The mRNA expression level of Bglap was normalized to that of β-actin, and the normalized level in undifferentiated cells was defined as 100. g Alizarin red staining of ATDC5 cells before (−) and after (+) 3 weeks of differentiation. Positive staining of differentiated cells indicated chondrocyte mineralization (Fig. S10). h As in f except that RNA was purified from ATDC5 cells shown in g, and the specified transcript levels were quantified. The results demonstrate high expression of chondrocyte hypertrophy and mineralization markers in the differentiated ATDC5 cells. Gel images are representative of three independent replicates. The bar graphs present the average of three independent replicates ± SEM. **P < 0.01; ****P < 0.000 1 using unpaired Student’s t test

Abundant levels of type I collagen are also secreted by mineralizing chondrocytes (Fig. 1d).1,3,23 To investigate whether the expression of Col1a1 and Col1a2 APA isoforms is specific to bone cells or whether these isoforms are also expressed in hypertrophic chondrocytes, we differentiated ATDC5 chondrogenic cells under conditions that induced chondrocyte hypertrophy and mineralization (Figs. 3g, h and S10a, b). We detected the same APA isoforms in hypertrophic ATDC5 cells (Fig. 3d, e). Notably, we sequenced the 3′ RACE-PCR products of Col1a1 and Col1a2 in all these experiments and confirmed that the 3′ end of each of the lAPA and sAPA isoforms is conserved among the studied cell lines and calli (Fig. S9a, b). Accordingly, APA-mediated 3′ UTR shortening in different mineralized tissues and cell lines removes ∼80% and 60% of the 3′ UTRs of the Col1a1 and Col1a2 mRNAs, respectively.

3′ UTR shortening reverses the inhibition of Col1a1 and Col1a2 expression

To investigate the biological significance of the APA-mediated shortening of the Col1a1 and Col1a2 3′ UTRs, we performed mechanistic studies using MC-3T3 cells, which express APA isoforms of Col1a1 and Col1a2 that are identical to those expressed in the callus (Fig. 3b–e). We constructed a reporter containing the 3′ UTR of the Col1a1 lAPA isoform (lAPA 3′ UTR), which represents the full-length 3′ UTR, downstream of the firefly luciferase (FLuc) coding region (FLuc-lAPA 3′ UTR; Fig. 4a). As expected, when introduced into MC3T3 cells, this reporter expressed 2 APA isoforms of FLuc mRNA (Fig. 4b). The 3′ end of each of the 2 isoforms was identical to that of endogenous Col1a1 (Fig. S9a). Similarly, a reporter containing the lAPA 3′ UTR of Col1a2 mRNA expressed 2 APA isoforms of FLuc mRNA that were identical to those of endogenous Col1a2 (Figs. 4c and S9b). Accordingly, the FLuc-lAPA 3′ UTR reporters recapitulated the alternative cleavage and polyadenylation of the corresponding endogenous transcripts. We next examined the 3′ UTRs of both Col1a1 and Col1a2 mRNAs to define putative PAs that might modulate proximal cleavage and polyadenylation. In each transcript, we identified an A-rich simple sequence repeat (SSR) that contained overlapping canonical (AAUAAA) PAs that were located <50 nt upstream of the proximal CPS (Fig. S9a, b). Deletion of the SSR (FLuc-ΔpPA 3′ UTR reporter; Fig. 4a) almost completely abolished the 3′ UTR shortening of both the Col1a1 and Col1a2 reporters, resulting in the expression of each reporter as a single APA isoform with a full-length 3′ UTR (Fig. 4b, c). These results indicate that SSR-embedded pPAs are the main drivers of 3′ UTR shortening in both the Col1a1 and Col1a2 genes.

Fig. 4

Shortening of the 3′ UTR promotes the expression of Col1a1 and Col1a2. a Schematic of FLuc reporters that contain either the 3′ UTR of the lAPA isoform (i.e., full-length 3′ UTR) of Col1a1 or Col1a2 (FLuc-lAPA 3′ UTR) or the deletion mutant in which the SSR was deleted (FLuc-ΔpPA 3′ UTR; X denotes the position of the deleted SSR). CDS: coding region. b 3′ RACE-PCR of FLuc mRNA was performed on RNA that was isolated from MC-3T3 cells transfected with Col1a1 3′ UTR reporters shown in a. PCR was performed using an FLuc-specific forward primer. The results indicate that deletion of the SSR in the FLuc-ΔpPA 3′ UTR reporter abrogated 3′ UTR shortening, resulting in the expression of FLuc mRNA as a single, long APA isoform. c As in b, except Col1a2 3′ UTR reporters were used. d Schematic of FLuc reporters that contain different isoforms or regions of the Col1a1 or Col1a2 3′ UTR. e Bar graph of the results of the dual luciferase assays that were performed on MC3T3 cells transfected with FLuc-Col1a1 3′ UTR reporters shown in d. FLuc activity was normalized to Renilla luciferase (RLuc) activity, and the FLuc/RLuc ratio in cells transfected with the “empty” FLuc (FLuc-) reporter was defined as 100. f As in e except that FLuc-Col1a2 3′ UTR reporters were used. n = 4–5 independent transfections, each measured in triplicate. g Bar graph of the results of the pulse-chase experiments used to measure RNA decay. MC3T3 cells were incubated with 5-ethynyl uridine (EU), which is an analog of uridine, for 2 h to label the nascent RNA during the “pulse step”. Cells were then either collected to determine the total level of EU-labeled mRNA (EU mRNA) or grown for 24 h in EU-free medium during the “chase step”. At each timepoint, EU RNA was purified (see Materials and methods), and the level of each of the specified EU mRNA was determined using RT‒qPCR. Primers that bind to the aUTR were used to specifically quantify the expression of the lAPA isoform of Col1a1 or Col1a2 (Table S14). The level of EU mRNA measured after the pulse step was defined as 100%, and the graph shows the % remaining after the 24-h chase step. The results demonstrate faster loss of the Col1a1 and Col1a2 lAPA isoforms compared to the total mRNA. β-actin mRNA was used as a positive control for EU mRNA decay, and consistent with its reported half-life (∼6.6 h),50 <10% remained after the 24-h chase. The gel images are representative of three independent replicates. The bar graphs present the average ± SEM. (ns) P > 0.05; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.000 1 using one-way analysis of variance (ANOVA) followed by Tukey’s post hoc test for e and f, and unpaired Student’s t test for g

To investigate how 3′ UTR shortening impacts gene expression, we constructed an additional FLuc reporter containing the sAPA 3′ UTR of Col1a1 or Col1a2 (FLuc-sAPA 3′ UTR; Fig. 4d). We compared the reporter activities of the FLuc-ΔpPA 3′ UTR and FLuc-sAPA 3′ UTR constructs, and the results indicated strong inhibition of FLuc activity (>70%) by the lAPA 3′ UTR of either Col1a1 or Col1a2 (Fig. 4e, f), while the repression that was caused by the sAPA 3′ UTR of Col1a1 and Col1a2 was very moderate (∼25%) (Fig. 4e) and not significant (Fig. 4f), respectively. These results demonstrate that the aUTR region, which is removed during 3′ UTR shortening, exerts a potent inhibitory effect on the expression of Col1a1 and Col1a2. Consistent with this observation, a reporter that contains the aUTR region of Col1a1 or Col1a2 (FLuc-aUTR; Fig. 4d) repressed Fluc activity to the same extent as the full-length 3′ UTR (Fig. 4e, f). To corroborate the reporter results and investigate the impact of 3′ UTR shortening on endogenous transcripts, we performed pulse-chase experiments and assessed the stability of endogenous Col1a1 and Col1a2 APA isoforms. The results indicated substantially faster degradation of the lAPA isoforms than the sAPA isoforms (Fig. 4g). Thus, 3′ UTR shortening enhances Col1a1 and Col1a2 expression via, at least in part, the stabilization of the transcripts.

Shortening of the 3′ UTR enhances Col1a1 and Col1a2 expression by removing miR-29a-3p binding sites

To define the inhibitory elements that are embedded within the aUTR of Col1a1 mRNA, we searched for putative miRNA response elements (MREs) among other factors. Among the top scoring miRNAs with conserved MREs in Col1a1 mRNA, four were identified by both TargetScan31 and miRDB;32 these four miRNAs included three members of the miR-29 family and miR-6980-5p (Fig. 5a and Tables S8 and S9). In fact, the miR-29 family members have been reported to inhibit the expression of Col1a1 with variable efficiency in different tissues.33,34,35,36. miRNA sequencing (miR-seq) of small RNAs that were isolated from d14 calli identified miR-29a-3p as one of the 15 most highly expressed miRNAs, generating ∼20 000 counts per million (CPM) (Table S10). The expression level of miR-29a-3p was ∼8-fold and 34-fold higher than that of miR-29c-3p and miR-29b-3p, respectively (Table S10). On the other hand, miR-6980-5p was expressed at very low levels, generating <10 CPM (Table S10). Thus, we focused our subsequent studies on miR-29a-3p. In situ hybridization confirmed the high expression of miR-29a-3p in woven bone osteoblasts and lining cells (Fig. 5b). Interestingly, miR-29a-3p has three MREs in the 3′ UTR of Col1a1 mRNA, all of which are in the aUTR (Figs. 5c and S11a and Tables S8 and S9). miR-29a-3p also has a single MRE in the aUTR of Col1a2 mRNA (Figs. 5d and S11b and Tables S8 and S9), suggesting that miR-29a-3p suppresses the expression of both Col1a1 and Col1a2 via MREs embedded in the aUTR. To further test this hypothesis, we used the same FLuc reporters that we constructed to express different isoforms and to contain different regions of the Col1a1 or Col1a2 3′ UTR (Fig. 5e, f). A miR-29a-3p mimic effectively repressed the reporters that contained the ΔpPA 3′ UTR or aUTR of either Col1a1 or Col1a2 by ≥50% (Fig. 5e, f). In contrast, the miR-29a-3p mimic had no effect on the sAPA reporters (Fig. 5e, f). Mutating the single MRE in the Col1a2 aUTR (FLuc-ΔMRE) abolished the inhibitory effect of miR-29a-3p (Fig. 5f). These data confirm that miR-29a-3p effectively inhibits the expression of Col1a1 and Col1a2 by binding to the aUTR. To corroborate the FLuc reporter data, we examined the effect of the miR-29a-3p mimic on endogenous transcripts. The mimic downregulated the total level of endogenous Col1a1 mRNA by ∼30% but exerted substantially stronger inhibitory effects on the lAPA isoform, decreasing its expression by ∼75% (Fig. 5g). Comparable results were obtained with endogenous Col1a2 mRNA (Fig. 5h). Consistent with this, although a miR-29a-3p inhibitor did not significantly change the total level of endogenous Col1a1 or Col1a2 mRNA, it increased the expression level of the lAPA isoforms of Col1a1 and Col1a2 by ∼1.8- and 2.3-fold, respectively (Fig. S12a, b). Taken together, these results highlight the role of 3′ UTR shortening in reversing the inhibition of Col1a1 and Col1a2 expression in biological contexts where miR-29a-3p, and possibly other members of the miR-29 family, are highly expressed.

Fig. 5

miR-29a-3p mediates the degradation of the lAPA, but not sAPA, isoforms of Col1a1 and Col1a2 mRNAs. a List of miRNAs that have conserved MREs in the 3′ UTR of mouse Col1a1 mRNA (identified by both TargetScan and miRDB) (Tables S8 and S9). miR-29a-3p (highlighted in red) exhibited the highest expression among the 4 miRNAs in the d14 callus (Table S10). b ISH using either a control (left) or miR-29a-3p (right) probe. The results show intense staining (blue) of miR-29a-3p in the woven bone area of the d14 callus. c Schematic of the lAPA isoform of Col1a1 mRNA showing the putative miR-29a-3p MREs in the aUTR. d As in c except Col1a2 mRNA is shown. e Bar graph of the results of dual luciferase assays that were performed in MC3T3 cells transfected with the indicated Col1a1 FLuc reporters along with either a control (blue) or miR-29a-3p (red) mimic. FLuc activity was normalized to RLuc activity, and the FLuc/RLuc ratio in cells transfected with the control mimic was defined as 100. n = 4–5 independent transfections, each measured in triplicate. f As in e except the cells were transfected with Col1a2 FLuc reporters. In the FLuc-ΔMRE reporter, we replaced the binding region of the miR-29a-3p seed sequence with a random sequence (denoted by “X”). g RT‒qPCR quantification of the expression of either total mRNA or the lAPA isoform of Col1a1 in MC3T3 cells transfected with control (blue) or miR-29a-3p (red) mimic. Primers that bind within the aUTR region were used to specifically amplify and quantify the levels of lAPA isoform (Table S14). The expression level was normalized to that of β-actin mRNA, and the normalized level in control mimic-transfected cells was defined as 100. n = 5. h As in g except that the expression levels of total mRNA and the lAPA isoform of Col1a2 were quantified. The bar graphs present the average ± SEM. (ns) P > 0.05; **P < 0.01; ***P < 0.001; ****P < 0.000 1. The significance of the difference was calculated using unpaired Student’s t test in e and f and two-way ANOVA in g and h

Endochondral ossification is accompanied by prevalent 3′ UTR shortening

Our RNA-seq data indicate that the progression of endochondral ossification is associated with increased cell proliferation (Figs. 1i and S5 and S6 and Tables S2 and S3). As rapidly proliferating cells exhibit widespread 3′ UTR shortening,