RNA sequencing-based identification of microRNAs in the antler cartilage of Gansu red deer (Cervus elaphus kansuensis)

Ethics approval statement

All experimental protocols were approved by the Institutional Animal Care and Use Committee of Qinghai University (Xining, China), and all methods were carried out in accordance with approved guidelines and regulations (code: SL-2022024). All procedures involving animals were conducted in accordance with the US National Institutes of Health Guide for the Care and Use of Laboratory Animals (National Academies Council, 2011). This study was carried out in compliance with the ARRIVE guidelines for reporting animal research (Percie du Sert et al., 2020).

Animals and samples collection

Samples were collected from a C. elaphus kansuensis population of 200 individuals, raised in a semi-wild setting and fed under the same conditions (Shandan County, Gansu Province, China). Cartilage tissues in the antler tips at different growth stages (30, 60, and 90 d) were collected from three healthy adult individuals (4–5 years old, male) for RNA extraction. Three antler samples from three different individuals were prepared for RNA-seq and qRT-PCR. Antlers were collected after anaesthetising deer with special Mian Naining (anaesthetic, No.9812, People’s Liberation Army (PLA) Military Supplies University Research Institute, China), and hemostasis measures were taken immediately after velvet cutting. After recovering from anaesthesia, deer were returned to the herd.

RNA extraction, library preparation, and sequencing

Total RNA was extracted using TRIzol™ Reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. Extracted RNA was treated with DNase I (Promega, Madison, WI, USA) to remove contamination. To ensure that the RNA samples were qualified for sequencing, concentration and integrity were assessed using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and RNA Nano 6000 Assay Kit with the Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA), respectively.

In this study, sample collection, high-throughput sequencing and data collection were conducted as previously described in Chen et al. (2022). Specifically, a total of 2.5 ng RNA per sample was used as input material for library preparation. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (New England Biolab, Ipswich, MA, USA) following manufacturer’s protocol, and index codes were added to attribute sequences. Library quality was assessed using Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).

Pooled miRNA libraries were sequenced for 50 cycles using HiSeq® Rapid SBS Kit v2 and HiSeq® 2500 platform (Illumina, San Diego, CA, USA). Sequencing data are available at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) (BioProject: PRJNA731682).

Data analysis

Raw reads in FASTQ format were preprocessed using in-house Perl scripts. In this step, clean reads were obtained via trimming adapters and removing reads containing poly(A) tails, low-quality bases, and sequences <18 or >30 nucleotides long. Additionally, Q20 and Q30 scores, GC content (%), and sequence duplication level of cleaned data were determined.

