Analysis of conserved miRNAs in cynomolgus macaque genome using small RNA sequencing and homology searching

Sample collection

Peripheral whole blood samples were taken from four unrelated Vietnamese-origin cynomolgus macaques (2–4 years old, one male, three females). The monkeys were healthy and housed at Guangdong Landau Biotechnology Co., Ltd. (Guangzhou China). All experiments were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Guangdong Landau Biotechnology Co., Ltd. (project number: IACUC-003). Peripheral whole blood samples (2 ml) were used for small RNA extraction. Firstly, red blood cells were selectively lysed using an EZNA Blood RNA Kit (Omega, USA) and white blood cells were collected by centrifugation (450 g, 10 min, 4 °C). To avoid interference of 18S and 28S rRNA, small RNAs (18–200 nt) were isolated from white blood cells using RNAiso (TaKaRa, China). The samples were dissolved in 30 µl of diethylpyrocarbonate-treated water and stored at −80 °C prior to small RNA sequencing at the Beijing Genomics Institute (BGI, Shenzhen, China).

Small RNA library construction and sequencing

Four small RNA libraries were prepared and sequenced using standard methods at BGI. Before library construction, the integrity and concentration of samples were checked using an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit, USA). The RNAs were size-fractionated on a 15% PAGE gel, and small fragments (range of 18–30 nt) were collected and then amplified by real-time polymerase chain reaction (RT-PCR) using adaptor primers for 17 cycles. The fragments (140–160 bp in length) were then fractionated with a 15% PAGE gel. The final library was quantified by determining average sequence length using the Agilent 2100 Bioanalyzer (Agilent DNA 1000 Reagents, USA) and quantifying the library using quantitative RT-PCR (StepOnePlus Real-Time PCR System, Applied Biosystems, Foster, USA). The fragment sizes of the four libraries were all between 140 bp and 160 bp and the concentration of each library was more than 8.5 nM, thus meeting the minimum requirements for the construction of a sequencing library. Finally, the libraries were used for sequencing analysis on an HiSeq 4000 System (Illumina, USA). The primary data were then subjected to low-quality filtering. Reads that satisfied any of the following parameters were removed: (1) reads aligned to adaptors or primers with less than three mismatches, (2) reads with >10% unknown bases (N bases), (3) reads with >50% low-quality bases in one read, and (4) reads <18 nt. The remaining clean reads were used in the following analyses. The clean reads of the four samples were submitted to the European Nucleotide Archive (ENA) under study accession number PRJEB38080.

Analysis of sequencing data

