Cross-Species Analysis Reveals Co-Expressed Genes Regulating Antler Development in Cervidae

Antlers constitute an interesting model for basic research in regenerative biology. Despite decades of being studied, much is still unknown about the genes related to antler development. Here, we utilized both the genome and antlerogenic periosteum (AP) transcriptome data of four deer species to reveal antler-related genes through cross-species comparative analysis. The results showed that the global gene expression pattern matches the status of antler phenotypes, supporting the fact that the genes expressed in the AP may be related to antler phenotypes. The upregulated genes of the AP in three-antlered deer showed evidence of co-expression, and their protein sequences were highly conserved. These genes were growth related and likely participated in antler development. In contrast, the upregulated genes in antler-less deer (Chinese water deer) were involved mainly in organismal death and growth failure, possibly related to the loss of antlers during evolution. Overall, this study demonstrates that the co-expressed genes in antlered deer may regulate antler development.

Introduction

An ability to regenerate complex structures is widespread among lower organisms and is retained in some vertebrate species such as urodele amphibians. However, adult mammalian examples of epimorphic regeneration are extremely rare with the most dramatic example being the annual replacement of antlers in deer (Goss, 1983; Bubenik et al., 1990). Therefore, antlers constitute an interesting model for basic research in regenerative biology. On the other hand, despite decades of being studied, much is still unknown about the molecular regulation of antler development. Although being called cranial appendages, deer antlers do not grow directly from the head; instead, they generate and regenerate from the fully grown pedicles (antecedents of antlers). Pedicles and first antlers are both derived exclusively from the periosteum overlying the frontal crest of the deer head, known as the antlerogenic periosteum (AP) (Hartwig and Schrudde, 1974; Goss and Powel, 1985; Li and Suttie, 2000). Removal of the AP prior to pedicle initiation stops the pedicle and antler formation, and transplantation of the AP autologously to other sites of the body causes ectopic pedicle and antler to grow (Goss and Powel, 1985; Li et al., 2009a), thus enabling investigation at the molecular level. Each year during spring, fully calcified antlers are cast from the pedicles, which trigger the initiation of new antler regeneration from the pedicle stump (Li et al., 2002). In late spring and summer, antlers enter into the most rapid growth period (up to 2 cm/day) and are covered by a special type of skin, called velvet skin. In autumn, antlers are intensively calcified and shed the velvet skin to expose the bony antlers.

Antlers exhibit polymorphism/polyphenism ranging from the large and complex structures grown by large species (Figure 1A) such as reindeer (Rangifer tarandus, antler length: ∼84 cm and multiple times) and sika deer (Cervus nippon, ∼78 cm and 3–5 times), to the small, simple antlers grown by small species such as the muntjac (Reeves’s muntjac, ∼5.4 cm and 1-2 times) and to the antler-less small deer species such as the Chinese water deer (Hydropotes inermis). Given the range of phenotypes across the deer species and their common evolutionary origin from an antlered ancestor (Goss, 1983; Bubenik et al., 1990; Rössner et al., 2020), here, we selected four deer species (Cervinae subfamily: sika deer and muntjac; Capreolinae subfamily: reindeer and Chinese water deer) that span all two deer subfamilies in Cervidae (Heckeberg, 2020), namely, three-antlered deer (ATD) and the antler-less. We utilized both genome and AP transcriptome data (or the presumptive AP tissue in the case of Chinese water deer) of these four species to reveal antler-related genes through cross-species comparative analysis.

FIGURE 1

FIGURE 1. (A) Antlers of three-antlered deer (ATD) species including reindeer (RD), sika deer (SD), and muntjac (MJ); and one antler-less deer species, Chinese water deer (CWD). (B) Principal component analysis (PCA) of the gene expression of 12 samples across four deer species (three samples/species). (C) Left: hierarchical cluster analysis based on the gene expression matrix of 12 samples among four deer species; Right: evolutionary topology of four deer species obtained from the previous study (Chen et al., 2019).

Materials and Methods

Sample Collection, RNA Preparation, and Sequencing

