miRNAome expression profiles in the gonads of adult Melopsittacus undulatus

Ethics statement

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 v1.2.1.1 (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 ddH₂O. 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).