Skip to main content

Advertisement

Transcriptomic analysis of different tissue layers in antler growth Center in Sika Deer (Cervus nippon)

Abstract

Background

With the unprecedented rapid growth rate (up to 2.75 cm/day), velvet antler is an invaluable model for the identification of potent growth factors and signaling networks for extremely fast growing tissues, mainly cartilage. Antler growth center (AGC) locates in its tip and consists of five tissue layers: reserve mesenchyme (RM), precartilage (PC), transition zone (TZ), cartilage (CA) and mineralized cartilage (MC). The aim of this study was to investigate the transcription dynamics in the AGC using RNA-seq technology.

Results

Five tissue layers in the AGC were collected from three 3-year-old male sika deer using our previously reported sampling method (morphologically distinguishable). After sequencing (15 samples; triplicates/tissue layer), we assembled a reference transcriptome de novo and used RNA-seq to measure gene expression profiles across these five layers. Nine differentially expressed genes (DEGs) were selected from our data and subsequently verified using qRT-PCR. The results showed a high consistency with the RNA-seq results (R2 = 0.80). Nine modules were constructed based on co-expression network analysis, and these modules contained 370 hub genes. These genes were found to be mainly involved in mesenchymal progenitor cell proliferation, chondrogenesis, osteogenesis and angiogenesis. Combination of our own results with the previously published reports, we found that Wnt signaling likely plays a key role not only in stimulating the antler stem cells or their immediate progeny, but also in promoting chondrogenesis and osteogenesis during antler development.

Conclusion

We have successfully assembled a reference transcriptome, generated gene expression profiling across the five tissue layers in the AGC, and identified nine co-expressed modules that contain 370 hub genes and genes predorminantly expressed in and highly relevant to each tissue layer. We believe our findings have laid the foundation for the identification of novel genes for rapid proliferation and chondrogenic differentiation of antler cells.

Background

A growth system, where normal cells exhibit a rapid proliferation and differentiation without becoming cancerous, would be desirable in identification of potent growth factors, unique signal transduction pathways and novel regulation systems. In this regard, deer antler is an invaluable model to meet these requirements. Antlers are male secondary sexual characteristics and each year form anew from the permanent frontal bony protuberances, called pedicles [1, 2]. During the growth phase, elongation of antlers in some large deer species (such as north America wapiti) can exceed 2 cm/day [3], nonetheless with a well-organized tissue structure [4]. The antler growth center (AGC) is located in the antler tip [5] and histologically consists of five tissue layers from distal to proximal: reserve mesenchyme (RM), pre-cartilage (PC), transition zone (TZ), cartilage (CA) and mineralized cartilage (MC) (Fig. 1) [6, 7].

Fig. 1
figure1

Schematic drawing to show the five tissue layers in an antler growth center. RM: reserve mesenchyme; PC: pre-cartilage; TZ: transition zone; CA: cartilage; MC: mineralized cartilage. This figure also presented in the Fuzzy c-means Clustering Analysis of gene expression patterns in seven clusters (see Results)

Antler is an organ of cartilage/bone, and its elongation is achieved through appositional growth [1, 5]. Results from histology, histochemistry, ultrastructure, and BrdU labeling [6, 8, 9] show that the cells in the RM layer are either in mitotically quiescent state (outer sublayer) or at mitotic state (inner sublayer). Below the RM layer, the mesenchymal cells start to differentiate towards a cartilaginous lineage (PC layer), where pre-chondroblasts and newly-formed-isolated vascular streaks reside. Further below the PC layer, pre-chondroblasts gradually maturate toward chondroblasts and chondrocytes; at the same time, the rod-like vascular streaks progressively extend from the two ends in parallel with the long axis of the antler, to form the TZ layer. Within the TZ layer, different stages of cartilaginous cells including prechondroblasts, chondroblasts and chondrocytes co-exist; with different stages of the vascular system (broken and continuous vascular channels co-exist) formation. Finally, continuous cartilage columns (which mainly contain chondroblasts in the periphery and chondrocytes in the center) alternate with continuous vascular channels (functional vascular system) to form the CA layer, the vascularized cartilage (unique feature of antler cartilage). At the end of this proximal layer, chondroclasia, osteogenesis and osteoclasia take place simultaneously to form the MC layer, within which the smooth-surfaced the osteo-cartilage columns are converted into irregular and broken trabeculae [4].

In order to facilitate the discovery of novel genes and/or regulatory systems for rapid antler growth and chondrogenesis using molecular techniques, we established a standardized method to allow rapid and precise sampling of each of these five tissue layers in the AGC of a fresh-cut antler based on morphologically distinguishable markers, and without having to let them undergo histological processes in order to do so [7, 9]. Unfortunately, thus far not much progress has been made in the discovery of the novel genes from the AGC since publication of the method more than a decade ago. This undesirable situation would have been at least partially attributed to the reason that the relevant studies neither used the high-throughput RNA sequencing (RNA-seq) technique, hence a large number of the genes were not detected due to heterologous microarray (deer cDNA versus mouse template) [10, 11]; nor applied the tissue layer sampling method in the study, hence failed to put those identified genes in the biological context [12, 13].

The aim of present study was to use the novel antler model again, but this time to combine the RNA-seq technique with the published standardized tissue layer sampling method, to seek to identify novel genes and regulatory pathways that underpin the system where the fastest growth and chondrogenesis are executed.

Results and discussion

Sequencing, de novo assembly and transcript abundance

