Transcription profiling of the developing humerus during early ossification and joint formation (TS23)

RNA sequencing of control humeri at TS23 gives an insight into the transcriptome at this key stage of development when the rudiment is undergoing early stages of ossification and tissue zones in the joint are being defined (Figure 1). Transcripts from 15,214 individual genes (n = 3 biological replicates, Figure 2; data submitted to EMBL-EBI ArrayExpress repository: [E-MTAB-1745]) were detected in this tissue at this stage. A minimum of 5 transcript reads (on average), with at least one read from each replicate sample, was chosen as the cut off point to reliably indicate reproducible expression across biological replicates[24]. Expressed genes were divided into groups according to their relative expression level (Figure 2A). 787 genes showed the lowest level of expression represented by between 5 and 10 read counts, while the majority of genes (71%; 10,789/15,214) showed between 100 and 5000 read counts. Only 732 genes are in the most highly expressed categories (≥5,000 read counts) (Figure 2A), representing 4.8% of expressed genes. Expression levels of selected individual genes are represented in Figure 2B. The most abundantly transcribed gene is Col2a1 (collagen 2, alpha 1) with a read count of 452,576, and among the 8 genes with read counts of more than 100,000, there are 4 other collagen encoding genes (Figure 2A and B). In total 41 collagen subtype genes are expressed (Figure 2B).

