Gene identification and analysis of transcripts differentially regulated in fracture healing by EST sequencing in the domestic sheep

ESTs and assembly into contigs and creation of an Ovis aries sequence database

In addition to the 47209 sheep ESTs generated from three libraries in the current project (Table 1), we extracted 10497 ESTs and 2391 mRNA entries from the EMBL database (release 82, Sept 2005). 344 sequences were rejected because less than 40 bp remained unmasked following quality control procedures. The lengths of the ESTs ranged from 40 to 928 with an average length of 555 nucleotides (Fig. 1a). All bases beyond 750 bp were masked. Using our sequences as well as the previously available sequences, 19087 contigs were assembled with an average length of 603 nucleotides (Fig. 1b). 397 contigs were shorter than 100 nucleotides and were excluded from further analysis. On average, each contig was assembled from 3.1 ESTs from all libraries, with a range of 0 to 2392 ESTs per contig (some contigs consisted only of previously available mRNA data and contained no ESTs). If the analysis is restricted to ESTs from the libraries described in this work, then each contig contained an average of 2.5 ESTs (Range: 0 – 2387). The contigs comprised of the largest number of ESTs were identified as collagen Iα 1 (2387 ESTs), collagen Iα 2 (1094 ESTs), and osteonectin (1086 ESTs). A total of 8965 contigs contained only a single EST clone.

Library	ESTs	Number of contigs
PF7	14884	5859
PF10	17148	6217
Pooled	15177	5421

We used the GeneNest package [12, 13] to analyze all available mRNA and EST sequences and to present the sequences in a user-friendly internet site. The main features of the GeneNest Sheep database are summarized in Fig. 2.

BLAST analysis of Ovis aries contigs to assign putative gene identities based on homology

BLAST searches were used to attempt to assign putative identities to the contig sequences identified in this study. Analysis against the NCBI RefSeq protein sequence database was performed, and 13149 contigs demonstrated a best e-value of 0.01 or lower (Fig. 3a). We used an arbitrary e-value cutoff of 10^-25 to identify significant homology to a known gene and allow a probable assignment of a sheep contig. With this cutoff, a total of 10387 contigs demonstrated significant homology to at least one other protein in the RefSeq protein database and could be assigned a putative identity, consisting of a total of 8789 distinct proteins. In cases where hits to a proteins from several species were identified, the most significant match was used to assign a putative identity. 292 contigs blasted to an Ovis aries RefSeq. A total of 2863 contigs (15%) had e-values of less than 10^-100 and are thus likely to represent true orthologs. Of the 19087 contigs, 11591 were assembled from ESTs only from the current project. 6577 were assigned to 5839 distinct RefSeq sequences by virtue of blastx e-values of 10^-25 or less. These sequences represent previously unidentified sheep genes. The sequences for these genes are available in FASTA format with the online supplemental material as [Additional file 1]

Gene-Ontology analysis of the sheep fracture libraries

GO analysis was performed by assigning annotations to the sheep contigs from the MGI annotations of putative murine homologs (e-value threshold 10^-25) in order to get an overview of the subset of genes expressed in healing fracture tissue of the sheep. Fig. 3b provides an overview of the results of this BLAST analysis. Putative homologs were found for 5637 distinct sheep genes, some of which were defined from multiple contig sequences. MGI annotations were available for 4600 of the 5637 genes. The results of the GO analysis using these annotations are summarized in Fig. 4. Many of the most frequently annotated terms in the subontology biological process, such as development, signal transduction, transcription, biosynthesis, and phosphate metabolism, reflect processes known to be involved in fracture healing and bone development. The distribution of annotations to terms in the subontology cellular component shows that the gene products corresponding to the ESTs identified in this work are active in many different cellular and extracellular locations.

A comparison of the sheep fracture libraries with normal human bone libraries

During development and in related processes such as fracture healing, genes exhibit complex patterns of spatial and temporal regulation. The identification of genes that are differentially expressed during fracture healing therefore represents an important initial investigation to discover genes that are essential for fracture healing.