A total of 654 million (86.58 Gbp) of clean paired-end reads from 15 libraries (triplicates/tissue layer) passed the quality filters (Additional file 1: Table S1). Then, a de novo assembly pipeline was applied to these reads to generate reference transcriptome, because of lacking of a proper sika deer genome at present time (Additional file 2: Figure S1). In total 88,369 non-redundant transcripts (≥ 300 bp) with FPKM ≥0.5 were generated through this process. Based on three different coding sequence prediction methods, a total of 44,177 high-quality coding transcripts were subsequently obtained (Additional file 2: Figure S1, Additional file 3: Table S2). Of these 44,177 high-quality transcripts, the length of N50 (defined as the sequence length of the shortest contig at 50% of the total genome length) was 2533 bp, which is longer than the previously published transcriptomes of the deer [12, 14, 15]. Length distribution results of the assembled transcripts showed that the coding transcripts were more abundant than the non-coding sequences in the bins of long transcripts (Additional file 4: Figure S2A). Transcript abundance was ranged from 3 to 4 orders of magnitude (Additional file 4: Figure S2B), and the coding transcripts (40–50%) were more abundant in the upper ranks of the distribution (Q4) than the non-coding sequences (20–30%) (Additional file 4: Figure S2C). Core Eukaryotic Genes Mapping Approach (CEGMA) [16] was applied to evaluate completeness of our transcript assembly, and found that a high percentage of core genes (94.35% had complete sequences; only 5.25% were partially sequenced) presented in our transcriptome data.

As a gene can have multiple transcript isoforms with different length, the longest transcript of a gene in this study was selected to represent that gene. A total of 13,203 genes were generated and then annotated using the top 1 orthologous of all currently known species based on UniProt database (May 2017), including human, mouse, cattle and sheep. We further screened the characteristics of the assembled 13,203 gene sequences by comparing them with protein sequences (not incorporated in UniProt) derived from genomes of closely related species, i.e., red deer [17], white tailed deer and cattle. For this analysis, BLASTX v2.5.0+ with default search parameters was carried out. The results showed that 13,131 (99.5%) and 13,110 (99.3%) out of 13,203 genes passed BLAST matched thresholds of E-value ≤10− 5 for cattle and white tailed deer respectively. But, the matched ratio with red deer proteins (12,145, 92.0%) was relatively low (Additional file 5: Table S3). The possible reason for this is that the current red deer genome is still not a completed version (1.96 Gbp in total). Our results also suggested that ~ 0.5–0.7% gene sequences could present the low conservativeness between sika deer and these closely related species. Interestingly, the slight redundancy sequences (~ 3.2%, Additional file 5: Table S3) could reside in these 13,203 genes although our transcriptome dataset was clustered using cd-hit-est tool (see methods). In addition, our results showed that ~ 60% out of deer genes is likely to be expressed in the AGC (Additional file 5: Table S3), which is compatible with the number of genes expressed in other mammalian tissues (Fig. 1 refer to [18]). Overall, a high-quality reference gene dataset (13,203 genes) was successfully generated, and was further used for the downstream bioinformatics analysis in this study (Additional file 6: Figure S3).

Gene expression profiling matches developmental state of tissue layers

To assess the global variation in gene expression across the AGC tissue layers, a matrix to the normalized expression values (mapping reads, 13,203 genes) were conducted. The results of pairwise Pearson Correlation indicated that the RM layer showed a stand lone group; whereas the rest of other four layers (PC, TZ, CA and MC) clustered together as a composite group, although the first three layers (PC, TZ and CA) in the composite group were closer compared to the last layer (MC) (Fig. 2A). This finding may reflect the use of distinct regulatory genes in the RM layer compared to the other layers, possibly because the cells in the RM layer have stem cell attributes while cells in the other layers are in the different differentiating states. Principal Component Analysis were also performed, and the results showed that PC1 explained 22.9% of the overall variation, which was found too low to separate the five tissue layers but could separate the RM layer (outmost) from the MC layer (innermost). PC2 explained 19.7% of the overall variation, which was found to allow either the RM layer or the MC layer to be distinguishable from the composite group (PC, TZ and CA) (Fig. 2B). Interestingly, the PC1 variation revealed a gradient of layer samples that perfectly fits their position in a proximo-distal axis of the AGC. Overall, our results suggest that the variation trajectory follows the similar pattern of the tissue layers in the AGC.

Fig. 2
figure2

Comparison of gene expression across five tissue layers in the antler growth center. a) Pairwise Pearson Correlations of gene expression across 15 samples (triplicates/each tissue layer). b) Principal Component Analysis of gene expression across 15 samples. c) Ratio of up-regulated to down-regulated genes across five tissue layers on different thresholds (|log2FoldChange|1, 2, 3 and 4)

In order to detect the changes in ratio of up-regulated to down-regulated genes across the five tissue layers inversely (i.e. a proximal layer vs its every distal layer), we used the four pre-set levels of thresholds (|log2FoldChange|1, 2, 3 and 4). The results showed that the higher the level the bigger the ratio (Fig. 2C), suggesting that the elevated level of up-regulated genes is positively correlated with the degree of cell differentiation (mesenchymal cells to chondrocytes) in the AGC.

Extensive changes detected in functional categories across tissue layers

To investigate the large-scale patterns of gene expression across the five tissue layers, Fuzzy c-means clustering model was used to group 5585 DEGs (selected based on this criteria (|log2foldchange| ≥ 1.5, adjusted Pvalue ≤0.001)), and these DEGs were grouped into seven clusters. Of these 5585 DEGs, 2740 were selected based on their member score (MS) within the seven clusters (MS 0.5) and used for further analysis (Fig. 3, also refer to Fig. 1). Next, we performed GO enrichment analyses using genes in each cluster to identify key biological process (BP) categories (Additional file 7: Table S4). The genes in Cluster 1 mainly included those with expression levels steadily increasing over sequential differentiation processes in the AGC (from the RM layer to MC layer). As expected, these genes were mainly related to key BP categories, e.g., extracellular matrix disassembly (EASE score = 0.002, e.g., CTSK, CD44, MMP13 and MMP15), osteoclast differentiation (EASE score = 0.007, e.g., CSF1 and TNFRSF11A), angiogenesis (EASE score = 0.015, e.g., TNFSF12 and SOX18) and wound healing (EASE score = 0.049, e.g., TIMP1). The genes in Cluster 2 included those showing an opposite trend in expression level to Cluster 1; these genes were mainly related to Wnt signaling and cell polarity (EASE score = 0.033, e.g., FZD1, FZD2, PSMA1 and PSMA4), and osteoblast proliferation (EASE score = 0.002, e.g., OSR2). These findings are consistent with our histological observations in that the transition from undifferentiated mesenchymal cells to chondrocytes distoproximally across the tissue layers is gradual, along with vascularization, chondroblast maturation and chondrification [9]. In addition, these genes may also play a role in the formation of vascularized cartilage, a unique structure which is thought to be required to meet the metabolic demands for rapidly growing antler tissue [6, 9] and as a conduit for hemopoietically derived chondroclast and osteogenic progenitors [19].

