Enrichment and verication of differentially expressed miRNAs in the ovarian tissue of Leizhou black ducks

Ovary is an important reproductive organ for poultry. MicroRNA (miRNA) is a highly conserved class of small non-coding RNA that function in a specic manner to post-transcriptionally regulate gene expression in organisms. Currently, miRNA has been studied extensively but research on duck ovarian tissue is relatively rare. Thus, in this study, we performed the rst miRNA analysis of ovarian tissues in Leizhou black ducks with low and high rates of egg production. Method: Using high-throughput sequencing technology, miRNA library was constructed to obtain miRNA expression prole; differentially expressed miRNAs were screened, to predict miRNA target genes.


Conclusion
The research results enrich the data resources of duck miRNAs, and provide a theoretical basis for further research on the laying performance of poultry and molecular-assisted breeding of poultry in the future.

Background
The duck industry occupies a pivotal position in China's agriculture, especially the livestock industry. The export volume of duck eggs even ranks rst in the world [1]. Therefore, the level of duck egg production has a great impact on China's economy. Leizhou Black Duck is a local high-quality small duck species found in the Leizhou Peninsula of China. The female duck has an early birth period (of 90-120 days), a long egg production period (of 2 years), and a high egg production rate (of 200/500 days). The egg production rate can reach more than 90% in the peak period, showing a good market prospect [2], which has been studied by our research group for many years [3][4][5][6]. However, due to several factors, egg production is extremely uneven. The egg production performance of poultry is associated with the ovary [7]. The egg production process of poultry is a complex developmental biological process, which is affected by the degree of ovarian development. The ovary is the most important organ in the reproductive system of female poultry and its development is very important as it is the direct factor affecting egg production in poultry [8]. Studying the genes related to the differences between the high and low-yield ovarian development of poultry will reveal the causes of the differences in high and low egg-laying ducks.
Like other poultry, the laying performance of Leizhou black duck may be regulated by many factors. The ovaries also have the exocrine function of releasing eggs, which affects the egg production performance of poultry [9]. A study found that the signi cant differences in the ovarian structure and follicular development of Gaoyou ducks in the double yolk high and low production groups affected the egg production characteristics of Gaoyou ducks [10]. Therefore, the ovary is used as the research material to nd different-expressed genes related to egg production traits, providing a new direction for improving duck egg production performance. It is of great signi cance to analyze the molecular genetic basis and molecular markers of duck egg production traits. Exploring the ovarian transcriptome is one of the necessary methods to improve animal fertility. Several studies have identi ed a large number of miRNAs in the ovarian tissues of different species, including mice, cattle, sheep, pigs, chickens, geese and other species through high-throughput sequencing technology, and found that miRNAs in these species play an important role in the ovarian tissue [11][12][13].
MicroRNAs (miRNAs) are non-coding RNAs that play an important role in animal post-transcriptional gene regulation by degrading target messenger RNA (mRNA) or inhibiting the translation of target genes [14,15]. In recent years, with the continuous development of various miRNA detection technologies including high-throughput sequencing technology, the research on miRNA in many livestock and poultry has increased. Studies have shown that mature miRNAs regulate downstream target genes in two main ways: inhibit translation and cut transcripts, to regulate the expression of target genes [16]. By highthroughput sequencing, recent studies have illustrated that miRNAs play a role in different organs of ducks, including embryonic breast muscle, ovary, feather follicle, and skin [17][18][19]. Many studies also have found that miRNAs are widely involved in animal cell proliferation, differentiation, apoptosis, and follicular development, including follicular growth, ovulation, and atresia processes [20,21]. Based on miRNA expression pro le data and functional analysis and through molecular biology techniques, the function and regulatory mechanism of key candidate miRNAs in follicular development are veri ed [14].
This provides a scienti c basis for improving the performance of poultry egg production and future molecular breeding research. Therefore, this research used high-throughput sequencing technology analysis to compare miRNAs expression in ovarian tissue of high yield and low yield Leizhou black ducks. The differences between the expressions of miRNAs are analyzed, to explore the different physical conditions of miRNAs in ovarian development and reproduction regulation mechanism to improve egg production traits.

Materials And Methods
Sample collection and RNA isolation The Leizhou black duck used in this experiment was provided by Zhanjiang Hengcheng Cultivation Cooperative. All ducks had the same nutritional standards and environmental conditions. Food and water were available ad libitum. 130 female ducks of the same generation, size, and health status were randomly selected in a single cage. The number of eggs laid in 16-43 weeks (w) was counted to group the ducks into high and low egg production groups. Those which produced eggs >150 were de ned as the high production group (HG) whereas those with <90 eggs were de ned as the low production group (LG).
Six (6) 43-week-old high and low yield Leizhou black ducks were randomly selected. The medulla layer of the ovary was collected from the ducks and quickly stored in liquid nitrogen. The total RNA of each sample was extracted to construct HG and LG gene pools for transcriptome sequencing analysis (HiPure Total RNA Plus Mini Kit, Magen). The remaining ovarian samples were stored in a refrigerator at -80°C for later use.
All procedures of the entire experiment were approved by the Animal Ethics Committee of Guangdong Ocean University.

