IMM and MAT were obtained from two chondral bones, the humerus and ceratobranchial (Figure 1), which are derived from lateral plate mesoderm and cranial neural crest, respectively (Noden 1982; Couly et al., 1993; Le Douarin et al., 1999; Burke and Nowicki 2001; Tani et al., 2020). In chick, the epiphyseal growth plate of such long bones as the humerus and ceratobranchial should contain IMM and MAT at HH36 (E10; Milz et al., 2002; Conen et al., 2009; Lui et al., 2014). Alcian blue identified cartilage, whereas Alizarin Red identified perichondral bone in whole-mount stains of both skeletal elements at HH36 (Figures 1A,B,J). Perichondral bone is often associated with underlying MAT, and Safranin O staining of histological sections demonstrated that MAT, underneath perichondral bone, have undergone hypertrophy in the mid-diaphyseal region (i.e., shaft of a long bone) of the humerus and ceratobranchial (Figures 1C–E,K–M). LCM was used to isolate IMM and MAT from the humerus and ceratobranchial (Figures 1F–I,N–Q). At HH36, vascular invasion had not yet occurred, and bone matrix was only detected in perichondral bone in both the humerus and ceratobranchial, as shown by Trichrome staining (Figures 1R,S), suggesting that no transdifferentiation of chondrocytes to osteoblasts had yet occurred (Conen et al., 2009; Zhou et al., 2014; Qin et al., 2020). RNA-seq was then carried out on RNA isolated from IMM and MAT.

To test the hypothesis that a core GRN is expressed during chondrocyte differentiation in the limb and head, chondrocyte transcriptomes from the HH36 chick humerus and ceratobranchial were compared (Figure 2; Supplementary Figure S1). To identify similarities and differences among limb and head chondrocyte transcriptomes, a principal component analysis (PCA) was performed (Supplementary Figure S1). The variation in the samples was captured with two components (39% variance explained by PC1 and PC2; Supplementary Figure S1). The limb IMM and MAT transcriptomes were separated from other samples in PC1/PC2 with 95% confidence, while head IMM and MAT transcriptomes overlapped in PC1/PC2 (Supplementary Figure S1).

When IMM and MAT datasets from limb and head were combined before normalization (i.e., IMM and MAT were considered as the same cell type) to reflect generally the chick “chondrocyte”, limb and head chondrocytes shared 77% of genes expressed above threshold (6190/8034 genes; Figure 2A; see normalization techniques in Methods). Limb and head chondrocytes expressed typical IMM genes, such as SOX9, SOX5, SOX6, ACAN, and COL2A1, and typical MAT genes, including RUNX2, COL10A1, MEF2C, MMP13, IBSP, and SPP1 at high levels (Tables 1, 2; Vortkamp et al., 1996; Zhao et al., 1997; Bridgewater et al., 1998; Watanabe et al., 1998; Bi et al., 1999; Smits et al., 2001; Smits et al., 2004; Arnold et al., 2007; Nicolae et al., 2007; Dy et al., 2012; Lu et al., 2014; Nakatani and Partridge 2017).

When IMM and MAT were each compared separately between the limb and head before normalization (i.e., IMM and MAT were considered as different cell types), each cell type still shared the vast majority of genes expressed above threshold in both the limb and head (for IMM: 5249/7648 = 69%; for MAT: 8209/10,217 = 80%; Figure 2B). Overlapping genes in the limb and head again included many typical cartilage genes, such as SOX9, COL2A1, ACAN, and COL9A1 for IMM, and RUNX2, COL10A1, SPP1, MMP13, and IBSP for MAT (Tables 1, 2). Gene ontogeny analyses on IMM and MAT from both limb and head demonstrated that cartilage-specific processes were enriched and conserved between the limb and head, even though genes associated exclusively with these processes only comprise approximately 2% of the GO term-associated genes (Figure 2C). The most enriched biological processes (>60% genes expressed above threshold) were related to basic cellular processes, such as cell proliferation, cell differentiation, transcription, and translation. Collectively, these data suggest that a core GRN driving chondrocyte differentiation is expressed in different regions of the skeleton.

