Small RNA sequencing of cryopreserved semen from single bull revealed altered miRNAs and piRNAs expression between High- and Low-motile sperm populations

Background Small RNAs present in bovine ejaculate can be linked to sperm abnormalities and fertility disorders. At present, quality parameters routinely used in semen evaluation are not fully reliable to predict bull fertility. In order to provide additional quality measurements for cryopreserved semen used for breeding, a method based on deep sequencing of sperm microRNA (miRNA) and Piwi-interacting RNA (piRNA) from individual bulls was developed. To validate our method, two populations of spermatozoa isolated from high and low motile fractions separated by Percoll were sequenced, and their small RNAs content characterized. Results Sperm cells from frozen thawed semen samples of 4 bulls were successfully separated in two fractions. We identified 83 miRNAs and 79 putative piRNAs clusters that were differentially expressed in both fractions. Gene pathways targeted by 40 known differentially expressed miRNAs were related to apoptosis. Dysregulation of miR-17-5p, miR-26a-5p, miR-486-5p, miR-122-5p, miR-184 and miR-20a-5p was found to target three pathways (PTEN, PI3K/AKT and STAT). Conclusions Small RNAs sequencing data obtained from single bulls are consistent with previous findings. Specific miRNAs are differentially represented in low versus high motile sperm, suggesting an alteration of cell functions and increased germ cell apoptosis in the low motile fraction. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-3394-7) contains supplementary material, which is available to authorized users.