Fig. 3
figure3

Fuzzy c-means clustering identifies general patterns of gene expression across the five tissue layers in the antler growth center (also refer to Fig. 1). The average FPKM values of DEGs (|log2FoldChange| ≥ 1.5, adjusted Pvalue ≤0.001) are used. The MS of a given gene within a cluster is represented in color, with red (MS = 1) indicating high association. The number of genes with MS ≥ 0.5 in clusters is also shown. For each cluster, the key genes involved in significantly enriched GO BP categories (Additional file 7: Table S4) were shown

For clusters 3, 4, 5, 6 and 7, up-regulated layer-specific genes were found to be with expression peaks in each of the five tissue layers (RM, PC, TZ, CA and MC layers). The genes in Cluster 3 had expression peaks in the RM layer and were found to be involved in cell proliferation (EASE score = 0.009, e.g., IGF1, EGFR), stem cell population maintenance (EASE score = 0.044, e.g., NOTCH2 and KIT), regulation of apoptotic process (EASE score = 0.02, e.g., CASP8, CASP8AP2 and TP53BP2) and DNA repair (EASE score = 0). The rate of antler growth can reach 2 cm per day [3], and this rapid growth is achieved mainly through fast cell proliferation in the RM layer [6, 9]. The rapid rate of cell proliferation in antlers would not only require factors that stimulate fast cell proliferation, but also those that properly control the cell cycle and protect genomic stability. It is also notable that genes that were involved in epigenetic regulation, including gene silencing by RNA (EASE score = 0.001, e.g., AGO1 and DICER1) and histone modification (EASE score = 0.01, e.g., CTR9 and EP300) were found in the RM layer, strongly implying early imprinting of the stem cell populations in this layer.

The genes in Cluster 4 showed an expression peak in the PC layer, and were found to be related to positive regulation of Wnt signaling (EASE score = 0, e.g., FGFR2 and WNT10B). The PC layer mainly consists of newly differentiated prechondroblasts, which actively form and secrete extracellular matrix [8]. Wnt signaling was detected in this Cluster (also in Clusters 2 and 3), suggesting this signaling pathway could be involved in early antler cell differentiation toward chondrogenesis.

The genes in Cluster 5 showed an expression peak in the TZ layer, and were found to be related to cartilage condensation (EASE score = 0.02, e.g., BMP1, SOX9 and COL2A1) and response to hypoxia (EASE score = 0.03, e.g., AKT1 and ANGPTL4). In the TZ layer, cells are in the differentiation transition from prechondroblasts to chondroblasts and chondrocytes, with vascular channels progressing from broken (unfunctional) to continuous (functional) [9]. Cartilaginous cells at different stages of differentiation in the vicinity of unfunctional vascular channels (not yet linked to the body vascular system), or further away from the functional vascular system (linked to the body vascular system) would be subjected to low oxygen tension/hypoxia. It is not surprising that, in order to cope with this low oxygen tension, these cells would have to highly express genes that are responsive to hypoxia.

The genes in Clusters 6 showed expression peaks in the CA layers. As expected, these genes were found to be related to extracellular matrix organization (EASE score = 0, e.g., COL10A1 and ACAN), angiogenesis (EASE score = 0.004, e.g., MMP14), regulation of chondrocyte differentiation (EASE score = 0.015, e.g., CTGF), ossification (EASE score = 0.03, e.g., MMP9 and TNFSF11) and bone mineralization (EASE score = 0.006, e.g., PTH1R), and response to hormones (EASE score = 0.03, e.g. LYN). The genes in Clusters 7 had expression peaks in the MC layers, and were found to regulate osteoclast differentiation (EASE score = 0.044, e.g., BGLAP) and negative regulation of blood coagulation (EASE score = 0.044, e.g., CD34). These findings support the previous ultrastructural observations, in which almost all hypertrophic chondrocytes were eventually subject to degeneration and apoptosis [8]. Programmed cell death of hypertrophic chondrocytes in both the CA and MC layers would create more space for the brought-in-osteogenic progenitor cells to build up bone tissue [20]. The seasonal high levels of androgen hormones would initiate the final antler mineralization processes, further result in the interruption of the blood flow into antlers from their bases, and cause eventual demise of antlers [21, 22]. Overall, gene expression profiling across the five tissue layers within an AGC matched the developmental states in each corresponding tissue layer.

Hub genes detected using co-expression network analysis and associated with Chondrogenesis and angiogenesis