Construction of small RNA libraries and solexa
Total RNA was extracted from HG and LG ovarian tissue samples using TRIZOL reagent following the manufacturer's protocol (Magen, Guangzhou, China) to construct the sequencing library. The overall process was as follows: through agarose gel electrophoresis, small RNA of 18-30nt was selected and using the special structure of sRNA that 5' has the inte grated phosphate group and 3' has hydroxyl, the total sRNA was ligated to 3' and 5' adapters using T4 RNA ligase. Reverse transcription and PCR ampli cation were performed on the small RNAs connected with the adapters on both sides. The PCR products were puri ed by gel electrophoresis (PAGE) and the separated DNA fragments were the cDNA libraries. Then the library preparations were sequenced on Illumina Hiseq 2500 platform.

Analysis of sequence data
The original image data obtained by sequencing is converted into sequence data by base calling; they were called the raw reads. The raw reads obtained by sequencing were processed for low quality, decontamination, and joint removal, and the total number, species number, and length of the sequence were counted, and all valid sequences are obtained for subsequent analysis. After trimming the 3' adaptors and removing the reads which had poly A, poly T, poly C, and poly G, the clean reads were mapped to the duck genome database (BGI duck 1.0 reference Annotation Release 101) and compared with the distribution of the known miRNAs. Using the tag sequence to align to the position of the genome, the special secondary structure of miRNA was predicted and possible new miRNAs were identi ed using miRdeep2 software [22,23]. Analysis of differential miRNA was performed by edgeR software and edgeR default parameters.
The screening criteria used for differential miRNAs were, the expression level changes more than 2 times and the P-value less than 0.05. The TPM (tags per million) expression of each miRNA was calculated to obtain all miRNA expression pro les and the target genes of differentially expressed miRNAs were predicted using RNA hybrid + SVM-light, Miranda, TargetScan to predict [24]. The DAVID gene annotation tool was used to perform Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.genome.jp/kegg/) was used to analyze candidate target genes of differentially expressed miRNAs to obtain the types and numbers of candidate target genes.
Veri cation of differently expressed miRNAs using real-time polymerase chain reaction The total RNA was reverse transcribed into rst-strand cDNA according to the manufacturer's instructions of miRNA reverse transcription kit (Vazyme, Nanjing, China). Twelve differentially expressed miRNAs were randomly selected for real-time PCR veri cation. The primer sequence of miRNA is shown in the table below (Table 1). MiRNA Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) was used to determine the expression level of miRNA, which was carried out using an iCycler IQ5 Multi color real-time PCR detection system (Bio-Rad, Hercules, CA, USA). The procedures included 1 min of pre-denaturation at 95°C, followed by 40 cycles of denaturation at 95°C for 10 s and 60°C for 30s. Real-time PCR analysis of each sample was done in triplicate using U6 as the housekeeping gene. The data were calculated by the normalized relative quanti cation method followed by 2 -ΔΔCT [25]. All data were calculated according to the normalization of log2 Fold-change (log 2 [JCH/NH-P]) to compare the predicted value and experimental value to verify the accuracy of the sequence.

Data analyses
All data were processed using the MS Excel program. A t-test was employed to analyze the signi cance of the differences between the groups at a level of p < 0.05 using SPSS 13.0 software [26].

