Integrated analysis of miRNA and mRNA expression profiles in testes of Duroc and Meishan boars

DOI: https://doi.org/10.21203/rs.2.15327/v2

Abstract

Background: MicroRNAs (miRNAs) are small non-coding RNAs playing vital roles in regulating posttranscriptional gene expression. Elucidating the expression regulation of miRNAs underlying pig testis development will contribute to a better understanding of boar fertility and spermatogenesis.

Results: In this study, miRNA expression profile was investigated in testes of Duroc and Meishan boars at 20, 75, and 270 days of age by high-throughput sequencing. Forty-five differentially expressed miRNAs were identified from testes of Duroc and Meishan boars before and after puberty. Integrated analysis of miRNA and mRNA profiles predicted many miRNA-mRNA pairs. Gene ontology and biological pathway analyses revealed that predicted target genes of ssc-mir-423-5p, ssc-mir-34c, ssc-mir-107, ssc-165 mir-196b-5p, ssc-mir-92a, ssc-mir-320, ssc-mir-10a-5p, and ssc-mir-181b were involved in sexual reproduction, male gamete generation, and spermatogenesis, and GnRH, Wnt, and MAPK signaling pathway. Four significantly differentially expressed miRNAs and their predicted target genes were validated by quantitative real-time polymerase chain reaction, and phospholipase C beta 1 (PLCβ1) gene was verified to be a target of ssc-mir-423-5p.

Conclusions: This study provides an insight into the functional roles of miRNAs in testis development and spermatogenesis and offers useful resources for understanding differences in sexual function development caused by the change in miRNAs expression between Duroc and Meishan boars.

Background

Testis is an important male reproductive and endocrine organ which is a critical tissue for spermatogenesis. Spermatogenesis is a complicated process including mitotic cell division, meiosis, and the process of spermiogenesis [1]. Spermatogenesis is strictly regulated by the expression of stage-specific genes in testis at both transcription and post-transcription levels [2]. Identifying key regulators in testis development and spermatogenesis will provide valuable insights into the mechanism of sexual function development [3].

As a class of endogenous small (~22 nucleotides) noncoding RNAs, miRNAs mediate post-transcriptional gene expression in animals and fine-tune the expression of approximately 30% of all mammalian protein-coding genes [4-6]. In addition, miRNAs regulate gene expressions not only at post-transcriptional levels but also at the transcriptional level by RNA-RNA interactions [7]. A large number of miRNAs were found to have been involved in many biological processes including cell growth and differentiation, embryo development and sperm morphology and mobility [8, 9]. miRNA expression patterns were significantly different between immature and mature mouse testes and miR-449a/b and miR-34b/c function redundantly in the regulation of male germ cell development [10, 11]. Comparative profiling of miRNAs expressed in the newborn, young adult, and aged human epididymides showed that 127 miRNAs were exclusively or preferentially expressed in the newborn epididymis, but only 3 and 2 miRNAs abundantly expressed in the adult and aged epididymides, respectively [12]. miRNAs were also involved in spermatogenesis, and their presence or absence in mature sperm was highly related to aberrant development, function and/or fertility [8, 13, 14]. Adult porcine miRNAs in ovary and testis have been identified and co-expression patterns of X-linked miRNAs in adult porcine gonads were found [15]. miRNA expression patterns between sexually immature (60-day) and mature (180-day) pig testes also have been evaluated and indicated that miRNAs had an important role in regulating spermatogenesis [16]. However, little has been reported about miRNA expression in the testes at various stages of development across pig breeds.

One of our previous studies of testis transcriptional profile revealed numerous differentially expressed genes (DEGs) and important biological pathways were significantly correlated to mammalian reproduction between Meishan and Duroc boars at 20, 75, and 270 days of age [17]. It is interesting to further investigate how miRNAs are involved in regulating sexual function development by fine-turning gene expression [18]. Chinese Meishan pigs are one of the most prolific pig breeds in the world and reach puberty at a relatively younger age (56-84 days) than conventional boars (120-180 days) [19, 20]. Duroc sires are utilized most frequently as a Terminal/Paternal sire in a terminal cross-breeding programe. Our previous work also showed spermatogenesis occurred prior to 75 days in Meishan boars and their spermatogenesis came earlier than Duroc boars, but the number of spermatogonia and Sertoli cells in Meishan boars are less than that in Duroc boars at adulthood [17]. Testis size in Chinese Meishan boars is only half that of conventional boars at maturity. In addition, Meishan boars accumulate Sertoli cells and seminiferous tubules at a more rapid rate compared with white composite boar during the first month after birth [17, 21, 22]. The diameter and number of seminiferous tubules determine the onset of puberty in males [23]. The physiological attributes mentioned above of Chinese Meishan boars make them highly prolific, which render them a valuable animal model for examining the mechanism of sexual function development and sperm production of boars.

Based on miRNA-mRNA pairwise correlations and computational target prediction of miRNA, the miRNA and mRNA expression profiles were integrated to construct miRNA-mRNA regulatory networks which could potentially affect testicular development and spermatogenesis. This study revealed a large number of miRNAs that potentially regulates pig testis development and spermatogenesis and provides a better understanding of differences in sexual function between Meishan and Duroc boars.

Results

