MicroRNAs (miRNAs) are all non-coding RNAs that play vital role in post-transcriptional regulation of various animals and plants (Bartel, 2009). Almost thirty thousand entries have been released in the latest miRBase database (v21) (Kozomara & Griffiths-Jones, 2014) where each miRNA can target hundreds of messenger RNAs (mRNAs) in diverse binding sites, resulting in the enhancement or suppression of gene expression (Lim et al., 2005). The importance of miRNA is evident from their evolutionary conservation and by the various biological processes in which they are involved, including development and physiology (Friedman & Burge, 2014; Skalsky et al., 2014). Animal miRNA are involved in neuronal cell fate, cell proliferation and differentiation, metabolism, aging, apoptosis and organ morphogenesis, suggest that miRNAs are particularly critical in the development, health, and aging of animals (Ameres & Zamore, 2013).
Sexual dimorphism is a universal trait among animals where morphological and behavioral differences between genders play an important role in their sexual selection. Gonads are the principal reproductive organs that are involved in sexual differentiation wherein they are involved in the production of sex hormones and gametes. Although sexual dimorphism is most visual in birds (example peacock and peahen, hummingbirds, songbirds, paradise fly-catcher), the miRNAs involved in gender differentiation among Psittaciformes are not known. Most recently, few miRNAs that are involved in gender differentiation have been reported in many non-avian species, including human (Ali et al., 2016), fruit flies (Marco, 2014), sheep (Torley et al., 2011), pigs (Mai et al., 2016), marine bivalves (Rosani, Pallavicini & Venier, 2016), Schistosoma mansoni (Marco et al., 2013) and the yellow catfish (Jing et al., 2014). In birds, previous miRNA studies have mainly focused on miRNAs that regulate germ cells in various stages of chicken embryo development and breeder cock testes (Table 1). While the latest release (v21) of the miRBase database contains a total of 1328 mature miRNAs in birds, these sequences are based exclusively on domesticated poultry species and do not represent the vast diversity of Class Aves.
|1||Identified 84 miRNAs in 0.5–5-day-old chicken embryos||Darnell et al. (2006)|
|2||Discovered 449 new miRNAs in chicken embryos||Glazov et al. (2008)|
|3||Detected 160 miRNAs (14 novel) in the embryonic chicken||Hicks, Tembhurne & Liu (2008)|
|4||Found 663 miRNAs in breeder cock testes||Wu et al. (2016)|
|5||Suggested that the MIR202* of chicken might play a critical role in regulating testicular development||Bannister et al. (2009)|
|6||Observed that miR-101, miR-31 and miR-202-5p of chicken had roles in testicular and ovarian development||Cutting et al. (2012)|
|7||Found 55 differentially expressed miRNAs between the ovaries of laying and non-laying ducks||Yu et al. (2013)|
|8||Detected 353 differentially expressed miRNAs between the ovaries of laying and broody geese||Xu et al. (2014)|
|9||Identified 93 differentially expressed miRNAs between the ovaries of mature and immature chickens||Kang et al. (2013)|
Parrots, comprising the parakeets, cockatiels, macaws and cockatoos, are members of the Psittaciformes order. Currently, within the order Psittaciformes at least 397 species are recognized across 94 genera (Gill & Donsker, 2017). This order is classified into four families: Strigopidae (New Zealand parrots), Cacatuidae (Cockatoos), Psittacidae (African & New world parrots), and Psittaculidae (Old world parrots). Parrots are generally recognized as the most remarkable intelligent animals along with corvids, chimps, dolphins, and humans. The budgerigar (Melopsittacus undulatus) is a small parrot from Australia and a popular domestic pet throughout the world (Del Hoyo et al., 2015). It is also the most widely studied parrot species, frequently used as a model organism for behavioral studies, feather pigmentation and neuroscience, specifically in the fields of behavioral and neurosciences due to their easy availability, small size and easy breeding in controlled conditions. Budgerigars were also well-fitted for genetic mapping experiments, the polyketide synthase was abolished by the R644W substitution, which blocked the synthesis of yellow pigmentation in the budgerigars (Cooke et al., 2017). Several studies have shown that vocal learning in adult budgerigars is gender-biased. In males, right-sided dominance of molecular neuronal activation was found in answer to mate calls in the CMM male (Eda-Fujiwara et al., 2016), several studies found that vocal learning in adult budgerigars is sex-biased (Striedter et al., 2003). These differences in learning have been well documented at various developmental stages in males and females as well (Hoeschele & Bowling, 2016). Nevertheless, the biological functions of gonadal miRNAs in budgerigars are largely unknown and the sexual differential expression (DE) profile of gonadal miRNAs in budgerigars has not been reported.
In the current study, miRNA expression profiles were collected in two groups of budgerigars (male testes and female ovaries) using Illumina paired-end sequencing technology. Based on the new data generated for M. undulatus and the existing information from miRBase, we achieved the following goals: (1) identification of known miRNAs and novel miRNAs in M. undulatus gonads; (2) nucleotide bias in these miRNAs; (3) identification of those that are differentially expressed.
Animals in the current study were authorized by the Ethics Committee of Anhui Normal University (Anhui, China) with authorization number #20150612.
Tissue collection, RNA preparation and sequencing
Six adult 1.5-year-old budgerigars (three males and three females) were obtained from Wuhu (31°33′N, 118°37′E, southeast of China) in 2015. Total RNAs from six gonadal samples were extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The quality and integrity of the RNAs were examined using an ND-8000 spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA) and a 2100-Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), which having a RNA integrity number >7.0 (Fig. S1). To reduce the individual differences between samples, two RNA pools for deep sequencing were prepared using equal amounts of the extracted RNA from three ovaries and three testes. The Truseq Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) was used to isolate small RNA (sRNA) to construct a miRNA library according to manufacturer’s specifications. The resulting cDNA products were sequenced by the Illumina Hiseq 2500 sequencer (Illumina Inc, San Diego, CA, USA). The resulting sequence data have been submitted to the Short Read Archive at NCBI and are available through accession SRR5664259 and SRR5664260.
Analyses of sequencing data
The low-quality reads and adaptor contamination were identified by FastQC v0.11.5 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) using the PHRED algorithm and removed using cutadapt v1.14 (Martin, 2011), resulting in fragments corresponding to RNAs of 14-41 nt in length. Subsequently, analyses of the length distribution and clustering of sRNA reads revealed the characteristics of the sRNAs. Next, the resulting sRNA sequences were aligned against expressed sequence tags (ESTs) stored in NCBI (https://blast.ncbi.nlm.nih.gov), Rfam 11.0 (Burge et al., 2012), and RepBase (Jurka et al., 2005). The aligned sequences were removed, and the small RNAs were annotated into different categories, such as rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolus RNA (snoRNA), and reads mapped to those were removed as well. The Bowtie v220.127.116.11 (Langmead et al., 2009) and the SOAP2 (Li et al., 2009b) software were used, allowing for ≤1 mismatch, to map with the budgerigar genome (v6.3) (Accession: NW_004848282.1) (Ganapathy et al., 2014). Conserved miRNAs were identified in the budgerigar by matching the unannotated data to precursor and mature miRNA sequences in miRBase v21. Venn diagram drawn using BioVenn (Hulsen, Vlieg & Alkema, 2008).
Prediction of potential miRNAs and analysis of miRNA function
miRDeep2 software was used to discover known and novel miRNAs from sequence data (Friedländer et al., 2012). The Dicer-binding sites and the free energies were combined in evaluating these candidate miRNAs to assign them score numbers that correspond to the reliability if the miRDeep2 predictions (Gong et al., 2017), then, the miRDeep2 predictions were used to map the remaining results with the genome (mismatch < 1). The criteria we used to distinguish miRNAs from other classes of small RNAs amounted to a score of total reads >5, and a true positive prediction ≥90%. Mfold program based on the free energy minimization, was used to predict their propensity to form hairpin loops as potential pre-miRNAs (Zuker, 2003).
To investigate differentially expressed miRNAs, each library was normalized to transcripts per million (TPM) by DESeq (Anders & Huber, 2010). Three miRNAs that were differentially expressed between ovaries and testes were determined to be statistically significant using p-value ≤0.05 (q < 0.01) and a minimum fold change of 2 (Benjamini & Hochberg, 1995). Samtools (Li et al., 2009a) were used to extract 3UTR sequence. Using 3UTR regions of budgerigar for targeting. The potential target genes for the seed sequence of three significantly differentially expressed miRNAs and eighteen miRNAs were predicted using miRanda (score ≥150, MFE (minimum free energy) <−20 Kcal/mol) (Enright et al., 2004), seedVicious (Marco, 2017) and TargetScan (Yatsenko et al., 2005) software. The basic functional targets were classified using Gene Ontology (GO) annotations (Ashburner et al., 2000) and the KEGG pathway database (Kanehisa et al., 2017), false discovery rate (FDR) <0.05 was defined as statistical significance.
Quantitative PCR validation
Five randomly selected miRNAs were used to validate their expression profiles by real-time quantitative PCR (Q-PCR). All reactions were duplicated three times to validate the reliability of the predicted miRNAs. Primers to detect miRNAs were designed using miRprimer2, as previously described (Busk, 2014). The cDNA was synthesized using Rayscript cDNA Synthesis KIT (GENEray, GK8030). Based on the protocol of the AceQ TM Q-PCR Probe Master Mix (Vazyme, Q112-02), the PCR reaction and temperature conditions were performed in a two-step Q-PCR method using the ABI 7500 qPCR instruments (Applied Biosystems Inc., Foster City, CA, USA). The housekeeping gene 5S ribosomal RNA was used as an internal normalization control. The Q-PCR reaction mixture (20 µL) contained 5 µL cDNA product, 0.5 µL primer forward and 0.5 µL primer reverse, 10 µL Taqman Mix with 0.4 µL ROX, and 3.6 µL ddH2O. PCR cycles were as follows: 2 min at 95 °C, followed by 40 cycles of 10 s at 95 °C, and 49s at 60 °C. All reactions were duplicated three times, and after amplification, melting curves were analyzed for all reactions. We calculated relative quantification expression results by the 2−ΔΔCT method (Livak & Schmittgen, 2001). Analyses were performed in R (http://www.R-project.org/) using the t-test function with a p-value ≤0.05, which was used to detect the significantly differentially expressed miRNAs between ovaries and testes.
Small RNA sequence profile
Two small RNA libraries were sequenced using Illumina Hiseq 2500 in a single lane. A total of 9,191,836 and 6,383,927 raw reads were obtained from the ovarian and the testicular samples, respectively. After filtering for adaptors and low quality reads, 82.03% and 84.43% of the total reads were recovered (Table 2). The proportion of clean and unique reads that matched to the M. undulatus genome were 6,927,475 (91.88%) / 448,234 (78.34%) and 4,646,370 (86.21%) / 2,007,263 (84.68%) in the ovarian and testicular libraries, respectively. Fewer ovarian reads mapped to CDS, intron, and 3UTR upstream and downstream, compared to those of the testis; however, more ovarian reads mapped to 5UTR sequences (Fig. 1). We also mapped unique reads to the Rfam, RepBase, EST database and miRBase (v21). The sequences that mapped perfectly to the pre-miRNAs and mature miRNAs in miRBase (v21), were considered to be mappable reads. To avoid influencing miRNA identification, the unmappable reads were removed by searching against noncoding RNAs (tRNA, rRNA, snoRNA, and snRNA) deposited in the RepBase and Rfam databases (Fig. 2). Finally, 7,540,057 cleaned sequences of ovaries representing 572,201 unique reads and 5,389,781 cleaned sequences of testes representing 2,370,340 unique reads were used for subsequent analysis.
Identification of new members of known miRNAs and novel miRNAs
The unique reads from the two deep-sequencing libraries were mapped to the budgerigar genome. The secondary stable hairpin structures were predicted for identification of new miRNAs (see Fig. S2). The abundances of known miRNAs that mapped to the budgerigar genome and miRBase v21. Predicted pre-miRNAs were analyzed and represent 101 miRNA families, including 254 miRNAs previously identified in Aves (see Table S1). Differences in expression frequencies among miRNAs reads were detected.
MiRNA family represents sequences that evolved from a common ancestor. The let-7 family is a conserved miRNA family, both in sequence and function and plays a vital role in animal development. To date, more than thirteen types of let-7 family miRNAs have been identified in animals, and the seven detected in this study, mun-let-7, mun-let-7a, mun-let-7c, mun-let-7e, mun-let-7f, mun-let-7g, mun-let-7i have not been previously described. These members of let-7 family in budgerigar are highly differentially expressed.
We identified several putative novel miRNAs using the reads that did not map to known miRNAs. A total of 141 novel miRNA candidates were predicted (Table S2). In total, 395 miRNAs were identified including 254 known and 141 novel miRNAs from the two libraries. Furthermore, 282 of these were co-expressed in male and female gonads, while 113 were gender-specific: 88 female and 25 male (Table S3). Compared with the Aves database (zebra finch and chicken) in miRBase v21, we used Venn diagrams to compare three species (zebra finch, chicken and budgerigar), 135 miRNAs were found co-expressed in the three species, 76 miRNAs were found in zebra finch and in budgerigar, and 43 miRNAs were found in chicken and in budgerigar (Fig. 3).
miRNA nucleotide bias
The first 5′ end nucleotide of the budgerigar miRNAs, of any length, was not frequently uridine, as seen in 42.8% of gonadal samples, various nucleotides are detected at the 5′ end of miRNA sequences. In general, U is the predominant nucleotide at 5′ end of all miRNAs of all lengths. Adenine accounts for 23.8%, Cytosine accounts for 23.3%, Guanidine accounts for 10.1% across all miRNA nucleotide positions. Averaged across both gonad samples, 66.6% of the nucleotides consisted of A+U in the first bias (Fig. 4). In the seed region of the miRNAs, i.e., the second to eighth nucleotide positions, A-U was the most frequently observed pair.
Differentially expressed and sex-biased miRNAs identification
Notably, three of the miRNAs (mun-215-5p, mun-novel24-5p and mun-novel10-3p) were expressed differentially between the gonadal samples with statistical significance, differentially expressed miRNAs were detected using both the absolute value of fold change >2 and q-value <0.01 (Table 3). Consequently, 11 of these were among the 141 novel miRNAs found in the testes but not found in the ovaries. Likewise, mun-miR-203-5p was expressed in the testes, exclusively. Of the remaining detected miRNAs expressed in both gonads. For instance, mun-miR-194-5p and mun-miR-375-3p in ovaries, and mun-miR-2954-3p, mun-novel101-3p, mun-novel102-3p and mun-novel8-3p in testes were expressed more than average. Due to these, 21 miRNAs were selected for following functional analysis.
Quantitative-PCR analysis of gonadal miRNAs
To confirm the reliability of the miRNA-seq data and bioinformatic predictions, stem-loop Q-PCR assays were conducted. Mun-miR-215-5p, mun-miR-18b-5p, mun-miR-2954-3p, mun-novel10-5p, and mun-novel24-3p were randomly selected for validation. All results were concordant with their relative expression trends for our miRNA-seq analyses basically, and confirmed that these miRNAs exist in the budgerigar. The expression level of mun-miR-18b-5p, mun-miR-2954-3p, mun-novel10-5p, and mun-novel24-3p in the male samples were higher than in the female samples, whereas mun-miR-215-5p expression was higher in the ovaries. There were a minor differences between the miRNA-seq analysis and Q-PCR assays results. For example, mun-novel10-5p and mun-novel24-3p expression were not detected in the sRNA sequence data of the ovaries (Table 3), but weak expression was detected in the validation experiments (Fig. 5).
Putative target prediction for known and novel miRNAs
The putative target genes for these 21 miRNAs were predicted by miRanda, seedVicious and Targetscan. One thousand, five hundred and sixty-seven genes were determined to be feasible targets and had complete complementarity to the seed sequence of 21 miRNAs. These predictions suggest that a single miRNA might target more than one mRNA, such as mun-miR-2954-3p, which is predicted to target 129 budgerigar genes (Table S4). Similarly, one gene can be controlled by one or more miRNAs. For instance, mun-novel102-3p could target FOXG1, and FOXP1 has three miRNAs target sites (mun-novel45-3p, mun-novel102-3p and mun-miR-2954-3p), and mun-novel7-3p, mun-novel36-3p, mun-novel102-3p and mun-miR-203-3p could target FOXP2 (Table S4). The miRNAs identified in this study can target multiple transcription factor genes, such as GATA4, GTF2E1, GTF3C3, JUN, LZTFL1, MBTPS1, MTF1, NFYA, NFYC, PBX3, RFX4, and SOX10. Moreover, computational analysis of the 21 miRNA sequences with known, functional miRNA–mRNA regulatory modules suggested these miRNAs may target genes encoding KPNA3, PKP4, CHN1, SMG7, FAM53A, and calmodulin, which are included in a series of important physiological processes and metabolic networks.
GO annotation and KEGG pathway analyses
The genes found to be potentially regulated by miRNAs from this study were annotated using GO annotation and KEGG pathway analyses. Gene Ontogeny annotations were classified as cellular component, biological process and molecular function, using GO rank 2 with p-value ≤ 0.05 (FDR < 0.05). We found that many of the miRNAs detected in this study were involved in the organogenesis. Seven subcategories within “cellular components” were found, with “cell part” and “membrane-bounded organelle” being most represented. Sixteen subcategories of “biological processes” were also identified, with “single-organism cellular process” being most abundant. Furthermore, many genes were assigned to five subcategories in “molecular function,” with the largest proportion in “protein binding” (Fig. 6). Notably, mun-miR-215-5p, which demonstrated significantly biased expression in the ovaries, is likely a primary modulator of a protein complex involved in signaling and/or catalytic enzyme activity, like RABGAP1, PRR5, FGF13, SDCBP, TIAL1, PRDX4 (Table S4). From the GO term, we found nine GO annotations were related to female, which the detected miRNAs could regulated some genes, like FKBP4, C1QBP, VMP1, AKT1, BMPR1B, ADAMTS1, PHB2, IF2B2, MED1, DACH2, CHD7, AK8, TBP, TAF4, DACH1, LHX9, TYRO3, KIT, EIF2B2, DACH2 and PDGFRA, thirteen GO terms were found related to male and six GO annotations wete detected as gonadal function as well (Table S5). In addition, FOXP1 (mun-miR-2954-3p, mun-novel102-3p and mun-novel45-3p) and FOXP2 (mun-miR-203-5p, mun-novel102-3p, mun-novel36-3p and mun-novel7-3p) were detected to regulate vocal learning, moreover, VDAC1, CNTN2, GRIN1, FOXP1, FOXP2, HIF1A, VDAC3, CAMK4, PRKAR1B, ITGA8, GRIA1, PLK2,HMGCR, KCNAB1, CHST10, RELN, FGF13, ABI1, APP, SRF, LIS1, AK8, CTNND2 and ATAD1 were illustrated to regulate learning ability (Table S5). The KEGG pathway analysis demonstrated that the target genes were related to significantly expressed miRNAs. According to the KEGG pathway analysis (FDR < 0.05), three pathways were significantly enriched, such as the cell communication (FDR = 0.00317745), excretory system (FDR = 0.00317745) and signal transduction (FDR = 0.00446146) (Table S6).
MiRNA research is needed in more phylogenetically disparate avian species to obtain a more accurate concept of the miRNAome. Furthermore, some important avian miRNAs involved in regulating gonadal sex differentiation and development were demonstrated (Table 1). However, avian studies of sex-biased miRNA expression between ovaries and testes are limited. The role of non-coding RNAs in the gonads is an area of active research. Here we used Illumina sequencing to investigate the differentially expressed miRNAs in the male and female gonads of budgerigars. We observed numerous miRNA families in our data that may potentially work as key regulators of gene expression. For example, the let-7 family has been shown to function as a heterochronic switch, and loss of these could cause periods of cell fate reiteration in adults. In contrast, increasing the gene dosage would led to premature expression in adult fates (Reinhart et al., 2000). Consequently, the let-7 family miRNAs are considered highly conserved in Animalia (Hertel et al., 2012). Let-7 and miR-125 are associated with polycistronic transcripts and work as two key regulators of development in Bilateria.
MiR2954, which is known in chicken (Garcia-Riart et al., 2017; Liu et al., 2017) and zebra finch (Lin, Balakrishnan & Clayton, 2014), is an bird-specific gene (absent in the mammalian lineage) and is encoded on the Z chromosome which is known to result in its higher expression in males than females (Lin, Balakrishnan & Clayton, 2014; Luo et al., 2012). It has been proposed that this might affect the neurogenomic mechanisms that lead to sexually dimorphic bird song habituation (Lin, Balakrishnan & Clayton, 2014). Based on bioinformatical and experimental analysis from chicken and zebra finch, miR-2954-3p is male-biased and Z-linked miRNA, which targets across a range of bird species, could help the study field about the evolutionary dynamics of partial dosage compensation and the genetic architecture underlying gonadal characteristics (Warnefors et al., 2017). The bioinformatical prediction and the qPCR validation in the present study has confirmed that there is a male-biased expression of miR2954 in M. undulatus, further corroborating its involvement in sexual dimorphism. Putative genes targeted by miR2954 include a TLE4 transcription factor family, which might be associated with nervous system function, including Ca2+/ calmodulin-dependent protein kinase IV (CaMKIV), SCAMP1, and SMARCA2. Mun-miR-2954 is also related to development, environmental adaptation, the nervous system, signaling molecules and interaction, and substance dependence as indicated in the KEGG analysis. Previous study in zebra finch has elaborated miR2954 could target FOXP2 to regulate vocal learning and detected higher male expression in many tissues (Fu et al., 2014), however, we haven’t detected mun-miR2954-3p target FOXP2, whereas mun-miR-203-3p, mun-novel102-3p, mun-novel36-3p and mun-novel7-3p might regulate FOXP2, mun-miR-2954-3p, mun-novel102-3p and mun-novel45-3p might regulate FOXP1, which related to vocal learning in the budgerigars and we also measured higher male expression in gonads.
Further investigations of GO analysis, miRNA in current study might regulate several sex-related genes. For instance, FKBP4 was considered to be markers of hypospermatogenic testis (D’Aurora et al., 2017), we found higher male expression of mun-novel7-3p and mun-novel72-3p could target FKBP4. mun-miR-194-5p and mun-novel31-3p might regulate SALL1 which the transcriptional regulators of adipose-specific sex-different genes (Karastergiou & Fried, 2017). LHX9 is needed for ovarian function (Workman, 2017), mun-novel72-3p and mun-novel8-3p could regulate it. Mun-miR-215-5p which significantly expressed in ovaries could target PRDX4 (sex-linked gene) (Tippabathani et al., 2017). In mice gonadal development, miR-181a suppressed granulosa cell proliferation by targeting ACVR2A, we detected that mun-novel7-3p could regulated ACVR2A as well (Zhang et al., 2013). Currently, mun-novel7-3p could also target AKT1 (with prior report for sex differences (Seney et al., 2013)). And mun-novel31-3p could target RNF2 (associated with regulation of genetic imprinting (Li, Zhang & Wu, 2017)). mun-novel68-3p might regulate HSF2 which related to sex-determining (Literman et al., 2018), mun-miR-2954-3p and mun-novel72-3p could target KITLG (sex development) (Hersmus et al., 2017) (Table S4). Of all the miRNAs tested, we found three gonad-enriched miRNAs (mun-miR-215-5p, mun-novel10-5p and mun-novel24-3p). In the present study, the presence of these miRNAs in the gonads suggests that they might serve a similar function in the budgerigar. These miRNAs may produce sex-specific responses to potential biological mechanisms that have not yet been described. Although the physiological and biochemical functions of mun-novel10-5p and mun-novel24-3p remain unclear, their differentially expressed patterns indicate that they might play important roles in sexual differentiation and development.
In this study, the whole gonadal miRNAome of budgerigars was sequenced, consisting of a total of 12,929,838 clean reads and 2,942,541 unique reads. Moreover, differential expression of 254 known miRNAs and 141 novel miRNAs were analyzed in the gonadal tissues of budgerigars. The majority of these miRNAs were evolutionarily conserved within chordates while some of them were budgerigar- or avian-specific. In conclusion, this work describes the characteristics of sex-biased miRNAs of M. undulatus and adds new sequences to the avian miRNAome database to facilitate further functional genomic research.