Results
Sequence analysis of small RNA in Leizhou black duck ovary The speci c sequencing results are shown in Table 2. It can be seen from the The detailed information of the statistical analysis of small RNA length on the obtained high-quality sequences is shown in Figure 1. Combining the results of the 5 groups, it is found that the total readings and the length with the largest proportion was 22 nt.
Identi cation of conserved and novel microRNAs After miRNA ltering, the results showed that among known miRNAs, miR-143-Y had the highest expression, followed by miR-26-x, miR-125-x, miR-99-x, and other miRNAs. Among the newly predicted miRNAs, novel-m0121-3p expressed the highest amount, followed by novel-m0193-5p, novel-m0180-5p, and novel-m0181-5p (Table 3). In all miRNA expression pro les, miR-143-y, miR-125-x, and let-7-x were highly expressed. The differential expression of miRNAs was enriched by using two standards that were p-value<0.05 and fold-change log2 (LG/HG)>0.05. The speci c results are shown in Figure 2 and Table 4. The volcano plots ( Figure 2) of LG vs HG concluded that the distribution in differentially expressed miRNAs. There were 29 signi cant differentially expressed miRNAs enriched between the two libraries. Among the differentially expressed miRNAs of LG and HG, the up-regulated differentially expressed miRNAs included miR-1175-y (P <0.05), miR-12-z (P <0.05), miR-305-x. (P <0.05), and novel-m0304-3p (P <0.01), and the down-regulated differentially expressed miRNAs, included Let-7-z (P <0.05) and miR-10014-x (P <0.05), miR-1307-y (P <0.01), and miR-34-x (P <0.01), among which miR-192-x was the most signi cant (p-value = 0.000139). Gene ontology analysis of miRNAs target genes Target gene prediction for differentially expressed miRNAs showed that the number of known and newly predicted miRNA target genes is 54 619. The results of GO analysis of miRNAs target genes are shown in Figure 3. The signi cantly enriched GO terms were mainly distributed in biological processes and cellular components. In addition, directed acyclic graph (DAG) was used to display the detailed relationship of the en riched GO terms; the deeper the color, the higher the level of enrichment in DAG dates. In terms of biological processes, target genes were signi cantly enriched in cellular processes, single-organism processes, metabolic processes, and reproductive processes. In terms of cellular components, target genes were signi cantly enriched in cells, cell part, organelle, membrane, and macromolecular complexes. In terms of molecular functions, target genes were signi cantly enriched in binding, catalytic activity, molecular transducer activity, and nucleic acid binding transcription factor activity ( Figure 4).