Overview of small RNA libraries

In order to identify differentially expressed miRNAs during the process of testicular development of Duroc and Meishan boars, six small RNA (sRNA) libraries of testis tissues of 20-, 75-, 270-day-old Duroc and Meishan boars (D20, D75, D270, M20, M75, and M270) were constructed and sequenced by the Illumina HiSeqTM 2000 platform. In total, 13,335,120, 13,051,493, 13,352,724, 13,606,755, 12,695,970, and 13,721,075 raw reads were generated in D20, D75, D270, M20, M75, and M270, respectively. After removing the low-quality sequences and adaptors, and then discarding the sequences shorter than 18 nt, 13,133,806, 12,892,394, 13,167,240, 13,441,877, 12,561,392, and 13,522,792 clean reads were obtained and used for further analysis (Table 1). A total of 1,078,105 and 878,097 unique sRNA from Duroc and Meishan boar testes were mapped to the porcine reference genome (Sscrofa10.2), respectively (Table 2). Read length distribution analyses of the six sRNA libraries showed that the dominant length of sRNAs was 22 nt, accounting for at least 36.24% of the population (Fig. 1). More 22 nt sRNAs were found in D20, D75, and M20 than in D270, M75, and M270. While few sRNAs with the length of 18 to 19 nt and 25 to 30 nt were detected, and sRNAs with the length of 25 to 30 nt may mainly represent Piwi-interacting RNAs (piRNA). These results are similar to those of previous study of pig [24].

Differentially expressed miRNAs across Duroc and Meishan boars

The expression profiles of known miRNAs from the six samples were analyzed and 36-139 significantly differentially expressed miRNAs (DE miRNAs) (P ≤ 0.05, |log2Ratio|≥1) were filtered in each pairwise comparison (Table 3; Additional file 9: Table S6). For example, the comparisons between the two breeds indicated that the number of DE miRNAs in the 75- and 270- day time points were substantially larger than that in the 20- day time point. These results suggested that significant differences existed during testicular development between Duroc and Meishan boars at the age of 75 and 270 days. More significantly DE miRNAs were detected in Meishan boars at age of 20 to 75 days than those in Duroc boars, which suggested that the testis developed at a faster rate in Meishan boars than in Duroc boars from 20 to 75 days. These findings were consistent with those in previous study of mRNA expression data of matched samples [17].

Venn diagram showed 45 significantly DE miRNAs were filtered from the four pairwise comparisons (D20 vs D270, D75 vs D270, M20 vs M75, and M20 vs M270) (Fig. 2a). The four pairs represented the comparisons before and after puberty since D270, M75, and M270 have reached puberty. Figure 2b showed two main sample branches (D20, D75, and M20 versus D270, M75, and M270), which indicated that the expression pattern of M20 was similar to those of D20 and D75, and that the expression pattern of D270 was similar to those of M75 and M270. The results were consistent with those of previous study of expression pattern of mRNA [17], demonstrating large differences existed in the process of testicular development between Meishan and Duroc boars. It could be concluded that miRNAs were pivotal factors regulating sexual function development.

Integrated analysis between differentially expressed miRNAs and target mRNAs in Duroc and Meishan boars at different stages

The mRNA expression data of six samples from our previous study were used for a pairwise integrated analysis [17]. Through the Trinity de novo assembly method, 20,525 non-redundant genes were obtained from the six samples (Additional file 1: Table S1), 19,310 (94.08%) and 18,241(88.87%) genes were mapped against Kyoto Encyclopedia of Genes and Genomes (KEGG) and (Gene Ontology) GO databases, respectively (Additional file 2: Table S2).

In order to better understand the potential roles of the miRNAs during the process of testicular development, computational target prediction was performed using Targetscan and miRanda. Then we performed integrated analyses of differentially expressed miRNAs and target mRNAs at the expression levels. A large number of correlated miRNA-mRNA pairs were detected in each pairwise comparisons (Fig. 3). The number of miRNA/mRNA-negative pairs between Duroc and Meishan boars at 75-day time point was obviously higher than that at 20- and 270-day time points, and more negative pairs were detected from 20 to 75 days in Meishan boars than those in Duroc boars (Fig. 3a). The previous study demonstrated that Meishan boars attained puberty and their testes generated sperms prior to 75 days earlier than Duroc boars [17]. These findings indicate that miRNAs as negative gene expression regulators significantly control the expression of genes involved in regulating the process of testicular development. Meanwhile, the number of miRNA/mRNA-positive pairs was also compared between Duroc and Meishan boars (Fig. 3b) with a similar tendency found in Duroc and Meishan boars at different ages. These results reveal that the subset of miRNAs may function as enhancers activating the transcription of genes which play an important role during the process of testicular development.

