Comparative transcriptomics reveals the gene regulatory network underlying the differentiation and evolution of skeletal cells

Abstract

Deciphering the gene regulatory network (GRN) underlying skeletal cells is key to understanding the differentiation and even origins of cartilage and bone, but current knowledge is limited to a few candidate genes. The vertebrate skeleton is mostly composed of three specific cell types: immature chondrocytes (IMM), mature (hypertrophic) chondrocytes (MAT), and osteoblasts (OST). Importantly, mature chondrocytes share the expression of many genes with both immature chondrocytes and osteoblasts. Currently little is known about mechanisms of gene regulation in mature chondrocytes, but previous studies suggest that overlapping actions between portions of the GRN active in IMM and OST direct the differentiation of MAT. While one GRN might give rise to distinct cell fates, studies analyzing the interaction between distinct portions of a GRN during cell differentiation are lacking. Immature chondrocytes, mature chondrocytes, and osteoblasts can have distinct embryological origins, mesoderm and neural crest, but the similarities in gene expression among skeletal cells of these two embryonic lineages remain controversial. Here, we aim to test the hypothesis that the skeletal cell GRN is generally conserved throughout the body and across vertebrate clades, but the molecular relationship between portions of this GRN active in IMM and OST has been modified during evolution (i.e. positive in earlier diverged vertebrates vs antagonistic in later diverged vertebrates). To test this hypothesis, laser- capture microdissection (LCM) coupled to RNA-seq was performed on skeletal cells (i.e. immature chondrocytes, mature chondrocytes, and osteoblasts) of mouse, chick, and gar. The hypothesis was partially supported by our results. First, transcriptomic analyses on a single species, the mouse embryo, suggested that one GRN drives the differentiation of IMM, MAT, and OST, and the overlapping actions of portions of this GRN active in IMM and OST regulate MAT differentiation. Second, to test whether this skeletal cell GRN is conserved throughout the body, limb (humerus, mesoderm-derived) and head cartilage (ceratobranchial, neural-crest derived) isolated from the chick embryo were compared and their transcriptomes showed a high degree of similarity. Third, regarding bone evolution, previous molecular studies suggest that gar OST express higher levels of chondrocyte genes (e.g. Sox9, Col2a1, and Col10a1) compared to mouse and chick OST, suggesting that the osteoblast might have evolved from a chondrocyte. To further test this hypothesis, several bioinformatic approaches including differential gene expression, model-based clustering, and gene ontology analyses were performed. Pairwise differential gene expression analyses revealed that gar OST indeed expressed higher levels of some hallmark chondrogenic markers including Col2a1, Sox6, Col10a1, and Acan compared to mouse or chick OST. Moreover, model-based clustering analysis revealed one cluster that showed higher expression in the gar OST and included the hallmark mature chondrocyte gene Col10a1. Finally, with the goal of understanding how changes in GRNs underlie skeletal cell evolution, gene co-expression network (GCN) analysis was used to estimate skeletal cell GRNs. In mouse and chick, two portions of the GRN driving IMM and OST exhibit an antagonistic relationship, and they only interacted positively in MAT. In contrast, GCN analysis in gar skeletal cells showed that positive interactions between IMM and OST increased, supporting the hypothesis that the relationship between portions of the GRN active in IMM and OST has been modified over time. To provide insight into how portions of the GRN might regulate gene expression in MAT, model-based cluster analysis was performed, and specific categories of gene expression were identified. These categories of gene expression were conserved to some degree in skeletal cell transcriptomes of all vertebrate clades, but when gar skeletal cell transcriptomes were analyzed more clusters showed downregulation of gene expression levels in MAT compared to IMM and OST, while IMM and OST exhibited similar expression levels. This result in gar was particularly unexpected since results in mouse and chick showed that IMM and OST generally show opposite gene expression patterns, and MAT is mostly a combination of gene expression levels between both cell types. Adding more phylogenetic clades into these evolutionary comparisons will provide more insight into skeletal cell GRN regulation and structure. These results highlight the complexity of skeletal cell GRN organization and propose a novel unbiased approach through which to understand differentiation and evolutionary origins of cartilage and bone.