Analysis of genetic information from the antlers of Rangifer tarandus (reindeer) at the rapid growth stage

Reindeer is the only deer species in which both males and females regularly grow antlers, providing an excellent model for studying the rapid growth and annual regeneration of antlers. The study of genetic information from reindeer is the basis for revealing the unique mechanism of antler growth. In the present study, we obtained 18.86 GB of clean reads, which were assembled to obtain 94,575 unigenes (average length: 704.69). Among these reads, 30,980 sequences were identified by searching a database of known proteins and then annotated with Gene Ontology (GO) terms, Clusters of Orthologous Groups (COG) classifications and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. All 7,480 simple sequence repeats (SSRs) were detected. A total of 84,435 and 82,226 high-quality single-nucleotide polymorphisms (SNPs) were identified in male and female reindeer, respectively. We identified 31 genes that were highly expressed in reindeer antlers. These genes regulate cell activities that are closely associated with the process of rapid tissue growth. Our results provide a basis for studying reindeer antlers and for further studying the molecular genetics, population genetics, and functional genomics of reindeer.


Introduction
Deer antlers are appendages that undergo annual regeneration [1]. These structures are located on the fore frontal bone of animals of the Cervidae family. In general, deer antlers regenerate annually in spring, when the hard horn is cast, and then grow rapidly in summer and undergo ossification in autumn, after which the antler skin is shed, and the antlers regenerate completely in the following spring. The entire process of antler growth can be divided into the rapid growth stage and the ossification stage [2]. The rapid growth stage occurs approximately 60 days after the antler starts to grow. The extremely high rate of antler growth can reach 2 cm/day, which is one of the fastest rates of organ development in the animal kingdom [3,4].
The growth center responsible for antler growth, located in the antler tip region [5], determines the rapid rate of antler growth and is elegantly regulated without becoming cancerous. This region was referred to be as the proliferation zone [4]. The proliferation zone consists of three layers. From distal to the proximal, these layers are the reserve mesenchyme, precartilage, and cartilage. The reserve mesenchyme is enriched with tightly packed, spindle-shaped, actively dividing cells [6,7]. Rapid antler growth is mainly achieved through the activity of the cells residing in the proliferation zone [8]. Prince et al. found that cells in the reserve mesenchyme are in a higher undifferentiated state than those in other proximal regions [9,10]. These cells will differentiate into chondroblasts toward the proximal region (precartilage and cartilage) [11]. Finally, bone tissues are formed through intramembranous ossification and endochondral ossification [12]. At present, little is known about the molecular mechanisms of antler regeneration and rapid growth. According to previous studies, some factors that are likely responsible for antler growth include insulin-like growth factor (IGF)-I, IGF-II and their receptors [13]; fibroblast growth factor 2 (FGF-2); transforming growth factor beta 1 (TGFβ-1); and bone morphogenetic protein (BMP)-2, BMP-4 and BMP-14 and their receptors [4,14,15]. IGF-I, IGF-II and the IGF receptors were first found to stimulate mesenchymal cell proliferation and were hypothesized to act as the antler growth 'stimulus' [13]. Members of the FGF family and their receptors have also been identified in the primary antler (first antler growth). Fibroblast growth factor 2 has been found to stimulate the proliferation of mesenchymal cells from regenerating antlers. Barling and colleagues identified BMP-2, BMP-4, BMP-14 and their receptors in the primary antler. Faucheux et al. performed immunolocalization of TGF-β in regenerating antlers. Some studies have shown that BMP-2 and TGFβ-1 inhibit the proliferation of mesenchymal cells, but other functionally significant components require further study. A previous study presented several lines of evidence showing that retinoic acid (RA) also plays a role in controlling mesenchymal cell growth and differentiation; sequence information related to antler growth needs to be further examined.
Reindeer are a unique cervid species in which both males and females annually regenerate antlers from the permanent cranial bony protuberance (pedicle). In reindeer antlers, the histologically classified zones are similar to those of other major deer species, including white-tailed deer, mule deer, fallow deer, red deer, American elk, caribou, and sika deer [11]. However, the reindeer antler growth process is poorly understood, and there has been no correlational research on the molecular mechanism underlying the rapid growth of reindeer antlers. Sequencing the transcriptomes of tissues is an effective way to generate useful sequence information to further reveal the molecular mechanisms underlying antler growth and development. In this study, we used RNA-Seq technology to obtain genetic information on reindeer antlers in the rapid growth stage. This investigation on the antlers of reindeer will greatly enhance our understanding of the molecular mechanisms underlying antler growth.