All AP tissues of the four deer species were collected from male deer and approved by the Animal Ethics Committee of the Institute of Antler Science and Product Technology, Changchun Sci-Tech University (CKARI202002). Sika deer and muntjac will start to initiate pedicle growth from the AP around the 6th month onward (approaching puberty), so AP was sampled from these two species at the age of 6 months. Chinese water deer do not grow antlers; thus, their AP-equivalent tissue can be sampled at any time. We collected AP-equivalent tissue from Chinese water deer at the age of 6 months in order to make the age of these species be consistent. Reindeer grow their pedicles at about 2 weeks of their age; thus, we sampled their AP tissues at the age of 2 weeks. A total of 12 pieces of AP tissues were collected (three biological replicates/species; four species), and these fresh tissues were immediately frozen in liquid nitrogen and then stored at −80°C for RNA extraction. Around 100 mg/tissue sample was rapidly ground into a fine powder using liquid nitrogen and a Freezer/Mill 6770 (SPEX CertiPrep Ltd., United States). Total RNA from each sample was isolated using a TRIzol reagent (Invitrogen Inc., Camarillo, CA) according to the manufacturer’s procedure. RNA quality was confirmed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., United States), 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). The library quality was assessed on the Agilent Bioanalyzer 2100 system, and molecular fragments with the length of 200–300 bp were selected and sequenced with 150 bp paired-end reads using an Illumina HiSeq 4000 instrument at Beijing Genomics Institute (Shenzhen, China). Information regarding the clean reads from the 12 samples that passed the quality filtering is summarized in Supplementary Table S1.

Retrieval of the Currently Available Published Deer Genome Sequences

Genome sequences of four deer species were retrieved and downloaded from the published genome databases, including sika deer (CNGB: GWHANOY00000000), reindeer (GigaDB: DOI:10.5524/100370) (Li et al., 2017), muntjac (GenBank: GCA_008787405.2) (Mudd et al., 2020), and Chinese water deer (GenBank: GCA_006459105.1) (Chen et al., 2019). These genome sequences were of high-quality, such as both sika deer and muntjac have been assembled and annotated at the chromosome level.

Obtaining Deer-Orthologs

In order to obtain high-quality deer-orthologs among four deer species, we developed a strict pipeline (Supplementary Figure S1) that combined their genome and transcriptome data to obtain protein-coding deer-orthologs based on the previous published methods (Fushan et al., 2015; Wang et al., 2016; Darbellay and Necsulea, 2020). First, we used Trinity v2.4.1 (Haas et al., 2013) to carry out de novo assembly from the quality-filtered reads for four deer species, respectively. The high-quality reads were also mapped on their corresponding reference genomes to generate reference-based transcripts by using HISAT2/StringTie (Pertea et al., 2016). Then, both the de novo and reference-based transcripts were merged using a CD-HIT-EST tool in CD-HIT Suite v3.0.3 (Li et al., 2006). Subsequently, the protein-coding sequences of the four species were predicted by TransDecoder v2.0.1 (Haas et al., 2013), respectively, and called transcriptome-based protein-coding datasets (Tset). We also derived the protein-coding sequences of the four species from their corresponding genomes based on gene annotation file, called genome-based protein-coding datasets (Gset). For the two datasets, we used cattle as the reference species to determine protein-coding orthologs, respectively, based on a best-bidirectional BLAST hit criterion (Overbeek et al., 1999). After joining these orthologs from both the Tset and Gset, we used InParanoid v4.1 (Ostlund et al., 2010) to obtain the deer-orthologs based on cattle orthologs (CT) as an outgroup. Finally, we obtained a total of 11,006 deer-orthologs among these four species.

Principal Component Analysis and Hierarchical Cluster Analysis

To obtain the gene expression matrix (count and FPKM) of the 12 samples, the high-quality reads were mapped to their corresponding 11,006 deer-orthologs, respectively, with RSEM v1.3.0 (Li and Dewey, 2011). Both PCA and hierarchical cluster analyses were performed based on these 11,006 deer-orthologs by “prcomp” and “hclust” R function, respectively. Among these deer-orthologs, 9,540 (86.7%) had FPKM > 0.5 in all biological replicates of any one of the four species.

Generation of Differentially Expressed Genes

Considering the different gene lengths between these deer-orthologs of the four species, the SCBN v1.10.0R package (Zhou et al., 2019) was used to search for the optimal scaling factors from the count matrix of the 11,006 deer-orthologs. The derived scaling factors were manually input to the DESeq2 v 2.1.18R package (Anders and Huber, 2010) to produce the DEGs. All p values were adjusted for multiple testing using the Benjamini–Hochberg method as implemented in DESeq2.