Second, cleaned reads were aligned against the Silva (Quast et al., 2013) (https://www.arb-silva.de/), Genomic tRNA (Chan & Lowe, 2016) (GtRNAdb; http://lowelab.ucsc.edu/GtRNAdb/), Rfam (Kalvari et al., 2018) (http://rfam.xfam.org/), and Repbase (Bao, Kojima & Kohany, 2015) (http://www.girinst.org/repbase) databases using Bowtie (Langmead et al., 2009) to remove sequences that matched to rRNAs, tRNAs, snRNAs, snoRNAs, other ncRNAs, and repeats. Then clean small RNAs with a length range from 18 to 35 nucleotides were further analyzed. Third, Randfold tools (Bonnet et al., 2004) (http://www.aquafold.com) was used for predicting novel miRNA secondary structure. Fourth, the miRBase v21.0 database (Kozomara & Griffiths-Jones, 2014) (http://www.mirbase.org/) with 0 or 1 mismatch and miRDeep2 (Friedländer et al., 2012) with default parameters were used to identify known miRNAs and predict novel miRNAs. The miRNA expression profiles were determined based on TPM values using the normalized formula (Zhao, Ye & Stanton, 2020): TPM = 10⁶ ∗ [Reads mapped to transcript/Transcript length]/[Sum (Reads mapped to transcript/Transcript length)]. (Note: “Reads mapped to the transcript” refers to read number of one miRNA, “Transcript length” refers to length of that one miRNA).

Based on the normalization counts, the significantly differentially expressed miRNAs in three stages were determined by their fold change. Based on a power analysis using the IDEG6 (Romualdi et al., 2003), we determined that our design had over 99% power to detect differentially expressed miRNAs (DEMs) at | log2(FC) | > = 1 and FDR < = 0.01 (Hou et al., 2019). We then screened for DEMs in the 30 d vs. 60 d and 60 d vs. 90 d pairwise comparisons.

MiRNA target prediction and functional analysis of candidate target genes

We put more focus on the 25 selected DEMs and the 30 specifically expressed miRNAs in 60 d library. Potential target genes of novel miRNAs and specifically expressed miRNAs were predicted using Miranda v3.3a (https://anaconda.org/bioconda/miranda) and TargetScan 7.2 (http://www.targetscan.org/vert_72/) with default parameters. Genes simultaneously predicted by both methods were considered as candidate targets of selected miRNAs. Functional annotation of these candidates were performed using GO enrichment (Ashburner et al., 2000) and KEGG pathway (Kanehisa & Goto, 2000) analyses.

Verification of miRNA expression using qRT-PCR

Stem-loop qRT-PCR was used to verify the accuracy of the sequencing data. Twelve of the selected 25 DEMs had similar names (except species), sequences and expression levels, and the remaining 13 DEMs were used for qRT-PCR. Total RNA was extracted as previously described. Approximately 1 µg of total RNA was reverse transcribed using Revert Aid First Strand cDNA Synthesis Kit (Fermentas, Thermo Fisher Scientific, Waltham, MA, USA). Reactions were incubated at 42 °C for 15 min and then terminated at 80 °C for 5 s. The qRT-PCR assay was performed using TB Green® Premix Ex Taq™ II (Tli RNaseH Plus) Kit (Takara Bio, Shiga, Japan) and CFX96 Real-time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA). The miRNA expression level was detected using stem-loop qRT-PCR and miRNA-specific stem-loop primers. Gene-specific forward primers and universal reverse primers (Sangon Biotech, Shanghai, China) were used for qRT-PCR, with U6 RNA as internal control. Cycling profiles were as follows: 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s, and 60 °C for 30 s. Table 1 lists all primers used in this study. Melting curve analysis was also performed to determine product specificity. All statistical analyses were performed in triplicate. All data are presented as the mean ± standard deviation. Relative gene expression was subsequently analyzed using the comparative 2^−ΔΔCt method (Livak & Schmittgen, 2001; Rao et al., 2013).

MiRNA sequencing

To determine miRNA function during velvet antler RGS, we constructed three small RNA libraries for antler growth centers at 30, 60, and 90 d, generating 14,971,388, 16,985,192, and 15,732,412 raw reads, respectively.

After data preprocessing, we obtained 9,836,383, 13,092,918, and 10,769,127 clean reads for 30-, 60-, and 90-d antler tips, respectively. We then aligned clean reads with Silva, GtRNAdb, Rfam, and Repbase databases using Bowtie (Langmead et al., 2009), a comparison software for short sequences, especially reads generated from highlight sequencing. After filtering out sequences that matched ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), other ncRNAs, and repeats, we obtained 8,646,311 (30 d), 11,699,869 (60 d), and 9,461,208 (90 d) unannotated reads that only contained miRNAs.