Total RNA preparation
Three samples of antler mesenchyme were pooled into one sample for extraction of total RNA using Trizol reagent (Invitrogen) according to the manufacturer's protocol. RNA concentration was measured using the Qubit TM RNA Assay Kit in a Qubit 1 2.0 fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit and an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

PLOS ONE
Analysis of genetic information from the antlers of Rangifer tarandus (reindeer) at the rapid growth stage

Library preparation and Illumina sequencing
Sequencing libraries were generated using the NEBNext 1 Ultra™ RNA Library Prep Kit for Illumina 1 (NEB, USA) following the manufacturer's recommendations. First, mRNA was purified from total RNA using oligo(dT) magnetic beads and fragmented into small pieces in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H-). Secondstrand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. NEBNext Adaptor with hairpin loop structures were ligated to prepare for hybridization. The library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). PCR was then performed with Phusion high-Fidelity DNA polymerase, universal PCR primers and the Index (X) primer. Finally, the PCR products were purified (AMPure XP system), and library quality was assessed in an Agilent Bioanalyzer 2100 system. The library preparations were sequenced on the Illumina HiSeq platform, and paired-end reads were generated.

Data quality control
In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N sequences and low-quality reads from the raw data. All downstream analyses were based on clean data with high quality.

Transcriptome assembly
Clean short reads were used for the de novo assembly of reindeer unigenes. De novo assembly was conducted using Trinity software. The longest transcript from the potential alternative splicing transcripts was selected as the sample unigene. The unigenes of two samples were combined to create a unigene database.

Functional annotation of unigenes
The general unigene sequences were subjected to BLAST searches against the NCBI nonredundant (Nr) protein database and the SwissProt, Gene Ontology (GO), Clusters of Orthologous Groups (COG), Eukaryotic Orthologous Groups (KOG), EggNOG, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Using KOBAS2.0 software, annotation of unigenes according to KEGG Orthology was performed. After amino sequence prediction, annotation of unigenes was performed using HMMER software against the Pfam database.

Calculation of gene expression levels
The reads obtained from sequencing were compared with the unigene library using Bowtie [16], and expression levels were estimated using RSEM [17] according to the comparison results. Gene expression levels were expressed as FPKM (fragments per kilobase of transcript per million mapped reads) values using the formula FPKM = 10 9 C/NL. In this equation, C is the number of mappable fragments that were uniquely aligned to a unigene; N is the total number of mappable fragments that were uniquely aligned to all unigenes; and L is the length of a unigene in base pairs.

Quantitative real-time PCR analysis
Real-time PCR (qPCR) assays were performed by using the One-Step SYBR method in a 7500 Real-Time PCR system to test the gene expression levels. The expression levels of target genes were normalized against the internal reference gene β-actin. Relative gene expression was calculated using the 2 -ΔΔCt method [18].

Illumina sequencing, assembly, and sequence analysis
cDNA samples were extracted from a male reindeer antler (SRR9051566) and a female reindeer antler (SRR9051567). Using the Illumina sequencing platform, a number of raw reads were generated from the two samples (S1 Table). Clean reads were obtained by removing lowquality reads and adaptor sequences. The amount of total clean data for the assembly of unigenes was 18.86 GB. These clean data were assembled to obtain 124,139 transcripts and 94,575 unigenes. The mean lengths of the transcripts and unigenes were 911.37 nt and 704.69 nt, respectively ( Table 1).

Annotation of predicted proteins
These sequences were first searched using BLASTx against the NCBI nonredundant (Nr), SwissProt, GO, COG, KOG, EggNOG, and KEGG protein databases with a cutoff E value of 10-5. A total of 30,980 unigene annotations were obtained ( Table 2). The Nr database queries revealed that the highest percentages of matched sequences were 18.63% for Bos taurus, 12.83% for Bubalus bubalis, 8.55% for Bos mutus, 8.13% for Pantholops hodgsonii, 6.47% for Capra hircus, 6.14% for Bison bison and 6.08% for Ovis aries. The percentages of sequences matched to other species did not total 35% (Fig 2).

Gene functional annotation
Based on the GO classifications, 14,618 sequences were categorized into 60 functional groups.
We also classified the 7,648 annotated sequences into 25 COG categories. The cluster for "general function prediction only (2,509 members)" was the largest group, followed by "transcription (938 members)", "replication, recombination and repair (909 members)", and "translation, ribosomal structure and biogenesis (721 members)". The categories representing extracellular structures did not contain any members, and nuclear structure (4 members) was the smallest group (S2 Fig). We also annotated 30,980 sequences to the reference pathways in the KEGG database. There were 15,167 sequences mapped to 292 KEGG pathways (S2 Table).