Protein Sequence Divergence Analysis of Deer-Orthologs

The orthologous alignments between the ATDs and Chinese water deer were further conducted by using GBLOCKS (Castresana, 2000) to remove gaps and unreliable alignment columns. We then calculated the non-synonymous mutation rate (Ka) and the synonymous mutation rate (Ks) of these alignments by using the KaKs-calculator (Zhang et al., 2006) with the “MA” method. We excluded genes with Ks > 2 and Ka/Ks > 3 because high estimates of Ks may indicate saturation in synonymous sites or alignment errors (Uebbing et al., 2016).

Construction of the Protein–Protein Interaction Network

The online database, STRINGdb (https://string-db.org/), was used to construct the protein–protein interaction (PPI) network, with all interaction sources and a minimum required interaction score being set at ≥ 0.4 for our genes. The Cytoscape v3.6 (Shannon et al., 2003) was used to visualize the protein–protein network. Network statistics were performed through in-house commands in the Cytoscape. Key hub nodes in the network were defined by their connective degrees with other nodes.

Ingenuity Pathway Analysis

Functional enrichment of the data set was carried out using the IPA package (release date: 8 Feb 2019). Gene expression as defined by average fold change of all comparison groups in each gene set that we defined was used as the core analysis type with the Ingenuity Knowledge Base as the reference set. When analyses were performed, an adjust p value with the Benjamini–Hochberg method was set at < 0.01 and z-score (absolute value) > 2.

Statistical Analysis

Correlation coefficients and the Wilcox test were performed using the “cor.test” and “wilcox.test” R functions, respectively. Significant differences (at least p value < 0.05) between two groups were determined.

Results and Discussion

Global Gene Expression Patterns of the AP Perfectly Match the Status of Antler Phenotype Variation

We analyzed gene expression data of the AP for 11,006 deer-orthologs to ascertain whether the global gene expression pattern matches the status of antler phenotypes across the four Cervidae species. Much of the variation in gene expression was evident in the comparison between ATDs and Chinese water deer, which separated on the first PCA axis (36.2% of the variation) (Figure 1B). The second PCA axis (26.7%) reflects the differences among the ATDs. This variation was also confirmed by a hierarchical cluster analysis among the four species (Figure 1C). The two large ATDs were closest to one another, followed by the small ATD (muntjac). The correlation of the ATDs and Chinese water deer was the lowest. These results showed that it is the global gene expression pattern in the AP tissues rather than the species evolutionary topology (Chen et al., 2019) that matches the status of antler phenotypes (Figure 1D), supporting the fact that the genes expressed in the AP may be related to antler phenotypes.

Upregulated Genes in the ATDs Show Evidence of Co-Expression and Co-Variation

First, we analyzed pair-wise correlation coefficients between the DEG sets (|log₂FoldChange| ≥ 1 and adjust p value < 0.01; Supplementary Table S2). The results showed that the correlation coefficients (0.86–0.88) between three ATDs relative to Chinese water deer were higher than those of other groups (Figure 2), suggesting that these DEGs were likely involved in regulation of antler development. However, these high correlations were found to be mainly contributed by the upregulated genes (Figure 2, red dotted boxes).

FIGURE 2

FIGURE 2. Correlation of the DEG sets (|log₂FoldChange| ≥ 1 and adjusted p value < 0.01) of the ATDs relative to Chinese water deer: (A) SD/CWD vs RD/CWD, (B) SD/CWD vs MJ/CWD, and (C) RD/CWD vs MJ/CWD. Note that the high correlations were due mainly to the upregulated genes (red-dotted boxes). RD: reindeer, SD: sika deer, MJ: muntjac, and CWD: Chinese water deer.

Next, we asked whether these DEGs related to antler development among the ATDs would show evidence of co-expression. We divided the DEGs (|log₂FoldChange| ≥ 1 and adjusted p value < 0.01) in each of the ATDs vs. Chinese water deer, respectively, into six gene sets based on the overlap level of these DEGs by applying Venn analysis (Figure 3A). Of these six gene sets, three were upregulated (uSet1 (400), uSet2 (1,001), and uSet3 (1,699)) and three downregulated (dSet1 (372), dSet2 (523), and dSet3 (644); Supplementary Table S3). Genes in Set1 were up/downregulated in all three comparative groups, those in Set2 were up/downregulated in any two of the three groups, and those in Set3 were up/downregulated in only one of the three groups. Expression changes of the upregulated genes in the uSet1, 2, and 3 of each comparative group decreased successively and significantly (Wilcox test, p value < 0.01; Figure 3B). In contrast, expression patterns in the downregulated genes did not show a similar trend. We further investigated divergence of their protein sequences among the DEGs in the six gene sets. We applied the Ka value (non-synonymous mutation rate) of protein orthologs to measure the degree of protein divergence. The results showed that the degree of protein divergence increased successively and significantly from uSet1 to uSet2 to uSet3 (Wilcox test, p value < 0.01) (Figure 3C). In contrast, there was no such trend detected for the downregulated genes. Taken together, these findings indicated the higher the co-expression degree of upregulated genes among the ATDs, the more conserved their protein sequences.

FIGURE 3

FIGURE 3. (A) Venn diagram of the upregulated (red) and downregulated (green) genes of three comparative groups respectively based on the overlap level. Note that the higher the overlap level, the higher co-expressed degree of the genes. (B) Boxplot of the changes in gene expression in the upregulated (red) and downregulated (green) gene sets of each of three comparative groups, respectively. Wilcox test **p < 0.01 and *p < 0.05. (C) Boxplot of the degree of protein divergence (non-synonymous mutation rate (Ka)) of deer-orthologs in the upregulated (red) and downregulated (green) gene sets of each of three comparative groups, respectively. Wilcox test **p < 0.01 and *p < 0.05.

Upregulated Genes in the ATDs Are Growth Related and Likely Participate in Antler Development

As these upregulated genes (400) in the uSet1 were the most co-expressed among the ATDs, these genes were used to construct the protein–protein interaction network using STRINGdb, and visualized by using Cytoscape. The network analysis results showed that the number of nodes was 346 (86.5%) and average degree of connectivity was 7.4 (Figure 4A), indicating this interaction network is robust. Furthermore, both the degree of connectivity and log₂FoldChange of these genes was correlated (p = 0.035; Figure 4B), we focused on the role of key hub genes with the degree of connectivity ≥ 15 in the network, for example, CTNNB1 (61), MAPK3 (41), JUN (38), ITGB1 (21), THBS1 (20), and BCL2L1 (19). The CTNNB1 is a key downstream component of the canonical Wnt signaling pathway, which plays an important role in appendage regeneration (Stoick-Cooper et al., 2007; Ba et al., 2019). The protein THBS1 is reported to promote angiogenesis through interactions with a number of integrin heterodimers (e.g., ITGA3-ITGB1 complex) (Chandrasekaran et al., 2000). This is consistent with the in vivo findings that blood vessels are richly distributed in antler lineage tissues, even including antler cartilage, which is avascular in its somatic counterpart. The anti-apoptotic protein BCL2L1, upregulated by the MAPK/c-Jun signaling pathway (Zhang et al., 2018), is a potent inhibitor of cell death by blocking the voltage-dependent anion channel (VDAC) to prevent the release of death proteases called caspase (CASP) activator, cytochrome c (CYC) (Tsujimoto and Shimizu, 2000). The anti-apoptotic factors could be prerequisite for the formation of an antler tissue mass of 20 kg or so within 60 days from around 3 million cells (Li et al., 2009b).

FIGURE 4

FIGURE 4. (A) Interactive network. This network consisted of 346 (86.5%) genes in the uSet. Average degree of connectivity was 7.4. Nodes with the degree of connectivity ≥ 15 are emphasized. Value of log₂FoldChange of genes is indicated by the grade of the color. The size of the node indicates the degree of connectivity; the larger the node, the higher the degree of connectivity. (B) Correlation of 346 (86.5%) genes in the uSet1 between their degree and log₂FoldChange. Gene names with the degree of connectivity > 15 are shown (red dot). Bubble plot of enrichment analysis results of (C) molecular/cellular functions and (D) canonical pathways in IPA software using all genes in the uSet1 and uSet2 and dSet1 and dSet2. The key upregulated (red) and downregulated (green) annotated terms are highlighted.

We further utilized a total of 2,296 co-expressed genes in the uSet1, uSet2, dSet1, and dSet2 to perform an enrichment analysis of molecular/cellular functions and canonical pathways by applying IPA software. The analysis of molecular/cellular functions showed that of the enriched genes, those related to growth and development of the embryo/organism/body and differentiation of bone cells were highly upregulated (z-score > 2), while those related to organismal death and growth failure (morbidity and mortality) were highly downregulated (z-score < -2) in the ATDs as compared with the Chinese water deer, which are possibly related to the disappearance of antlers in the Chinese water deer during the course of evolution (Randi et al., 1998; Wang et al., 2019) (Figure 4C and Supplementary Table S4). These findings indicate that these co-expressed genes are involved in antler development and rescue antler failure.

The enrichment analysis of canonical pathways showed that the most upregulated signaling pathways were growth related (e.g., synaptogenesis, estrogen receptor, ERK/MAPK, BMP, and PI3K/AKT, z-score > 2) and the most downregulated signaling pathways included mitochondrial dysfunction and oxidative phosphorylation (z-score < -2; Figure 4D and Supplementary Table S5). Mitochondrial dysfunctions (e.g., CASP9, CYC1, VDAC1, and VDAC2) are a group of genetic disorders that are characterized by defects in oxidative phosphorylation (Gorman et al., 2016), which may be related to the disappearance of antlers in the Chinese water deer.

Among the upregulated pathways, synaptogenesis signaling was the most significantly enriched, which further supports the neural crest origin of the AP cells (Price et al., 2005; Landete-Castillejos et al., 2019). Antlers are organs of bone and directly formed from the proliferation and differentiation of AP cells (Goss and Powel, 1985; Li and Suttie, 2000). The BMP pathway is required for chondrogenesis/osteogenesis and thus is closely associated with the antler development. In the BMP pathway, genes for BMP1/2/5, BMPR2, and CREB1 were upregulated (Supplementary Table S3). BMPs are involved in endochondral bone formation and embryogenesis (Sasano et al., 1993). These proteins transduce their signals through the formation of heteromeric complexes of BMPR1 and 2. It was reported that estrogen receptor mediates cell proliferation through the cAMP/PKA/CREB1 axis in murine bone marrow mesenchymal stem cells (Chuang et al., 2020).

The activated ERK/MAPK and PI3K/AKT pathways in the present study were also identified in our previous studies (Li et al., 2012; Li et al., 2018). These two pathways are important for maintaining rapid cell proliferation. In the PI3K/AKT signaling pathway, genes for BCL2L1, CCND1, CTNNB1, and RAF1 were upregulated (Supplementary Table S3). The gene CCND1 (encoding activity of cyclin-dependent kinase), a cell cycle-dependent factor, involved in cell cycle progression, may promote proliferation of the AP cells. Interestingly, it is reported that a newly identified binding motif of the androgen receptor evolved upstream of the CCND1 gene, and may result in female reindeer antler growth (Lin et al., 2019). The CDKN1A protein inhibits the activity of CCND1; it was downregulated in the present study (Supplementary Table S3), which thus may increase cellular proliferation. Given that a limited number of AP cells (around 3.3 million cells) can generate 10 kg or more of antler tissue mass within 60 days (Li and Suttie, 2001), CCND1 and CDKN1A, two highly expressed DEGs in the AP of ATDs, would be at least one of the factors contributing to this phenomenal growth rate.

Overall, this is the first study on the genes involved in antler development through cross-species comparative analysis based on their genome and AP transcriptome data in the Cervidae, and it further demonstrates that the co-expressed genes in antlered deer may regulate antler development.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/, PRJNA768490.

Ethics Statement

All AP tissues of the four deer species were collected from male deer and approved by the Animal Ethics Committee of the Institute of Antler Science and Product Technology, Changchun Sci-Tech University (CKARI202002).

Author Contributions

CL and HB conceived and designed the research. HB and MC collected samples. HB analyzed the data. CL and HB wrote and revised the manuscript. All authors have read and approved the final manuscript.

Funding

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA16030305), National Natural Science Foundation of China (U20A20403), Jilin Province Education Department Support Program (No. JJKH20221324KJ), and Changchun Science and Technology Development Funds (No. 21ZY51).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank Rudiger Brauning of Agresearch New Zealand for allowing them to access IPA software and for partly helping with data analysis. They wish to thank Peter Fennessy of AbacusBio Limited, Dunedin, New Zealand, for reading through the manuscript and giving valuable comments.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.878078/full#supplementary-material

