Cartilage-selective genes identified in genome-scale analysis of non-cartilage and cartilage gene expression

Skeletogenesis begins with condensation of mesenchymal chondroprogenitor cells to form the cartilage anlagen that pattern the early skeleton. Subsequently, for bones that grow by endochondral ossification, the chondrocytes differentiate further to establish the growth plates. At the joint surfaces, development of articular cartilage facilitates and maintains joint movement during fetal life. These multi-step processes require the coordinated expression of many genes, including genes encoding extracellular matrix proteins and morphogens, as well as proliferative, angiogenic, and apoptotic signals [1]. Most of our knowledge of the function of the genes involved has been derived from developmental studies in model systems and cell lines [2] as well as from the identification of disease genes in skeletal disorders.

Whole genome analysis of chondrocyte gene expression has the potential to reveal novel genes and gene expression programs which define the tissue. Although the complete set of genes expressed in human cartilage has not yet been described, analysis of human cartilage cDNA libraries has provided an initial in vivo picture of the cartilage transcriptome [3–6]. These investigations have also identified expression of both known and novel genes. Comparative microarray studies in rat cartilage [7] and several chondrocyte cell lines [8, 9] have provided a larger set of genes of potential importance in chondrocytes, including genes specific to the stages of chondrocyte differentiation. Wang et al. (2004) identified 92 genes with two-fold variation in expression between hypertrophic and proliferative growth plate chondrocytes. In this in vivo study, significant gene expression changes were principally associated with cell cycle, transcription, extracellular matrix structure, receptor and transporter functions. In microarray studies of mouse micromass cultures [8], 212 genes exhibited at least a ten-fold difference in gene expression as the cultures differentiated. Thus global characterization of gene expression is beginning to describe the identities of key regulatory molecules and their targets in chondrocytes.

Disrupting genes involved in the organization and maturation of the growth plate and/or the stability of articular cartilage results in inherited skeletal disorders that range from perinatal lethal phenotypes to mild disorders with early-onset osteoarthropathy as their major feature [10, 11]. Of the approximately 370 clinically distinguishable skeletal dysplasias [12], mutations in 115 genes have been associated with about 150 disorders. Many of these disease genes are expressed in a cartilage-selective pattern, and therefore identifying additional genes uniquely expressed in cartilage should yield new skeletal dysplasia candidate genes.

To identify a larger set of genes uniquely expressed in chondrocytes, this study describes a genome-scale gene expression profile for 18–22 week human fetal cartilage. There were 161 genes which appeared to be selectively expressed in fetal cartilage, comprising a variety of novel genes that may contribute to skeletal development. The data suggest a complex pattern of cartilage gene expression and indicate that the number of genes selectively expressed in cartilage has been greatly underestimated.

Using genome-scale microarrays, gene expression in human fetal cartilage was compared with a robust set of other normal tissues. Hierarchical clustering showed remarkable similarity among the 18–22 week fetal cartilage expression profiles and demonstrated that a subset of the cartilage transcriptome is composed of a unique gene set not generally expressed in the other tissues studied. Using SAM, 2446 probesets measured preferential expression of 1712 genes with at least three-fold higher expression in cartilage as compared with other tissues. 1028 (42%) of these probesets matched genes identified in a cartilage growth plate cDNA library [4] validating their expression in cartilage via an independent dataset. The identification of genes known to have restricted patterns of expression in cartilage confirmed the presence of RNA derived from the reserve (GREM1), hypertrophic (BMP6, COL10A1), and terminally differentiated (MMP13) chondrocytes, in addition to genes expressed throughout all zones of the growth plate. This analysis suggested that there is differential transcriptional regulation of many genes in fetal cartilage and that the data could be used to identify genes selectively expressed in the tissue.