PLOS ONE
Analysis of genetic information from the antlers of Rangifer tarandus (reindeer) at the rapid growth stage

Analysis and statistics of SSRs
Using MISA software, SSR analysis of unigenes with lengths greater than 1 kb was performed, and the results are shown in Table 3. We identified 7,480 SSRs distributed in 5,955 sequences, accounting for 5.91% of the total number of unigenes. Mononucleotide repeats were the most common microsatellite form present (58.70%), followed by trinucleotide repeats (19.30%) and dinucleotide repeats (14.30%). Among the mononucleotide repeats, the (A)n and (G)n types were the most common repeats. The (AC)n motif was the most frequently recovered dinucleotide repeat. The numbers of GGC/GCC/GCG/ CCG/CGC/ CGG-type repeats were greater than those of other trinucleotide repeats, and these motifs are rich in GC base repeats. The percentages of other types of repeats (compound formation, tetranucleotide, pentanucleotide, hexanucleotide) among all SSRs did not reach 1%.

SNP calling analysis
The sequencing reads were named according to the SNP locus allele numbers (Allele). SNP loci can be divided into homozygous-type SNP loci (only one allele) and heterozygous-type SNP loci (two or more alleles) ( Table 4). The GATK2 calls generated a total of 84,435 and

PLOS ONE
Analysis of genetic information from the antlers of Rangifer tarandus (reindeer) at the rapid growth stage 82,226 high-quality SNPs for the male and female reindeer, respectively. All 50,100 SNPs were shared by the male and female reindeer and were classified into 3 types: 4,689 homozygoustype SNP loci in both males and females, 25,110 that were homozygous in males and heterozygous in females, and 20,301 that were heterozygous in males and homozygous in females.

Detection of gene sequences involved in rapid antler growth
By calculating gene expression levels, we characterized a large number of genes that were highly expressed during the rapid growth stage. The results of the real-time PCR (qPCR) assays were consistent with the transcriptome sequencing results (S3 Fig). These highly expressed genes were related to rapid tissue growth and included genes encoding signaling molecules, transcription factors and extracellular matrix (ECM) components (Tables 5-7).

Discussion
Reindeer are the only species in which females regularly grow antlers. The Aoluguya reindeer represents one of the southernmost reindeer species in the world and the only reindeer