In this work, we generated two EST libraries from early fracture healing stages (PF7 and PF10) with the goal of identifying genes important in the earliest stages of fracture healing. Following fracture, initial hemorrhage into the fractured bone is followed by hemostasis and a phase of inflammation, which lasts up to two weeks in the sheep (Fig 5). Intramembranous ossification begins around day 10. Therefore, we reasoned that genes that are overexpressed in the PF7 and PF10 libraries could represent genes with potential importance for the initiation of new bone formation. To estimate gene expression levels in the two libraries, we compared EST sequence counts using a two-tailed Fisher exact test. This analysis revealed statistically significant differential expression for only three genes; tartrate resistant acid phosphatase 5 (ACP5) increased by a fold change of 10.60 (5 to 53 ESTs), cytochrome c oxidase subunit I (COX1) increased by a fold change of 2.1 (87 to 178 ESTs), and perlecan decreased by a fold change of 3.9 (31 to 8 ESTs). We did not perform such a comparison using the third EST library, because the latter was derived from a mix of time points using PCR amplification of cDNA inserts in order to get a broader overview of the genes involved in the many stages of fracture healing in sheep.

It is striking that 24 of the 39 genes (62%) had previously been known to be involved with developmental processes in some way (references are given in the table). Eleven of the 24 genes were previously known to be involved in bone development or fracture healing. Given that fracture healing is thought to recapitulate bone development on a molecular level [9], these findings are suggestive that our EST analysis captured some essential aspects of fracture healing.

Expression analysis of candidate fracture-healing related genes

We then reasoned that differential expression during the course of fracture healing could be another kind of indicator for genes with potential importance in fracture healing. Therefore, we performed RT-PCR expression analysis using samples from postfracture day 7, 10, 14, and 42. The RNA source for days 7 and 10 was identical to that used for the PF7 and PF10 EST libraries, and the source for day 14 was one of the RNA samples used to construct the pooled EST library. The day 42 RNA sample was not used for any of the EST libraries. Samples from normal sheep bone were not available.

Almost all genes showed a large drop in expression at the 42-day time point. Early stage callus tissue is a metabolically highly active tissue. In comparison, healthy bone is less metabolically active, which is not surprising given that cells make up only around 2% of the total mass of bone. It has previously been noted that the expression of many genes related to fracture healing drops to very low or undetectable levels at 6 weeks after fracture in a rat model [14]. In our experiments, 77 of the 78 genes (99%) tested by ANOVA over all four timepoints were found to be significantly differentially expressed, a proportion that dropped to 96% if only the first three time points were considered, and 31% for a comparison between only the first two time points. [Additional file 2] presents the normalized expression values determined in these experiments.

We then performed hierarchical clustering to produce a graphic representation of the the expression patterns for the genes investigated in RT-PCR. We used the gap-stat method [15] to estimate the number of clusters in the data (Fig. 6). The estimated number of clusters is two. However, the gap function rises again with peaks at 4 and 6 clusters, suggesting that there are less well defined groups of four and of six clusters in the data. For the purposes of discussion and visualization, we chose to partition the genes into six clusters (Fig. 7). Clusters A, B, E and F demonstrated high frequencies of genes of specific functional classes (see below).

Genes differentially expressed in fracture healing

Many of the differentially expressed genes we identified were previously known to be involved in various biological processes involved in bone development or fracture healing (see Tables 2 and 3 for references), suggesting that our approach has captured some essential aspects of fracture healing. Additionally, we have demonstrated differential expression for many genes not previously known to be involved in fracture healing.

The differentially expressed genes can be divided into several different classes, as will be discussed below. Three of these classes, extracellular matrix proteins, resorption/remodeling/inflammation, and transcriptional regulation/signal transduction, have been well studied in the context of fracture healing, and our results have identified both previously known and novel differentially regulated genes with these functional roles.

In the following, we describe three other classes with numerous differentially regulated genes, angiogenesis, free-radical control, and ribosomal genes. These functional classes have received less attention to date, although they are known to be involved in the process of fracture healing. Fracture-related differential expression had not previously been shown for the great majority of genes in these classes (cf. Tables 2 and 3). These observations thus represent starting points for future studies on the role of processes such as angiogenesis and free-radical control in fracture healing.

Fracture model