A representative miRNA-mRNA regulatory network of biological pathways was shown in Fig. 4. Eight DE miRNAs (ssc-mir-423-5p, ssc-mir-34c, ssc-mir-107, ssc-mir-196b-5p, ssc-mir-92a, ssc-mir-320, ssc-mir-10a-5p, and ssc-mir-181b) were selected from 45 miRNAs deriving from Fig. 2 serving as functional miRNAs in testicular development of Meishan and Duroc boars according to their annotations and the potential relationship between miRNAs and spermatogenesis and gonad development. GO analysis and KEGG functional annotation of potential target genes of eight DE miRNAs were performed to detect the functional characteristics of miRNAs. A large number of target genes were assigned to the functional categories related to sexual reproduction, male gamete generation, spermatogenesis, sperm development as well as meiosis, indicating that the eight miRNAs were highly involved in spermatogenesis and testis development. Five important pathways including GnRH, Wnt, p53, mTOR, and MAPK signaling pathway related to the regulation of male sexual function were enriched by functional genes including phospholipase C beta 1 (PLCβ1) involved in GnRH and Wnt signaling pathway with PLCβ1 being the target of ssc-mir-423-5p and ssc-mir-34c, serine/threonine/tyrosine interacting protein (STYX) involved in MAPK signaling pathway with STYX being the target of ssc-mir-320, ssc-mir-10a-5p, ssc-mir-92a and ssc-mir-107; cyclin D2 (CCND2), phosphatase and tensin homolog (PTEN), and cyclin B1 (CCNB1) with the latter 3 genes involved in the p53 signaling pathway, and so on. The representative miRNA-mRNA regulatory networks contained biological pathways regulating male sexual function, which illustrated a complex relationship and interaction between the two biomolecular types.

Verification of DE miRNAs and target verification of ssc-mir-423-5p

To evaluate our DE miRNAs library, the expression profiles of 3 DE miRNAs (ssc-mir-181b, ssc-mir-423-5p, ssc-mir-196b-5p), four DEGs from the miRNA-mRNA interaction networks, and ssc-mir-4334-3p from the Venn diagram in Fig. 2a, all of which were highly related to boar sexual function and reproduction, were further analyzed by quantitative real-time PCR (qRT-PCR) with specific primers. As shown in Additional file 3: Figure S1, the results of RNA-seq data and qRT-PCR data were identical. CYLD did not show consistent expression between RNA-seq and qRT-PCR data from Duroc and Meishan boars at age of 20 and 270 days, which was probably caused by the sensitivity of the different methods. In general, the results of qRT-PCR validated the RNA-seq results and demonstrated the reliability of our data.

ssc-mir-423-5p was one of the differentially expressed miRNAs and was selected as a candidate miRNA for analyzing male sexual function. PLCβ1 was predicted to be a target of ssc-mir-423-5p (Fig. 5a). The dual-luciferase reporter assay system analyzed the interaction between ssc-mir-423-5p and PLCβ1 gene. The analysis results indicated that Luciferase activity was significantly suppressed after we co-transfected ssc-mir-423-5p mimic (Additional file 4: Table S3). However, luciferase activity was not significantly changed when we co-transfected ssc-mir-423-5p mimic and pmirGLO- PLCβ1-3´-UTR -mut into Swine Testis (ST) cells (Fig. 5b). Meanwhile, ssc-mir-423-5p inhibitor significantly promoted luciferase activity after we co-transfected ssc-mir-423-5p inhibitor and luciferase reporter vector (pmirGLO- PLCβ1-3´-UTR), and luciferase activity was unchanged after we co-transfected ssc-mir-423-5p inhibitor and pmirGLO-PLCβ1-3´-UTR-mut into ST cells (Fig. 5c). qRT-PCR and western blotting analyses revealed that PLCβ1 mRNA and protein expression levels were significantly reduced after ssc-mir-423-5p mimic was transfected into ST cells, whereas the inhibition of ssc-mir-423-5p increased the expression of PLCβ1 mRNA and protein in ST cells (Fig. 5d-g, Additional file 5: Figure S2, Additional file 6: Figure S3). These results suggested that PLCβ1 was the target gene of ssc-mir-423-5p and the regulation of ssc-mir-423-5p in the process of spermatogenesis may be mediated by PLCβ1.

Discussion

miRNA regulation is critical and effective mechanism underlying the development of testis and spermatogenesis. The sRNA-seq data in this study elucidated the differences in testis development at different stages between Duroc boars and Meishan boars. The present study obtained approximately 99% clean reads from raw reads in each sample, the percentage of high quality reads has reached nearly 100%. Most of the clean reads (80.56–84.88%) identified in this study could match the S. scrofa genome. This result is similar to that found in the study of the pig muscle and ovary transcriptome (78.7%) [25], which indicates that our sRNA-seq data are of high quality. Read length distributions of six libraries demonstrated that 20 to 24 nt represented the length of most sRNAs, of which, 22-nt accounting for the highest percentages. This finding is consistent with the normal size of miRNAs reported in previous study [26]. The present study also indicates that miRNAs are abundantly expressed in testes before puberty, which is consistent with the result that the immature testis had a higher expression level of miRNA than the mature testis [27]. piRNAs are most abundantly expressed in male germ cells, especially during spermatogenesis and these piRNAs are reported to be abundantly expressed in cells at prophase of meiosis I and to get lost at some point before the production of mature sperm [24, 27-29]. piRNAs were abundantly expressed in 75-day-old Meishan boars and 270-day-old Duroc boars, both of which reached puberty and begun to produce sperms [17], while lowly expressed in immature testes (D20, D75, and M20). However, lower expression of piRNAs were also detected in 270-day-old Meishan boars at adulthood than in 75-day-old Meishan boars and 270-day-old Duroc boars, indicating that piRNAs may be the crucial factors causing the differences in sexual function development between Meishan and Duroc boars.