Sequences of the obtained 13,203 genes were also analyzed to investigate the association with the chondrogenesis and angiogenesis in the AGC using Weighted Gene Co-expression Network Analysis (WGCNA). The results showed that majority of these genes (13,103 genes; 99.24%) were assigned to 34 modules (39–3329 genes per module). Of these 34 modules, nine (10,267 genes: 78.36%, marked with asterisk in Fig. 4A) were found to have significantly negative/positive module-trait correlations with the corresponding tissue layers when the |Pearson Correlation| ≥ 0.6 and Pvalue ≤0.01 criterion was applied, and to be enriched in the BP categories (Fig. 4B, Additional file 8: Table S5). These BP categories were concordant with the counterparts generated from our Fuzzy c-means clustering analysis. For example, the MEfloralwhite module showed a positive module-trait correlation with the CA layer and these enriched BP categories were consistent with those in Cluster 6, and partially in Cluster 1. Likewise, the MElightcyan module showed a positive module-trait correlation with the PC layer, such as Wnt signaling pathway, which is consistent with that in Cluster 4.

Fig. 4
figure4

Co-expression gene networks. a) Genes (13,103) are assigned to 34 different modules (39 to 3329 genes/module). Nine modules marked with asterisks have significant BP categories identified by a hypergeometric test after adjusted Pvalue ≤0.05 and show significant negative/positive module-trait correlations with the corresponding tissue layers (|Pearson Correlation| ≥ 0.6 and Pvalue ≤0.01). b) For these nine modules, significantly enriched GO biological process categories are shown (Additional file 8 Table S5). For clarity, only the selected key categories are shown in the figure

Next, we focused on highly connected intra-module hub genes in each module, and these genes had high module membership (MM) values (MM 0.7) in their respective modules. We detected 370 hub genes that might be involved in rapid antler growth, and in regulation of chondrogenesis and angiogenesis (Fig. 5). Of these hub genes, 187 (50.5%) were involved in seven Fuzzy c-means clusters, and more specifically, these associated hub genes from each module belong to one or two clusters (except for 26 genes in MEgreenyellow module), further suggesting that the results between WGCNA and Fuzzy c-means analysis are consistent. Twenty nine hub genes are also shown in the previously published studies using other molecular technologies (e.g., western blot, qRT-PCR, in-situ hybridization or cDNA microarrays) (Table 1); and were found to be involved in extracellular matrix organization, mineralization and degradation (e.g., COL1A1, COL2A1, COL10A, MMP9, MMP13, SPARC, IBSP, BGLAP, CTSK, ALPL MGP and MATN1), signaling molecules (e.g., IGF1, VEGF, IHH, WNT4, PTHLH, CSF1, TNFSF11 and TGFB1), receptors (e.g., FGFR1, FGFR3, PTH1R, CSF1R, CALCR and TNFRSF11A), binding (e.g., APOD) and transcriptional factors (e.g. RUNX3 and SOX9). Some of the molecules and their receptors (i.e. FGFR1, FGFR3, IHH, SOX9, PTHLH and PTH1R) in our hub gene pool are found to be expressed during embryogenesis, and reportedly regulate long bone formation and growth plate development [23, 24]. This finding supports a previous hypothesis that there is evolutionary conservation of the developmental signaling pathways that occur during embryogenesis and postnatal cartilage/bone regeneration [25].

Fig. 5
figure5

Expression of 370 highly connected intra-module hub genes across the five tissue layers. All hub genes are derived from nine modules marked with asterisks in Additional file 6: Figure S3A with high MM values (MM ≥ 0.7). Each line represents the scaled log2 transition of average FPKM value for each individual gene. Among them, 29 genes being highlighted with yellow are reported previously and 187 genes that emphasized with blue font are involved in seven Fuzzy c-means clusters. Particularly, the Wnt related genes and their receptors (Frizzled genes) are emphasized with bold italic. ‘@’: MM value ranged from 0.7 to 1. ‘#’: genes related to angiogenesis. ‘$’: genes related to cartilage/bone development. ‘C’: Cluster ID related to Fuzzy c-mean Clustering Analysis

Table 1 Summary of data for 29 hub genes reported in the previous studies of antler tip using molecular technologies

Six Wnt genes (WNT2B, WNT5A, WNT5B, WNT4, WNT10B and WNT11), and four their receptor genes (Frizzled genes; FZD1, FZD4, FZD5 and FZD9) in our hub gene pool were overrepresented. It is known that the Wnt signaling pathway plays an essential role in cartilage/bone development during embryogenesis [26,27,28]. Furthermore, the canonical Wnt signaling pathway is known to be involved in the establishment of AGC in the early antler regeneration, more specificallyβ-catenin is an important factor controlling survival and lineage specification of the mesenchymal progenitor cells toward chondrogenesis [29]. In our results, both Wnt genes and their receptors were identified and fell in seven of nine co-expressed modules, suggesting that this pathway may not only play an important role in promoting the proliferation of the stem cells or their immediately differentiated progeny in the AGC, but also be involved in chondrogenesis during antler development. Therefore, we think that further studies should be directed at elucidation of the mechanisms underlying the cross-interactions between Wnt signaling pathway and other local signaling pathways, such as TGF-beta and Hedgehog signaling pathways, as well as hormonal stimuli, for the regulation of antler cell proliferation, survival and chondrogenesis.

Four Fox genes (FOXC1, FOXC2, FOXO1 and FOXO4) and four Sox genes (SOX8, SOX9, SOX17 and SOX18) were found to be overrepresented in our hub gene pool. FOXC2 is reported to augment tumor propagation and metastasis in osteosarcoma [30]. In contrast, FOXO1 and FOXO4 are regarded as tumor suppressor genes for certain cancers including osteosarcoma through diversified mechanisms, such as initiating apoptosis [31]. Antler growth involves fast cell proliferation that is elegantly regulated without becoming cancerous, suggesting FOXO genes play a role in the maintenance of normal antler tissue growth.

In order to verify our RNA-seq results, nine hub genes (i.e., DLX3, FOXC2, FRZB, JUNB, SMAD6, SOX18, SRPX2, TNMD and TWIST2) were selected from our hub gene pool (370 hub genes in total) based on the criteria of |log2foldchange| 2 and adjusted Pvalue ≤0.001, and validated using qRT-PCR. Results of 25 pairwise PCR reactions showed high consistency with that of RNA-seq data (R2 = 0.80) (Fig. 6).