All animal experiments were carried out according to the policies and principles established by the Animal Welfare Act, the NIH Guide for Care and Use of Laboratory Animals and the national animal welfare guidelines. The study was approved by the local legal representative (Landesamt für Arbeitsschutz, Gesundheitsschutz und technische Sicherheit, Berlin: G 0224/01 and G 0172/04).

The experimental system has been described in detail elsewhere [48, 49, 38]. Briefly, twelve two-year old female Merino mix sheep with a mean weight of 66 ± 12 kg underwent a standardized midshaft osteotomy of the right tibia under general anesthesia. The osteotomy was distracted to a gap of 3 mm and stabilized with a monolateral external fixator. The fixator consisted of 6 Schanz screws (5 mm, 3 inserted proximally, 3 inserted distally of the osteotomy) and 2 steel tubes (10 mm) and was mounted medially to stabilize the defect.

The osteotomy wound was sutured and the leg, together with the external fixator, was covered with a tube bandage. The animals received an analgetic (Finadyne, Essex, Germany) for 3 days postoperatively. Daily pin care involved cleaning of the insertions of the Schanz screws with ethacridin lactate (Rivanol, 5%, Chinosol, Germany).

After removal of the fixator under general anesthesia the callus tissue was removed and sheep were sacrificed. The tissue was frozen in liquid nitrogen and stored at -80°C until further use.

Isolation of RNA and library construction

For generation of ESTs, three separate libraries were constructed. The first library (PF7) was isolated from four animals 7 days post-fracture, the second library (PF10) was isolated from three animals 10 days post-fracture, and the third library was isolated from a pool of one animal each at 2 weeks, 3 weeks, and 4.5 weeks post-fracture.

PF7 consisted of 3.6 × 10⁵ clones with an average insert length of 1.6 kb. PF10 consisted of 3.4 × 10⁵ clones with an average insert length of 1.4 kb. The pooled library consisted of 4.5 × 10⁶ clones with an average insert size of 1.4 kb. Average insert lengths were determined by screening at least 300 inserts per library.

For isolation of RNA the tissue was ground in a liquid-nitrogen-cooled stainless steel mortar. RNA was isolated with peqGOLD TriFast (peqLab) and subsequently purified with RNeasy Mini columns (Qiagen) followed by poly(A)-RNA isolation using the polyATtract mRNA Isolation System III (Promega).

For construction of cDNA libraries the CloneMiner cDNA library construction Kit (Invitrogen) was used according to the manufacturer's protocol with the following modifications. All precipitations were carried out using PelletPaint Co-Precipitant (Novagen) instead of glycogen. The first strand reaction was carried out in five separate reactions with 1 μ g poly(A)-RNA each. Four reactions were precipitated and pooled with the fifth for a single second strand reaction. Glass fiber filters for radioactive determination of the cDNA yield were processed separately. Clones for further analysis were arrayed into 384-well-microtiter plates by robotic colony picking.

Sequencing and analysis of ESTs

Template DNA for sequencing was prepared by PCR amplification of cDNA inserts with M13 primers or by automated plasmid purification. PCR product quality control was done by 384 well agarose gels, positive products where rearrayed and diluted to appropriate concentrations. For plasmid preparation, cultures were grown in 384 deepwell plates containing 200 μ l of 2YT Media. Incubation took 17–18 h at 1100 rpm.

Plasmids were isolated and purified automatically by a modified alkaline lysis procedure. Quality control of plasmids was done by randomly testing eight samples of each plate on agarose gels. Sequences were determined using M13 -40 primer, ABI Big Dye Terminator 3 chemistry and ABI 3730XL capillary sequencer systems (ABI, Weiterstadt, Germany). All raw sequences were processed by PHRED [50, 51] and screened for vector or Escherichia coli contamination.

Before clustering, sheep EST sequences were subjected to an extensive quality control procedure. In a first step sequences were clipped at both ends as long as PHRED quality values stayed below 20 in a window of size 15 bp. Secondly, for all sequences (ESTs + mRNAs) putative vector and repeat sequences were masked. Clustering and assembly of the clipped sequences was performed using the GeneNest software tools [12, 13].

An Ovis aries EST database

GeneNest is designed to generate and visualize gene indices based on ESTs and mRNA sequences [12, 13]. We have used the GeneNest software to create a database with 47209 sheep ESTs generated in the current project, as well as all previously available Ovis aries sequences. The database is available at http://genenest.molgen.mpg.de/cgi-bin/search_db?db=Oa5.