Many mammalian miRNAs play an important role in development and other processes and the expression patterns of miRNAs are tissue-specific or developmental stage-specific [30]. Bioinformatics analyses of miRNAs deriving from the four pairwise comparisons before and after puberty (D20 vs D270, D75 vs D270, M20 vs M75, and M20 vs M270) showed that M20, D20 and D75 were clustered and that M75, M270 and D270 were clustered together (Fig. 2a). These findings are in accordance with our previous mRNA-seq analyses [17]. The expression pattern of miRNAs of M75 was similar to that of M270 and D270, which were in adulthood. These results are consistent with the characteristics of early sexual maturity of Meishan boars reaching their puberty prior to 75 days. A large number of miRNAs were up-regulated in M20, D20, and D75, but down-regulated in M75, M270, and D270. This finding agrees with the analysis result of length distribution that miRNAs were abundantly expressed in testes before puberty. These miRNAs may serve potential roles in regulating testis development, and their specific expressions may induce the differences in sexual maturity and spermatogenesis. A growing number of reports have revealed that miRNAs play important roles in the complex processes of animal testis development and spermatogenesis via their regulations of cell proliferation, apoptosis, and differentiation [31-35].

The miRNAs usually negatively regulate gene expressions by binding specific mRNAs in their 3´-(Untranslated regions) UTRs based on sequence complementation to promote degradation of target mRNAs or to inhibit their translation. Recent studies indicate that miRNAs are also involved in positively-regulating gene transcription by targeting promoter elements [36-38]. Emerging evidence shows that miRNAs are also present in the nucleus in addition to their functions in the cytoplasm. Furthermore, the miRNAs in the nucleus exhibit gene-activation function by activating enhancer RNA (eRNA) expression, altering histone modification [39, 40]. The analysis results of negative or positive miRNA/mRNA pairs were consistent with those of the differential expression analyses of miRNAs and mRNAs, suggesting that Meishan boars have an earlier puberty than Duroc boars and that there was a significant difference in sexual function between two breeds. This study verifies PLCβ1 to be the target gene of ssc-mir-423-5p (Fig. 5) and this gene was also reported to be involved in regulating fertilization rate and embryo development in mice in previous study [41]. In addition, PLCβ1 was reported to have mapped to GnRH [42] and Wnt [43] signaling pathways regulating testis development and spermatogenesis and sexual reproduction. SOX30 predicted to be the target gene of ssc-mir-423-5p was also associated with testis development [44]. Thus, it can be inferred that ssc-mir-423-5p may be a potential regulator in the process of testis development and fertilization. To validate the reliability of the data, we selected important DE miRNAs and DEGs related to reproduction and performed qRT-PCR experiments. The qRT-PCR results matched the sRNA-seq and mRNA-seq data.

Conclusions

In this study, we obtained comprehensive sRNA-seq data and modelled the miRNA-mRNA regulatory networks related to testis development and spermatogenesis through the integrated analysis of differentially expressed miRNAs and target mRNAs between Duroc and Meishan boars. These complex networks identified in our study may participate in regulating boar fertility. In summary, the present study provides an insight into the functional roles of sRNAs in testis development and offers a useful resource for understanding the differences in sexual function and testicular development caused by the change in miRNAs expression between Meishan and Duroc boars.

Methods

Tissue collection 

The experimental animals in this study were raised under the same conditions. Three 20-, 75-, and 270-day-old clinically healthy Duroc and Chinese Meishan boars were selected from the Fine Farm of Hua Zhong Agricultural University. Boars were given general anesthesia (Zoletil 50, Virbac Co., France) before sampling. Testis samples were collected by castration and immediately snap-frozen in liquid nitrogen and stored at -80 °C until RNA extraction. Boars were then fed normally for meat industry.

Illumina small RNA sequencing analysis

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. RNA quality and quantity were measured using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Then the equal amount of RNA from three individuals of the same age and breed were mixed to form 6 RNA pools, named D20, D75, D270, M20, M75, and M270. Sequencing libraries were generated using the TruSeq Small RNA Library Prep Kit (Illumina, San Diego, CA, USA), following the manufacturer’s recommendations. Briefly, the small RNA (sRNA) fragments (18-30 nt) were isolated from RNA pools by polyacrylamide gel electrophoresis (PAGE) and 3´ adaptor (TGGAATTCTCGGGTGCCAAGG) was first ligated to the RNA 3′ ends. Then the 5´end adaptor (GTTCAGAGTTCTACAGTCCGACGATC) was ligated to 5´end of the preparation. T4 RNA Ligase (Takara, Dalian, China) was used in the ligation reaction. The adaptor-ligated sRNA was then converted to cDNA using SuperScript II Reverse Transcriptase (Life Technologies, Carlsbad, CA, USA). The resulting cDNA was amplified on the PCR machine. The purified PCR products were recovered with QIAquick Gel Extraction Kit (Qiagen, Beijing, China) following the manufacturer’s instruction and assessed on an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Each sRNA library was sequenced individually using Illumina HiSeqTM 2000 platform (BGI, Shenzhen, China) and 50 nt single-end reads were generated. To predict the genes targeted by miRNAs, two target prediction algorithms (Targetscan 5.0 and miRanda 3.3a) were used to identify miRNA binding sites. Finally, the data obtained were combined and the overlaps were calculated. Integrated analysis between mRNA [17] and miRNA libraries was performed using miRTrail bioinformatics tool.

