Identifying deer antler uhrf1 proliferation and s100a10 mineralization genes using comparative RNA-seq

Background Deer antlers are bony structures that re-grow at very high rates, making them an attractive model for studying rapid bone regeneration. Methods To identify the genes that are involved in this fast pace of bone growth, an in vitro RNA-seq model that paralleled the sharp differences in bone growth between deer antlers and humans was established. Subsequently, RNA-seq (> 60 million reads per library) was used to compare transcriptomic profiles. Uniquely expressed deer antler proliferation as well as mineralization genes were identified via a combination of differential gene expression and subtraction analysis. Thereafter, the physiological relevance as well as contributions of these identified genes were determined by immunofluorescence, gene overexpression, and gene knockdown studies. Results Cell characterization studies showed that in vitro-cultured deer antler-derived reserve mesenchyme (RM) cells exhibited high osteogenic capabilities and cell surface markers similar to in vivo counterparts. Under identical culture conditions, deer antler RM cells proliferated faster (8.6–11.7-fold increase in cell numbers) and exhibited increased osteogenic differentiation (17.4-fold increase in calcium mineralization) compared to human mesenchymal stem cells (hMSCs), paralleling in vivo conditions. Comparative RNA-seq identified 40 and 91 previously unknown and uniquely expressed fallow deer (FD) proliferation and mineralization genes, respectively, including uhrf1 and s100a10. Immunofluorescence studies showed that uhrf1 and s100a10 were expressed in regenerating deer antlers while gene overexpression and gene knockdown studies demonstrated the proliferation contributions of uhrf1 and mineralization capabilities of s100a10. Conclusion Using a simple, in vitro comparative RNA-seq approach, novel genes pertinent to fast bony antler regeneration were identified and their proliferative/osteogenic function was verified via gene overexpression, knockdown, and immunostaining. This combinatorial approach may be applicable to discover unique gene contributions between any two organisms for a given phenomenon-of-interest. Electronic supplementary material The online version of this article (10.1186/s13287-018-1027-6) contains supplementary material, which is available to authorized users.