Fig. 6
figure6

Correlation between RNA-seq and qRT-PCR for nine highly differentially expressed hub genes. Each color dot represents a qRT-PCR pairwise reaction with the corresponding gene

Conclusions

Overall, we believe our results (such as 370 hub genes involved in nine co-expressed modules, particular genes that predominantly expressed in and highly revenant to each tissue layer) provide a foundation for future studies of more detailed molecular analyses to the development of the AGC. Besides, our transcriptome data would be valuable for other functional genomic research in sika deer or closely-related species. In the long term, establishing the molecular and cellular mechanisms involved in regulating chondrogenesis may lead to the development of strategies for enhancing cartilage/bone repair and regeneration in other mammals including humans.

Methods

Sample collection

Antler tissues were collected from three 3-year-old healthy sika deer (Cervus nippon) at about 30 days after casting of the previous hard antlers, and processed as previously described [9]. Briefly, the distal 8 cm of the growing tip was removed and sectioned sagittally along the longitudinal axis. Five layers of tissues of the tip were immediately dissected and then cut into 4–6 mm pieces, which were then frozen in liquid nitrogen and stored at − 70 °C for RNA preparation and sequencing.

RNA preparation and sequencing

The tissue samples were rapidly grinded into fine powder by using Freezer/Mill 6770 (SPEX CertiPrep Ltd., USA). Total RNA was extracted from the sample powder using a Trizol reagent (Invitrogen Inc., Camarillo, CA) according to the manufacturer’s procedure. RNA quality was confirmed using Bioanalyzer with a minimum RNA integrity number of 7.0. Six micrograms of total RNA were used to construct libraries according to the manufacturer’s instructions (Illumina TruSeq Library Preparation Kit v3). Libraries were sequenced using an Illumina HiSeq X Ten at BGI (Shenzheng, China). We sequenced three biological replicates of each layer of tissue with 150 bp paired-end sequencing.

Transcriptome assembly, annotation and differential expression

We used Trinity v2.4.1 [32] with fixed default k-mer size of 25 to carry out de novo assembly of a reference transcriptome from the quality-filtered reads. The paired-end reads were mapped to the assembled transcriptome with Bowtie2 v2.0.5 (−no-mixed --no-discordant --gbar 1000 --end-to-end -k 200) [33], and abundance estimation (FPKM, fragments per kilobase of transcript per million mapped reads) was performed using RSEM v1.3.0 [34]. Next, we developed a strict pipeline to filter the assembly errors and background sequences (Additional file 2: Figure S1). Briefly, 1) Removal of short sequences (length ≤ 300 bp); 2) Removal of the background sequences (FPKM of all replicates in any one tissue ≥0.5); 3) Removal of the redundancy transcripts (identity ≥95%) using cd-hit-est v3.0.3 [35]; 4) Removal of non-coding transcripts. The coding transcripts (≥ 100 amino acid) were predicted by ESTScan v2.2.1 [36] with human model and TransDecoder v2.0.1 [32], and further annotated by searching against UniProt database using BLASTX (E-value ≤10− 5). Finally, differential gene expression analysis was conducted based on the mapped counts using DESeq2 v2.1.18 R package [37] at an adjusted Pvalue of 0.001.

Fuzzy c-means Custer analysis

The mean FPKM values were clustered using Fuzzy c-means clustering from Mfuzz v2.42 R package [38]. Only genes with significant differences in expression between at least two layers (|log2FoldChange| ≥ 1.5, adjusted p-value ≤0.001) were used as input for this clustering analysis. The optimal number of clusters was set to 7 and the fuzzifier coefficient set to 2.01. The number of clusters was determined at which the minimum centroid distance plateau was achieved using the Dmin function (Additional file 9: Figure S4). Genes with a membership score (MS) of at least 0.5 were plotted and used as input for categorical enrichment analysis.

Weighted gene co-expression network analysis

A co-expression network was constructed using WGCNA v1.48 R package [39]. In brief, the gene FPKM matrix was subjected to a variance stabilizing transformation using DESeq2 v2.1.18 R package [37]. A soft-threshold power value of 16 was selected in this analysis, which corresponds to R2 (> 0.9). Co-expression modules were identified as clusters from the dendrogram using the cutreeDynamic function with a minimum module size of 30 genes. Modules with eigengene correlations no less than 0.75 were subsequently merged using the mergeCloseModules function with a height cutoff of 0.25. We then conducted module-trait correlations between the module eigengenes and libraries corresponding to each of the developmental layers.

Gene ontology enrichment analysis

We used two bioinformatics tools (DAVID version 6.8 website [40] and GOstats v2.44 R package [41]) to perform GO enrichment analysis for the gene set from Fuzzy c-means cluster and WGCNA analysis. The DAVID version 6.8 website was used to obtain overrepresented GO BP categories with an adjusted Fisher exact P value (EASE score). The GOstats v2.44 R package was performed with Benjamini and Hochberg correction using p.adjust program in R package. The human orthologs of the corresponding deer genes were used in the GO enrichment tests to take advantage of the more complete GO annotation available for human genes.

Quantitative real-time PCR (qRT-PCR) analyses