Quantitative real-time PCR (qRT-PCR) of miRNA and mRNA expression.

To validate the sequencing data, miRNAs and mRNAs were reversely transcribed using RevertAid First Strand cDNA Synthesis Kit (Thermo, Wuhan, k1622) in accordance with the manufacturer’s instructions. The qRT-PCR of the miRNAs and mRNAs was performed by using a standard UltraSYBR Mixture (CWBIO, Beijing, China) in the Roche LightCyler 480 system (Roche, Mannheinm, Germany) according to the manufacturer’s instructions. U6 and β-actin were used as endogenous control genes of miRNA and mRNA, respectively. The primer sequences used for the qRT-PCR were listed in Additional file 7: Table S4. At least three independent biological replicates were used for each of the miRNAs and mRNA. The qRT-PCR data were analyzed by using the 2-∆∆CT method, as previously described [45].

GO annotation and the KEGG pathway

To further investigate the biological processes and functions which DE miRNAs are involved in through the nine pairwise comparisons, we conducted GO and KEGG pathway analyses. GO enrichment analysis was performed with software Blast2GO, and genes were classified in terms of cellular component, molecular function and biological process using GO annotation. Pathway enrichment analysis was based on KEGG database (http://www.genome.jp/kegg/) and P value was used to determine the threshold of significance in multiple tests and analyses with P value <0.05 defined as significant enrichment.

Plasmid construction

The porcine PLCβ1 3´UTR was amplified with primers PLCβ1-3´UTR-F/PLCβ1-3´UTR-R in Additional file 8: Table S5, and their products were double-digested with Sac I (Takara, Dalian, China) and Xho I (Takara, Dalian, China), and then was cloned into the pmirGLO vector (Promega, Madison, WI, USA) as a wild-type plasmid. And the wild-type plasmid was used as template to construct mutated plasmid with primers PLCβ1-3´UTR-F/PLCβ1-3´UTR-MR in Additional file 8: Table S5, and the primer PLCβ1-3´UTR-R contained the binding site of ssc-mir-423-5p.

Cell culture, transfection and dual-luciferase reporter assays  

Swine testis (ST) cells [46] (ATCC, Shanghai, China) were cultured with Dulbecco’s modified Eagle medium (DMEM) (Hyclone, Logan, UT, USA) containing 10% fetal bovine serum (CLARK, Worcester, MA, USA) at 37 °C in a humidified atmosphere with 5% CO2. The cells were plated into 24 wells. When these cells grew until they reached 70-80% confluent, it was time of transfection. miRNAs were co-transfected with plasmids into the cells using lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). For luciferase assays, wild-type or mutated plasmids at 200 ng together with 3 µL/well of miRNA mimic, NC, miRNA inhibitor or inhibitor NC (Genepharma, Shanghai, China) were transfected into the cells. 24 h after transfection, the luciferase activities were measured with a PerkinElmer 2030 Multilabel Reader (PerkinElmer, Waltham, MA, USA).

Western blotting analysis 

Cellular protein lysates were generated using RIPA Lysis Buffer (Beyotime, Jiangsu, China). Cellular proteins were extracted 72 h after transfection. Proteins were heated to 99 °C for 10 min in 5×SDS buffer, separated by SDS-PAGE, and a Mini Trans-Blot Cell (Bio-Rad, Hercules, CA, USA) was used to transfer protein onto polyvinylidene fluoride membranes (Millipore, Billerica, MA, USA). Then the membranes were blocked with 5% non-fat dried milk and incubated overnight with primary antibodies specific for PLCβ1 (1:1000; ABclonal, USA, A1971). β-actin (1:1000; Wuhan, Hubei, China, GB13001-1) used as a loading control. An Image Quant LAS4000 mini (GE Healthcare Life Sciences, Piscataway, NJ, USA) was used to detect protein expression.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               Graphics and heatmaps 

A heatmap of miRNAs expression was generated using R heatmap package [47]. All the histograms and graphs were generated with GraphPad Prism and Adobe Photoshop CS5 software, respectively. Molecular interaction network was analyzed by Cytoscape_3.3.0.

Statistical analysis 

All experiments are performed at least in triplicate and the corresponding results are presented in terms of the mean ± SD. Two-tailed t-test was used to reveal the differences between two groups. Significant differences were evaluated through an independent-sample t-test. P-value < 0.05 (*P < 0.05; **P < 0.01; ***P < 0.001) was considered to be significant.

Abbreviations

DEGs: differentially expressed genes; DE miRNAs: Differentially expressed miRNAs; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; miRNAs: microRNAs; NC: negative control; PLCβ1: phospholipase C beta 1; piRNA: Piwi-interacting RNA; qRT-PCR: Quantitative real-time PCR; sRNA: small RNA; ST: Swine testis; UTR: Untranslated Region;

Declarations

Acknowledgments

Not Applicable

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files. The transcriptome data is also freely available.

Funding