The relative levels of expression of genes associated with signalling pathways involved in regulating skeletal development reviewed in[25, 26] are highlighted in Figure 2B. This shows the potential components that can contribute to these signalling pathways at this stage of skeletal development. For example the hedgehog (HH) pathway is known to play an important role in ossification through the action of Ihh (640 reads), binding to its receptor Ptch1 (1,342 reads) activating Smo (5,154 reads). Similarly, 14 FGF ligand encoding genes and 5 FGF receptor encoding genes were detected, highlighting the potential for multiple FGF signalling interactions. The BMP signalling pathway genes also reveal potential for multiple signalling interactions with 9 BMP encoding genes expressed. Bmp1 is by far the most highly expressed (5719 reads); although it is not previously reported in this tissue at this stage in gene expression databases (http://www.informatics.jax.org/). The relative expression levels of Wnt ligand (Wnt), Wnt receptor (Fzd, Lrp, Kremen), extracellular Wnt interactor (Sfrp, Dkk, Rspo, Sost), intracellular Wnt pathway component (Dvl, Axin, Ctnnb) and Wnt pathway transcription factor (Tcf, Lef) encoding genes are represented. The detection of previously unreported Wnt gene expression in the humerus (e.g. Wnt2 and Wnt2b) opens up new considerations for functional roles, especially as both genes are up-regulated in muscle-less rudiments (described below). The low density lipoprotein receptor-related genes Lrp5 and Lrp6, which are Wnt co-receptors, are most highly expressed among the Lrp gene family (4,727 and 3,310 reads respectively). Interestingly eight Fzd receptor encoding genes are detected. Other genes known to be involved in skeletal development are highly expressed; Sox9 (4,550 reads), Runx2 (1,285reads), Spp1 (8,155 reads) and Mmp13 (1,103 reads).

Identification of differentially expressed genes in muscle-less versus control developing humeri and associated joints

The Venn diagrams in Figure 3 represent overlap of DE genes identified by the two platforms. In total 1,132 independent genes were identified as DE in the muscle-less mutant compared to its phenotypic control (Figure 3, Table 1). Of the 374 genes identified by microarray (Table 1), 73.2% of these genes were also represented in the RNA-seq data. RNA-seq detected a greater total number of DE genes than the microarray (1,037 versus 374, respectively, using the standard cut off criteria). In addition if the stringency for DE gene selection is weakened for the microarray data by moving the cut off point to corrected p-value of ≤0.08, the number of genes detected as DE in common by the two platforms increases to 426 from 274 (not shown). This suggests that the RNA-seq approach was more sensitive in detecting differential expression. Grouping genes based on the degree of change showed that the spread of fold change values was similar across platforms although there were differences in the number of genes in different categories (Additional file1: Table S1). Additional file1: Table S2 and S3 list all DE genes and fold change values estimated by Microarray and RNA-seq.

Biological Interpretation of differentially expressed genes: Down-regulated genes are associated with development and differentiation, cytoskeletal architecture and cell signalling

110 genes associated with the Cytoskeleton were down-regulated in Pax3 ^Spd/Spd humeri, including those encoding microfilament, microtubule and intermediate filament components (Figure  4); 33 are directly associated with microfilaments (e.g. Acta1, Ablim2, Pdlim3), 13 with microtubules (e.g. Rassf5, Tubb2b, Tppp) and 4 with intermediate filaments (Des, Krt33b, Sync, Tchh); other DE genes associated with the cytoskeleton encode proteins that interact with myosin, or the extracellular matrix (ECM), including integrin and cadherin encoding genes (e.g. Lamb3, Cdh4, Itga4 and Pdgfa (not shown)).

Up-regulated genes are associated with cell adhesion, cell signalling and development and differentiation

Genes up-regulated in muscle-less humeri revealed by microarray and RNA-seq were similarly analysed for enrichment of genes associated with particular biological processes or cellular components, using associated GO terms (Table 4). For terms within the subontology biological process, the strongest enrichment was for cell adhesion and ECM associated genes (total of 50/447 genes, with p-values down to 1.99 × 10^-10). The genes identified in this category include ECM glycoproteins (Lsamp, Svep1), ECM structural constituents (Col8a1, Col8a2, Frem3), cell-adhesion molecules (Cntn1, Cntn4, Cntn5) and calcium-dependent cell adhesion proteins (Pcdh8, Pcdhb2, Pcdhga1). This grouping in addition included genes involved in signalling pathways which overlap with the next most enriched terms; cell signalling (GO:0007154) and cell-cell communication (GO:0007267). The signalling pathway components identified in this category include: Hedgehog (Hhip, Ptch1), fibroblast growth factor (Fgf10, Fgf2), transforming growth factor (Frem2, Bmp3), Notch (Ctn1, Nrg1) and Wnt signalling (Dkk2, Cpz Rspo2, Rspo3, Sfrp2, Wnt2, Wnt2b, Wnt4 and Wnt16) and others, including receptors (Epha3, Epha4, Epha5, Grin2a, Grm4, Grm7, Grm8, Gfra2, and Pdgfra). Other signals identified as up-regulated included c-fos induced gowth factor (Figf), hepatocyte growth factor (Hgf) and Insulin-like growth factor (Igf1). The gene lists in the next most enriched set, Development and Differentiation, similarly show large overlap due again to the presence of the signalling pathway genes mentioned above, and also transcription factors Foxc2, Foxo3, Lmx1a, Lmx1b. z.

Under the subontology Cellular Component there was also striking enrichment of extracellular (total 81/447; with p-values down to 2.34 x10^-10) and membrane associated gene products (total 227/447 p-values down to 2.46 × 10^-16), including cell adhesion molecules (Cadm2, Cntn1), receptors (Fat2, Fat4, Grin2a, Grm4, Gfra2, Pdgfra, Robo2, Sorcs1, Sstr4), cell surface molecules (Cd83, Cd96), cadherins (Cdh20, Cdh8), trans-membrane proteins (Tmem26, Tmem28), voltage gated channels (Kcna1, Kcna2, Kcnc2) and cell adhesion and extracellular components: Alcam, Cntn4, Epha4, Col8a1, Col8a2, Pappa, Pcdh8.

Signalling pathway analysis of differentially expressed genes

Given the strong enrichment of genes associated with Signalling Pathways and Development and Differentiation functions (Tables 3 and4), we sorted DE genes according to participation in major developmental regulatory pathways (Figure 5 and Tables 5 and6). By far the most strongly impacted cell communication pathway is Wnt signalling with 34 DE genes encoding signalling molecules, receptors, pathway antagonists, known targets or potential targets of the pathway (Figure 5, Table 5). From the diagrammatic representation of pathway components shown in Figure 5, it is apparent that the encoded products of DE genes act either at the cell surface in Wnt signal generation/modulation/interpretation or are targets of the pathway. The genes listed include known targets of the pathway and seven potential target genes (not verified), included here due to their similarity to known targets; for example Sall1 is included because the orthologous gene Sall4 is a known direct target of the pathway[36]. In general, genes encoding cell surface components of the pathway are up-regulated, including signalling ligands and agonists (Wnt2, Wnt2b, Wnt4, Wnt16, Rspo2 and Rspo3), and extracellular antagonists (Dkk2, Sfrp2), while down-regulated genes identified are more commonly targets of the pathway (e.g. Fgf4, Cacng1, Pitx2, Dll1, Prg4, Lrrn1, Met, and Cd44). Interestingly nine known Wnt target genes are up-regulated, including Dkk2, Rspo2, Rspo3, Cldn1, Grem2, Kcnd1, Epha4 and Sfrp2, which encode membrane associated proteins, some of which regulate the Wnt pathway (Dkk2, Sfrp2, Rspo2, Rspo3).

Table 6 lists DE genes associated with other developmentally relevant signalling pathways, including the BMP, Hedgehog, Fibroblast growth factor, Hippo and Notch signalling pathways.

Spatial alteration of gene expression patterns

While Microarray and RNA-seq analysis provides data on quantitative changes in gene expression levels across the whole developing rudiment, it does not reveal alterations in the spatial distribution of transcripts or give clues to the specific developmental events affected. We therefore performed in situ hybridisation on control and muscle-less mutant (Pax3 ^Spd/Spd ) limb sections at TS23 for a selected subset of genes. Three of the genes selected encode components of Wnt signalling pathways known to be important in skeletal development: Cd44, Sfrp2 and Wnt4. Spp1 encodes an ECM protein Osteopontin which is a prominent component of mineralised matrices of bone and teeth[76]. Cd44 is a cell-surface glycoprotein involved in cell-cell interactions, cell adhesion and migration. It is a receptor for hyaluronic acid and can interact with other extracellular proteins, such as osteopontin, collagens and matrix metalloprteinases (MMPs) reviewed in[77]. It is a target gene of the Wnt signalling pathway[50]. Cd44 gene expression is down-regulated 2.28 fold in muscle-less humeri RNA (Table 5). The in situ hybridisation analysis reflects this down regulation dramatically with Cd44 transcripts hardly detectable in either shoulder or elbow joints of muscle-less mutant embryos, compared to the clear joint line restricted expression seen in controls (Figure 6A-D). Sfrp2 encodes a secreted protein that acts as a modulator of the Wnt signalling pathway, in particular during normal skeletal patterning in developing limbs[78]. Normal expression at TS23 (Figure 6E, G) can be detected in the elbow and shoulder joints. The level of up-regulation from microarray and RNA-seq analysis is 2.62 (p-value 0.01) and 2.09 (p-value 1.13E-06) fold respectively (Table 5). This up-regulation was reflected in the intensity of expression seen in mutant sections following in situ hybridisation (Figure 6F, H). In addition to the increased level of expression a change in the spatial pattern is also evident. In both the elbow and the shoulder joints expression is expanded, particularly on the ventral aspect of the joint and the staining is unevenly distributed (Figure 6F, H). A similar expression increase was seen in the phalangeal and carpal joints of the handplate, compared to control (not shown). Wnt4 is a signalling ligand of the Wnt signalling pathway. Expression of this gene has previously indicated its role in joint formation during limb development[79–81]. Control Wnt4 expression is seen at a low level in the ventral portion of the elbow joint (Figure 6I) at TS23; there is also distinct expression in the epidermis, as previously detected[82]. In the mutant there is a high level of expression in the elbow joint region; expression is spread across the whole joint line from ventral to dorsal, more extensive on the dorsal side (Figure 6J).

The Spp1 gene is normally expressed in the hypertrophic zone and adjacent perichondrium[83]. Despite the appearance of hypertrophic chondrocytes at the mid-diaphysis of immobile Spd embryos (Figure 6N arrowhead), no Spp1 gene expression is detected in these cells. There is weak but detectable staining in the perichondrium at the site of the hypertrophic region, but again apparently lower than in control tissue (Figure 6M, N). Expression is also detected in the hypertrophic zone of the ulna (Figure 6O) and this is reduced but not absent in the muscle-less mutant ulna (Figure 6P), reflecting the reduced phenotypic effect seen in this rudiment[9]. Expression is seen only in the perichondrium of the radius in both the control and mutant (Figure 6O, P).

The transcriptome of a developing skeletal rudiment

Utilising RNA-sequencing technology to reveal the transcriptome in the normal developing humerus and associated joints at TS23 provides an insight to the processes that are occurring during this stage of skeletal development when chondrocytes are undergoing hypertrophy, the cartilaginous rudiment is beginning to ossify at the mid diaphysis reviewed in[84] and specific zones within the joint are differentiating reviewed in[85]. This adds a valuable resource to a growing set of data that can be combined to explore skeletal development. Previous transcriptome profiling studies have examined cartilage condensation in the tibia-fibula from E11.5 to E13.5[86], ossification of metatarsals[87], in different zones of the growth plate in 2 week old mice[88] and in Runx2^-/- mutant mice[89]. Earlier stages of limb development, although not specific to skeletal development, are informed by whole limb bud profiles between E9.5 and E11.5[32]. Another important resource that has recently become available examines the spatial expression patterns of 18,000 coding genes and over 400 microRNAs within the whole mouse embryo at the same developmental stage examined here[90]. This digital transcriptome atlas is a powerful tool (Eurexpress; http://www.eurexpress.org/ee/intro.html) that can be used to examine spatial analysis of specific genes, exploring possible functional associations. Combining these resources gives information on quantitative and spatial expression of individual genes providing the basis to explore regulatory networks active during the development of skeletal rudiments.

Several of the findings of the transcriptome analysis are as expected; the most highly expressed genes (8 with ≥100,000 read counts) include 5 collagen encoding genes (Figure 2B); collagens have been shown to be the most abundant structural proteins in cartilage and show characteristic distribution patterns as skeletal rudiments develop[31, 91]. The Insulin-like growth factor genes Igf2 and Igf1 and their associated receptors Ifg2r and Igf1r are also highly expressed (4,194 to 221,194 read counts); these are reported to play a prominent regulatory role in skeletal development (reviewed in[92]). Similarly, aggrecan (Acan) and osteopontin (Spp1), both involved in skeletal development, are highly expressed (Figure 2B). Although much is known about the regulatory network that controls early chondrogenesis and joint formation (Kronenberg[25]; Provot and Schipani[26]), open-ended whole transcriptome studies are required to add new information. Centred on regulatory signalling pathways known to be involved in skeletal development; HH, FGF, TGFβ (including BMP) and Wnt; a complete list of components of these pathways expressed in the humerus as TS23 was extracted (Figure 2B) indicating the potential role players in each of these pathways at this specific stage. The importance of Indian hedgehog (Ihh) expression in the early prehypertrophic and hypertrophic chondrocytes of cartilage condensations is well established[89, 93, 94]. The full spectrum of possible interacting molecules in HH signalling (Figure 2B) include the receptors Ptch1, Ptch2, Smo and transcription factors Gli1, Gli2, and Gli3. In addition to Ihh, Desert hedgehog (Dhh) expression was also detected (read count 70) and, examining the data presented by Cameron et al. confirms that Dhh is up-regulated (≥3 fold) in E13.5 fibual and tibual cartilage[86]. Dhh has not previously been functionally linked to skeletogenesis but this opens the possibility of regulatory contribution, perhaps co-operating with Ihh. No expression of Shh was detected. Similarly, novel components of the FGF, TGBβ and Wnt pathways were identified (Figure 2B).

The Wnt signalling pathway plays a central role during embryonic development and is known to be an important regulator of bone formation and bone remodelling reviewed in[95]. It also plays a pivotal role in joint formation and maintenance, shown through gain and loss of function experiments[80, 96]. The key intracellular mediator of canonical Wnt signalling, β-catenin (Ctnnb1; 28,608 read counts) is the most highly expressed Wnt signalling component in the TS23 humerus and associated joints. The most highly expressed Wnt ligand is Wnt5a, previously associated with expression in joints and perichondrium[79] and proliferating chondrocytes[97]. Other highly expressed ligands include Wnt9a, Wnt5b, Wnt11, and Wnt4. Wnt5b and Wnt11 expression has been shown in the pre-hypertrophic chondrocytes and Wnt4, Wnt9 and Wnt16 in the developing joints[79–81, 98]. High expression of Wnt9a could be due to its role in the temporal and spatial regulation of Ihh[96].

Numerous extracellular modulators of the pathway were detected; all five secreted frizzled related protein (Sfrp) genes, Dickoff (Dkk) 1, 2 and 3 genes and four R-spondins (Rspos), indicating a huge potential for pathway modulation. The most highly expressed antagonist modulators of the pathway were Dkk3 and Sfrp2 both of with are detected in joint cells at E13.5 and E15.5[82]. The most highly expressed R-spondin agonist of the pathway was Rspo3, previously detected in phalanges[55].

Identification of Mechanoresponsive genes

The identification of differentially expressed genes between humeri from control and muscle-less embryonic limbs allows an investigation of the biological processes and the developmental regulatory signalling pathways that are affected by the removal of mechanical stimulation on skeletogenesis in vivo. We previously reported that muscleless (Splotch) mutants display abnormal ossification in the humerus, altered humeral morphogenesis and altered elbow and shoulder joint formation and these effects were first observed at TS23[9]. This was chosen as the point of analysis for differential gene expression because, although limb muscles begin to contract from approximately E12.5[99] it is uncertain how much stimulation is transmitted to the skeletal rudiments when the developing tendons are at early stages of morphogenesis[100]; the sole indication that the force is functionally transmitted is the mutant phenotype seen at TS23. Although analysis at TS23 may miss some of the earliest effects, it is relatively early in the response and the earliest time at which it is certain that the system is disturbed. Alteration in expression pattern of some selected candidate genes and pathways was previously revealed[8, 10], but here we carry out the first genome wide study identifying a total of 1,132 independent genes as differentially expressed: with approximately 60% down-regulated and 40% up-regulated. The finding of more genes being down-regulated than up-regulated and to a greater extent is consistent with the proposal that mechanical stimuli support the correct differentiation of cells, as observed in the ossification phenotype[9], and for the maintenance of tissue patterning, as seen in the developing joint[16]. GO annotation analysis identified specific biological processes that are affected when mechanical stimuli are removed. This type of analysis has been used previously to interpret biological processes associated with developing skeletal tissue[87, 89, 101]. Analysis of the down-regulated DE gene set identified genes associated with development and differentiation as the most highly enriched categories, including developmental regulatory signalling pathway molecules and transcription factors. Similarly, analysis of up-regulated DE gene sets indicated genes associated with cell signalling and development and differentiation.

DE genes were also highly enriched for genes associated with the cytoskeleton. The cytoskeleton controls cell shape, organelle transport, cell motility and division, and connects the extracellular matrix to internal cell processes reviewed in[102]. It maintains the mechanical integrity of cells and has been implicated in relaying mechanical signals to downstream biochemical responses[7, 103]. This was seen in the embryonic lung where cytoskeletal network inhibitors resulted in altered tissue morphogenesis and conversely when cytoskeletal tension was activated lung development was accelerated reviewed in[104], indicating the dynamic role the cytoskeleton has in morphogenesis.

In chondrocytes the actin microfilaments are predominantly located at the periphery of the cytoplasm[105], tubulin microtubules are uniformly distributed throughout the cytoplasm[106] as are intermediate filaments, connecting the nuclear membrane with the cell periphery[107]. In this study 84 genes annotated as cytoskeletal were down-regulated when mechanical stimulation was removed. These include 33 genes directly associated with actin microfilaments, 13 with microtubules and 4 with intermediate filaments (Figure 4). The most highly affected group, the Filamentous-actin cytoskeleton, has been shown to be involved in articular cartilage chondrocyte mechanotransduction, converting a mechanical stimulus into a biochemical response[103, 108–110]. Other studies have confirmed the involvement of the actin cytoskeleton in cartilage chondrocyte mechanotransduction via manipulation of the actin accessory proteins[108, 111], but there are few reports on the affect of mechanical stimulation on microtubule and intermediate filaments[112]. Among the DE genes is an actin-binding protein, cofilin2 (Cfl2); cofilin was previously shown to be increased following cyclic mechanical loading of chondrocytes[111].

The identification of cytoskeletal genes down-regulated following the removal of mechanical stimulation indicates that the cytoskeleton is affected, but is this because the mechanical integrity of the cell is altered or because mechanotransduction from the ECM is affected, or perhaps a combination of both? The finding that ECM and cell adhesion associated genes are also affected further supports changes in mechanotransduction pathways. The cell-adhesion associated integrins (Itga4, Itga7, Itgb6) and cadherins (Cdh4, Cdh15, Cdhr3) are down-regulated and these proteins potentially function to physically couple cells to the ECM and play a role in mechanical signal transduction[103, 104]. Articular chondrocytes have been shown to express both integrin[113] and non-integrin[114] ECM receptors. Another actin-associated protein identified to be down-regulated is actinin-α2 (Actn2); this protein also couples the cytoskeleton to the ECM and may be involved in transducing mechanical stimulation.

Secreted phosphoprotein 1 (Spp1), previously known as Osteopontin (Opn), is one of the abundant non-collagenous proteins in bone matrix produced by osteoblasts and osteoclasts reviewed in[115]. Spp1 binds to hydroxyapatite and is a potent inhibitor of the mineralisation process, inhibiting the growth of bone matrix crystals[116]. Spp1 is expressed early in bone development, however it was concluded not to be required for normal development of bones as null mice (OPN^-/-) have no apparent effect on the structure or distribution of cells within bone tissue[117]. However, Spp1 expression has been shown to be regulated by mechanical stimulation both in vitro and in vivo[35, 118–120]. We found Spp1 to be down-regulated in the developing humerus at TS23 in muscle-less embryos and in situ hybridisation showed a dramatic absence of detectable Spp1 expression in hypertrophic chondrocytes whereas it is still detectable in the perichondrium (Figure 6N), indicating a specific effect on expression in hypertrophic chondrocytes and not a delay in the onset of normal expression. It was previously shown that OPN^-/- mice did not suffer bone loss in response to mechanical unloading[120], suggesting that mice lacking Spp1 could not sense the changes in mechanical stress, thus indicating its potential role in the signal transduction of mechanical stimulation. It has been suggested that mechanotransduction through Spp1 is dependent on microfilament integrity, as mechanically stimulated increases in Spp1 expression was blocked by disruption with cytochalasin-D in osteoblasts[119]. This again highlights the link between an ECM component and the cytoskeleton in a mechanoresponse implicating these components in signal transduction, either directly through the cytoskeleton or through cell adhesion complexes via the cytoskeleton.

An example of a non-integrin ECM component that is down-regulated in the absence of mechanical stimulation is Cd44, a target gene of the Wnt signalling pathway[50], encoding a single-pass membrane glycoprotein that binds proteoglycan and hyaluronan (HA) to produce a pericellular matrix surrounding chondrocytes reviewed in[121]. Cd44 has been implicated in joint cavitation through interaction with HA[122]. We previously showed loss of expression of the Cd44 gene in the interzone of the forming knee joints of immobilised chick embryos[16], one of a number of gene expression patterns reflecting a general loss of organisation of differentiating tissue territories, and here we show a similar effect on Cd44 expression in the elbow and shoulder joint of muscle-less mouse limbs, where we previously showed a similar loss of tissue organisation[9]. The very restricted expression of Cd44 in the interzone of forming joints in control embryos at TS23 is barely detectable in muscle-less mutants (Figure 6A-D). As well as the gene being sensitive to mechanical stimulation, as an integral part of the ECM and a regulator of joint formation, the gene product may also be an important mediator of mechanical stimuli.

The link between Cd44 and the Wnt signalling pathway highlights perhaps the most striking finding of this analysis; the altered expression of 34 genes implicated in the Wnt signalling pathway (Figure 5). Canonical Wnt signalling has been shown to be involved in maintaining joint integrity[80, 96] and is disturbed in the joints of muscle-less mouse embryos[8]. Wnt signalling might also be involved in co-ordinating ossification and joint development; both processes altered in muscle-less embryos[5]. Non canonical signalling has also been implicated in planar cell polarity during growth plate regulation[123, 124]. There are also previous indications that the Wnt pathway is responsive to mechanical stimulation in mesenchymal stem cells[125] in mature bone in vivo[126] and in response to injury of articular cartilage[127]. Here, genes encoding four Wnt ligands are up-regulated in muscle-less embryos and in the case of Wnt4 we show specific up-regulation in the elbow and shoulder joint region (Figure 6I-L). Two of the up-regulated genes (Wnt2 and Wnt2b) have not previously been associated with skeletal development. The Wnt signalling antagonist Sfrp2 is also up-regulated specifically in the joint region (Figure 6E-H). The majority of known Wnt target genes affected are down-regulated (Figure 5), perhaps due to increased expression of negative regulators Sfrp2 and Dkk2 and down-regulation of the Fzd10 receptor; however ligands and agonists Rspo2 and Rspo3 are up-regulated as are some target genes indicating effects at multiple levels of regulation of the pathway. It is interesting that a number of the up-regulated targets feedback as negative regulators of the pathway (for example: Dkk2, Sfrp2). It is now important to functionally test the mechanisms linking mechanical stimulation with Wnt signalling. This work provides sets of candidate genes to use in functional assays to excavate this important link. Understanding how mechanical stimuli influence the Wnt signalling pathway would be a major step forward in understanding how mechanical cues work together with classical molecular positional information to guide spatially appropriate tissue differentiation and provide indications of how conditions can be effectively recreated in vitro to guide stem cell differentiation.

In situ hybridisation analysis showed altered gene expression of 3 Wnt pathway genes in the developing shoulder and elbow joints; one down regulated (Cd44) and two upregulated (Sfrp2 and Wnt4) (Figure 6). Changes in the spatial restriction of Sfrp2 and Wnt4 expression were also seen. This does not represent a delay in normal expression in the mutant because the altered pattern is not reminiscent of earlier stages[82] and the changes are consistent with altered patterning of the tissue territories in the forming joint and the fusion phenotype seen in the mutant, with cartilage forming across the joint at later stages[8, 9].

In this study differential expression in developing skeletal rudiments is documented in the absence of limb muscle; this will include genes that respond to lack of mechanical stimulation but perhaps also as a paracrine response to adjacent muscle cells. We know that phenotypic effects on ossification and joint formation are due to the lack of mechanical stimulation rather than physical absence of muscle cells because phenotypic analysis of a range of mouse mutants where muscle is immobile[8], reduced[9] or absent[8, 128] have similar effects and we see similar effects in immobilised chick embryos[15, 16]. Therefore, although some of the genes identified here may respond to lack of adjacent muscle tissue, many must underlie the phenotypic effects seen in response to lack of mechanical stimulation. This is further supported by the overlap of some of the genes identified here and in skeletal cells in culture or adult tissues, in response to mechanical stimulation[18, 19].

Another important limitation in this work is the possibility that a proportion of the down-regulated genes may be due to contamination of the dissected control humeri with adjacent mesenchyme/muscle, since this is being compared to tissue from muscle-less embryos. Although care was taken with the dissections, it is impossible to be sure that all muscle tissue was eliminated from the control. The down-regulated gene set also showed enrichment for muscle associated genes (Additional file1: Table S4), consistent with possible contamination of the dissected control humeri by neighbouring muscle. To inform this we also sequenced the transcriptome of mesenchyme adjacent to the humerus of control embryos at TS23 and compared it to the transcriptome of control humeri. We then cross referenced this to the down-regulated gene set in control versus muscle-less humeri (Additional file1: Table S2), noting any genes enriched more than 3 fold in mesenchyme compared to control humeri; these are indicated in column 2 of Additional file1: Table S2. It is possible that these genes are involved in both cartilage and muscle development so no genes have been removed from the data set, however, DE genes also showing higher expression in mesenchyme compared to control humeri must be treated with caution with respect to a skeletal specific response to mechanical stimulation. Such genes have not been prioritised in any of our subsequent exploration of candidate mechanosensitive genes.

The developing humerus at TS23 constitutes different cell and tissue populations at different stages of differentiation including the joint region, the perichondrium and the organised zones within the cartilage rudiment. Therefore the experimental design employed here will capture genes associated with different cells types at different stages of differentiation. It will now be important to sort out which cells and tissues have altered expression of specific genes. This can be addressed for a sub set of genes by in situ hybridisation, with an initial analysis of 4 genes presented in Figure 6. It can be addressed in a high throughput manner by isolating specific cell populations using laser microdissection from tissue sections (laser capture), purification of RNA and quantitative RT-PCR gene expression profiling, comparing control and mutant tissue from,, for example the hypertrophic, prehypertrophic or the elbow joint region alone.

We used both RNA-sequencing (Illumina) and Microarray (Agilent) technologies in parallel to determine differential expression. Microarray technology has been utilised to determine expression of chondrogenic and osteogenic genes from developing whole tissues[32, 87, 89], and from in vitro differentiation procedures[19, 34, 129–131]. The use of RNA-seq technology to describe the transcriptome is more recent[132–134]. Previous direct comparisons between microarray and RNA sequencing-based approaches to reveal alterations in gene expression between tissues reported that RNA-seq identified more DE genes[23, 24, 135]. We also found that RNA-seq is more sensitive in reproducibly detecting alterations in gene expression, detecting more genes altered at lower quantitative levels (Additional file1: Table S1; 5 > 2 fold). This was further emphasised by reducing the stringency of the statistical analysis to p ≤ 0.08, which increased the number of genes detected by microarray specifically (not shown). An example of the importance of the increased sensitivity and reproducibility of RNA-seq is shown by the Spp1 gene which did not show statistical significance by microarray but has been verified by qRT-PCR and in situ hybridisation (Table 2 and Figure 6). The larger dynamic range[24] and higher reproducibility across replicates[135] has also been found in other studies.

RNA preparation

Heterozygous Splotch-delayed (Pax3 ^Spd/+ )[12] mice were purchased from Jackson Laboratories (Jax®). All animal work was carried out under the guidelines of Trinity College Dublin Bioresources Unit and Bioethics Committee. The generation of homozygous Pax3 ^Spd/Spd mutant embryos was achieved by crossing heterozygous Pax3 ^Spd/+ males and females. Embryonic material was collected from timed pregnancies on the afternoon of the 14th day (E14.5). Individual embryos were dissected and the developmental stage according to Theiler criteria[21], and the phenotype were recorded. All embryos were genotyped following PCR amplification as described in[136]. The humeri, including the associated joint regions, were finely dissected from control and mutant embryos at stage TS23 (Figure 1A). Tissue was mechanically homogenised and total RNA extracted (SV Total RNA Isolation System: Promega, UK). Pooling of rudiment tissue from multiple embryos of the same genotype (2–4) was performed. RNA integrity was assessed on a 2100 Bioanalyser (Agilent Technologies); RNA samples with RIN (RNA integrity number) values of 8.2-9.6 were used for Microarray and RNA-seq analysis.

Microarray

Four independent pooled sets of samples were used for microarray (n = 4 biological replicates) analysis. All microarrays were processed at IMGM® Laboratories (Martinsried, Germany). 100 ng of total RNA per sample was reverse transcribed into cDNA and then converted into labelled cRNA by in vitro transcription incorporating cyanine-3-CTP (Low input quick-amp labelling kit one-colour, Agilent Technologies). Genome wide expression profiling was carried out using the Agilent Mouse GE v2 Microarrays (4x44K format) (G4846A, Agilent Technologies) which contains 39,485 coding and non-coding sequences of the mouse genome (Figure 1B). A one-colour based hybridisation protocol was performed at 65°C for 17 hours on separate mouse GE v2 microarray platforms. Microarrays were then washed with increased stringency using Gene Expression Wash Buffers (Agilent Technologies) followed by drying with Acetonitrile (Sigma). Fluorescent signal intensities were detected with Scan Control A.8.4.1 software (Agilent Technologies) in the Agilent DNA microarray scanner and extracted from the images using Feature Extraction 10.7.3.1 software (Agilent Technologies). The software tools Feature Extraction 10.7.3.1, GeneSpring GX 11.5.1 and Spotfire Decision Site 9.1.2 (TIBCO) were used for quality control and statistical data analysis. Quantile normalisation was applied to each data set in order to impose the same distribution of probe signal intensities for each array[137], thus adjusting them to a uniform level that can allow for comparable downstream analysis. Welch’s approximate t-test (“unpaired unequal variance”, parametric) was applied to compare the control and mutant groups. A corrected p-value was calculated based on the algorithm of Benjamini and Hochberg[27], based on control of the False Discovery Rate (FDR). A fold change of ≥ 2 and FDR-adjusted p-value of ≤0.05 were used as criteria to indicate differential expression between the two groups.

RNA-sequencing: alignment and differential expression analysis

Three independent pooled sets of samples were used for RNA-seq (n = 3 biological replicates) analysis. The DNase-treated RNA (3 μg) was used to prepare RNA-Seq libraries with the TruSeq RNA Sample Prep kit. A total of six cDNA libraries were constructed, representing triplicate biological replicates for each group. 40 bp single end reads were obtained from an Illumina GAII in FASTQ format, one sample per sequencing lane. The Tophat aligner (http://tophat.cbcb.umd.edu/) was used to align the reads to the mouse reference genome (mm9). After alignment the read counts for each gene were extracted using htseq-count (http://www-huber.embl.de/users/anders/HTSeq/) based on an mm9 Refseq gff file. Differential expression in our two groups was evaluated using DESeq version 1.4.1, implemented in R 2.14.1. DESeq uses a negative binomial distribution to model genic read counts following normalisation based on size factors and variance. As for the microarray analysis, p-values were adjusted by the procedure of Benjamini and Hochberg to control the type I error rate, and a cut off of p ≤ 0.05, and a fold change of ≥ 2 were used as a threshold to define differential expression.

Quantitative real-time reverse transcription-polymerase chain reaction

Gene ontology annotation analysis

Gene Ontology (GO) terms were utilised to reveal significant enrichment of groups of genes among the DE datasets from the microarray and the RNA-seq analysis using the Database for Annotation, Visualisation and Integrated Discovery, DAVID (http://david.abcc.ncifcrf.gov/)[28], and GOstat (http://gostat.wehi.edu.au/) software. Analysis of GO terms associated with biological process, molecular function and cellular component was performed on all data-sets independently and combined to identify significantly (p ≤ 0.05) enriched gene sets. The strength of the enrichment of any GO term-associated gene set is reflected in the calculated p-values, comparing the proportion of genes in the data-set and the proportion of genes in the genome bearing that annotation.

In situ Hybridisation

Availability of supporting data

The data sets supporting the results of this article are available in the EMBL-EBI ArrayExpress repository (http://www.ebi.ac.uk/arrayexpress/). The differential expressed data set from the Microarray [E-MTAB-1744], the differential expressed data set from the RNA-sequencing [E-MTAB-1746] and for the transcriptome [E-MTAB-1745]. Lists of differentially expressed genes are available in Additional file1.