To confirm the DEGs identified from RNA-seq assay, nine highly expressed genes (see results) with great alteration expression levels were chosen and validated using qRT-PCR. The specific primers located in gene coding regions were designed using Primer 5 software and listed in Additional file 10: Table S6. Actin, cytoplasmic 1 (ACTB) was used as a standard control according to our in-house selection standard. Total RNA was firstly treated with DNase I before reverse transcription by superscript III double-stranded cDNA synthesis kit (Invitrogen Inc., Camarillo, CA). The qRT-PCR was then performed using the SYBR Kit (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s protocol using an Applied Biosystems 7500 detection system. The melting curves for verification of amplification specificity by a thermal denaturing step. The relative quantitative method (2-ΔΔCT) was used to calculate the fold change in the expression levels of target genes [42]. All reactions were performed in three biological replicates using the independent RNA samples. Linear regression analysis and loess smooth plot were performed by ggplot2 R package [43].

Abbreviations

AGC:

Antler growth center

BP:

Biological process

CA:

Cartilage

CEGMA:

Core Eukaryotic Genes Mapping Approach

DAVID:

Database for Annotation, Visualization and Integrated Discovery

DEGs:

Differentially expressed genes GO: gene ontology

FPKM:

Fragments per kilobase of transcript per million mapped reads

MC:

Mineralized cartilage

MM:

Module membership

MS:

Membership score

PC:

Pre-cartilage

qRT-PCR:

Quantitative reverse transcription poly chain reaction

RM:

Reserve mesenchyme

RNA-Seq:

RNA sequencing

TZ:

Transition zone

WGCNA:

Weighted correlation network analysis

References

  1. 1.

    Goss RJ: Deer Antlers. Regeneration, Function and Evolution. New York, NY: Academic Press; 1983.

  2. 2.

    Li C, Suttie JM, Clark DE. Histological examination of antler regeneration in red deer (Cervus elaphus). The anatomical record Part A, Discoveries in molecular, cellular, and. Evol Biol. 2005;282(2):163–74.

  3. 3.

    Goss RJ. Problems of antlerogenesis. Clin Orthop. 1970;69:227–38.

  4. 4.

    Li C. Histogenetic aspects of deer antler development. Front Biosci. 2013;5:479–89.

  5. 5.

    Chapman DI. Antlers-bones of contention. Mammal Rev. 1975;5(4):121–72.

  6. 6.

    Banks WJ, Newbrey JW. Light microscopic studies of the ossification process in developing antlers. In: Antler development in Cervidae: 1982; Caesar Kleberg Wildl. Res. Inst., Kingsville, Texas; 1982. p. 231–60.

  7. 7.

    Li C, Suttie JM. Tissue collection methods for antler research. Eur J Morphol. 2003;41(1):23–30.

  8. 8.

    Szuwart T, Kierdorf H, Kierdorf U, Clemen G. Ultrastructural aspects of cartilage formation, mineralization, and degeneration during primary antler growth in fallow deer (Dama dama). Annals of anatomy = Anatomischer Anzeiger : official organ of the Anatomische Gesellschaft. 1998;180(6):501–10.

  9. 9.

    Li C, Clark DE, Lord EA, Stanton JA, Suttie JM. Sampling technique to discriminate the different tissue layers of growing antler tips for gene discovery. Anat Rec. 2002;268(2):125–30.

  10. 10.

    Molnar A, Gyurjan I, Korpos E, Borsy A, Steger V, Buzas Z, Kiss I, Zomborszky Z, Papp P, Deak F, et al. Identification of differentially expressed genes in the developing antler of red deer Cervus elaphus. Mol Gen Genomics. 2007;277(3):237–48.

  11. 11.

    Gyurjan I Jr, Molnar A, Borsy A, Steger V, Hackler L Jr, Zomborszky Z, Papp P, Duda E, Deak F, Lakatos P, et al. Gene expression dynamics in deer antler: mesenchymal differentiation toward chondrogenesis. Mol Gen Genomics. 2007;277(3):221–35.

  12. 12.

    Yao B, Zhao Y, Zhang H, Zhang M, Liu M, Liu H, Li J. Sequencing and de novo analysis of the Chinese sika deer antler-tip transcriptome during the ossification stage using Illumina RNA-Seq technology. Biotechnol Lett. 2012;34(5):813–22.

  13. 13.

    Yao B, Zhao Y, Wang Q, Zhang M, Liu M, Liu H, Li J. De novo characterization of the antler tip of Chinese sika deer transcriptome and analysis of gene expression related to rapid growth. Mol Cell Biochem. 2012;364(1–2):93–100.

  14. 14.

    Jia BY, Ba HX, Wang GW, Yang Y, Cui XZ, Peng YH, Zheng JJ, Xing XM, Yang FH. Transcriptome analysis of sika deer in China. Mol Gen Genomics. 2016;291(5):1941–53.

  15. 15.

    Yang XG, Chen YX, Liu XD, Liu QZ, Pi XM, Liu YH, Zheng D. De novo characterization of velvet skin transcriptome at the antlers tips of red deer (Cervus elaphus) and analysis of growth factors and their receptors related to regeneration. Pakistan Journal of Zoology. 2016;48(1):151–7.

  16. 16.

    Parra G, Bradnam K, Korf I. CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics. 2007;23(9):1061–7.

  17. 17.

    Bana NA, Nyiri A, Nagy J, Frank K, Nagy T, Steger V, Schiller M, Lakatos P, Sugar L, Horn P, et al. The red deer Cervus elaphus genome CerEla1.0: sequencing, annotating, genes, and chromosomes. Molecular genetics and genomics : MGG. 2018.

  18. 18.

    Yu Y, Fuscoe JC, Zhao C, Guo C, Jia M, Qing T, Bannon DI, Lancashire L, Bao W, Du T, et al. A rat RNA-Seq transcriptomic BodyMap across 11 organs and 4 developmental stages. Nat Commun. 2014;5:3230.

  19. 19.

    Faucheux C, Nesbitt SA, Horton MA, Price JS. Cells in regenerating deer antler cartilage provide a microenvironment that supports osteoclast differentiation. J Exp Biol. 2001;204(Pt 3):443–55.

  20. 20.

    Colitti M, Allen SP, Price JS. Programmed cell death in the regenerating deer antler. J Anat. 2005;207(4):339–51.

  21. 21.

    Bubenik GA, Bubenik AB, Stevens ED, Binnington AG. The effect of neurogenic stimulation on the development and growth of bony tissues. J Exp Zool. 1982;219(2):205–16.

  22. 22.

    Suttie JM, Fennessy PF, Lapwood KR, Corson ID. Role of steroids in antler growth of red deer stags. J Exp Zool. 1995;271(2):120–30.

  23. 23.

    Kronenberg HM. Developmental regulation of the growth plate. Nature. 2003;423(6937):332–6.

  24. 24.

    Provot S, Schipani E. Molecular mechanisms of endochondral bone development. Biochem Biophys Res Commun. 2005;328(3):658–65.

  25. 25.

    Faucheux C, Nicholls BM, Allen S, Danks JA, Horton MA, Price JS. Recapitulation of the parathyroid hormone-related peptide-Indian hedgehog pathway in the regenerating deer antler. Developmental dynamics : an official publication of the American Association of Anatomists. 2004;231(1):88–97.

  26. 26.

    Poss KD, Shen J, Keating MT. Induction of lef1 during zebrafish fin regeneration. Developmental dynamics : an official publication of the American Association of Anatomists. 2000;219(2):282–6.

  27. 27.

    Caubit X, Nicolas S, Le Parco Y. Possible roles for Wnt genes in growth and axial patterning during regeneration of the tail in urodele amphibians. Developmental dynamics : an official publication of the American Association of Anatomists. 1997;210(1):1–10.

  28. 28.

    Stoick-Cooper CL, Weidinger G, Riehle KJ, Hubbert C, Major MB, Fausto N, Moon RT. Distinct Wnt signaling pathways have opposing roles in appendage regeneration. Development. 2007;134(3):479–89.

  29. 29.

    Mount JG, Muzylak M, Allen S, Althnaian T, McGonnell IM, Price JS. Evidence that the canonical Wnt signalling pathway regulates deer antler regeneration. Dev Dynam. 2006;235(5):1390–9.

  30. 30.

    Gozo MC, Jia DY, Aspuria PJ, Cheon DJ, Miura N, Walts AE, Karlan BY, Orsulic S. FOXC2 augments tumor propagation and metastasis in osteosarcoma. Oncotarget. 2016;7(42):68792–802.

  31. 31.

    Fu Z, Tindall DJ. FOXOs, cancer and regulation of apoptosis. Oncogene. 2008;27(16):2312–9.

  32. 32.

    Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.

  33. 33.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

  34. 34.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics. 2011;12.

  35. 35.

    Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.

  36. 36.

    Iseli C, Jongeneel CV, Bucher P. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proceedings International Conference on Intelligent Systems for. Mol Biol. 1999:138–48.

  37. 37.

    Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10).

  38. 38.

    Kumar L. M EF: Mfuzz: a software package for soft clustering of microarray data. Bioinformation. 2007;2(1):5–7.

  39. 39.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics. 2008;9:559.

  40. 40.

    Huang d W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

  41. 41.

    Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2007;23(2):257–8.

  42. 42.

    Schefe JH, Lehmann KE, Buschmann IR, Unger T, Funke-Kaiser H. Quantitative real-time RT-PCR data analysis: current concepts and the novel "gene expression's CT difference" formula. J Mol Med. 2006;84(11):901–10.

  43. 43.

    H W. ggplot2: elegant graphics for data analysis. New York: Springer Publishing Company; 2009.

  44. 44.

    Sadighi M, Haines SR, Skottner A, Harris AJ, Suttie JM. Effects of insulin-like growth factor-I (IGF-I) and IGF- II on the growth of antler cells in vitro. The Journal of endocrinology. 1994;143(3):461–9.

  45. 45.

    Suttie JM, Gluckman PD, Butler JH, Fennessy PF, Corson ID, Laas FJ. Insulin-like growth factor 1 (IGF-1) antler-stimulating hormone? Endocrinology. 1985;116(2):846–8.

  46. 46.

    Elliott JL, Oldham JM, Ambler GR, Bass JJ, Spencer GS, Hodgkinson SC, Breier BH, Gluckman PD, Suttie JM. Presence of insulin-like growth factor-I receptors and absence of growth hormone receptors in the antler tip. Endocrinology. 1992;130(5):2513–20.

  47. 47.

    Price JS, Oyajobi BO, Nalin AM, Frazer A, Russell RG, Sandell LJ. Chondrogenesis in the regenerating antler tip in red deer: expression of collagen types I, IIA, IIB, and X demonstrated by in situ nucleic acid hybridization and immunocytochemistry. Developmental dynamics : an official publication of the American Association of Anatomists. 1996;205(3):332–47.

  48. 48.

    Li C, Harper A, Puddick J, Wang W, McMahon C. Proteomes and signalling pathways of antler stem cells. PloS one. 2012;7(1):e30026.

  49. 49.

    Wang ST, Gao YJ, Duan CC, Li DD, Tian XC, Zhang QL, Guo B, Yue ZP. Effects of PTHrP on expression of MMP9 and MMP13 in sika deer antler chondrocytes. Cell biology international. 2013;37(12):1300–7.

  50. 50.

    Faucheux C, Nesbitt SA, Horton MA, Price JS. Cells in regenerating deer antler cartilage provide a microenvironment that supports osteoclast differentiation. The Journal of experimental biology. 2001;204(Pt 3):443–55.

  51. 51.

    Korpos E, Molnar A, Papp P, Kiss I, Orosz L, Deak F. Expression pattern of matrilins and other extracellular matrix proteins characterize distinct stages of cell differentiation during antler development. Matrix biology. journal of the International Society for Matrix Biology. 2005;24(2):124–35.

  52. 52.

    Clark DE, Lord EA, Suttie JM. Expression of VEGF and pleiotrophin in deer antler. The anatomical record Part A, Discoveries in molecular, cellular, and evolutionary biology. 2006;288(12):1281–93.

  53. 53.

    Guo B, Wang ST, Duan CC, Li DD, Tian XC, Wang QY, Yue ZP. Effects of PTHrP on chondrocytes of sika deer antler. Cell and tissue research. 2013;354(2):451–60.

  54. 54.

    Faucheux C, Horton MA, Price JS. Nuclear localization of type I parathyroid hormone/parathyroid hormone-related protein receptors in deer antler osteoclasts: evidence for parathyroid hormone-related protein and receptor activator of NF-kappaB-dependent effects on osteoclast formation in regenerating mammalian bone. Journal of bone and mineral research : the official journal of the American Society for Bone and Mineral Research. 2002;17(3):455–64.

  55. 55.

    Zhang HL, Yang ZQ, Duan CC, Geng S, Wang K, Yu HF, Yue ZP, Guo B. WNT4 acts downstream of BMP2 to mediate the regulation of ATRA signaling on RUNX1 expression: Implications for terminal differentiation of antler chondrocytes. Journal of cellular physiology. 2017.

  56. 56.

    Steger V, Molnar A, Borsy A, Gyurjan I, Szabolcsi Z, Dancs G, Molnar J, Papp P, Nagy J, Puskas L, et al. Antler development and coupled osteoporosis in the skeleton of red deer Cervus elaphus: expression dynamics for regulatory and effector genes. Molecular genetics and genomics : MGG. 2010;284(4):273–87.

  57. 57.

    Price J, Allen S. Exploring the mechanisms regulating regeneration of deer antlers. Philosophical transactions of the Royal Society of London Series B, Biological sciences. 2004;359(1445):809–22.

  58. 58.

    Lai AK, Hou WL, Verdon DJ, Nicholson LF, Barling PM. The distribution of the growth factors FGF-2 and VEGF, and their receptors, in growing red deer antler. Tissue & cell. 2007;39(1):35–46.

  59. 59.

    Zhang HL, Cao H, Yang ZQ, Geng S, Wang K, Yu HF, Guo B, Yue ZP. 13cRA regulates the differentiation of antler chondrocytes through targeting Runx3. Cell biology international. 2017;41(3):296–308.

  60. 60.

    Molnar A, Gyurjan I, Korpos E, Borsy A, Steger V, Buzas Z, Kiss I, Zomborszky Z, Papp P, Deak F, et al. Identification of differentially expressed genes in the developing antler of red deer Cervus elaphus. Molecular Genetics And Genomics. 2007;277(3):237–48.

  61. 61.

    Yao B, Zhang M, Liu M, Wang Q, Liu M, Zhao Y. Sox9 Functions as a Master Regulator of Antler Growth by Controlling Multiple Cell Lineages. DNA Cell Biol. 2018;37(1):15–22.