The project was supported by National Science and Technology Major Project of China (2016ZX08006002), Agricultural Science and Technology Innovation Action of Hubei Province (2018skjcx01), National Natural Science Foundation of China (31201764, 31172179), State Key Development Program for Basic Research of China (2014CB138504), Natural Science Foundation of Hubei Province (2015CFA101), Fundamental Research Funds for the Central Universities (2662015PY192, 2662015PY014, 2662015PY108), Agricultural Science Innovation Foundation of Hubei Province, China. The funding agencies played no further roles in the design of the study, data collection, analysis, and interpretation or in writing the manuscript.

Authors´ contributions

ML, DQX, HSD, and YBY conceived this study; HSD and CFZ performed the experiments. HSD, CFZ, XBY, and TS conducted the data analysis and prepared figures and tables. ML, DQX, HSD and YBY wrote the manuscript. All of the authors reviewed and approved the manuscript.

Ethics approval and consent to participate

All animal experiments were conducted according to the Guidelines of the Care and Use of Laboratory Animals in China. The verbal informed consent was obtained from all participants (veterinary officers and farmers) owing to distance and traffic. All animals did not suffer unnecessarily at any stages of this experiment. The study was also approved by the Ethics Committee of Experimental Animals of Huazhong Agricultural University.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References

  1. de Kretser DM, Loveland KL, Meinhardt A, Simorangkir D, Wreford N. Spermatogenesis. Hum Reprod. 1998;13(suppl 1):1–8.
  2. Weng B, Ran M, Chen B, He C, Dong L, Peng F. Genome-wide analysis of long non-coding RNAs and their role in postnatal porcine testis development. Genomics. 2017;109(5-6):446–456.
  3. Zhang X, Li L, Jiang H, Ma J, Li J, Chen J. Identification and differential expression of microRNAs in testis and ovary of Amur sturgeon (Acipenser schrenckii). Gene. 2018; 658:36–46.
  4. Sevignani C, Calin GA, Siracusa LD, Croce CM. Mammalian microRNAs: a small world for fine-tuning gene expression. Mamm genome. 2006;17(3):189–202.
  5. Kasimanickam V, Kastelic J. MicroRNA in sperm from Duroc, Landrace and Yorkshire boars. Sci Rep. 2016;6:32954.
  6. Wu S, Ren X, Li Y, Guo W, Lei X, Yao J, et al. Effect of dietary Astragalus Polysaccharide supplements on testicular miRNA expression profiles and enzymatic changes of breeder cocks. Sci Rep. 2017;7:38864.
  7. Guil S, Esteller M. RNA-RNA interactions in gene regulation: the coding and noncoding players. Trends biochem sci. 2015;40(5):248–256.
  8. Curry E, Safranski TJ, Pratt SL. Differential expression of porcine sperm microRNAs and their association with sperm morphology and motility. Theriogenology. 2011;76(8):1532–1539.
  9. Bushati N, Cohen SM. microRNA functions. Annu Rev Cell Dev Biol. 2007;23:175–205.
  10. Yan N, Lu Y, Sun H, Tao D, Zhang S, Liu W, et al. A microarray for microRNA profiling in mouse testis tissues. Reproduction. 2007;134(1):73–79.
  11. Bao J, Li D, Wang L, Wu J, Hu Y, Wang Z, et al. MicroRNA-449 and microRNA-34b/c function redundantly in murine testes by targeting E2F transcription factor-retinoblastoma protein (E2F-pRb) pathway. J Biol Chem. 2012; 287(26):21686–21698.
  12. Zhang J, Liu Q, Zhang W, Li J, Li Z, Tang Z, et al. Comparative profiling of genes and miRNAs expressed in the newborn, young adult, and aged human epididymides. Acta Biochim Biophys Sin. 2010;42(2):145–153.
  13. Amanai M, Brahmajosyula M, Perry ACF. A Restricted Role for Sperm-Borne MicroRNAs in Mammalian Fertilization. Biol Reprod. 2001;75(6):877–884.
  14. Curry E, Ellis SE, Pratt SL. Detection of porcine sperm microRNAs using a heterologous microRNA microarray and reverse transcriptase polymerase chain reaction. Mol Reprod Dev. 2009;76(3):218–219.
  15. Li M, Liu Y, Wang T, Guan J, Luo Z, Chen H, et al. Repertoire of porcine microRNAs in adult ovary and testis by deep sequencing. Int J Biol Sci. 2011;7(7):1045–1055.
  16. Luo L, Ye L, Liu G, Shao G, Zheng R, Ren Z, et al. Microarray-based approach identifies differentially expressed microRNAs in porcine sexually immature and mature testes. PloS one. 2010;5(8):e11744.
  17. Ding H, Luo Y, Liu M, Huang J, Xu D. Histological and transcriptome analyses of testes from Duroc and Meishan boars. Sci Rep. 2016;6:20758.
  18. Qiu W, Zhu Y, Wu Y, Yuan C, Chen K, Li M. Identification and expression analysis of microRNAs in medaka gonads. Gene. 2018;646:210–216.
  19. Lunstra DD, Ford JJ, Klindt J, Wise TH. Physiology of the Meishan boar. J Reprod Fertil Suppl. 1997;52(1):181–193.
  20. Ford JJ, Wise TH. Sertoli cell differentiation in pubertal boars. J Anim Sci. 2009;87(8):2536–2543.
  21. Borg KE, Lunstra DD, Christenson RK. Semen characteristics, testicular size, and reproductive hormone concentrations in mature duroc, meishan, fengjing, and minzhu boars. Biol Reprod. 1993; 49(3):515–521.
  22. McCoard SA, Wise TH, Lunstra DD, Ford JJ. Stereological evaluation of Sertoli cell ontogeny during fetal and neonatal life in two diverse breeds of swine. J Endocrinol. 2003;178(3):395–403.
  23. Ford JJ, McCoard SA, Wise TH, Lunstra DD, Rohrer GA. Genetic variation in sperm production. Soc Reprod Fertil Suppl. 2006;62:99–112.
  24. Lau NC, Seto AG, Kim J, Kuramochi-Miyagawa S, Nakano T, Bartel DP, et al. Characterization of the piRNA complex from rat testes. Science. 2006;313(5785):363–367.
  25. Nie Q, Fang M, Jia X, Zhang W, Zhou X, He X, et al. Analysis of muscle and ovary transcriptome of Sus scrofa: assembly, annotation and marker discovery. DNA res. 2011;18(5):343–351.
  26. Ma J, Li M, Li X. MicroRNA Transcriptomes Relate Intermuscular Adipose Tissue to Metabolic Risk. Int J Mol Sci. 2013;14(4):8611.
  27. Lian C, Sun B, Niu S, Yang R, Liu B, Lu C, et al. A comparative profile of the microRNA transcriptome in immature and mature porcine testes using Solexa deep sequencing. FEBS J. 2012;279(6):964–975.
  28. Hannon GJ. A germline-specific class of small RNAs binds mammalian Piwi proteins. Nature 2006, 442(7099):199–202.
  29. Hirano T, Iwasaki YW, Lin ZY, Imamura M, Seki NM, Sasaki E, et al. Small RNA profiling and characterization of piRNA clusters in the adult testes of the common marmoset, a model primate. Rna. 2014; 20(8):1223–1237.
  30. Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of Mammalian MicroRNA Targets. Cell. 2003;115(7):787–798.
  31. Hu P, Guan K, Feng Y, Ma C, Song H, Li Y, et al. miR-638 Inhibits immature Sertoli cell growth by indirectly inactivating PI3K/AKT pathway via SPAG1 gene. Cell Cycle. 2017;16(23):2290-2300.
  32. Ran M, Weng B, Cao R, Li Z, Peng F, Luo H, et al. miR-26a inhibits proliferation and promotes apoptosis in porcine immature Sertoli cells by targeting the PAK2 gene. Reprod Domest Anim. 2018;53(6): 1375-1385.
  33. Gao H, Ran M, Luo H, Weng B, Tang X, Chen Y, et al. miR-499 promotes immature porcine Sertoli cell growth by the PI3K/AKT pathway by targeting the PTEN gene. Reproduction. 2019; pii: REP-19-0303.R1.
  34. Ran M, Weng B, Cao R, Peng F, Luo H, Gao H, et al. miR-34c inhibits proliferation and enhances apoptosis in immature porcine Sertoli cells by targeting the SMAD7 gene. J Integr Agr. 2019;18(2): 449-459.
  35. Ran M, Luo H, Gao H, Tang X, Chen Y, Zeng X, et al. miR-362 knock-down promotes proliferation and inhibits apoptosis in porcine immature Sertoli cells by targeting the RMI1 gene. Reprod Domest Anim. 2020; doi: 10.1111/rda.13626.
  36. Xiao M, Li J, Li W, Wang Y, Wu FZ, Xi YP, et al. MicroRNAs activate gene transcription epigenetically as an enhancer trigger. Rna Biol. 2017;14(10):1326–1334.
  37. Majid S, Dar AA, Saini S, Yamamura S, Hirata H, Tanaka Y, et al. MicroRNA-205-directed transcriptional activation of tumor suppressor genes in prostate cancer. Cancer. 2010;116(24):5637–5649.
  38. Matsui M, Chu Y, Zhang H, Gagnon KT, Shaikh S, Kuchimanchi S, et al. Promoter RNA links transcriptional regulation of inflammatory pathway genes. Nucleic Acids Res. 2013;41(22):10086–10109.
  39. Place RF, Li LC, Pookot D, Noonan EJ, Dahiya R. MicroRNA-373 induces expression of genes with complementary promoter sequences. Proc Natl Acad Sci USA. 2008;105(5):1608–1613.
  40. Li G, Ruan X, Auerbach RK, Sandhu KS, Zheng M, Wang P, et al. Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation. Cell. 2012;148(1-2):84–98.
  41. Ballester M, Molist J, Lopez-Bejar M, Sánchez A, Santaló J, Folch JM, et al. Disruption of the mouse phospholipase C-beta1 gene in a beta-lactoglobulin transgenic line affects viability, growth, and fertility in mice. Gene. 2004;341(1):279–289.
  42. Naor Z, Benard O, Seger R。 Activation of MAPK Cascades by G-protein-coupled Receptors: The Case of Gonadotropin-releasing Hormone Receptor. Trends Endocrinol Metab. 2000;11(3):91–99.
  43. Dong WL, Tan FQ, Yang WX. Wnt signaling in testis development: Unnecessary or essential? Gene. 2015;565(2):155–165.
  44. Han F, Dong Y, Liu W, Ma X, Shi R, Chen H, et al. Epigenetic regulation of sox30 is associated with testis development in mice. PLoS One. 2014;9(5):e97203.
  45. Ren Q, Xu ZL, Wang XW, Zhao XF, Wang JX. Clip domain serine protease and its homolog respond to Vibrio challenge in Chinese white shrimp, Fenneropenaeus chinensis. Fish Shellfish Immun. 2009;26(5):787–798.
  46. Ma C, Song H, Guan K, Zhou J, Xia X, Li F. Characterization of swine testicular cell line as immature porcine Sertoli cell line. In Vitro Cell Dev-An. 2016; 52(4):427–433.
  47. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010.11(10):R106.