Using published work on mouse, chick, frog, and fish, regulatory interactions among important genes from chick chondrocyte transcriptomes were summarized into a GRN (Figure 2D; Longabaugh et al., 2005). As master regulators of IMM and MAT, respectively, SOX9 and RUNX2 were placed at the top of the GRN hierarchy (Figure 2D; Komori et al., 1997; Bi et al., 1999; Lian and Stein 2003; Eames et al., 2004). In the IMM portion of the GRN (highlighted in pink in Figure 2D), SOX9 binds to SOX5 and SOX6 during early stages of chondrocyte differentiation and activates the expression of typical IMM markers, such as COL2A1, COL9A1, and ACAN (Lefebvre et al., 2001; Akiyama et al., 2002; Liu and Lefebvre 2015). In the MAT portion of the GRN (highlighted in yellow in Figure 2D), RUNX2 activates the expression of typical MAT markers, including COL10A1, MMP13, SPP1, and IBSP (Ducy et al., 1997; Komori et al., 1997; Inada et al., 1999; Inada et al., 1999; Lee et al., 2000; Wang et al., 2009; Li et al., 2011; Peacock et al., 2011; Lu et al., 2014). Since SOX9 and RUNX2 generally exhibit an antagonistic relationship, SOX9 inhibits the expression of MAT markers, such as RUNX2, SPP1, and IBSP, while RUNX2 inhibits the expression of SOX9 and thus likely other IMM markers (Figure 2D; Zhou et al., 2006; Cheng and Genever 2010; Peacock et al., 2011; Lui et al., 2019).

Genes expressed above threshold only in limb or head might reflect the embryonic origin of the cells. For genes located in the humerus portion of the Venn diagrams, such as TBX5, DLX6, SALL4, HOXA10, and HOXD10, embryonic limb/forelimb morphogenesis was an enriched process, and these genes are all known regulators of limb development (Figure 2B; Wahba et al., 2001; Robledo et al., 2002; Zakany and Duboule 2007; Vieux-Rochas et al., 2013; Neufeld et al., 2014; Akiyama et al., 2015). More general skeletal processes, such as skeletal system development and ossification, were also enriched in limb chondrocytes. Many neural crest-dependent biological processes were enriched specifically in ceratobranchial IMM and MAT transcriptomes, such as cranial skeleton morphogenesis, cell migration, neuron differentiation, middle ear morphogenesis, heart development, and melanocyte differentiation (Figure 2B). Many orthologs of a proposed GRN driving neural crest-derived cartilage, such as PAX7, SIX1, and ID1, were also identified (Figure 2B; Meulemans and Bronner-Fraser 2005; Meulemans and Bronner-Fraser 2007; Betancur et al., 2010; Murdoch et al., 2012; Wu et al., 2019).

GRNs rely upon functional data to verify regulatory interactions, but such studies are limited (Figure 2D; Su et al., 2009; Peter and Davidson 2011). To infer regulatory interactions underlying chondrocyte differentiation, and to compare GRN organization in the limb and head, GRNs were estimated using Cytoscape to graph gene co-expression networks (GCNs) of chondrocyte transcriptomes of the HH36 humerus or ceratobranchial (McCall 2013; Khosravi et al., 2015). All genes expressed above threshold were used to construct these GCNs (>8,000 genes, Figures 2A,B). GRNs of both limb and head data were organized into two large groups of genes expressed during chondrocyte differentiation (Figures 3A,B). Expression within each group was positively correlated (red lines in Figures 3A,B), but expression between the two groups was negatively correlated (blue lines in Figures 3A,B). One group was enriched for genes that were differentially expressed in IMM, and the other group was enriched for genes that were differentially expressed in MAT (for ease of view, selected IMM and MAT differentially expressed genes of limb and head chondrocytes are depicted in Figures 3C,D, respectively).