Download references

Acknowledgements

We wish to thank Dr. Peter Fennessy for reading through the paper and giving valuable comments.

Funding

This work was funded by Natural Science Foundation of Jilin Province of China (No. 20170101003JC), the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDA16010403), and Central Public-interest Scientific Institution Basal Research Fund (No. 1610342016003). The funders had no role in study design, data collection, and analysis, decision to publish, or the preparation of the manuscript.

Availability of data and materials

This project has been deposited at NCBI under the BioProject accession PRJNA408029.

Author information

HXB and CYL conceived and designed the experiment. DTW participated in sample collection, RNA extraction, carried out the qRT-PCR. HXB and YDS performed bioinformatics and statistical analyses. HXB and CYL drafted the manuscript. TOY, HXB, DTW and CYL revised the manuscript. All authors read and approved the final manuscript.

Correspondence to Chunyi Li.

Ethics declarations

Ethics approval and consent to participate

This study complied with Chinese Velveting Code (for removal of growing antlers), and was approved by Animal Ethics Committee of Institute of Special Wild Economic Animals and Plants, Chinese Academy of Agricultural Sciences.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Summary of sequencing data of 15 samples across the five tissue layers. (XLSX 11 kb)

Additional file 2:

Figure S1. The flow chart illustrates the major steps involved in the de novo reference transcriptome assembly to filter out the low quality and error assembled transcripts. (JPG 5583 kb)