References

Anders, S., and Huber, W. (2010). Differential Expression Analysis for Sequence Count Data. Genome Biol. 11, R106. doi:10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

Ba, H., Wang, D., Yau, T. O., Shang, Y., and Li, C. (2019). Transcriptomic Analysis of Different Tissue Layers in Antler Growth Center in Sika Deer (Cervus Nippon). BMC genomics 20, 173. doi:10.1186/s12864-019-5560-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Bubenik, A. B. (1990). “Epigenetical, Morphological, Physiological, and Behavioral Aspects of Evolution of Horns, Pronghorns, and Antlers,” in Horns, Pronghorns, and Antlers: Evolution, Morphology, Physiology, and Social Significance. Editors G. A. Bubenik,, and A. B. Bubenik (New York, NY: Springer New York), 3–113. doi:10.1007/978-1-4613-8966-8_1

CrossRef Full Text | Google Scholar

Castresana, J. (2000). Selection of Conserved Blocks from Multiple Alignments for Their Use in Phylogenetic Analysis. Mol. Biol. Evol. 17, 540–552. doi:10.1093/oxfordjournals.molbev.a026334

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandrasekaran, L., He, C.-Z., Al-Barazi, H., Krutzsch, H. C., Iruela-Arispe, M. L., and Roberts, D. D. (2000). Cell Contact-dependent Activation of α3β1 Integrin Modulates Endothelial Cell Responses to Thrombospondin-1. MBoC 11, 2885–2900. doi:10.1091/mbc.11.9.2885

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Qiu, Q., Jiang, Y., Wang, K., Lin, Z., Li, Z., et al. (2019). Large-scale Ruminant Genome Sequencing Provides Insights into Their Evolution and Distinct Traits. Science 364, eaav6202. New York, N.Y. doi:10.1126/science.aav6202

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuang, S. C., Chen, C. H., Chou, Y. S., Ho, M. L., and Chang, J. K. (2020). G Protein-Coupled Estrogen Receptor Mediates Cell Proliferation through the cAMP/PKA/CREB Pathway in Murine Bone Marrow Mesenchymal Stem Cells. Int. J. Mol. Sci. 21, 6490. doi:10.3390/ijms21186490