Tissue-selective genes have been previously defined as genes with enriched expression in a particular tissue [14] and characterized with algorithms dependent on the degree of differential expression relative to other tissues, including t-test [15], SAM [16], fold change [14, 17, 18], and enrichment scores [14]. While these approaches successfully identify tissue-selective genes, the reliance on fold change reduces the significance of many selectively expressed genes with low fold change. To compensate for this and identify cartilage-selective genes expressed at lower levels, the approach presented here placed increased significance on the preferentially expressed genes that showed the least variation of expression in non-cartilage tissues. This was made possible by the use of publicly-released reference gene expression data performed on the same platform and led to the reliable identification of genes with lower fold changes, but high cartilage selectivity. The impact of the use of coefficient of variation on the ranked gene list is apparent in Tables 1 and 2. In the U133A dataset, nine of the top 25 genes were ranked higher than 100 in significance in the SAM ranking. The average fold change of the probes for these nine genes was 10.7, while the average fold change of the probes for the other 14 of the top 25 genes was 42.2. One of these probes, COL10A1 was among the top four cartilage-selective genes using the CV algorithm but ranked at 284 by SAM (Table 1). In the U133B dataset, which contains a higher percentage of unannotated genes, 4 of the top 50 probes had a SAM ranking below 100, and the average SAM ranking was 576. Overall, to identify only the most cartilage-selective genes, a threshold of 50% coefficient of variation was used across all three platforms, yielding 161 cartilage-selective genes. A subset of 13 of the 161 cartilage-selective genes was studied by quantitative RT-PCR in cartilage and eight non-cartilage tissues to independently assess tissue selectivity. The data confirmed the cartilage-selectivity of genes with less than 50% CV, validating the selection procedure and suggesting that the gene expression patterns determined by microarray analysis are representative.

The coefficient of variation selection approach could, in theory, equally select for three different patterns of expression: cartilage-specific genes; genes with a consistent level of baseline expression in non-cartilage tissues; and genes with significant but equal expression in all tissues (e.g. housekeeping genes). In this data analysis, however, the most highly ranked genes consistently demonstrated little or no expression in non-cartilage tissues. The data thus demonstrate that incorporating coefficient of variation preferentially selected for genes not significantly expressed in non-cartilage tissues, yielding genes likely to have important and perhaps unique roles in cartilage.

Regardless of expression level, a cartilage-selective expression pattern suggests that the product of each identified gene may have a functional role in the development of the skeleton. Concordant with this hypothesis, mutations in 25 of the 161 selected genes have been associated with skeletal phenotypes in humans and/or mice. Included among them were the products of the well characterized genes encoding aggrecan and the cartilage-specific collagens, gene products known to have a prominent role in skeletal development and endochondral ossification. By this measure, the remaining genes may be candidate genes for skeletal dysplasias in which the disease gene has yet to be identified. As new skeletal dysplasia loci are defined, coincidence between a locus and a cartilage-selective gene may promote rapid identification of the disease gene. Knockout of the orthologous genes in mice would also facilitate exploring the role of each gene in skeletal development.

Classification of the biological roles of the products of the cartilage-selective gene set reveals genes with diverse functions including structural proteins of the cartilage extracellular matrix, enzymes that modify them, and 41 gene products with unannotated function. There were 65 genes that are components of signaling pathways, and only 43% of these were identified by sequence analysis of a comparable fetal cartilage cDNA library [4]. Among the genes were elements of the nitric oxide, VEGF, TNF/RANK, and gp130 pathways, all of which have known roles in the growth plate [19–23]. Mutations in the genes encoding some of the molecules in these pathways, including RPS6SKA3, LIFR, TNFRSF11A and IKBKG, have been associated with human skeletal dysplasias [12], again suggesting that the remaining genes may also serve critical roles in endochondral ossification.