Background
Reproductive success is crucial for species' survival. Infertility is a disorder affecting humans as well as other animals. Concerning these, latter infertility is a major cause of economic losses and a major limitation to the achievement of optimum efficiency in the livestock production system. The causes of infertility can be numerous and complexes. In human, infertility is prevalently due to anatomical problems and endocrine disorders causing low sperm counts and poor sperm quality, and in part to genetic disorders [1]. In cattle, a number of bulls considered of high-merit based on their spermatozoa motility and morphology were reported to be unable to produce successful full-term pregnancies, according to extensive fertility data and progeny records [2,3], suggesting that molecular defects affect the ability of spermatozoa to fertilize and contribute to normal embryo development [4][5][6]. Individual bulls differ in their ability to fertilize oocytes in vitro depending on different sperm traits, like motility, membrane and acrosome integrity, and the ability to penetrate oocytes [7]. Cryopreserved semen is used worldwide in farm animal husbandry and for animal genetic resources conservation. Several advanced technologies can be used to examine quality of spermatozoa -as Computer-Assisted Semen Analysis (CASA) and flow cytometry (FCM) -which can provide accurate and unbiased evaluation of sperm functions. It is generally accepted that sperm motility is a determining factor in normal male fertility because of its essential role in reaching the site of fertilization [8], as a consequence, the evaluation of sperm motility is useful for the diagnosis and treatment of low fertility and infertility [9]. Despite their relevance, the molecular mechanisms controlling sperm motility are still partially unknown. The integration of several tests, from standard procedures for the evaluation of sperm motility and viability, to sperm molecular investigation, is a promising approach to achieve a better understanding of sperm functions as well as to evaluate semen quality and predict bull fertility.
During fertilization, besides the paternal genome, spermatozoa transport coding and non coding RNAs into the oocyte. Mammalian sperm contains an array of RNAs including messenger RNAs (mRNAs), ribosomal RNAs (rRNAs) and small RNAs (sRNAs), largely representing remnant transcripts produced during spermatogenesis [10][11][12]. RNA-Seq characterization of bovine spermatozoa revealed the presence of degraded and full-length nuclear-encoded transcripts involved in capacitation and fertilization, suggesting that RNA could be translated after spermatogenesis and potentially contribute to capacitation and early embryogenesis [13]. Furthermore, sperm transcripts retain information of the past events of spermatogenesis and probably contribute to egg fertilization and development. Comparisons between sperm from fertile and infertile males in different species indicate that sperm transcripts may have diagnostic value, and suggest a relationship between sperm transcripts composition and proper sperm functions [8,[14][15][16][17].
sRNAs are a class of short non-coding RNAs including different types of RNAs (i.e. microRNA (miRNA) and Piwi-interacting RNA (piRNA)), that play an essential regulatory role in spermatogenesis, such as maintenance and transposon silencing. piRNAs are known to be important to maintain fertility, as confirmed by the defects in fertility observed in mutants lacking Piwi in C. elegans [18], Danio rerio [19] and Mus musculus [20]. miRNAs were found to regulate spermatogonial stem cell (SSCs) renewal at the post-transcriptional level via targeting specific genes [21]. The testicular expressed miRNAs were reported to change depending on the stage of spermatogenesis [22,23]. miRNAs participate in the control of many functions, such as maintenance of spermatogonial stem cells (SSCs) status, regulation of SSCs differentiation, meitoic and post-meiotic processing and spermiogenesis [24]. Dysregulation in miRNAs' expression patterns is severely affected in different types of reproduction abnormalities [25][26][27]. Sperm miRNA profiling alteration was detected in bulls with high vs low fertility level, indicating a possible role of miRNAs in male infertility [28].
Since the first genome-wide miRNA and piRNA profiling in human testis was reported [29], the Next Generation Sequencing (NGS) technology was adopted to detect sRNAs dysregulation associated to sperm characteristic alterations. Recently, the bull sperm microRNAome was found to be altered in the "fescue toxicosis" syndrome, a disease related to consumption of alkaloids contaminated feed, which has negative effects on growth and reproduction in animals [30]. However, due to the low yields in miRNA recovery from frozen semen, analyses were conducted on RNAs from several pooled individuals.
Here, we propose the first integrated approach to compare miRNA and piRNA expression between high and low motility sperm populations isolated after Percoll gradient from cryopreserved spermatozoa collected from single bulls. Deep sequencing information from single animal was achieved to explore how miRNA and piRNA expression variations can potentially affect bovine sperm characteristics, such as motility and kinetic parameters. The development of a reliable method for small RNA profiling in bovine sperm isolated from frozen thawed sperm through NGS could be an important step in deciphering the contents of miRNA and piRNA sequences in animals that are well characterized for different traits such as fertility.

Isolation of spermatozoa through Percoll gradient
Frozen semen straws from four mature progeny tested Holstein bulls with satisfactory semen quality were obtained from an Artificial Insemination AI center (INSEME, Zorlesco, Lodi, Italy).
For each bull 12 frozen semen doses (0.5 mL, 20x10 6 cells per dose) were simultaneously thawed in a water bath at 37°C for 20 seconds and pooled. The pool (6 mL) was split in 3 aliquots of 2 mL that were overlaid on a dual-layer (90-45%) discontinuous Percoll gradient (Sigma-Aldrich, St. Louis, USA) in three 15 ml conical tubes and centrifuged at 700 × g for 30 min at 20°C. The Percoll layers were prepared by diluting Percoll solution as previously described [31]. The Percoll gradient is a colloidal suspension of silica particles coated with polyvinylpyrrolidone (PVP). By using two discontinuous layers (45% and 90%) by centrifugation it is possible to obtain a different sedimentation according to sperm motion. The two fractions obtained (High Motile = HM and Low Motile = LM) from each of the three tubes (replicates) were washed in Tyrode's albumin lactate pyruvate (TALP) buffer at 700 × g for 10 min at 20°C; the obtained pellets were re-suspended in 150 μl of TALP. For each bull an aliquot of semen of the High Motile and Low Motile fractions was evaluated immediately after Percoll density gradient centrifugation. Three technical replicates per bull were evaluated for sperm kinetic parameters by CASA, and sperm viability and acrosomal status by flow cytometer in both fractions. Aliquots from each replicate were kept at −80°C until RNA extraction (approximately 1 month later).

Flow cytometry analysis
Measurements were performed on a Guava Easycyte TM 5HT microcapillary flow cytometer (Merck KGaA Darmstadt) with the CytoSoft™ and IMV EasySoft software for semen analysis (IMV Technologies, France). The fluorescent probes were excited by an Argon ion blue laser (488 nm). A forward and side-scatter gate were used to separate sperm cells from debris. Non sperm events were excluded from further analysis. Detection of fluorescence was set with three photomultiplier tubes (green: 525/30 nm, orange/yellow: 586/26 nm, and red: 690/ 50 nm). Compensation for spectra overlap between fluorochromes was set (http://www.drmr.com/compensation). Calibration was carried out using standard beads with the Guava Easy Check Kit (Guava Technologies, Inc., Millipore). Acquisitions were performed using the CytoSoft™ software. A total of 5000 events per sample were analyzed with a flow rate of 200 cells/s. The assessment of sperm viability and acrosome integrity was performed by using EasyKit 5 (IMV Technologies, France). The percentage of cells with disrupted acrosome within viable or dead sperm fractions was measured. Each well of the ready-to-use 96well plate was filled with 200 μL of Embryo Holding solution (IMV Technologies, France), 40.000 sperm cells were added and incubated for 45 min at 37°C in the dark.
Spermatozoa with disrupted acrosomes were labeled with a green probe, dead spermatozoa with damaged plasma membrane were labeled with a red fluorochrome, consequently the percentages of alive and dead sperm fractions with intact or damaged acrosomal membrane were computed.

RNA extraction
For each bull, HM and LM sperm fractions obtained from three technical replicates (equivalent to approximately four frozen semen doses each) were used for RNA isolation. RNA was extracted using TRIzol ® (Invitrogen, Carlsbad, CA) according to Govindaraju et al. [28], with some modifications. Briefly, 400 μl of TRIzol were added into each sperm cell pellet and then homogenized at high speed for 30 s. Glycogen (3 μl of 20 mg/ml) was added to the tubes and another 400 μl of TRIzol ® were then added, mixed and incubated for 15 min at 65°C. Total RNA was then purified with the NucleoSpin ® miRNA kit (Macherey-Nagel, Germany), following the protocol in combination with TRIzol ® lysis with small and large RNA in one fraction (total RNA). RNA concentration and quality were determined by Agilent 2100 Bioanalyzer (Santa Clara, CA). The isolated RNAs were stored at −80°C until use.

Library preparation and sequencing
Six sperm RNA samples, representing three technical replicates for both HM and LM fractions, were obtained from each single bull. RNA extraction from semen straws typically resulted in few picograms of RNA: a quantity not compatible for single small RNA library sequencing. Therefore pool of sample has been usually used for semen small RNA sequencing. In order to avoid pooling samples, our approach provide a library preparation from each single RNA sample with proper index. Libraries from single samples were then combined, approximately fifteen-fold concentrated in volume and isolated. Small RNA libraries were generated using the Illumina Truseq Small RNA Preparation kit according to manufacturer's instructions with the following modifications: before size selection, libraries were pooled together and added with Agencourt ® AMPure ® XP (Beckman, Coulter, Brea, CA) (1 Vol. sample: 1.8 Vol. beads). Libraries were eluted in 1/15 volume of the initial pool solution (15X libraries pool). The libraries pool was purified on a Pippin Prep system (Sage Science, MA, USA) to recover the 125 to 167 nt fraction containing mature miRNAs (Additional file 1). The quality and yield after sample preparation was measured with an Agilent 2200 Tape Station, High Sensitivity D1000. Libraries were sequenced on a single lane of Illumina Hiseq 2000 (San Diego, CA).

piRNA analysis
Preliminary quality control of raw reads was carried out with FastQC (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/). Illumina raw sequences were then trimmed with Trimmomatic [32] to remove primers, Illumina adapters and low quality regions and sequences. A minimum average base quality of 15 over a 4 bases sliding window and a minimum length of 12 bases of the trimmed sequence were used as thresholds.
Small RNA sequences ranging from 26 to 33 nt in length after trimming were selected for piRNA detection. Sequences were collapsed to remove identical sequences but retain information on read counts using the collapse tool from the NGS toolbox [33]. Furthermore, low-complexity reads were removed using the duster tool from the NGS toolbox. The resulting sequences were mapped to the Bos taurus 3.1 (Bt3.1) genome assembly and to chromosome Y from the 4.6.1 assembly with sRNA mapper. Only the best-scoring alignments were taken into account, and up to two non-templated 3′ nucleotides were allowed in order to successfully map sequences that were subject to post-transcriptional 3′ editing [34]. After mapping, the program reallocate (http://www.smallrnagroup-mainz.de/software.html) was used to assign read counts of multiple mapping sequences according to estimated local transcription rates based on uniquely mapping sequences.
piRNA cluster detection was performed with pro-TRAC version 2.1 [35,36], imposing a piRNA length of 26 to 33 bp and a minimum cluster length of 5000 bp. Genes falling within the detected clusters were retrieved according to Bt3.1 NCBI annotation, repeats and transposable elements were also retrieved, according to the Repeat Masker annotation available at the NCBI. Overlaps among HM and LM clusters were assessed with BedTools Intersect (http://bedtools.readthedocs.org).
miRNA detection and analysis miRNA detection and discovery was carried out with Mirdeep2 on Illumina high quality trimmed sequences. Bos taurus miRNAs available at MirBase (http://www.mirbase. org/) were used to accomplish known miRNA detection on the trimmed sequences. Known miRNAs from related species (sheep, goat and horse) available at MirBase were also input into Mirdeep2 to support the individuation of novel miRNAs.
The Mirdeep2 quantifier module was used to quantify expression and retrieve counts for the detected known and novel miRNAs. Differential expression analyses between the HM and LM fractions were run with the Bioconductor edgeR package [37]. miRNA cluster analysis was performed with Genesis [38]. Box-plot graphic was generated with BoxPlotR [39]. miRNA target prediction and functional analysis were performed by Ingenuity Pathway Analysis (IPA, Ingenuity System, www.ingenuity. com). Human homologous miRNAs were analyzed with microRNA Target filter (IPA) to attribute (experimentally observed) target genes. Finally miRNA target mRNA and the corresponding experimental Log Ratios were used for pathway analysis.

Statistical analysis
Data obtained from CASA and flow cytometry measurements were analyzed using the SAS™ package v 9.4 (SAS Institute Inc., Cary, NC, USA). The General Linear Model procedure (PROC GLM) was used to analyze the effect of technical replicates on semen quality parameters in the two fractions. The model included as fixed effects the bull and the replicate nested in the sperm fractions (HM and LM).
A mixed model procedure (PROC MIXED) was used to perform analysis on sperm quality parameters in order to evaluate the efficiency of the sperm separation into the HM and LM sperm fractions. The mixed model included the fixed effect of the sperm fraction (HM and LM), and bull as random. Results are given as adjusted least squares means ± standard error means (LSM ± SEM).

Discussion
Due to limitations in sRNA recovery from frozen thawed sperm, usually miRNA sperm profiling was achieved exclusively by microarray experiments or Real Time PCR [28,[30][31][32][33][34][35][36][37][38][39][40]. sRNAs profiling through NGS sequencing provides different advantages vs microarray experiments, such as discrimination of miRNAs that are very similar in sequence (isomiRs), detection of novel miRNAs [41] and the simultaneous piRNAs detection. This work firstly provides a comprehensive description of small RNAs isolated from HM and LM fractions from individual bull cryopreserved spermatozoa obtained by NGS sequencing and analysis. HM and LM fractions obtained after Percoll showed a good reproducibility between technical replicates, confirming the efficiency of the isolation of sperm fraction by Percoll gradient.
This procedure attenuated samples variability [42,43]. The increase of the proportion of living spermatozoa with intact acrosome (VIA) in HM fraction reported in this study was in agreement with previously results [44][45][46] in frozen-thawed bull and buffalo semen, indicating that the major part of dead spermatozoa were retained in the  Table 2 Comparison between known DEmiRNAs found in our study and in other studies. miRNAs were obtained from A) high and low motile fractions isolated from bovine cryopreserved semen, B) adult testis tissue from sheep well or under fed, C) human spermatozoa isolated from patients with normal or abnormal semen, D) human spermatozoa isolated from patients with normal semen or spermatogenic impairments, E) human spermatozoa isolated from patients with normal or vasectomized epididymis This study Guan et al., [55] Liu et al., [26] Abu-Halima et al., [27] Belleanée et al., [56] Specie upper layers of the gradient. Moreover, the HM fraction obtained in this study was characterized by sperm with fast motion characteristics and membranes integrity, two aspects strictly related to the fertilizing capacity [47][48][49].
Since piRNAs were firstly observed to have a putative role in gametogenesis in developing mouse male germ cells, they have been thought to be absent from mature spermatozoa [50]. Later, a survey of small RNAs in human sperm revealed sequence reads aligned to piRNA clusters located on several chromosomes and speculated their possible role in early embryo development in sperm [51]. A recent study identified a panel of piRNAs presents in seminal plasma that can serve as markers to distinguish fertile from infertile males [52]. Here we present the first characterization of piRNAs in two sperm fractions obtained by Percoll fractionation, showing a high diversity of piRNA content between the HM an LM fractions and the presence of a higher number of Table 2 Comparison between known DEmiRNAs found in our study and in other studies. miRNAs were obtained from A) high and low motile fractions isolated from bovine cryopreserved semen, B) adult testis tissue from sheep well or under fed, C) human spermatozoa isolated from patients with normal or abnormal semen, D) human spermatozoa isolated from patients with normal semen or spermatogenic impairments, E) human spermatozoa isolated from patients with normal or vasectomized epididymis (Continued)  To our knowledge, this is the first study able to compare piRNA expression in different sperm samples. Because of the low level of piRNAs conservation between even closely related species [50][51][52][53], the functional role of piRNAs dysregulation in HM and LM fraction remain to be further understood. On the contrary, bovine miRNA content in sperm was explored in different studies. In this study the total number of known bovine miRNAs isolated from semen cryopreserved in straws was in agreement with previously reported data obtained from NGS miRNA profiling of sperm isolated from caudal epididymis [54] or frozen sperm pellet [30] in Bos taurus. Finally, several of the top expressed miRNAs in this study have been previously reported as the most abundant in bovine sperm [30]. Moreover, miRNA expression comparison between the two fractions showed that about 10% of the miRNA are differentailly expressed. A similar percentage of expressed miRNA variation was observed in the semen of infertile men with semen abnormalities analyzed by microarray [25]. The high level of miRNA conservation among species supports a direct comparison of our data with data presented in previous studies on different species. About half of the known miRNAs found in this study have been previously reported in sperm or testis tissue from other species: Ovis aries [55] and Homo sapiens [26,27,56]. The relative expression in HM and LM fractions of several of these known miRNAs, including bta-miR-103, bta-miR-30b-5p, bta-miR-17-5p, bta-miR-106b, bta-miR-142-3p, bta-miR-34b, bta-miR-18a, bta-miR-34c, bta-miR-455-5p, bta-miR-10b, bta-miR-99b, bta-miR-1246, bta-miR-99a-5p, and bta-miR-1388-5p, was consistent with the relative abundance of their homologous miRNAs, observed in the normal vs abnormal sperm. Conversely, bta-miR-26a, bta-miR-24 and bta-miR-100 showed an opposite expression in our samples with respect to what previously described in literature [26,27,55]. bta-miR-122 and bta-miR-574 expression in the HM or LM fraction was only partially in agreement with what reported in previous studies [26,27,56]. Different miRNAs, including miR10b, miR-26, miR-34c and miR-99a, were also seen to change their expression level in underfed animals, and these variations were postulated to cause reduction in spermatozoa quality by disruption to Sertoli cell function and to increase germ cell apoptosis [55]. Functional analysis of the known DEmiRNAs showed targeting to mRNAs involved in different pathways, in particular "STAT3 Pathway", "PI3K/AKT Signaling" and "PTEN Signaling". The PTEN pathway is a crucial mediator of mitochondria-dependent apoptosis [57]. The role of PTEN in mammalian spermatogenesis under normal physiological conditions, consists in suppressing AKT activity to maintain activation of the RAF1/ERK signaling, which in turn maintains the normal function of the initial segment and, therefore, normal sperm maturation [58]. PTEN function is linked to its capacity of antagonizing the PI3K/AKT signaling. Akt1 and Akt2 knockout was seen to increase PTEN activity, probably inducing sperm apoptosis, decreasing spermatogenesis, sperm maturation and fertilization in male mice [59]. The inhibition of the STAT pathway in spermatozoa was reported to increase ROS production and calcium levels, and to decrease cellular ATP levels and mitochondrial membrane potential, that is consistent with cells undergoing apoptosis [60].
We postulate that the simultaneous low expression and up-regulation of different miRNAs could dysregulate PTEN, PI3K/AKT and STAT signaling and influence the apoptosis, vitality and motility in spermatozoa (Fig. 4). PTEN could be targeted by the simultaneous action of miR-17-5p, miR-26a-5p, miR-486-5p. According to previous results, miR-17-5p, miR-26a-5p up-regulation enhances AKT pathway activation by PTEN suppression and promotes cancer [61,62]. On the contrary, miR-486 plays a pro-apoptotic tumor-suppressor role [63], and its high expression was associated with a good prognosis in gastric adenocarcinoma [64]; however, miR-486 over expression in dystrophin-deficient mice was also observed to reduce PTEN expression [65]. In agreement with our results, two miRNAs found to be highly expressed in the LM fraction were reported to target the AKT pathway and promote apoptosis. miR-122 was reported to play a pivotal role as tumor suppressor by decreasing AKT3 levels, inhibiting cell migration and proliferation and inducing apoptosis [66], whereas miR-184 was found to be involved in suppressing cell survival and growth by targeting AKT2 in neuroblastoma cells [67]. Finally, miR-17-5p and miR-20a-5p, that we found to be under-expressed in the LM fraction and potentially target PTEN and STAT signaling, if down-regulated were proved to trigger cell apoptosis [68].

Conclusion
In conclusion we provide a protocol, based on small RNA sequencing, enabling to characterize miRNA and piRNA contents in cryopreserved bovine spermatozoa from single animals. We also provide a dataset of novel bovine miR-NAs and a first description of piRNA genomic clusters expressed in bovine spermatozoa in high and low motile sperm population. Small RNAs were seen to differ between HM and LM sperm fractions. Furthermore, some miRNAs differentially expressed in HM and LM fraction targeted genes associated with cell apoptosis, mitochondrial membrane potential and spermatogenesis alteration, indicating a functional redundancy, which might influence sperm motility and thus bull fertility.