KEGG pathways enrichment of differential expressed miRNAs target genes
The KEGG database annotation analysis of predicted target gene pathways showed that differentially expressed miRNA target genes were signi cantly enriched in axon guidance (360 target genes predicted), ErbB signaling pathway (209 target genes predicted), TRP channel In ammatory mediator regulation of TRP channels (276 target genes predicted), Proteoglycans in cancer (400 target genes predicted) and other signaling pathways. Signi cantly enriched in oxytocin signaling pathway (360 predicted target genes), GnRH signaling pathway (240 predicted target genes), progesterone-mediated oocyte maturation (147 predicted target genes), and AMPK signaling pathway (196 predicted target genes), and other signaling pathways related to ovarian development ( Figure 5, Table 5). Veri cation of the differential expressed miRNAs To verify the accuracy of sequencing data, we randomly selected 12 miRNAs with known differential expression for veri cation. The selected miRNAs include the up-regulated miRNAs, such as miR-305-X, miR-1175-Y, miR-306-X, novel-m0306-3p, and the down-regulated miRNAs which were miR-192-X, miR-216-X, miR-25-Y, miR-361-X, miR-34-X, miR-34-Y, miR-484-X, and novel-m0007-3p. Through calculating the results that the fold change log 2 (LG/HG), the results showed that had the same variation trend compared with the results (Figure 6).
Discussion miRNAs generally bind to target gene mRNA at the post-transcriptional level and play a regulatory role through translation inhibition or target gene degradation [14,27]. In recent years, the large-scale application of high-throughput sequencing technology and the popularization of bioinformatics technology have greatly promoted the development of animal genomics [28]. There are more and more researches on small RNA sequencing, in animals and plants. For poultry, the level of egg production e ciency is closely related to economic bene ts and is mostly determined by the development of ovaries and follicles [29]. At present, a lot of work has been carried out on the sequencing analysis and identi cation of miRNA, especially the research related to egg production traits [30]. However, in the current research on small RNA sequencing of ducks, no one has performed any research on small RNA sequencing by comparing the ovarian tissues of high-and low-laying ducks.
In this study, high-throughput sequencing methods were used to perform small RNA sequencing analysis on the ovarian tissues of Leizhou black ducks in the high-and low-laying groups. By analyzing the differential expression between the LG and HG groups, 7 up-regulated and 22 down-regulated differentially expressed miRNAs were obtained. Among the differentially expressed miRNAs, miR-34-x was signi cantly down-regulated. A study used high-throughput sequencing technology to analyze the differences in the expression of miRNAs in the ovarian tissues of low-and high-production layers, and screened 17 differentially expressed miRNAs of which gga-miR-34b-5p was signi cantly up-regulated [31]. This result is contrary to that found in this study which may be attributed to the differences in the species. The differentially expressed miRNAs obtained in the sequencing of small RNAs at different developmental stages of Hu sheep's ovary [32] are different from those obtained in the current study. The results obtained by Chen Rong et al. [33] are similar to this study, however, they used the lower samples of Liancheng White Duck and Cherry Valley Duck. Sequencing analysis of thalamic tissue showed that the differentially expressed miRNAs obtained in the results were different from those in this experiment [33]. This difference may be due to differences in the tissues used.
For the veri cation of small RNA sequencing data and results, most of the current researches use uorescence quantitative PCR methods by randomly selecting the differentially expressed miRNAs obtained by sequencing for quantitative detection. The qPCR result is compared to the original sequencing results to judge the accuracy of the sequencing results. A previous study used transcriptome sequencing technology to identify miRNAs in the ovarian tissues of 1 and 8 months old Hu sheep and randomly selected 3 differentially expressed miRNAs for qPCR veri cation [32]. The results were consistent with the differentially expressed miRNAs in the sequencing results. To ensure the accuracy of the small RNA sequencing results, 12 differentially expressed miRNAs were randomly selected for qPCR veri cation. It was found that the qPCR results were consistent with the up-regulation and downregulation trends of the differentially expressed miRNA in the sequencing results. The veri cation method of this test is consistent with the veri cation method used in the above-mentioned research [32]. In this study, the statistical analysis of the sequence length in the small RNA sequencing library was 20-24nt, which accounts for the vast majority with 22 nt as the largest proportion was. This is consistent with the sequence length obtained by Ji et al. [34] on small RNA sequencing of goat mammary tissue. This proved the accuracy of the small RNA sequencing data and results in this study.
In this study, let-7-x was found to be one of the miRNAs highly expressed in the ovarian tissues of the HG and LG groups. Yu et al. [13] used high-throughput projection to analyze the expression pro le of nest goose ellipse and pre-follicle miRNA, and found that gga-let-7b was expressed in high abundance in the library. Zhu Long [35] conducted a small RNA sequencing analysis on goat ovaries in follicular and luteal phases and found that let-7f was highly abundant miRNA, which is similar to the results obtained in this experiment. Early studies found that the let-7 miRNA has a regulatory effect on the development of the ovarian follicles as it is one of the most abundant miRNAs in the ovary [36]. Let-7 can participate in the regulation of follicular development, cumulus-oocyte information exchange [37], ovarian hormone synthesis and release, ovarian cancer cell proliferation, and apoptosis [38]. Therefore, this study concluded that let-7 may be related to the development of ovaries and follicles in Leizhou black duck.
The results also showed that miR-34-x, miR-361-x, and miR-484-x are miRNAs differentially expressed in high-and low-yield ovarian tissues. Studies have found that miR-34a promotes human granulosa cell apoptosis by regulating the expression of apoptosis-related proteins and genes [39]. It has also been shown that in porcine granulosa cells, miR-361-5p can directly act on the 3′UTR of the VEGFA gene, providing a molecular basis for further studies on the regulation of VEGFA expression and the role of miR-361-5p in ovarian follicle development [40]. Again, miR-484 is involved in the growth of mouse granulosa cells [41]. Therefore, it is speculated that these differentially expressed miRNAs may play an important role in the ovarian follicle development in Leizhou black duck.
In this study, the KEGG pathway analysis of the differentially expressed miRNA target genes found that the oxytocin signaling pathway, GnRH signaling pathway, insulin signaling pathway, AMPK signaling pathway, and TGF-β signaling pathways were signi cantly enriched. Studies have shown that these pathways [42,43] are all related to ovarian development. BMPs (bone morphogenetic protein), an important member of the TGF-β signaling pathway, widely act on all levels of ovarian oocytes, granulosa cells, follicles, etc., and regulate the development of follicles through paracrine [44][45][46]. The insulin signaling pathway plays an important role in the physiological and reproductive activities of the bovine ovary and participates in the cell differentiation, development, and apoptosis of the ovary, as well as the integration of reproductive activities [47]. It is inferred that miRNA may have a certain in uence on the development of ovaries by regulating target genes.

Conclusions
In conclusion, the miRNA pro les revealed that different miRNA regulation mechanisms might exist in the ovarian tissues of HG and LG ducks. In this study, the ovarian tissues of high-and low-yield Leizhou black ducks were analyzed for the rst time through small RNA sequencing, and a total of 29 differentially expressed miRNAs were identi ed, of which 7 were up-regulated and 22 were down-regulated. The target genes of differentially expressed miRNAs were predicted to obtain a total of 54 619 target genes. The GO and KEGG function analysis of the target genes found that the targets of differential expression miRNAs were mainly involved in ovarian development. The results of this study enrich the miRNA data resources of ducks and provide a theoretical basis for the follow-up study of miRNA's regulation of ovarian development and poultry reproduction; it also lays a foundation for marker-assisted selection of poultry in the future.

Availability of data and materials
All data generated or analyzed during this study are included in this published article.
Ethics approval and consent to participate All the animals were maintained and studied following the National Institute of Health (NIH) guidelines for care and use of laboratory animals, and all protocols were approved in advance by the Animal Care and Ethics Committee of Guangdong Ocean University of China (No. NXY20160172).

Consent for publication
Not applicable.  Fluorescence quantitative veri cation of differentially expressed miRNA results