PLOS ONE
Analysis of genetic information from the antlers of Rangifer tarandus (reindeer) at the rapid growth stage population in China; this species is limited to a small region in the northeastern part of the Greater Khingan Mountains and numbers less than 1000 individuals [19]. Such small groups are threatened with extinction. It is necessary to carry out research on reindeer in China to protect and utilize reindeer populations. However, little is known about the underlying genetic basis of some traits of reindeer. In 2017, the reindeer genome size was published (2.64 G), and the clustering results showed that reindeer, domestic cattle and sheep share a common ancestor [20]. Our results show that the gene sequences of reindeer share the highest homology with those of B. taurus. There are a large number of SSRs and SNPs distributed in reindeer genes, and the numbers of SNPs distributed in male and female reindeer genes are different. These SSRs and SNPs are expected to serve as molecular markers for studying the genetic diversity and genetic structure of reindeer populations, which will provide guidelines for the conservation of reindeer populations. Our sequence annotation showed that 14,618 sequences were categorized into 60 GO classifications; 7,648 sequences were categorized into 25 COG classifications; and 15,167 sequences were categorized into 292 KEGG pathways. There are 4 pathways related to cellular activities, including pathways associated with cancer (88 members), the PI3K-Akt signaling pathway (427 numbers), focal adhesion (338 members) and the MAPK signaling pathway (322 members). The annotation information provides a reference for further study of the genes that play a role in the rapid growth of reindeer antlers.
Recently, comparative genomic studies have shown that there are gene mutants that are specifically expressed in the reindeer genome [21]. Research shows that tumor suppressor genes and proto-oncogenes are strongly positive for selection in cervids and exhibit strong specific expression in the antlers (e.g., ADAMTS18, FOS, REL, and FAM83A), also find that the fast-growing antlers present a more osteosarcoma-like profile than normal bone tissue [22]. In this study, we observed that a number of genes were highly expressed in the reserve mesenchyme of reindeer antlers, including nine genes encoding signaling proteins ( Table 5). Five of these signaling proteins (CDKN1B, CDKN1C, PTN, GJA1, and MORF4L2) can act as tumor suppressors. The CDKN1B gene exhibited the highest expression level. The major functions of CDKN1B and CDKN1C are to stop or slow down the cell division cycle and control cell cycle progression at G1 [23]. The PTN gene is highly expressed in several tumor cell types [24] involved in tumor angiogenesis and presents mitogenic activity for fibroblasts [25]. GJA1 can regulate cell death, proliferation, and differentiation [26]. MORF4L2 may be required for replicative senescence, apoptosis, and DNA repair. The other four genes (IHH, MK or MDK, ENG, and TGF-βs) perform many cellular functions, including controlling cell growth, cell proliferation, cell differentiation, apoptosis, migration and adhesion [27][28][29]. Five transcription factor genes (ATF4, AP-1, RUNX2, DLX-5, and SOX-9) were identified ( Table 6). ATF4 was the most highly expressed transcription factor, similar to the findings of a study on the antler tips of Chinese sika deer [30]. Wang W et al. found that ATF4 is expressed in growth plate chondrocytes and controls chondrocyte proliferation and differentiation. AP-1 has been shown to control cellular processes, including differentiation, proliferation, and apoptosis [31][32][33]. RUNX2 plays a cell proliferation-related regulatory role in cell cycle entry and exit in osteoblasts. This protein suppresses preosteoblast proliferation by affecting cell cycle progression in the G1 phase and inhibits osteoblastic differentiation via DNA binding [34]. DLX-5 is necessary for osteoblast differentiation as an early BMP-responsive transcriptional activator [35]. The transcription factor SOX-9 is crucial for chondrocyte differentiation [36].
We also observed seventeen ECM protein genes that were highly expressed (Table 7), mainly including TN-N, COLA, Gal-1, WISP-1, OPN, and BSPII. TN-N was the most highly expressed ECM protein gene. Previous research has shown that TN-N is able to induce the migration of tumor cells as well as endothelial cells [37]. Gal-1 is a lectin with a broad range of biological activities [38]. Recently, a study suggested that total absence of Gal-1 in the bone microenvironment allows rapid development of bone tumors [39]. Increased expression of Gal-1 has been correlated with a variety of processes in cancer progression, including cellular aggregation, tumor formation, the metastatic spread of cancer, angiogenesis, and apoptosis [40]. OPN and BSPII belong to a family of glycoproteins that have been linked to cancer metastasis and progression and are overexpressed in a variety of cancers, including lung cancer, breast cancer, colorectal cancer, stomach cancer, ovarian cancer, papillary thyroid carcinoma, melanoma and pleural mesothelioma [41,42]. WISP-1 promotes mesenchymal cell proliferation and osteoblastic differentiation and represses chondrocytic differentiation [43,44]. Among these ECM protein genes, 10 types of COLA genes were highly expressed, namely, type II, III, IV, V, VI, IX, XIII, XV, XVIII, and XXVII genes. In our laboratory, using the differential display reverse transcription PCR method, Zhai Jiancheng et al. previously found that COL6A3 exhibits a high expression level in the mesenchymal layer of male and female reindeer antlers [45]. In the present study, COL6A3 was found to be widely expressed in both male and female antlers. Notably, collagen II, which is a marker of mature cartilage, was shown to be highly expressed in the reindeer antler during the period of rapid antler growth (60 days) in the present study. Rucklidge GJ et al. previously demonstrated that collagen type II is expressed in young antlers (6 weeks) in small quantities and is not expressed in old antlers (5 months) [46]. The collagen II expression pattern in reindeer antlers needs to be explored in further detail. The composition and function of the ECM play an important role in anticancer effects in naked mole rats [47]. The antler mesenchymal layer expresses a large number of ECM genes, the functions of which need to be further studied.

Conclusions
Our results are identical to those of a previous study showing that the gene sequences of reindeer share higher homology with those of B. taurus than with those of other species. The number and distribution of SNPs are different in male and female reindeer. In the mesenchymal layer of the antler tip in male and female reindeer, a large number of highly expressed genes may be related to rapid antler growth. Perhaps these genes are commonly involved in activities such as cell proliferation and differentiation control. Some tumor suppressors are highly expressed during rapid antler growth. The study of the correlation between rapid antler growth and tumor progression may provide a new perspective on the mechanism of antler development.