Additional file 3:

Table S2. Statistics of transcript assembly in the pipeline. (DOCX 16 kb)

Additional file 4:

Figure S2. De novo assembly of the reference transcriptome. A) Length distribution of coding and noncoding transcripts across the length bins. Triangles and dots represent the percentages of coding and noncoding transcripts in the length bins, respectively. B) Expression levels displayed differences of 3–4 orders of magnitude. For each tissue layer, transcripts with FPKM ≥0.5 are all replicates. C) In all cases, the coding transcripts (40–50%) are more abundant in the upper ranks of the distribution (Q4) than the non-coding sequences (20–30%). (PDF 468 kb)

Additional file 5:

Table S3. Summary of 13,203 gene sequences matched with non-redundancy proteins of three closely related species using BLASTX. (DOCX 14 kb)

Additional file 6:

Figure S3. The flow chart shows the steps involved in down-stream bioinformatics analysis, including Principal Component Analysis, Fuzzy c-mean Clustering Analysis and WGCNA. The blue boxes represent executive bioinformatics tools with the corresponding threshold values and the white boxes represent the statistical results from tools. (JPG 5070 kb)

Additional file 7:

Table S4. GO BP enrichment of the DEGs of seven Fuzzy c-means clusters. (XLSX 72 kb)

Additional file 8:

Table S5. GO BP enrichment of the genes of nine WGCNA modules. (XLSX 508 kb)

Additional file 9:

Figure S4. Plot of the minimum centroid distance across a range of cluster number (JPG 3924 kb)

Additional file 10:

Table S6. Primers of the nine hub genes were used for qRT-PCR. (XLSX 9 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Antler
  • Antler growth center
  • Tissue layer
  • Chondrogenesis
  • Transcriptome
  • RNA-seq