PubMed Abstract | CrossRef Full Text | Google Scholar

Darbellay, F., and Necsulea, A. (2020). Comparative Transcriptomics Analyses across Species, Organs, and Developmental Stages Reveal Functionally Constrained lncRNAs. Mol. Biol. Evol. 37, 240–259. doi:10.1093/molbev/msz212

PubMed Abstract | CrossRef Full Text | Google Scholar

Fushan, A. A., Turanov, A. A., Lee, S. G., Kim, E. B., Lobanov, A. V., Yim, S. H., et al. (2015). Gene Expression Defines Natural Changes in Mammalian Lifespan. Aging Cell 14, 352–365. doi:10.1111/acel.12283

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorman, G. S., Chinnery, P. F., DiMauro, S., Hirano, M., Koga, Y., McFarland, R., et al. (2016). Mitochondrial Diseases. Nat. Rev. Dis. Primers 2, 16080. doi:10.1038/nrdp.2016.80

PubMed Abstract | CrossRef Full Text | Google Scholar

Goss, R. J. (1983). Deer Antlers: Regeneration, Function, and Evolution. New York, NY: Academic Press.

Google Scholar

Goss, R. J., and Powel, R. S. (1985). Induction of Deer antlers by Transplanted Periosteum I. Graft Size and Shape. J. Exp. Zool. 235, 359–373. doi:10.1002/jez.1402350307