In the limb, these two portions of the GRN included 859 genes [absolute log2 fold change greater than 2 (p < 0.01)] that were differentially expressed between IMM and MAT of the HH36 chick humerus (see Supplementary Tables S1, S2 for a full list of genes). A total of 263 genes were upregulated in IMM, whereas 596 genes were upregulated in MAT (Figure 3A; Supplementary Tables S1, S2). Upregulated genes in IMM included MATN3, GLI2, GLI3, DDR2, and NOG, and also some HOX genes that have a role during chondrocyte differentiation in the limb (Figure 3A, labelled pink in GRN; Koziel et al., 2005; Kruger and Kappen 2010; Tan et al., 2018; Yamamoto et al., 2019). Typical maturation genes were upregulated in MAT, including COL10A1, MMP13, SPP1, MEF2C, IBSP, and SPARC (Figure 3A, labelled yellow in GRN; Bianco et al., 1991; D'Angelo et al., 2000; Arnold et al., 2007; Peacock et al., 2011; Lu et al., 2014; Rosset and Bradshaw 2016). These differences in gene expression patterns were also demonstrated by unsupervised model-based clustering analysis (Supplementary Table S2). Some clusters showed enriched expression of genes in IMM including hallmark cartilage genes, while others showed enhanced expression in MAT including several important maturation markers (Supplementary Table S2). Typical chondrocyte differentiation genes, such as SOX9, SOX5, SOX6, COL2A1, and ACAN showed high expression levels in both IMM and MAT of the HH36 chick humerus (Figure 3A, labelled blue in GRN; Tables 1, 2).

In the head, the two portions of the GRN included 118 genes that were differentially expressed between IMM and MAT of the HH36 chick ceratobranchial (see Supplementary Tables S3, S4 for a full list of genes). A total of 70 genes were upregulated in IMM, whereas 48 genes were upregulated in MAT (Figure 3B; Supplementary Tables S3, S4). Upregulated genes in IMM included COMP and MMP16 (Figure 3B, labelled pink in GRN; Hecht et al., 2005). Typical maturation genes were upregulated in MAT, including RUNX2, COL10A1, IBSP, and SPP1 (Figure 3B, labelled yellow in GRN; Bianco et al., 1991; Komori and Kishimoto 1998; Arnold et al., 2007; Peacock et al., 2011; Lu et al., 2014).

To directly visualize an estimate of the general chick chondrocyte GRN, limb and head data were combined before normalization and graphed as a GCN. The overall organization of this GRN was similar to that for just limb or head chondrocyte data. Two portions of positively correlated genes that were differentially expressed in each type of chondrocyte (i.e., IMM or MAT) were negatively correlated with each other (Figures 4A,B). A total of 458 genes were differentially expressed between IMM and MAT (see Supplementary Tables S5, S6 for a full list of genes). Of these, 195 genes were upregulated in IMM, and 263 genes were upregulated in MAT (Supplementary Tables S5, S6). Upregulated genes in IMM and MAT included many of the same genes that were upregulated in either limb or head data alone (compare genes labelled pink in Figures 4A,B to Figures 3A,B). Other upregulated IMM genes were not previously implicated in chondrocyte differentiation, such as those encoding the transcription factors POU3F3 and ZFP703, growth factors CSF1 and PTN, and the Netrin receptor UNC5C (Figures 4A,B, labelled pink; Cecchini et al., 1997; Tare et al., 2002; Srivatsa et al., 2014; Kumar A. et al., 2016; Kumar S. et al., 2016). Upregulated genes in MAT included IFITM5, whose role had only been studied previously in osteoblast differentiation, and many other signalling pathway genes previously implicated in cartilage maturation, including BMP4, FGF9, HES5, TEK, TGFBR2, TGFBR3, WNT5B, and WNT11 (Figures 4A,B, labelled yellow; Hoffmann and Gross 2001; Yang et al., 2003; Karlsson et al., 2010; Shen et al., 2013; Usami et al., 2016; Zhang et al., 2021). Additional genes upregulated in MAT had never been associated previously with chondrocyte maturation, including those encoding the kinase FAM20C and transcription factor TCF7L2 (Hirose et al., 2020; Li et al., 2021).

Although typical chondrocyte differentiation genes, such as SOX9 and COL2A1, were not differentially expressed between limb and head (Figures 3, 4), differences in gene expression were obvious when the normalized counts were compared (Tables 1, 2). In the limb, for instance, the SOX trio (SOX9, SOX5, and SOX6) showed higher expression levels in IMM, and expression levels decreased in MAT, as reported by others (Ikegami et al., 2011; Lui et al., 2019). In the head, however, expression levels of SOX9 and SOX6 remained at comparable levels in IMM and MAT, and MAT even showed higher levels of SOX5 expression compared to IMM. A similar situation is observed when RUNX2 levels are compared between limb and head. In the head, RUNX2 expression strikingly increases during IMM to MAT transition (Tables 1, 2; Figure 3B), while in the limb this increase in RUNX2 expression is not as as dramatic (Tables 1, 2). These differences in the expression might be the result of variation in the timing of cartilage maturation among skeletal elements, as well as differences related to growth and shape of endochondral bones in distinct locations of the body (Chiba et al., 1995; Patton and Kaufman 1995; Mitgutsch et al., 2011). Together, these data demonstrated that the GRN driving chondrocyte differentiation is organized similarly throughout the body and provided novel genes that might regulate IMM and MAT differentiation.

To identify genes that might influence which type of chondrocyte differentiates in specific embryonic regions, IMM or MAT data from the limb and head was used to estimate a GRN underlying IMM or MAT differentiation by graphing a GCN. Two portions of each GRN were identified, this time enriched for genes that were differentially expressed in chondrocytes of the limb or the head (Figures 5A,B). Many of these region-specific genes were positively correlated to each other, while negatively correlated to genes enriched in chondrocytes from the other embryonic region (Figures 5A–D). Many of these genes were patterning genes of the limb or head.

Genes upregulated in limb chondrocytes were very similar for the IMM and MAT data, while genes upregulated in head chondrocytes varied among IMM and MAT data. The two portions of the IMM GRN included 1223 genes that were differentially expressed between limb and head chondrocytes (Figure 5A; see Supplementary Tables S7, S8 for a full list of genes). A total of 216 genes were upregulated in limb IMM, while 1,007 genes were upregulated in head IMM (Supplementary Tables S7, S8). The two portions of the MAT GRN included 1,215 genes that were differentially expressed between limb and head chondrocytes (Figure 5B; Supplementary Tables S9, S10). A total of 664 genes were upregulated in limb MAT, whereas 551 genes were upregulated in head MAT (Supplementary Tables S9, S10). Upregulated genes in limb IMM and MAT included several HOX genes, DLX5, TBX5, and SHOX2, all of which are known to have a role during limb morphogenesis (Figures 5A,B, labelled green in GRN; Agarwal et al., 2003; Koziel et al., 2005; Yu et al., 2007; Kruger and Kappen 2010; Vieux-Rochas et al., 2013; Tan et al., 2018; Yamamoto et al., 2019). The other portion of the IMM GRN included genes upregulated in the head that are involved in neural crest related processes, such as ZIC1, SIX1, SIX4, PAX7, ID1, GATA2, GATA3, and MSX2, as identified by gene ontology analyses and previous studies (Figure 5A, labelled orange in the GRN; Han et al., 2007; Murdoch et al., 2012; Simões-Costa and Bronner 2013; Garcez et al., 2014; Plouhinec et al., 2014; Simões-Costa and Bronner 2015; Wu et al., 2019; Seal and Monsoro-Burq 2020). Similar to the IMM GRN, SIX1, and SIX4 were also upregulated in the head in one portion of the MAT GRN. Also upregulated in head MAT were TBX1 and PRRX1, which actually can regulate both limb and cranial neural crest (Figure 5B, labelled orange in GRN; Moraes et al., 2005; Balic et al., 2009; Simões-Costa and Bronner 2015). Similar results were obtained when IMM and MAT data from the limb and head were considered as four separate datasets before normalization, and included in the same GRN (Supplementary Figure S3). Unsupervised model-based clustering analyses also identified specific categories of gene expression when comparing limb and head chondrocyte data (Supplementary Figure S2). Some clusters showed enriched expression of genes in the limb, including several classic limb patterning genes, while other clusters showed increased expression in the head, including many genes involved in cranial neural crest differentiation (Supplementary Figure S2). Again, comparisons between limb and head transcriptomes revealed that typical IMM and MAT genes, including SOX9, COL2A1, RUNX2, COL10A1, and SPP1, were not differentially expressed between head and limb, supporting the hypothesis that a core transcriptional program driving chondrocyte differentiation is expressed wherever cartilage forms in the body.

All animal procedures were performed according to guidelines approved by the University of Saskatchewan Animal Care and Use Committee. White leghorn chicken eggs were incubated in a humified incubator at a constant temperature of 37°C. Embryos were harvested at Hamburger-Hamilton stage 36 (∼E10; Hamburger and Hamilton 1951). Each embryo was decapitated, and the forelimbs and lower jaws were dissected and either fixed in 4% paraformaldehyde overnight or immediately placed in 1X PBS/DEPC, followed by embedding in OCT (Tissue Tek, Torrance, CA, United States), and immediately flash-frozen using liquid N2 and 2-Methylbutane (isopentane).

Chick HH36 embryos were stained with Alcian blue and Alizarin red using an acid-free solution that included MgCl₂ to differentiate staining, and then cleared in glycerol/KOH as described elsewhere (Eames et al., 2011). Importantly, in our hands, various sources of Alcian blue do not work with the acid-free protocol, but one that does is from Acros Organics (Alcian Blue 8GX). Safranin O/Fast Green staining was performed on 10 μm thick frozen sections of the HH36 chick humeri and ceratobranchial, as described previously (Ferguson et al., 1998). Trichrome staining was performed on 10 μm thick frozen sections of the HH36 chick mandible, as described elsewhere (Ashique et al., 2022).

LCM was performed on a Laser Microdissection—Molecular Machines & Industries (MMI) CellCut apparatus. Immature and mature chondrocytes were captured from developing chick HH36 humeri (IMM, n = 4; MAT, n = 5) and ceratobranchial (IMM, n = 3; MAT, n = 3). At this early stage of development, perichondral bone was evident in both the humerus and ceratobranchial, but no osteoblasts or other cell types were present in the cartilage template, and invasion by the vasculature had not yet occurred. Tissue slices were processed and sequenced individually (not pooled at any stage), and the captured cells were collected onto the inner lid of 0.5 ml MMI IsolationCaps (either Diffuser caps (Prod#50202) or Transparent caps (Prod#50204; MMI Molecular Machines & Industries).

RNA was isolated from laser-captured tissue using the ARCTURUS PicoPure RNA Isolation Kit (ThermoFisher Scientific; Cat# KIT0204), according to the manufacturer’s instructions, and DNase treatment was done using RNase-Free DNase (Qiagen; Cat#79254). RNA was then amplified with one round using MessageAmp II aRNA Kit (ThermoFisher Scientific; Cat# AM1751). RNA integrity was evaluated on the observation of a signature eletropherogram pattern (Bioanalyzer).

RNA-seq libraries were prepared by the National Research Council (NRC, Saskatoon) using the Illumina TruSeq RNA Sample Prep Kit v2 with the following modification: the protocol was started at the Elute, Prime, Fragment step using 5 µl amplified mRNA [minimum amount was 5–10 ng mRNA as determined using Quant-iT RiboGreen RNA Assay Kit (Invitrogen)]. The quality of each cDNA library was checked on a DNA 1,000 chip using the 2,100 Bioanalyzer (Agilent Technologies Inc.). In average, the sequencing depth was 21-million reads per sample (min. 14-million reads per sample; max. 31- million reads per sample).

The paired-end Illumina reads were trimmed using a Java -based tool, Trimmomatic v0.30 (Bolger et al., 2014), and the reads were then mapped to the chicken genome on Ensembl using STAR v 2.5.2 (Dobin et al., 2012). The location of each read was matched to genome annotation using HTSeq-count (Anders et al., 2015). The distribution of average log₂ expression across three replicates of each tissue produced three bimodal distributions, which were used to set the count thresholds to 142 and 23 for IMM and MAT isolated from the ceratobranchial, and 37 and 58 for IMM and MAT isolated from the humerus. When head and limb data were combined before normalization, thresholds were set to 64 and 35 for IMM and MAT, respectively. Including different sets of samples before normalization affects the number of genes expressed above threshold, because the gene counts will change depending on the exact samples they are normalized against. Differential expression analysis was performed using EdgeR after excluding genes with zero or very low counts (less than three counts for all cell types) across the cell type. Pairwise comparisons between tissues were made with Fisher’s exact test, and a gene was considered differentially expressed if it had an absolute log₂ fold change greater than 2 (p < 0.01). Venn diagrams were constructed using gplots v3.0.1 for isoforms and RNA-seq expression data.

To evaluate similarities and differences among the IMM and MAT transcriptomes obtained from limb and head chondrocytes, a principal component analysis (PCA) was performed. PCA was performed on the data using prcomp from the stats library in R to determine if the biological replicates of each cell type separated into distinct groups based on gene expression variance. The 95% confidence ellipses were included using R package factoextra version 1.0.7. The variation in the samples was captured with two components (39% variance explained by PC1 and PC2; Supplementary Figure S1).

The algorithms from MBCluster.Seq 1.0 package in R were used to cluster the genes from our RNA-seq data (Si et al., 2014). Genes were assigned to 10 clusters based on expression profiles across all IMM and MAT isolated from limb and head (Supplementary Figure S2).

The skeletal cell GRN was constructed using BioTapestry version 7.1.2 (www.BioTapestry.org/) following developer’s protocol (Longabaugh et al., 2005). Regulatory interactions were validated based on published studies including genetic molecular experiments and cis-regulatory analyses performed in bones of mouse, chick, frog, and zebrafish (Ducy et al., 1997; Komori et al., 1997; Lecanda et al., 1997; Zhao et al., 1997; Bridgewater et al., 1998; Drissi et al., 2000; Sekiya et al., 2000; Akiyama et al., 2002; Zhang et al., 2003; Zheng et al., 2003; Akiyama et al., 2004; Bastepe et al., 2004; Stock et al., 2004; Yang et al., 2004; Hong et al., 2005; Magee et al., 2005; Meech et al., 2005; Yagi et al., 2005; Arnold et al., 2007; Holleville et al., 2007; Liu et al., 2007; Yun and Im 2007; Grogan et al., 2008; Higashihori et al., 2008; Vincourt et al., 2008; Teplyuk et al., 2009; Wang et al., 2009; Dao et al., 2010; Fazenda et al., 2010; Higashiyama et al., 2010; Leung et al., 2011; Nagy et al., 2011; Peacock et al., 2011; Gu et al., 2012; Ionescu et al., 2012; McGee-Lawrence et al., 2013; Oh et al., 2014; Liu and Lefebvre 2015; Ohba et al., 2015; Heilig et al., 2016; Takegami et al., 2016; Watanabe et al., 2016; Komori 2017; Yao et al., 2017; Kawane et al., 2018; Komori 2018; Liu et al., 2018; Qin et al., 2018; Tan et al., 2018; Xu et al., 2018; Komori 2019; Kurakazu et al., 2019; Mokuda et al., 2019; Yamashita et al., 2019; Wuelling et al., 2020).

DAVID v6.8 (http://david.abcc.ncifcrf.gov/home.jsp) functional annotation analysis was performed on genes expressed above threshold in head and limb. The GO term biological process (BP) in DAVID was used to perform the gene-annotation enrichment analysis.

For constructing GCNs, the Pearson correlation between genes was calculated using the TMM normalized gene expression data. Thresholding the edge weights (+/−0.85) was then performed to remove potentially irrelevant edge weights before visualization. All processing and normalization of the RNA-seq counts were performed using the edgeR package using R version 4.0.0. GCNs were visualized in Cytoscape version 3.8.2 and all color coding of edges and nodes was performed using Cytoscape.