Multiple genes encoding members of the LIF/gp130 signaling pathway met the definition of cartilage-selective genes. LIF is a cytokine that is expressed in terminally-differentiated growth plate chondrocytes [24] and signals through the gp130/LIFR complex. Homozygosity for loss of function mutations in the LIF receptor produces the recessively inherited skeletal dysplasia, Stuve-Wiedemann syndrome [25]. In addition to their skeletal features, these patients have cardiovascular, pulmonary, gastrointestinal, neurologic and metabolic abnormalities, likely attributable to the role that LIFR plays in embryonic or fetal development. Genes on the cartilage-selective gene list upstream of the receptor include RELA and RELB, NF-KB survival transcription factors that increase transcription of LIF [26], as well as the LIF gene itself. Through the LIFR/gp130 complex, LIF can regulate both the JAK/STAT and ERK MAP kinase pathways. Pathway components downstream of the receptor include ATF1, part of the ATF1/CREB transcription factor complex that participates in ERK MAP kinase signaling [27, 28]. The ATF1/CREB complex is also regulated by phosphorylation by the product of the RPS6SKA3 gene [29, 30], another gene in the MAPK/ERK pathway that is associated with a skeletal phenotype. The gene encoding RPS, a phosphorylation target of RPS6SKA3 [30], was also cartilage-selective, but the role of this protein in growth plate differentiation has yet to be determined. Finally, the gene encoding FOSL1, a FOS-like transcription factor activated by the ERK/MAPK pathway which binds cJUN to form a transcription complex [31, 32], was among the cartilage-selective genes identified. Thus comparative microarray analysis has identified multiple components of a regulatory pathway that can be explored to further evaluate their importance in growth plate differentiation and endochondral ossification.

While this study has provided a deep set of genes that exhibit a cartilage-selective expression pattern, there are some limitations to the analysis. First, the study focused on total cartilage RNA, including all types of growth plate chondrocytes, as well as articular cartilage. As a result, it cannot be determined if the selected genes are expressed in all types of chondrocytes or only a subset of cells. In this context, nine of the cartilage selective genes have been shown to be more highly expressed in hypertrophic cells relative to proliferating chondrocytes in the rat and/or mouse [7, 8]. Second, the cartilage samples were derived from a single anatomic site and a narrow window of fetal development, so it is unclear to what extent the observed gene expression pattern can be generalized. Third, neither all possible non-cartilage tissues nor each type of cell within each tissue were studied, so cartilage-selectivity could be affected if additional fetal and/or adult tissues that express the identified genes were found. This may be particularly important for other connective tissues such as bone, tendon and ligament which contain cells known to express some of the cartilage-selective genes identified here.

Not all genes selectively expressed in developing cartilage will necessarily be identified using this approach. For instance, the COL2A1 gene fell just below the rigorous 50% CV standard set to define cartilage-selectivity. The underlying reasons for this are complex. Probe performance as well as the known expression of COL2A1 in fetal liver and heart, are likely to have had an effect, as both factors could have contributed to the variation in expression in non-cartilage tissues. In addition, the approach presented here treated the three expression platforms, U133A, U133B and U133 2.0 equally from the viewpoint of the threshold for cartilage-selectivity. Because the comparative dataset of normal tissues was both broader and deeper for the 133A platform, additional genes from this platform, albeit with greater than 50% CV, could be considered to be tissue-selective (e.g. COL2A1). Thus, a platform independent threshold would likely yield additional genes of interest within the U133A dataset. Finally, tissue-specific genes were identified using only microarrays and a single generalized algorithm. Additional genes selectively expressed in cartilage could be identified by less stringent criteria or other methods.

Genome-scale comparative expression analysis using human fetal cartilage and a broad set of normal human tissues has identified 161 cartilage-selective genes, including 27 uncharacterized genes. The data identify novel gene products that may provide essential roles in normal skeletogenesis and suggest new candidates for the over 100 inherited skeletal disorders in which the disease gene has not been identified. The results demonstrate that fetal cartilage is a complex and transcriptionally active tissue, and that the set of genes selectively expressed in cartilage has been greatly underestimated.

A flow chart outlining methods and results as well as other supplemental information is provided in additional data file 1.