PubMed Abstract | CrossRef Full Text | Google Scholar

Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2013). De Novo transcript Sequence Reconstruction from RNA-Seq Using the Trinity Platform for Reference Generation and Analysis. Nat. Protoc. 8, 1494–1512. doi:10.1038/nprot.2013.084

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartwig, H., and Schrudde, J. (1974). Experimentelle Untersuchungen zur Bildung der primären Stirnauswüchse beim Reh (Capreolus capreolus L.). Z. für Jagdwissenschaft 20, 1–13. doi:10.1007/bf01901843

CrossRef Full Text | Google Scholar

Heckeberg, N. S. (2020). The Systematics of the Cervidae: a Total Evidence Approach. PeerJ 8, e8114. doi:10.7717/peerj.8114

PubMed Abstract | CrossRef Full Text | Google Scholar

Landete-Castillejos, T., Kierdorf, H., Gomez, S., Luna, S., García, A. J., Cappelli, J., et al. (2019). Antlers - Evolution, Development, Structure, Composition, and Biomechanics of an Outstanding Type of Bone. Bone 128, 115046. doi:10.1016/j.bone.2019.115046

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: Accurate Transcript Quantification from RNA-Seq Data with or without a Reference Genome. BMC bioinformatics 12, 323. doi:10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Clark, D. E., Lord, E. A., Stanton, J.-A. L., and Suttie, J. M. (2002). Sampling Technique to Discriminate the Different Tissue Layers of Growing Antler Tips for Gene Discovery. Anat. Rec. 268, 125–130. doi:10.1002/ar.10120

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Gao, X., Yang, F., Martin, S. K., Haines, S. R., Deng, X., et al. (2009). Development of a Nude Mouse Model for the Study of Antlerogenesis-Mechanism of Tissue Interactions and Ossification Pathway. J. Exp. Zool. 312, 118–135. doi:10.1002/jez.b.21252

