Sequence analysis of small RNA in Leizhou black duck ovary
The specific sequencing results are shown in Table 2. It can be seen from the table that the total reads obtained in the LG group were 11,785,116 and 12,651,078, among which high-quality sequences were 11,609,326 and 12,469,633, respectively accounting for 98.5084% and 98.5658% of the original small RNA sequences. The total reads in the HG group were 13,930,857, 12,573,457, and 12,223,767 original small RNA sequences, of which high-quality sequences were 13,756,770, 12,390,224, and 12,061,537, accounting for 98.7503%, 98.5427% and 98.6728% respectively of the original small RNA sequences. The library is of high quality and meets the test requirements. After a series of screenings, the small RNA reads obtained were: LG group; 11,262,740 and 12,166,225 accounting for 97.0146% and 97.5668% of high-quality sequence respectively and HG group; 13,500,147, 12,117,644, and 11,713,803 accounting for 98.1346%, 97.8000% and 97.1170% of high-quality sequence respectively.
Table 2. High and low yield Leizhou black duck ovarian tissue sequencing data statistics table
Id
|
L1
|
L2
|
H1
|
H2
|
H3
|
Clean reads
|
11785116
(100%)
|
12651078
(100%)
|
13930857
(100%)
|
12573457
(100%)
|
12223767
(100%)
|
High quality
|
11609326
(98.5084%)
|
12469633
(98.5658%)
|
13756770
(98.7503%)
|
12390224
(98.5427%)
|
12061537
(98.6728%)
|
3'adapter_null
|
33664
(0.2900%)
|
31112
(0.2495%)
|
31617
(0.2298%)
|
27912
(0.2253%)
|
38702
(0.3209%)
|
Insert null
|
17892
(0.1541%)
|
11991
(0.0962%)
|
15480
(0.1125%)
|
11316
(0.0913%)
|
12074
(0.1001%)
|
5'adapter_contaminants
|
4371
(0.0377%)
|
4986
(0.0400%)
|
8450
(0.0614%)
|
3364
(0.0272%)
|
4286
(0.0355%)
|
Smaller than_18nt
|
140535
(1.2105%)
|
107521
(0.8623%)
|
80093
(0.5822%)
|
102332
(0.8259%)
|
118281
(0.9806%)
|
Poly A
|
57
(0.0005%)
|
49
(0.0004%)
|
32
(0.0002%)
|
50
(0.0004%)
|
42
(0.0003%)
|
low cutoff
|
150067
(1.2926%)
|
147749
(1.1849%)
|
120951
(0.8792%)
|
127606
(1.0299%)
|
174349
(1.4455%)
|
Clean tags
|
11262740
(97.0146%)
|
12166225
(97.5668%)
|
13500147
(98.1346%)
|
12117644
(97.8000%)
|
11713803
(97.1170%)
|
Small RNA sequence length distribution in Leizhou black duck ovary
The detailed information of the statistical analysis of small RNA length on the obtained high-quality sequences is shown in Figure 1. In L1, the 22 nt sequence length accounted for the largest proportion, 47.49%, followed by 21 nt (22.58%) and 23 nt (14.33%). In L2, the 22 nt sequence length accounted for 47.82%, followed by 21 nt (21.96%) and 23 nt (16.09%). In H1, the largest sequence length was 22 nt, and the specific proportion was 48.89%, followed by 21 nt (21.59%) and 23 nt (16.06%). In H2, the 22 nt sequence length accounted for 47.34%, followed by 21 nt (22.66%) and 23 nt (15.43%). In H3, the sequence length at 22 nt accounted for 47.45%, followed by 21 nt (21.77%) and 23 nt (15.99%). Combining the results of the 5 groups, it is found that the total readings and the length with the largest proportion was 22 nt.
Identification of conserved and novel microRNAs
After miRNA filtering, 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).
Table 3. Known miRNAs and newly predicted miRNAs with the highest expression in ovarian tissues (top 8)
miRNA
|
known miRNA
|
miRNA
|
newly predicted miRNA
|
total
|
TPM value
|
total
|
TPM value
|
miR-143-y
|
16 315 981
|
1 391 037
|
novel-m0121-3p
|
17 664
|
1 499.103
|
miR-26-x
|
8 805 902
|
747 761.8
|
novel-m0193-5p
|
2 223
|
189.3829
|
miR-125-x
|
3 973 385
|
339 412.5
|
novel-m0180-5p
|
1 891
|
160.6285
|
miR-99-x
|
3 458 884
|
295 616.6
|
novel-m0181-5p
|
1 891
|
160.6285
|
miR-199-x
|
2 617 097
|
222 298.8
|
novel-m0009-3p
|
1 742
|
153.7457
|
miR-100-x
|
2 560 542
|
218 603.3
|
novel-m0119-3p
|
1 570
|
131.7257
|
miR-10-x
|
2 089 460
|
177 720.7
|
novel-m0015-5p
|
643
|
54.9524
|
let-7-x
|
2 004 649
|
170 971.1
|
novel-m0294-5p
|
643
|
54.9524
|
Analysis of differential expressed miRNAs
In all miRNA expression profiles, 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 specific 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 significant 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 significant (p-value = 0.000139).
Table 4. Differentially expressed miRNA details
miRNA
|
Log2(fc)
|
up/down
|
(LG/HG)P-value
|
Let-7-z
|
-2.09752
|
down
|
0.042916
|
miR-10014-x
|
-5.77109
|
down
|
0.049729
|
miR-1307-y
|
-6.5809
|
down
|
0.008262
|
miR-1788-y
|
-3.47557
|
down
|
0.044594
|
miR-192-x
|
-3.28617
|
down
|
0.000139
|
miR-194-x
|
-2.70287
|
down
|
0.003833
|
miR-205-x
|
-1.66821
|
down
|
0.011228
|
miR-215-x
|
-1.88212
|
down
|
0.038472
|
miR-25-y
|
-2.67909
|
down
|
0.035621
|
miR-28-y
|
-5.61368
|
down
|
0.040224
|
miR-34-x
|
-3.27515
|
down
|
0.005329
|
miR-34-y
|
-2.74243
|
down
|
0.041979
|
miR-361-x
|
-3.76559
|
down
|
0.02982
|
miR-423-y
|
-3.96174
|
down
|
0.031725
|
miR-449-x
|
-2.82011
|
down
|
0.018527
|
miR-462-x
|
-2.82117
|
down
|
0.042008
|
miR-484-x
|
-4.504
|
down
|
0.043299
|
miR-722-y
|
-4.50772
|
down
|
0.00758
|
miR-7475-x
|
-5.63387
|
down
|
0.044794
|
miR-9344-y
|
-2.76445
|
down
|
0.031751
|
novel-m0007-3p
|
-2.70013
|
down
|
0.033385
|
novel-m0248-3p
|
-2.70013
|
down
|
0.032746
|
miR-1175-y
|
9.719355
|
up
|
0.011029
|
miR-12-z
|
4.808272
|
up
|
0.041664
|
miR-305-x
|
9.266302
|
up
|
0.016344
|
miR-306-x
|
7.099856
|
up
|
0.044834
|
novel-m0304-3p
|
6.873403
|
up
|
0.0087
|
novel-m0305-3p
|
6.873403
|
up
|
0.008807
|
novel-m0306-3p
|
6.873403
|
up
|
0.008877
|
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 significantly 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 enriched GO terms; the deeper the color, the higher the level of enrichment in DAG dates. In terms of biological processes, target genes were significantly enriched in cellular processes, single-organism processes, metabolic processes, and reproductive processes. In terms of cellular components, target genes were significantly enriched in cells, cell part, organelle, membrane, and macromolecular complexes. In terms of molecular functions, target genes were significantly 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 significantly enriched in axon guidance (360 target genes predicted), ErbB signaling pathway (209 target genes predicted), TRP channel Inflammatory mediator regulation of TRP channels (276 target genes predicted), Proteoglycans in cancer (400 target genes predicted) and other signaling pathways. Significantly 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).
Table 5. Pathway analysis of differentially expressed miRNA target genes KEGG (top 10)
KEGG_A_class
|
KEGG_B_class
|
Pathway
|
LG vs HG
|
P value
|
Organismal Systems
|
Development
|
Axon guidance
|
360
|
1.497404E-22
|
Environmental Information Processing
|
Signal transduction
|
ErbB signaling pathway
|
209
|
4.030866E-22
|
Organismal Systems
|
Sensory system
|
Inflammatory mediator regulation of TRP channels
|
276
|
1.300796E-20
|
Human Diseases
|
Cancers
|
Proteoglycans in cancer
|
400
|
2.932151E-20
|
Organismal Systems
|
Nervous system
|
Cholinergic synapse
|
273
|
1.528637E-19
|
Organismal Systems
|
Endocrine system
|
Insulin signaling pathway
|
227
|
6.608405E-18
|
Organismal Systems
|
Nervous system
|
Neurotrophin signaling pathway
|
232
|
3.596927E-16
|
Organismal Systems
|
Endocrine system
|
Oxytocin signaling pathway
|
360
|
4.312446E-15
|
Environmental Information Processing
|
Signal transduction
|
Phospholipase D signaling pathway
|
299
|
3.90E-14
|
Metabolism
|
Glycan biosynthesis and metabolism
|
Glycosphingolipid biosynthesis - lacto and neolacto series
|
87
|
4.904162E-14
|
Verification of the differential expressed miRNAs
To verify the accuracy of sequencing data, we randomly selected 12 miRNAs with known differential expression for verification. 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 log2 (LG/HG), the results showed that had the same variation trend compared with the results (Figure 6).