Tables

Table 1. Quality analyses of small RNA-seq data

Sample

Total reads

high_quality

Smaller_than_18nt

PolyA

Clean reads

D20

13,335,120

13,268,836(100%)

79,018(0.60%)

15

13,133,806(98.98%)

D75

13,051,493

12,987,508(100%)

41,877(0.32%)

21

12,892,394(99.27%)

D270

13,352,724

13,289,346(100%)

70,860(0.53%)

16

13,167,240(99.08%)

M20

13,606,755

13,537,863(100%)

42,915(0.32%)

22

13,441,877(99.29%)

M75

12,695,970

12,642,443(100%)

23,239(0.18%)

28

12,561,392(99.36%)

M270

13,721,075

13,658,694(100%)

83,811(0.61%)

6

13,522,792(99.01%)

 

Table 2. The reads mapping to reference genome from small RNA-seq data

Sample

Unique sRNAs

Total sRNAs

Total

mapped

Total 

mapped

D20

402,022

269,359 (67%)

13,133,806

10,689,788 (81.39%)

D75

433,787

293,153 (67.58%)

12,892,394

10,942,820 (84.88%)

D270

787,613

515,593 (64.46%)

13,167,240

10,849,401 (82.4%)

M20

417,327

269,558 (64.59%)

13,441,877