CrossRef Full Text | Google Scholar

Li, C., Harper, A., Puddick, J., Wang, W., and McMahon, C. (2012). Proteomes and Signalling Pathways of Antler Stem Cells. PLoS One 7, e30026. doi:10.1371/journal.pone.0030026

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., and Suttie, J. M. (2001). Deer Antlerogenic Periosteum: a Piece of Postnatally Retained Embryonic Tissue? Anat. Embryol. (Berl) 204, 375–388. doi:10.1007/s004290100204

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., and Suttie, J. M. (2000). Histological Studies of Pedicle Skin Formation and its Transformation to Antler Velvet in Red Deer (Cervus elaphus). Anat. Rec. 260, 62–71. doi:10.1002/1097-0185(20000901)260:1<62:aid-ar70>3.0.co;2-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Yang, F., and Sheppard, A. (2009). Adult Stem Cells and Mammalian Epimorphic Regeneration-Insights from Studying Annual Renewal of Deer antlers. Cscr 4, 237–251. doi:10.2174/157488809789057446

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Zhao, H., Wang, D., McMahon, C., and Li, C. (2018). Differential Effects of the PI3K AKT Pathway on Antler Stem Cells for Generation and Regeneration of antlers I In Vitro I. Front. Biosci. 23, 1848–1863. doi:10.2741/4676

CrossRef Full Text | Google Scholar

Li, W., and Godzik, A. (2006). Cd-hit: a Fast Program for Clustering and Comparing Large Sets of Protein or Nucleotide Sequences. Bioinformatics 22, 1658–1659. doi:10.1093/bioinformatics/btl158

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Lin, Z., Ba, H., Chen, L., Yang, Y., Wang, K., et al. (2017). Draft Genome of the Reindeer (Rangifer tarandus). GigaScience 6, 1–5. doi:10.1093/gigascience/gix102

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z., Chen, L., Chen, X., Zhong, Y., Yang, Y., Xia, W., et al. (2019). Biological Adaptations in the Arctic Cervid, the Reindeer (Rangifer tarandus). Science 364, eaav6312. New York, N.Y. doi:10.1126/science.aav6312

PubMed Abstract | CrossRef Full Text | Google Scholar

Mudd, A. B., Bredeson, J. V., Baum, R., Hockemeyer, D., and Rokhsar, D. S. (2020). Analysis of Muntjac Deer Genome and Chromatin Architecture Reveals Rapid Karyotype Evolution. Commun. Biol. 3, 480. doi:10.1038/s42003-020-1096-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostlund, G., Schmitt, T., Forslund, K., Kostler, T., Messina, D. N., Roopra, S., et al. (2010). InParanoid 7: New Algorithms and Tools for Eukaryotic Orthology Analysis. Nucleic Acids Res. 38, D196–D203. doi:10.1093/nar/gkp931