Accession numbers

All ESTs generated in this project were submitted to NCBI's GenBank and are available with the accession numbers [GenBank:DY475559–DY505707] and [GenBank:DY505714–DY522957].

Identification and annotation of genes

In order to identify and annotate novel and previously known sheep genes, similarities to known proteins were identified with a blastx [52] search against the RefSeq protein database [53] of 28 September, 2005 using default settings for blastx. Contigs with a best e-value of 10^-25 or lower were assigned a putative identity.

Gene Ontology

For the purpose of Gene Ontology [54] annotation, the contigs were mapped to mouse mRNAs using blastn with default settings. Putative identities were assigned using an e-value threshold of 10^-25. Mouse-Genome Informatics (MGI) annotations [55] were assigned to the sheep homologs of annotated mouse genes, and GO analysis was performed using the Ontologizer [56].

Differential expression of early fracture genes

In order to gain insight into the genes that are upregulated in the early stages of fracture healing in the sheep, we compared the distribution of ESTs in the PF7 and PF10 sheep libraries against the distribution of ESTs in the UniGene [6] human bone libraries ID.655 [57] and ID. 1124, which have a total of 6573 ESTs. Mapping of sheep contigs to human UniGene clusters was performed using blastn with an e-value cutoff of 10^-25 as described above. We then compared the EST counts for the PF7/PF10 libraries against those in the ID.655/ID.1124 libraries using a two-tailed Fisher-Exact test. The result of the Fisher-exact test is a p-value that reflects the probability of observing the given differences in EST counts in the respective libraries assuming no differential expression [58]. The expression ratio was calculated based on the total counts of 32032 ESTs in libraries PF7 and PF10 and 6573 ESTs in UniGene libraries ID.655/ID.1124. A similar test was applied to compare EST counts in PF7 against those in PF10. If the p-value of this association test is below a certain threshold, the gene is considered as putatively differentially expressed. Correction for multiple testing was performed using a false discovery rate (FDR)-control method [59]. The FDR represents the expected proportion of errors among the rejected null hypotheses and is thus a different kind of procedure from multiple testing correction procedures such as a Bonferroni correction that control the family-wise error rate, or the probability that one or more rejected null hypotheses are falsely rejected. The Benjamini-Hochberg-adjusted p-values are shown in the Tables for the EST comparisons.

RT-PCR

RNA was isolated for quantitative reverse-transcriptase polymerase chain reaction (RT-PCR) from 7-day postfracture (PF7, n = 4 animals), 10-day postfracture (PF10, n = 3), 14 day (PF14, n = 1) and 42-day (PF42, n = 1) tissue. PF7 and PF10 were identical to the RNA sources used for EST library construction, and PF14 was used for the pooled library. For semi-quantitative RT-PCR measurements, single-stranded cDNA was prepared from the different healing stages using random hexamer primers with the Taqman Reverse Transcription Reagents Kit (Applied Biosystems, Forster City, USA). Gene-specific primers were designed using the Primer Express software (Applied Biosystems) and quantitative RT-PCR was performed with SYBR green PCR Master Mix (Applied Biosystems) on an ABI Prism 7900 HT (Applied Biosystems). The mRNA levels were quantified in triplicate by the standard curve method normalized to a exogenous green fluorescent protein (GFP) transcript of which 1 ng was spiked into the sample per 1 μ g of total RNA. See [Additional file 2] for primer sequences and RT-PCR profiles.

Three measurements were made for each gene at each of the four time points. Statistical analysis of the RT-PCR results was performed by a two-sided Student's t test to compare expression between PF7 and PF10, and one-way ANOVA was performed to identify differential expression among the groups PF7, PF10, PF14 or among all four groups. A Bonferroni correction for multiple testing was performed.

Clustering and visualization of RT-PCR results

Complete-linkage hierarchical clustering was used to visualize the results of RT-PCR on the 78 putative fracture-related genes. For each gene, a vector of four values representing the four time points was normalized to have zero mean and unit variance. In order to estimate the optimal number of clusters, the gap statistic method [15] was used.