10,828,754 (80.56%)

M75

857,893

554,629 (64.65%)

12,561,392

10,279,770 (81.84%)

M270

117,259

53,910 (45.98%)

13,522,792

10,279,770 (81.84%)

 

Table 3. Analysis of differentially expressed miRNAs

Sample

Total

Up regulated miRNA

Down regulated miRNA

D20-vs-M20

39

22

17

D75-vs-M75

100

21

79

D270-vs-M270

120

104

16

D20-vs-D75

36

13

23

D20-vs-D270

138

19

119

D75-vs-D270

139

22

117

M20-vs-M75

111

16

95

M20-vs-M270

125

64

61

M75-vs-M270

116

93

23

 

Additional File Legends

Additional file 1: Table S1 The number of genes detected in the six samples.

Additional file 2: Table S2 The number of genes involved in KEGG pathway and GO term.

Additional file 3: Figure S1 Expression analyses of 4 DE miRNAs and 4 DEGs using qRT-PCR and RNA-Seq. Analyses of relative expression levels of miRNAs and genes in Duroc (black bars) and Meishan (gray bars) boars in the left; RPKM of miRNAs and genes in Duroc (blue bars) and Meishan (purple bars) boars a in the right. The X-axis indicates age of samples. The left Y-axis and right Y-axis show the relative expression level and RPKM of miRNAs and genes, respectively. a ssc-mir-181b, b ssc-mir-423-5p, c ssc-mir-196b-5p, d ssc-mir-4334-3p, e CYLD, f CDYL, g PLCβ1, h SOX30. **P< 0.01, *P< 0.05. DE miRNAs, differentially expressed miRNAs; DEGs, differentially expressed genes; qRT-PCR, quantitative real-time PCR; RNA-Seq, RNA sequencing; RPKM, Reads Per Kilobase per Million mapped reads; CYLD, Cylindromatosis; CDYL, chromodomain Y like; PLCβ1, phospholipase C beta 1; SOX30, SRY-box containing gene 30.

Additional file 4: Table S3 Synthesis of ssc-mir-423-5p sequences. NC, negative control.

Additional file 5: Figure S2 The expression levels of ssc-mir-423-5p were analyzed by qRT-PCR. a ssc-mir-423-5p expression level after transfecting ssc-mir-423-5p mimic, mimic NC. b ssc-mir-423-5p expression level after transfecting ssc-mir-423-5p inhibitor or inhibitor NC into ST cells. NC, negative control; ST, swine testis

Additional file 6: Figure S3 Full image of Western Blotting. a PLCβ1 protein expression levels at 72 h after ST cells were transfected with ssc-mir-423-5p mimic and mimic NC. b PLCβ1 protein expression levels at 72 h after ST cells were transfected with ssc-mir-423-5p inhibitor and inhibitor NC.

Additional file 7: Table S4 Primers of genes and miRNAs used for quantitative real-time PCR.

Additional file 8: Table S5 Primers for amplification of porcine PLCβ1 3´UTR.

Additional file 9: Table S6 The identified miRNA information