For miRNA prediction, the remaining reads were aligned to the deer genome (https://www.ncbi.nlm.nih.gov/genome/?term=red+deer), with no mismatch allowed between unannotated reads and the reference sequence. Respectively, 5,430,521 (62.81%), 7,081,757 (60.52%), and 5,621,148 (59.41%) reads from 30-, 60-, and 90-d antler cartilage tissues were successfully mapped and annotated (Table 2).

Identification of known and novel miRNAs in velvet antler

To identify miRNAs in the three libraries, all mapped reads that were 18–35 nucleotides long were compared with mature animal miRNAs in the miRBase (v21.0) database. Perfectly matched sequences were considered as known miRNAs, while the remaining sequences were designated as novel miRNAs. The analysis identified 1,073 known miRNAs (1,003, 1,052, and 1,021 from 30-, 60-, and 90-d antler tips) and 204 novel miRNAs (184, 190, and 183 from 30-, 60-, and 90-d antler tips). The size distribution of identified miRNAs was similar in all samples, and most of them changed from 20 to 24 nt, which was consistent with the typical size range of miRNAs (Fig. 1; Table S1). The most abundant size class was 22 nt, which accounted for approximately 38.92%, 39.29% and 39.53% in the three libraries. All miRNAs were classified into 188 miRNA families, with the mir-2284 family being the most represented (Table S2).

The expression levels of each miRNA were normalized to the expression of tags per million (TPM) (Table S3). The top 20 known miRNAs with the highest number of reads (Table 3) accounted for 82.36%, 74.23%, and 76.69% of total reads from 30-, 60-, and 90-d antler cartilage tissues, respectively. The top 20 novel miRNAs also gathered the most reads (Table 4), accounting for 69.36%, 70.27%, and 79.65% of total reads from 30-, 60-, and 90-d antler cartilage tissues, respectively. Novel miRNAs had considerably lower sequencing frequencies than known miRNAs (Table S3). The same pattern has also been reported in other species, suggesting that novel miRNAs are typically weakly expressed and known miRNAs are highly expressed in different organisms.

There were 1,128 unique miRNAs co-expressed in 30- (95.03%), 60- (90.82%), and 90-d (93.69%) antler cartilage tissues (Fig. 2). The most abundant miRNA was ggo-miR-148a, then aca-miR-148a-3p, bta-miR-21-5p, chi-miR-21-5p, and sha-miR-21 (Table 3). In addition, eight (0.67%), 30 (2.42%), and 11 (0.91%) miRNAs were specifically expressed (SEMs) in each sample (Fig. 2; Table 5).

Differential expression analysis of miRNAs

Differential expression analysis (| log2 fold change [FC] values | ≥ 1 and P ≤ 0.05) revealed 243 differentially expressed miRNAs (DEMs) across the three growth stages in this study (Tables S4–S5). Pairwise comparisons of 30 d vs. 60 d, 30 d vs. 90 d, and 60 d vs. 90 d revealed 177 (136 upregulated and 41 downregulated), 116 (97 upregulated and 19 downregulated), and 85 (45 upregulated and 40 downregulated) DEMs, respectively. Venn diagrams of DEMs in the three pairwise comparisons are presented in Fig. 3.

To identify key miRNAs controlling rapid antler growth, we selected 25 DEMs that had different expression patterns in the 30 d vs. 60 d and 60 d vs. 90 d comparisons only (Table 6). These are the miRNAs that likely play important roles in antler development during RGS at day 60. Among the 25 DEMs, cfa-miR-146b has the highest expression in 60 d, followed by cgr-miR-146b-5p, and bta-miR-146b, bta-miR-27a-5p, and cgr-miR-181a-3p. Two DEMs were specifically expressed in 60 d antler cartilage tissue, namely chi-miR-411b-3p and CM008039.1_315920. We were more interested in the potential biological function of novel miRNA CM008039.1_315920.

MiRNA target prediction and functional annotation of candidate target genes

MiRNAs usually act via translational repression and/or mRNA cleavage, but some evidence shows that they can also upregulate translation through diverse mechanisms. The rapid growth stage is a special period in the growth process of antler velvet; thus, we pay more attention to miRNAs in 60d. In this study, the target genes of 25 DEMs, 30 miRNAs specifically expressed in 60 d and novel miRNA CM008039.1_315920 were predicted. In particular, 726, 1,759 and 139 targets were predicted for 25 selected DEMs, specifically expressed miRNAs in 60 d and CM008039.1_315920, respectively. We then functionally annotated candidate target genes by using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses.

The GO analysis showed that 726 target genes of the 25 DEMs were classified into 52 functional subcategories. “cellular process” (431 genes), “cell part” (487 genes) and “cell” (486 genes), and “binding” (397 genes) were the top enriched GO terms under the biological process (BP), cellular component (CC), and molecular function (MF) categories, respectively (Fig. 4A). For the specifically expressed miRNAs in 60 d, 1,759 target genes were classified into 53 functional subcategories. Similarly, “cellular process” (1,028 genes) “cell” (1,169 genes) and “binding” (943 genes) were the top enriched GO terms under BP, CC, and MF (Fig. 4B). For novel miRNA CM008039.1_315920, 139 target genes were classified into 44 functional subcategories. “cellular process” (89 genes), “cell part” (98 genes) and “binding” (81 genes) were also the top enriched GO terms under BP, CC, and MF (Fig. 4C).

The KEGG classification showed that target genes of the 25 DEMs were mainly associated with the “Metabolic pathways”, followed by pathways for “Pathways in cancer”, and “PI3K-Akt signaling pathway” (Fig. 5A). Similarly, most target genes of the specifically expressed miRNAs in 60 d were related to the “Metabolic pathways”, followed by pathways for “PI3K-Akt signaling pathway” and “Pathways in cancer” (Fig. 5B). For novel miRNA CM008039.1_315920, most targets were classified to “Metabolic pathways”, as well as pathways of “Gap junction”, “Oocyte meiosis”, “PI3K-Akt signaling pathway”, “Rap1 signaling pathway”, “Ras signaling pathway”, “Proteoglycans in cancer” and “Progesterone-mediated oocyte maturation” (Fig. 5C). KEGG enrichment of miRNA target genes was also analyzed (Fig. S1). Notably, among predicted targets of the 25 DEMs, we discovered that CM008039.1_315920 miRNA was highly enriched in “NF-kappa B signaling pathway” (Fig. S1C), which was mainly associated with apoptosis and cell proliferation.

Quantitative real-time PCR (qRT-PCR) verification

To validate DEM expression profiles, we randomly selected 13 miRNAs, including 10 known and three novel miRNAs, for stem-loop qRT-PCR. We discovered that qRT-PCR results and high-throughput sequencing data were correlated (Fig. 6), validating DEM expression patterns and confirming the reliability and accuracy of our miRNA-seq data.