The cynomolgus macaque reference genome (v. macFas5) was downloaded from the UCSC Genome Browser (http://genome.ucsc.edu). Clean reads were aligned to the reference genome using BGI software SOAP2 (Li et al., 2009b). Sequences perfectly mapped to the genome along their entire length were further annotated by alignment against macaque genome annotations and known annotation documents from the NCBI/UCSC/Ensembl databases (Benson et al., 2013; Kalvari et al., 2018). Sequences matching rRNA, tRNA, snRNA, snoRNA, and repeat associated sequences or overlapping with exons were discarded, and the remaining sRNA tags were aligned to a reference set of all animal miRNAs from miRBase (release 22) using BLAST by allowing no more than two mismatches outside the seed region. To distinguish authentic miRNAs from false annotations, candidates corresponding to known miRNAs from miRBase were considered conserved miRNA orthologs in the cynomolgus macaques when they were supported by at least 10 reads of mature sequences from the four samples and detected in at least three individuals according to criteria described previously (Veeranagouda et al., 2015; Kozomara & Griffiths-Jones, 2014). The pre-miRNA sequences were predicted based on the presence of hairpin structures analyzed by MIREAP under default settings (https://sourceforge.net/projects/mireap/) (Tan et al., 2013).

Detection of miRNA genes by homology searching

BLAT and BLAST softwares were used for homology searching (Kent, 2002). All mammalian pre-miRNA sequences from miRBase (release 22) were aligned to the cynomolgus macaque genome. The candidates were assumed to be conserved pre-miRNA sequences in cynomolgus macaques only when they fulfilled three criteria: (1) similarity of full length pre-miRNA >93%; (2) gene location corresponding to only one miRNA (although can have multiple gene mapping information); and (3) in mature sequences, no more than two mismatches outside the seed region and no more than one mismatch in the seed region compared with existing miRNAs (Veeranagouda et al., 2015). Furthermore, 10 nt up and downstream of the conserved pre-miRNA sequences were analyzed by MIREAP (Tan et al., 2013) and RNAfold (Lorenz et al., 2011). Stem-loop hairpins were considered acceptable only when they fulfilled three criteria: (1) the mature miRNA present in one arm of the stem; (2) absent of large internal bulges (n < 6); and, (3) a minimal folding free energy (MFE) of predicted pre-miRNA secondary structures of no more than -15 kcal/mol. After mapping all pre-miRNA candidates with acceptable secondary structures to the genome, they were further annotated by alignment against known annotation documents from the NCBI/UCSC/Ensembl databases (Benson et al., 2013; Kalvari et al., 2018). Small RNAs were sorted by their chromosomal locations on the genome. Based on criteria described previously (Marco et al., 2013), if two neighboring miRNA loci were located within 10 kb and on the same strand, they were considered clustered miRNA genes. The miRNA orthologs were grouped into different miRNA families according to sequence similarity to known miRNA family members. ClustalW was used for the alignment of sequences. Data patch processing was implemented in python.

Target prediction

The target genes of miRNAs were predicted by TargetScan 7.2 (http://www.targetscan.org), miRanda (http://www.microrna.org), and PITA (https://tools4mirs.org/). Gene ontology (GO) was applied for functional annotation analysis of the predicted targets. Enriched pathways involved with the miRNAs were determined according to the Kyoto Encyclopedia of Genes and Genomes (KEGG). Both GO and KEGG analyses were performed using the clusterProfilers package in R (Yu et al., 2012), and an adjusted P-value (Benjamini–Hochberg false discovery rate correction) of less than 0.05 was considered significant.

Correlation between age and miRNA expression

Considering the variable RNA concentrations and biases in the RNA-seq procedures for each sample, the expression levels of mature miRNAs were normalized by transcripts per million, i.e., dividing the number of reads for each miRNA by the total number of clean reads for each sample and multiplying it by 10⁶ (Xu et al., 2013). The ages of the miRNAs were assessed based on the set of homologous sequences identified in other vertebrate species in miRBase (Meunier et al., 2013). If the miRNA homologous sequences were shared among vertebrates or mammals, these miRNAs were considered old and grouped into vertebrate or mammalian categories, respectively. If the miRNAs were specific to primates or Cercopithecidae, they were considered young and grouped into primate or Cercopithecidae categories, respectively. The average expression levels of miRNAs in different age groups of the four cynomolgus macaques were statistically analyzed by GraphPad Prism 8. Spearman correlation coefficients (R values) between the age of miRNAs and average expression levels of mature miRNAs were computed by matrix analysis in GraphPad and P-values were derived at a confidence interval of 95%.

Identification of miRNA genes in cynomolgus macaque genome by homology searching

To identify previously unknown miRNA sequences in macaques (Yang et al., 2014; Veeranagouda et al., 2015), we used all mammalian pre-miRNA sequences deposited in miRBase as a reference for homology searching. According to the criteria described above, a total of 1,697 homologous sequences were perfectly mapped to the cynomolgus macaque genome with nonoverlapping gene locations and were then used to analyze their secondary structures. 1,526 of the 1,697 homologous sequences were predicted to form acceptable stem-loop hairpins and were considered conserved miRNA candidates. Additionally, we annotated the 1,526 candidates in the cynomolgus macaque genome using the NCBI/UCSC/Ensembl databases and miRBase for screening of homologies in other species. A candidate identified by homology searching was deemed to be of high confidence if it met any of the following conditions: (1) annotated as miRNA in other species or homolog has experimental support; and, (2) reported as a Macaca fascicularis miRNA previously. In total, 1,455 of the 1,526 pre-miRNA candidates showed high confidence. The remaining 71 pre-miRNA candidates lacked support in the other species and were excluded from the following analysis. 672 of the 1,455 candidates were reported previously and the remaining 783 gene loci were firstly annotated in cynomolgus macaque genome (Yang et al., 2014; Veeranagouda et al., 2015). Compared the 1,455 pre-miRNAs with their homologs in other mammals, the stem-loop sequence similarities of 1,186 candidates (81.5% of 1,455) exceeded 95%. Furthermore, 1,074 candidates (73.8% of 1,455) contained 0 mismatches in the -5 p or -3 p mature sequences, with the remaining 299 and 82 having one and more than one mismatch, respectively. In comparison to human miRNAs, 1,083 of the 1,917 human pre-miRNAs deposited in miRBase had orthologs in the cynomolgus macaque genome. Thus, most pre-miRNA and mature miRNA sequences were highly conserved. Most of the 1,455 miRNA sequences in the cynomolgus macaque were single copies, although others had paralogs at several gene locations. For instance, mir-548, a poorly conserved primate-specific miRNA gene family, was represented by 53 paralogs distributed across 18 chromosomes, similar to that observed in humans (Liang, Guo & Liu, 2012). The names, stem-loop sequences, gene locations, homology similarities, mismatches in mature sequences, and MFE values of the predicted hairpin structures for the 1,455 pre-miRNA candidates are listed in Table S1.

Identification of miRNAs in cynomolgus macaque using small RNA sequencing

A major obstacle in research on primate miRNAs is the lack of experimentally verified miRNAs in nonhuman primates. To better characterize the identified miRNAs in cynomolgus macaques, small RNA-seq was performed using peripheral blood leucocytes from the four animals. On average, 14 million clean reads were obtained and about 9 million small RNAs (∼5,000 unique reads) were mapped to miRBase for each sample (Table 1). As described in the methods, sequences that matched against orthologs in miRBase were considered as known miRNAs, otherwise, they were considered novel. A total of 546 known mature miRNA sequences and 136 potentially novel mature miRNAs (Table S2) were identified. As the 136 potentially novel miRNAs need further confirmation, only the 546 known miRNAs were used in subsequent analyses. To understand whether these miRNAs resulted from one or both of arms of a given precursor, we assigned the 546 mature miRNAs based on pre-miRNA information. As shown in Table S3, 337 mature miRNAs were expressed from both the -5 p and -3 p arms of 170 pre-miRNAs, and the remaining 209 mature miRNAs were expressed from either the -5 p or -3 p arms of 213 pre-miRNAs. As seven mature miRNAs were derived from more than one pre-miRNA, the 546 mature miRNAs were processed from 383 stem-loop miRNA sequences. We identified 379 of the 383 pre-miRNAs by homology searching (as above), including 357 pre-miRNAs with mapping information for single gene and 22 with mapping information for 2–4 genes. We re-screened the remaining four pre-miRNAs, mir-383, mir-378c, mir-4536, and mir-5698, which were deposited in miRBase against the cynomolgus macaque genome. All showed sequence similarities below 60%, indicating that these four pre-miRNA sequences were distinct to cynomolgus macaques. Analysis of sequencing data indicated that 408 of the 1,455 miRNA loci were expressed in peripheral blood leucocytes. These 408 miRNA loci demonstrated high confidence as they were identified by both RNA sequencing and homology searching and showed homolog support in other species. In total, 335 have been reported in cynomolgus macaques previously (Yang et al., 2014; Veeranagouda et al., 2015), with the remaining 73 considered to be new pre-miRNA loci in cynomolgus macaques (Table 2).

Thus, in the current study, we confirmed 1,455 high-confidence pre-miRNA sequences by homology screening and RNA sequencing. According to their sequence similarities (Dostie et al., 2003; Yue, Sheng & Orwig, 2008; Bartel, 2009), 915 were grouped into 491 families (Table S1). Another 540 were not clustered into currently known families, suggesting new targeting properties. To assess the age of the miRNAs, their precursor sequences were searched in miRBase. The age of miRNAs was based on the group of species in which miRNA homologous sequences were identified. We classified the 1,455 miRNAs into four groups: i.e., 267 conserved in vertebrates, 332 conserved in mammals, 806 primate-specific, and 50 Cercopithecidae-specific. Among the 408 high-confidence pre-miRNA sequences, 179 were conserved in vertebrates, 140 were conserved in mammals, 77 were conserved in primates, and 12 were conserved in Cercopithecidae. We then examined the relationship between average expression level of miRNAs and their age and found a positive correlation (Spearman R: 0.27 to 0.42, P < 10⁻⁴; Fig. 1). Primate-specific miRNAs demonstrated lower average expression levels (10-fold and 160-fold, respectively) than ancient miRNAs in mammals and vertebrates. In addition, Cercopithecidae-specific miRNAs showed the lowest expression level among the four groups.

Characteristics of miRNA gene locations in cynomolgus macaques

According to the annotations of the 1,455 high-confidence pre-miRNA sequences (Table S1), the miRNA encoding genes in the cynomolgus macaque genome were also distributed in diverse locations, similar to that found in other mammals (Kim & Kim, 2007; Morin et al., 2008). Of the 1,455 pre-miRNAs, approximately 46% were located in intergenic regions, 49% were located in introns and 5% were located in RNA coding regions. It is worth noting that 64 pre-miRNA sequences were processed from exons of protein-coding transcripts or non-coding transcripts, including 46 derived from untranslated regions. Six pre-miRNA sequences were found overlapping intron-exon junctions, i.e., SO-miRNAs. Most miRNA gene is encoded by one strand of one locus as above, but in rare cases, two distinct pre-miRNAs are found transcribed and processed on both strands of one locus (Tyler et al., 2008; Picao-Osorio et al., 2015; Mohammed et al., 2018). Two pairs of these so-called antisense or mirror miRNAs were identified in our study, i.e., miR-3120/miR-214 and miR-24-1/miR-3074 (Table S1). Here, miR-3120 was located within an intron of the dynamin-3 gene and miR-214 was embedded on the antisense strand of the same locus. For the miR-24-1/miR-3074 pair, we found that miR-24-1 was embedded within an intron of the AOPEP gene and miR-3074 was processed from the antisense strand. The miR-24-1/miR-3074 mirror miRNAs were first reported in this study in cynomolgus macaques (Cordero et al., 2015).

MiRNA gene clusters in cynomolgus macaques

Previous studies have revealed that miRNAs are significantly enriched in clusters in introns or intergenic regions (Marco et al., 2013; Wang et al., 2016; Altuvia et al., 2005). To profile the clustering patterns of miRNAs in cynomolgus macaques, we grouped pre-miRNAs into distinct clusters following previous study (Marco et al., 2013). If two neighboring miRNA loci were located within 10 kb and on the same strand, they were considered to be clustered miRNA genes. Among the 1,455 pre-miRNAs, 285 miRNA genes were grouped into 78 distinct clusters. When comparing the 78 miRNA clusters to their human counterparts, if a human miRNA was not represented in the macaque clusters, it was screened cross the macaque genome in order to validate the existence of an ortholog in macaques (Wang et al., 2016). As such, 62 miRNA clusters were found to be similar to that in the human genome and six miRNA clusters showed minor changes (Table 3), indicating that most of the clustered miRNAs were evolutionarily conserved. Ten of the 78 clusters were newly emerged in the cynomolgus macaques. According to the similarity of seed regions of the miRNA sequences in each cluster, we grouped the 78 clusters into two sub-classes: 21 homo-seed clusters (miRNAs having identical sequences in the seed region) and 57 hetero-seed clusters (miRNAs having distinct sequences in the seed region). These results are consistent with previous studies suggesting that the composition of miRNA clusters is heterogeneous with a diversiform evolutionary mechanism (Wang et al., 2016). Homo-seed clusters are usually products of local duplications, and hetero-seed clusters are mainly shaped by a combination of de novo formations and local duplications (Wang et al., 2016).

For most miRNA clusters revealed in this study, the number of pre-miRNAs ranged from two to three. Nonetheless, 10 large clusters with no less than five miRNA precursors were also observed (Table 3), which were composed of members from different miRNA families but with the same expression state. Among the 10 large clusters, most miRNAs in the six of them were detected by RNA-seq in this study, but none was detected by RNA-seq in the rest of the four clusters (Table 3). These results indicate that each of the 10 clusters was transcribed as a polycistronic element (Ryazansky, Gvozdev & Berezikov, 2011; Saini, Griffiths-Jones & Enright, 2007). Comparing the compositions of the 10 miRNA-clusters in cynomolgus macaques to their orthologs in humans, eight showed identical composition, whereas the mir-379/mir-656 and mir-512/mir-517 clusters showed minor changes. The mir-379/mir-656 cluster in macaques was positioned on chromosome 7 and contained 42 miRNA genes. Compared to humans, 40 out of the 42 miRNA genes were the same as humans. Simultaneously, a mir-376a copy was added and a mir-300 ortholog was lost. In the macaques, the mir-512/mir-517 cluster was located on chromosome 19 and contained 40 miRNA genes. Most shared common seed regions of the miR-515 family and likely originated from a common ancestor (Rodriguez et al., 2004; Zhang, Wang & Su, 2008). Comparing the miRNA composition of the mir-512/mir-517 cluster to the C19MC in humans, 38 of the 40 miRNA genes had orthologs in humans, two (mir-518g and mir-523b) were unique to macaques, and eight, including mir-522, mir-527, and six mir-520 variants, were not found in the macaque genome. These data show that although most C19MC members are conserved throughout the primate kingdom, frequent miRNA gains and losses, such as the gain of mir-518 and mir-523 copies and loss of mir-520 copies in macaques, have occurred (Berezikov et al., 2006; Iwama et al., 2013; Awan et al., 2017). These findings are consistent with previous results suggesting that the C19MC region is under rapid evolution (Iwama et al., 2013; Bentwich et al., 2005). Specifically, this region is enriched by dispersed Alu elements, which are believed to have facilitated the expansion of C19MC (Zhang, Wang & Su, 2008; Lehnert et al., 2009). To understand the functions of miRNA clusters, the potential targets of the two very large miRNA clusters were predicted using TargetScan, followed by GO functional annotation. Among the target genes of the mir-379/mir-656 cluster, those relevant to cellular components (CC) were enriched in synapse, membrane, and cell projection; those relevant to molecular function (MF) were mainly clustered in binding and catalytic activity; those relevant to biological process (BP) were mainly enriched in cellular process, biological regulation, metabolic process, development process, and response to stimulus (Fig. 2). For the mir-512/mir-517 cluster, the functional annotations of their target genes were similar to those of the mir-379/mir-656 cluster. KEGG analysis of the two very large miRNA clusters suggested that both regulate normal growth, development, and pathological processes through many pathways, such as MAPK signaling and Axon guidance (Fig. 2), consistent with previous studies on humans (Haller et al., 2010; Lavon et al., 2010; Cortez & Calin, 2009; Grün et al., 2005; Donker et al., 2012; Noguer-Dance et al., 2010; Tsai et al., 2009; Li et al., 2009a; Kleinman et al., 2014; Ward et al., 2014). The other eight large miRNA clusters in macaques are with conserved compositions to that in humans, indicating that their functions are also evolutionarily conserved.