Characterization of FD RM cells and hMSCs for cell proliferation studies. The proliferative capacity of FD RM cells and hMSCs were determined by cell counting studies and cell cycle analysis.
For cell counting studies, cells were seeded into 48-well plates at a density of 0.26 x 10 4 cells/cm 2 overnight (Day 0). Three different media formulation were used -1) DMEM, 10 % FBS, 1 % P/S, 2) Mesenchymal Stem Cell Growth Media and 3) Mesenchymal Stem Cell Growth Media supplemented with 10 ng/mL fibroblast growth factor-2 (FGF-2; Peprotech, Rocky Hill, NJ). Media were changed every 48 h. At 2, 4 and 6 days, cells were counted using an automated cell counter (Beckman Coulter Z2 Particle Counter, Beckman Coulter, Brea, CA). Cell doubling times were calculated using R-studio (R Studio, Boston, MA, http://www.rstudio.com) by visually determining the exponential phase of growth, and plotting the log of the cell counts against time to determine the slope (ln(2)/slope = doubling time) for each sample.
For cell cycle analysis, cells were seeded into T75 tissue culture flasks at a density of 0.27 x 10 4 cells/cm 2 overnight. The following day (Day 0), media were changed to DMEM, 10 % FBS, 1 % P/S. After 4 days, cells were dissociated with 0.05 % trypsin, resuspended in PBS, pelleted by centrifugation (210 g at 4 °C for 5 min for FD RM cells and 500 g at 25 °C for 5 min for hMSCs), fixed in 4 % paraformaldehyde for 10 min, pelleted by centrifugation (210 g at 4 °C for 5 min for FD RM cells and 500 g at 25 °C for 5 min for hMSCs) and stored in PBS at 4 °C until analysis. On the day of analysis, cells were pelleted by centrifugation (210 g at 4 °C for 5 min for FD RM cells and 500 g at 25 °C for 5 min for hMSCs), permeabilized with 0.1 % Triton X-100 in PBS (Sigma Aldrich, St. Louis, MO) for 15 min and stained using propidium iodide/RNAse solution (Cell Signaling Technology, Danvers, MA) for 30 min. Analyses were performed on a BD Aria II flow cytometer and data were analyzed using Flowjo 9.7.5. For each cell cycle distribution, data were fitted to a Watson-Pragmatic model.

Characterization of FD cells and hMSCs for cell differentiation studies.
The mesenchyme lineage commitment of FD cells and hMSCs were determined by several cell differentiation assays. The ability of FD cells to differentiate into adipocytes and chondrocytes was confirmed by Oil Red O staining and Alcian Blue staining, respectively. The ability of FD cells to differentiate into osteoblasts was confirmed by osteogenic gene expression, ALP staining and Alizarin Red S staining whereas the ability of hMSCs to differentiate into osteoblasts was confirmed by Alizarin Red S staining.
For chondrogenic studies, cells were seeded into 24-well plates at a density of 8 x 10 4 cells/5 μL drops to generate micromass cultures. After 2 h (Day 0), media were changed to StemPro Chondrogenic Basal Media (Gibco, Thermo Fisher Scientific, Waltham, MA), 10 % FBS, 1 % P/S (Control media) or StemPro Chondrogenic Differentiation Kit Media (Gibco, Thermo Fisher Scientific, Waltham, MA; Chondrogenic media). Media were changed every 72 h. After 24 days, cells were washed with PBS, fixed with 10 % neutral buffered formalin for 30 min, stained with 1 % Alcian Blue (in 1N HCl; Electron Microscopy Sciences, Hatfield, PA) for 30 min, washed three times with distilled water and air-dried. Brightfield images were acquired using an inverted Zeiss AxioObserver Z1 microscope equipped with an Axiocam ICC 1 color camera.
For osteogenic studies involving ALP staining, cells were seeded into 24-well plates at a density of 1.57 x 10 4 cells/cm 2 overnight. The following day (Day 0), media were changed to DMEM, 10 % FBS, 1 % P/S (Without BMP-2) or DMEM, 10 % FBS, 1 % P/S, 100 ng/mL BMP-2 (With BMP-2). Media were changed every 48 h. After 6 days, cells were fixed for 1 min in 3.7 % formaldehyde. ALP activity (Kit 86C, Sigma Aldrich, St. Louis, MO) was detected according to the manufacturer's instructions. Brightfield images were acquired using an inverted Zeiss AxioObserver Z1 microscope equipped with an Axiocam ICC 1 color camera as well as a Nikon digital SLR camera (Nikon Digital Camera D70, Nikon Corp., Japan). Where necessary, the average pixel intensity was determined using the image histogram tool in Adobe Photoshop (Adobe Systems, San Jose, CA, http://www.adobe.com) as previously described (2,3).

Construction of FD RM cell and hMSC cDNA libraries for RNA-seq. RNA-seq studies were performed for FD RM cells (Isolate 2) and hMSCs (Isolate 24268) under proliferation and mineralization conditions to identify proliferation and mineralization genes, respectively.
For proliferation studies, both FD RM cells and hMSCs were seeded at a density of 1.26 x 10 4 cells/cm 2 overnight and media conditions included 1) DMEM, 0 % FBS, 1 % P/S (0 % serum) and 2) DMEM, 10 % FBS, 1 % P/S (10 % serum). In parallel, another set of cells was seeded at similar densities in 48-well plates to monitor cell growth daily using an automated cell counter (For 6 days). Media were changed every 48 h. After 2.5 days (when cells were observed to be in the exponential phase of growth), RNA was harvested.
For each condition, two replicate RNA samples were isolated. Cells were dissociated with 0.25 % trypsin, the RNA harvested (Qiagen RNeasy Plus Mini kit, Qiagen, Germany) and reverse-transcribed into cDNA (Ovation RNA-seq System V2 kit, NuGEN, San Carlos, CA) for RNA-seq library construction. First, the remainder of the cDNA was sheared (S2 focused-ultrasonicator, Covaris, Woburn, MA). Following this, end repair of the fragmented cDNA, dA-tailing of the end-repaired cDNA, adaptor ligation of dA-tailed cDNA (using custom primers), and PCR enrichment of adaptorligated cDNA were performed for 6 cycles to prepare the cDNA library for RNA-Seq (NEBNext DNA Library Prep Master Mix Set for Illumina, New England Biolabs, Ipswich, MA). Human and FD proliferation samples were sequenced to approximately 87,276,798 reads per library while human and FD mineralization samples were sequenced to approximately 74,904,447 reads per library ( Supplementary Table 1 and 2). RNA-seq and bioinformatics analysis. After cDNA library preparation, samples were sequenced using 100 base-pair, paired-end RNA-seq technology (HiSeq 2000, Illumina, San Diego, CA) and data were analyzed using several bioinformatics software (4,5). The raw sequencing data were concatenated as necessary and the adaptor sequences were removed from the reads to prepare them for analysis. In order to analyze the data, all reads were aligned to the appropriate genome using Spliced Transcripts Alignment to a Reference (STAR; Version 2.3.0, https://code.google.com/p/rna-star/) software (4). Both FD RM proliferation and mineralization reads were aligned to the bovine (Bos taurus) genome (bosTau7 from UCSC Genome Browser, indexed using STAR), as Bos taurus is the closest relative to the fallow deer whose genome has been sequenced (at the time of analysis). Since the human genome is readily available, both hMSC proliferation and mineralization reads were aligned to the human (Homo sapiens) genome (preindexed hg19 released by STAR authors). After converting the STAR output to the appropriate file type and sorting the files using SAMtools (Version 0.1.19 http://www.htslib.org/), the Cufflinks package (Version 2.1.1.1, https://github.com/cole-trapnell-lab/cufflinks) was used to assemble the gene transcripts and to determine differentially-expressed genes between control and treatment conditions (5). Initial attempts at running the Cufflinks package on both FD RM cell and hMSC datasets resulted in errors. These errors originated from Cufflinks' treatment of the soft-clipped regions of aligned reads from the STAR output for some scaffolds (Alexander Dobin, personal communications). When these error-causing soft-clipped regions were removed from both FD and human datasets, the Cufflinks package ran with no errors. The output from Cufflinks was subsequently processed in R-Studio using the cummerbund package (http://compbio.mit.edu/cummeRbund/) to visualize RNA-seq data as well as data quality. Differentially-expressed genes between control and treatment groups were identified based on cut-off values for probability (p ≤ 0.05) and false discovery rate (q ≤ 0.05). Microsoft Excel (Microsoft Corp., Redmond, WA; http://www.microsoftstore.com/), Ingenuity Pathway Analysis (Qiagen, Germany; http://www.ingenuity.com/products/ipa) and Gene Ontology Enrichment Analysis (http://geneontology.org/page/go-enrichment-analysis) were used to compare and analyze differentially-expressed genes in FD RM cell and hMSC datasets. Genes-of-interest for subsequent cloning were identified based on two arbitrary criteria -1) Differentially-expressed genes that exhibited more than 5-fold upregulation in control versus treatment conditions and 2) Genes that were uniquely-expressed in the FD RM cell dataset (i.e. not differentially-expressed in the hMSC dataset).
Cloning of FD genes. Uniquely-expressed, FD genes were PCR-cloned, ligated into a plasmid, transformed into bacteria and purified for subsequent overexpression studies in mammalian cells.
To perform PCR-based cloning, primers were designed to isolate the gene-of-interest from the cDNA of FD RM cells based on the sequences of homologous genes in the bovine genome (National Center for Biotechnology Information, http://www.ncbi.nlm.nih.gov/). These primers included complementary regions 15 -20 base pairs into the 5' and 3' ends of the bovine genes, and contained enzyme restriction sites necessary for eventual insertion into the pVitro2-MCS-Blast plasmid (InvivoGen, San Diego, CA), which contains 2 multiple cloning sites (MCS), a blasticidin resistance gene as well as necessary elements for bacterial and mammalian expression. Each gene was checked to ensure that the selected enzyme restriction sites were not present within the coding sequence of the published mRNA sequence. PCR (Platinum Blue PCR SuperMix, Invitrogen, Thermo Fisher Scientific, Waltham, MA) was performed according to the manufacturer's instructions. PCR products were purified (Wizard SV Gel and PCR Clean-up System, Promega, Sunnyvale, CA) and successful cloning was confirmed via gel electrophoresis and DNA sequencing (ElimBio, Hayward, CA).
To construct the plasmid containing the gene-of-interest, purified PCR products were digested using the appropriate restriction enzymes (New England Biolabs, Ipswich, MA), ligated into pVitro2-MCS-Blast plasmid according to the manufacturer's instructions (NEB Quick Ligation Kit, New England Biolabs, Ipswich, MA). Successful construction of the plasmid containing the gene-of-interest was confirmed by DNA sequencing.
To obtain large quantities of DNA for subsequent overexpression studies, plasmids were transformed into chemically competent DH5α E. coli (One Shot MAX Efficiency DH5α-T1 Competent Cells, Invitrogen, Thermo Fisher Scientific, Waltham, MA) and purified. Bacterial cultures were amplified on LB-agar-blasticidin (Invivogen, San Diego, CA) plates in a 5 % CO2 incubator at 37 °C overnight. Successful colonies were grown in TB-blasticidin liquid media at 37 °C, shaking at 225 -250 rpm overnight, before isolating the amplified plasmid DNA (Purelink Quick Plasmid Miniprep Kit or Purelink HiPure Plasmid Midiprep Kit, Invitrogen, Thermo Fisher Scientific, Waltham, MA). Several plasmids including pVitro2-MCS-Blast (empty plasmid) as well as pVitro2-mRuby2-Blast (6), a plasmid encoding the fluorescent protein mRuby2 (7), were similarly obtained to serve as appropriate transfection controls.

Confirmation of gene overexpression using non-competitive, semi-quantitative PCR.
To confirm gene overexpression, non-competitive, semi-quantitative PCR was employed. RNA was harvested from untransfected cells, cells stably transfected with the empty pVitro2-MCS-Blast plasmid and cells stably transfected with the pVitro2-MCS-Blast plasmids containing the gene-of-interest (Qiagen RNeasy Plus Mini kit, Qiagen, Germany) under appropriate culture conditions. The RNA samples were reverse-transcribed into cDNA according to the manufacturer's instructions. PCR was performed on the cDNA templates for 22 -35 cycles with primers used to originally clone the geneof-interest as well as Bovine gapdh as a reference gene for normalization. PCR-amplified DNA was separated by on a 1.0 % agarose gel (Sigma Aldrich, St. Louis, MO) at 80V for 1 -1.5 h (Bio-Rad Laboratories Inc., Hercules, CA). Images of gel electrophoresis samples were analyzed using the image histogram tool in Adobe Photoshop as previously described (2,3).

Proliferation capability of identified FD gene(s). The function of FD gene(s)-of-interest in cell proliferation was measured in stably-transfected C3H10T1/2 cells using cell counting studies and in FD RM cells using gene knockdown studies.
For cell counting studies, cells were seeded into 48-well plates at a density of 0.26 x 10 4 cells/cm 2 overnight (Day 0) in DMEM, 10 % FBS, 1 % P/S. Media were changed every 48 h. Cells were counted daily using an automated cell counter for 5 days. Cell doubling times were calculated using R-studio (R Studio, Boston, MA, http://www.rstudio.com) by visually determining the exponential phase of growth, and plotting the log of the cell counts against time to determine the slope (ln(2)/slope = doubling time) for each sample.
For gene knockdown studies, cells were seeded into 48-well plates at a density of 0.26 x 10 4 cells/cm 2 overnight in DMEM, 10 % FBS, 1 % P/S. The following day (Day 0), cells were transfected with 30 nM uhrf1 siRNA A and E (Custom-designed based on bovine uhrf1 sequence, Santa Cruz Biotechnology Inc., Dallas, TX) according to the manufacturer's instructions (Polyplus, France) for 72 h. No media change was performed. Confirmation of siRNA-mediated uhrf1 knockdown was determined using non-competitive, semi-quantitative PCR and gel electrophoresis. At 0 and 3 days, cells were counted using an automated cell counter.
For ALP staining, cells were seeded into 24-well plates at a density of 1.57 x 10 4 cells/cm 2 overnight. The following day (Day 0), media were changed to DMEM, 10 % FBS, 1 % P/S (Without BMP-2) or DMEM, 10 % FBS, 1 % P/S, 100 ng/mL BMP-2 (With BMP-2). Media were changed every 48 h. After 4 -12 days, cells were fixed for 1 min in 3.7 % formaldehyde. ALP activity was detected according to the manufacturer's instructions. Brightfield images were acquired using an inverted Zeiss AxioObserver Z1 microscope equipped with an Axiocam ICC 1 color camera. Where necessary, the average pixel intensity was determined using the image histogram tool in Adobe Photoshop as previously described (2,3).
For immunofluorescence staining, cells were fixed in 4 % paraformaldehyde, washed 3 times in PBS, permeabilized with 0.2 % Triton X-100 for 10 min and washed 3 times in PBS. Following this, antibody staining was performed. Cells were incubated in 10 % donkey serum for 20 min followed by incubation with 10 μg/mL rabbit anti-UHRF1 (Sc98704, Santa Cruz Biotechnology Inc., Dallas, TX) or 1 μg/mL mouse anti-S100A10 (Ab89438, Abcam Inc., Cambridge, MA) primary antibody overnight at 4 °C. The following day, cells were washed 3 times in wash buffer (5 min each), incubated in 15 μg/mL donkey anti-rabbit Alexa 647 (711-605-152, Jackson Immunoresearch, West Gove, PA) or 15 μg/mL donkey anti-mouse Alexa 647 (715-605-150, Jackson Immunoresearch, West Gove, PA) secondary antibody for 1 h at 25 °C and washed 5 times in wash buffer (5 min each). Fluorescence images were acquired using an inverted Zeiss AxioObserver Z1 microscope equipped with an X-Cite® Series 120Q metal halide lamp, appropriate filters and an AxioCam MRm camera. Where necessary, the average pixel intensity was determined using the image histogram tool in Adobe Photoshop as previously described (2,3).
Histological and immunofluorescence staining of antler tissue. To examine antler tissue regeneration and ascertain in vivo physiological relevance of in vitro comparative RNA-seq results, deer antler tissue were harvested from an independent fallow deer herd at another local deer ranch (Walking Beam Ranch, Santa Paula, CA) in accordance with the guidelines established by Stanford University's Administrative Panel on Laboratory Animal Care and subjected to histological and immunofluorescence staining. These deer were approximately 2 -3 years old and antler tissues were harvested during early stages of antler regeneration (4 -7 inches in height) as described previously. Tissue samples were fixed in 10 % formalin and stored in 70 % ethanol. Tissue samples were subjected to histological processing via a graded ethanol dehydration series at 4 °C (70 % ethanol overnight, 85 % ethanol overnight, 95 % ethanol overnight and 100 % ethanol overnight) followed by xylene infiltration at 25 °C (50 % xylene in ethanol for 1 h, 50 % xylene in ethanol overnight, 100 % xylene for 1 h and 100 % xylene overnight) and then paraffin infiltration at 60 °C (50 % paraffin in xylene for 2 h, 100 % paraffin for 2 h, 3 washes of 100 % paraffin for 20 min each, 100 % paraffin overnight and 100 % paraffin for 20 min). Subsequently, tissue samples were embedded in paraffin blocks and sectioned at 4 -8 μm intervals using a Leica rotary microtome (RM 2255, Leica Biosystems Inc., Buffalo Grove, IL). Prior to staining, tissue sections were deparaffinized (3 washes of 100 % xylene for 3 min each) and rehydrated (50 % ethanol in xylene for 3 min, 2 washes of 100 % ethanol for 3 min each, 2 washes of 95 % ethanol for 3 min each, 2 washes of 70 % ethanol for 3 min each and PBS for 3 min).
For histological staining, tissue sections were stained with Alcian Blue or Alizarin Red S. To stain for cartilage, tissue samples were incubated with 1 % Alcian Blue (in 1N HCl) for 30 min, washed three times with distilled water, counter-stained with Neutral Red for 10 min, washed three times with distilled water, dehydrated in ethanol and mounted. Brightfield images were acquired using an inverted Zeiss AxioObserver Z1 microscope equipped with an Axiocam ICC 1 color camera. To stain for mineralized bone, thick (1 -3 mm) tissue samples were incubated with 2 % Alizarin Red S for 30 min and washed five times with distilled water. Thick tissue samples were imaged using a Nikon digital SLR camera.
Statistical analysis. Statistical analyses involving RNA-seq were performed by Cufflinks and R-Studio (5). Statistical significance for differentially-expressed genes was established at p ≤ 0.05 and q ≤ 0.05. Statistical analyses not involving RNA-seq were performed using IBM SPSS Statistics for Windows 22 (IBM Corp., North Castle, NY, http://www.ibm.com). These experiments were performed with at least 3 replicates per condition. Sample sizes were estimated to detect a group mean difference of 50 % ± 1 to 2 standard deviations with a power (1 -β) of 0.8 and α = 0.05 (http://powerandsamplesize.com/Calculators/Compare-k-Means/1-Way-ANOVA-Pairwise). Quantitative data was presented as means ± standard error of mean (mean ± SEM) where appropriate. Relative fold changes for PCR data were log transformed in order to make the data distribution more symmetrical since gene expression data are often log normally distributed (8). To determine whether data were normally-distributed and whether there was equality of variances among groups, p values were computed via the Shapiro-Wilk test and the Levene test, respectively. For two mean comparisons, p values were computed via the t-test. If there was equal variance between groups, p values were calculated using pooled variance. Otherwise, p values were calculated using separate variance. For more than two mean comparisons, p values were computed via Analysis of Variance (ANOVA). If majority of data (two-thirds or more) were normally-distributed or there was equal variance among groups, p values were calculated using ANOVA followed by Tukey's Honest Significant Difference post-hoc multiple comparison test. This approach was based on the robustness of ANOVA under conditions of non-normality and heterogeneity of variance (9, 10). Otherwise, p values were calculated using Welch's ANOVA followed by Games-Howell post-hoc multiple comparison test. This approach enables improved control of Type I errors and greater power under conditions of non-normality and heterogeneity of variance (11). Statistical significance was established at p ≤ 0.05.

Supplementary Figures S1-S7
Supplementary Figure S1. Overview of in vitro comparative RNA-seq. The colored lines indicate sites where deer (FP, facial periosteum; PP, pedicle periosteum; RM, reserve mesenchyme) and human (hMSC, human mesenchymal stem cells) skeletal progenitor cells were harvested. Isolated cells were subjected to RNA-seq under proliferation and mineralization conditions independently and the resulting datasets were compared to identify genes that were highly-expressed (> 5-fold increase) and unique to deer. Subsequently, the physiological relevance and contribution of these genes to antler regeneration were determined.  Table S1. Data quality of FD RM cell (Isolate 2) RNA-Seq samples.