PubMed Abstract | CrossRef Full Text | Google Scholar

Overbeek, R., Fonstein, M., D’Souza, M., Pusch, G. D., and Maltsev, N. (1999). The Use of Gene Clusters to Infer Functional Coupling. Proc. Natl. Acad. Sci. U.S.A. 96, 2896–2901. doi:10.1073/pnas.96.6.2896

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level Expression Analysis of RNA-Seq Experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667. doi:10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, J. S., Allen, S., Faucheux, C., Althnaian, T., and Mount, J. G. (2005). Deer antlers: a Zoological Curiosity or the Key to Understanding Organ Regeneration in Mammals? J. Anat. 207, 603–618. doi:10.1111/j.1469-7580.2005.00478.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Randi, E., Mucci, N., Pierpaoli, M., and Douzery, E. (1998). New Phylogenetic Perspectives on the Cervidae (Artiodactyla) Are provided by the Mitochondrial Cytochrome B Gene. Proc. R. Soc. Lond. B 265, 793–801. doi:10.1098/rspb.1998.0362

PubMed Abstract | CrossRef Full Text | Google Scholar

Rössner, G. E., Costeur, L., and Scheyer, T. M. (2020). Antiquity and Fundamental Processes of the Antler Cycle in Cervidae (Mammalia). Sci. Nat. 108, 3. doi:10.1101/2020.07.17.208587

CrossRef Full Text | Google Scholar

Sasano, Y., Ohtani, E., Narita, K., Kagayama, M., Murata, M., Saito, T., et al. (1993). BMPs Induce Direct Bone Formation in Ectopic Sites Independent of the Endochondral Ossification In Vivo. Anat. Rec. 236, 373–380. doi:10.1002/ar.1092360211

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoick-Cooper, C. L., Weidinger, G., Riehle, K. J., Hubbert, C., Major, M. B., Fausto, N., et al. (2007). Distinct Wnt Signaling Pathways Have Opposing Roles in Appendage Regeneration. Development (Cambridge, England) 134, 479–489. doi:10.1242/dev.001123

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsujimoto, Y., and Shimizu, S. (2000). VDAC Regulation by the Bcl-2 Family of Proteins. Cell Death Differ 7, 1174–1181. doi:10.1038/sj.cdd.4400780

PubMed Abstract | CrossRef Full Text | Google Scholar

Uebbing, S., Künstner, A., Mäkinen, H., Backström, N., Bolivar, P., Burri, R., et al. (2016). Divergence in Gene Expression within and between Two Closely Related Flycatcher Species. Mol. Ecol. 25, 2015–2028. doi:10.1111/mec.13596

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Yang, Y., Wang, L., Ma, T., Shang, H., Ding, L., et al. (2016). Different Gene Expressions between Cattle and Yak Provide Insights into High-Altitude Adaptation. Anim. Genet. 47, 28–35. doi:10.1111/age.12377

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Zhang, C., Wang, N., Li, Z., Heller, R., Liu, R., et al. (2019). Genetic Basis of Ruminant Headgear and Rapid Antler Regeneration. Science 364, eaav6335. doi:10.1126/science.aav6335

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Xu, M., Zhang, X., Chu, F., and Zhou, T. (2018). MAPK/c-Jun Signaling Pathway Contributes to the Upregulation of the Anti-apoptotic Proteins Bcl-2 and Bcl-xL Induced by Epstein-Barr Virus-Encoded BARF1 in Gastric Carcinoma Cells. Oncol. Lett. 15, 7537–7544. doi:10.3892/ol.2018.8293

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Li, J., Zhao, X.-Q., Wang, J., Wong, G. K.-S., and Yu, J. (2006). KaKs_Calculator: Calculating Ka and Ks through Model Selection and Model Averaging. Genomics, Proteomics & Bioinformatics 4, 259–263. doi:10.1016/s1672-0229(07)60007-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhu, J., Tong, T., Wang, J., Lin, B., and Zhang, J. (2019). A Statistical Normalization Method and Differential Expression Analysis for RNA-Seq Data between Different Species. BMC bioinformatics 20, 163. doi:10.1186/s12859-019-2745-1

PubMed Abstract | CrossRef Full Text | Google Scholar