MicroRNA signature in spermatozoa and seminal plasma of proven fertile men and in testicular tissue of men with obstructive azoospermia

MicroRNAs (miRNAs) have recently received a significant amount of attention due to their remarkable influence on post‐transcriptional gene regulation. In this study, we aim to provide a catalogue of miRNAs present in spermatozoa, seminal plasma and testicular tissue. Expression profiles of miRNA in spermatozoa and seminal plasma of 16 proven fertile men and testicular tissue of eight men with morphologically and/or histologically confirmed obstructive azoospermia were determined by microarray and RT‐qPCR in combination with bioinformatics analyses. A total of 123, 156 and 133 miRNAs were consistently detected in spermatozoa, seminal plasma and testicular tissue respectively. Sixty‐four miRNAs were shared across all sample types. Based on miRNAs expression level present in each group, correlation analysis showed moderate‐to‐strong correlations within the spermatozoa and seminal plasma samples and a wider range of correlations within the testicular tissue samples. The target genes of known miRNAs appeared to be involved in a wide range of biological processes related to reproduction, development and differentiation of germ cells. Our results suggest that there is a certain similarity between spermatozoa and seminal plasma for the relative miRNA expression changes with respect to testicular tissue and provide an overview of the miRNAs present in each sample type.

regulators of gene expression (Ha & Kim, 2014). To date, accumulating studies have indicated that miRNAs are participating in the control of many, if not all, biological and pathological processes (Gurtan & Sharp, 2013). In mammalian spermatogenesis, evidence indicates that even a slight alteration in miRNA expression activities can be detrimental to germ cell function, ultimately leading to impaired spermatogenesis and compromised fertility (Chang et al., 2012;Hayashi et al., 2008;Maatouk, Loveland, McManus, Moore, & Harfe, 2008). Altered expression of miRNAs in men with reduced fertility has been investigated (Abu-Halima et al., 2019, 2016. Specifically, differentially expressed miRNAs were reported in spermatozoa of infertile men with normozoospermic versus. subfertile men with asthenozoospermic and oligoasthenozoospermic (Abu-Halima et al., 2013), in testes of men with normal versus. impaired spermatogenesis (in Sertoli cell only, mixed atrophy and germ cell arrest) (Abu-Halima,  and in extracellular microvesicles obtained from seminal plasma from normozoospermic versus oligoasthenozoospermic men (Abu-Halima et al., 2016). Although several miRNAs have been found altered in certain spermatogenic impairments, the exact cellular function and biological mechanism of these miRNAs have been revealed for a very limited number of miRNAs. With this study, by using miRNA microarrays along with the bioinformatics analysis, we aim to provide a catalogue of miRNAs present in spermatozoa and seminal plasma of proven fertile men and testicular tissue of men with obstructive azoospermia.

| Study population
Following Institutional Ethics Committee of Saarland University (Nr. Ha195/11) approval, eight semen samples were obtained from men aged between 26 and 34 years (29 ± 3.12, Mean ± SD) with proven fertility, that is men with normal semen parameters and parenthood after the investigation (Table S1). Semen samples were provided after at least three days of sexual abstinence, and semen parameters, including volume, pH, viscosity, sperm count (10 6 /ml), motility (motile %) and morphology (%) were performed based on the World Health Organization, 2010, (WHO, 2010) guidelines. Using Puresperm 45%-90% discontinuous density gradients, 1 ml of the semen sample was placed on the upper layer before centrifugation at 500 g for 15 min at room temperature (RT). After centrifugation, the upper seminal plasma layer was aspirated and transferred to a new sterile 1.5-mL Eppendorf tube for total RNA including miRNAs purification. In addition, purified RNA from spermatozoa samples (n = 8) were obtained from our previous study (Abu-Halima et al., 2013). These samples were collected from proven fertile men aged between 24 and 34 years (34 ± 5.6, Mean ± SD) (Table S1). Purified RNA from testicular samples (n = 8) were obtained from our previous study (Abu-Halima, . These samples were obtained from men aged between 24 and 30 (27 ± 2.1, Mean ± SD), with normal testicular volume, no evidence of inflammatory granuloma or malignancy, and the testis exhibited nearly normal spermatogenesis process with many mature sperms noted and populated by normal Sertoli cells and Leydig cells number.

| Purification and quality assessment
Purification of miRNA from seminal plasma was performed as previously described (Wu et al., 2012), with a slight modification. Briefly, 200 µl of the seminal plasma was used for total RNA purification.
The total RNA was isolated by mixing the seminal plasma with an equal volume of TRIzol reagent (Thermo Fisher Scientific), and the procedure was completed according to the manufacturer's protocol. Three steps of phenol/chloroform purification were added to get rid of proteins. The quantity and purity were tested using the NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific).
The integrity of the purified RNAs was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies). DNase treatment was carried out to remove the gDNA carry-over as described previously (Abu-Halima et al., 2016).

| Analysis of miRNAs by microarray
The miRNA expression profile in the spermatozoa and testicular tissue samples was obtained from our previously generated and published raw data using SurePrint™ 8X60K Human v16 miRNA platform (Agilent Technologies) (Abu-Halima, Abu-Halima et al., 2013). As for seminal plasma samples, miRNA profiling analysis was carried out on the purified RNA fraction from eight seminal plasma samples using the SurePrint™ 8X60K Human v16 miRNA platform (Agilent Technologies). These platforms contain probes for the detection of 1,205 human miRNAs.
One hundred nanograms (ng) of the purified RNA were labelled and subsequently hybridized to the miRNA microarray chip as previously described (Abu-Halima, Abu-Halima et al., 2013).

| Reverse Transcription and Quantitative Real-Time PCR (RT-qPCR) of miRNA
The expression level of certain miRNAs in each sample was determined using the miScript PCR System (Qiagen). One hundred nanograms (ng) purified RNA from each sample, that is spermatozoa (Abu-Halima et al., 2013), testicular tissue (Abu-Halima,  and seminal plasma, were reverse transcribed using the miScript RT II kit (Qiagen).
The cDNA was then diluted 1:10, and 1 µl of cDNA was mixed with 10 µl 2X QuantiTect SYBR Green, 2 µl 10X Universal Primer Assay, 2 µl 10X Specific Primer Assay and H 2 O in a total volume of 20 µl.
The miRNA primer assays were selected based on their abundance level changes in men with subfertility compared to men with proven fertility and based on their known associations with spermatogenesis.

| Bioinformatics analysis
Since the samples came from different sources, we cannot directly compare the expression level of miRNAs across the sample groups. Therefore, we focused on relative differences between the groups.
Arrays data were quantile normalised by using the preprocessCore package of the R software (www.R-proje ct.org). The resulting values were shifted by the absolute value of the minimal expression level (determined over the whole data set) and additionally by 0.1 to ensure that all values are positive. Finally, a log2 transformation of the data was performed ( Figure 1a). Additionally, a detection matrix was generated where a miRNA is considered to be detected in a sample of its total probe signal is bigger than 3*probe error.
Probe error is the robust average of all the processed green signal errors for each replicated probe multiplied by the total number of probe replicates: the effective feature size fraction, the nominal spot area and the weight. Therefore, it is an estimate of the probe signal error provided by the Feature Extraction software from Agilent. To classify a miRNA as absent, present or undetermined, we defined two thresholds for each sample group computed from the processed expression levels of miRNAs detected in 0% and 100% of all samples in the corresponding group. In our analysis, we concentrated on miRNAs classified as present with respect to the three sample groups. The present identified miRNA list was then converted from v16 to v21 according to the updated miRBase database (Release 22.1). The miRNAs which were not matched/ identified were subsequently removed from the database. Then, we determined the within-group sample correlation computing Spearman's correlation coefficient based on miRNAs present in the corresponding group. Furthermore, pairwise miRNA correlations with the corresponding p-values were calculated for present miRNAs for each group separately. P-values were corrected using the Bonferroni correction. Finally, we performed a quantile rank analysis (Smalheiser et al., 2012) applied to all miRNAs in our data set to identify miRNAs whose rank changed significantly between two sample groups. The 2 -ΔCq relative quantitative procedure (Livak & Schmittgen, 2001) was used to detect the expression changes of single miRNAs. The small nuclear RNA (snRNA) RNU6B was chosen as a reference endogenous for normalisation.

| Subjects characteristics
In total, we compared 24 miRNA expression profiles, including eight profiles from the seminal plasma of fertile men, that is men with normal semen parameters and proven fertility. Eight miRNA expression profiles from spermatozoa of men with proven fertility were included from our previous study (Abu-Halima et al., 2013). We had, however, problems in collecting testicular tissues from men with proven fertility. Therefore, eight miRNA expression profiles from testicular tissues of men with morphologically and/or histologically confirmed obstructive azoospermia were included as control tissue from our previous study (Abu-Halima, .

| Present and absent miRNAs
The data set was quantile normalised using the R-package pre-processCore (R Development Core Team 2010) ( Figure 1a). To classify miRNAs as present or absent, we defined two expression level thresholds (one for absence and one for presence) for each sample group. These thresholds were computed as follows: Assuming that not detected miRNAs demonstrate lower expression levels and detected miRNAs tend to have a high expression level, we estimated the density distribution of miRNA expression values for miRNAs with 0% and 100% detection rate in the considered sample group using R-package sm (Shah, Leidinger, Blin, & Meese, 2010). The 1 percentile of the 100% detection rate density distribution and the 99 percentile of the 0% detection rate density distribution were set as the thresholds for the absent and present miRNAs respectively ( Figure 1b): We defined a miRNA as present or absent under the following conditions: a miRNA is considered as present if its expression value is above the upper threshold (99 percentile) and absent if its expression value is below the lower threshold (1 percentile), and otherwise, the miRNA was defined as undetermined. Second, a miRNA is present/absent in a sample group if it is present or absent in all samples of the considered group, that is in spermatozoa, seminal plasma and testicular tissue.
Out of 1,205 investigated miRNAs, 697, 729 and 691 were not detected (absent) and 126, 161 and 137 miRNAs were detected (present) in spermatozoa, seminal plasma and testicular tissue respectively (Tables S2 and S3). Using the miRCarta web-service tool (Backes et al., 2018), we converted the present miRNAs from v16 to v21 according to the updated miRBase database (Release 22.1, http://www.mirba se.org/). Three miRNAs (miR-1274b, miR-720 and miR-1280), five miRNAs (miR-1274b, miR-3676, miR-1274a, miR-720 and miR-1280) and four miRNAs (miR-1274b, miR-1274a, miR-720 and miR-1280) were not detected in the updated version of miRBase v22.1 and have subsequently been removed from the present miRNA list of each group, that is spermatozoa, seminal plasma and testicular tissue respectively. A detailed list of miRNAs that cannot be found in the destination annotation (in a more recent miRBase 22.1 version) is reported in bold in Table   S3. Therefore, in total, 123, 156 and 133 miRNAs were present in spermatozoa, seminal plasma and testicular tissue, respectively, according to miRBase v22.1. The Venn diagrams in Figure 2a,b show the overlap of present miRNAs sets between the three sample groups. As shown in Table 1, 18 miRNAs present only in spermatozoa, and 50 miRNAs in seminal plasma and 38 miRNAs in testicular tissue samples were identified. F I G U R E 1 (a) Box plots of shifted (by |minimal value| and 0.1) raw expression levels per sample on log2-scale and (b) estimated density distributions based on processed expression levels of miRNAs with a detection rate of 0% (red), 50% (green) and 100% (blue) in the sample groups. The two thresholds are visualised by two dashed vertical lines F I G U R E 2 Venn diagram of (a) absent miRNAs and (b) present miRNAs

| Correlation analysis results
To determine the degree of correlation across individual samples within each group, we computed the pairwise sample correlations (Spearman's correlation coefficient) based on miRNAs expression level present in the considered sample group (Table S4). Figure 3a shows a histogram of the resulting correlation values for each sample group separately. We found the strongest correlation be-

| Quantile rank analysis of miRNA expression data
To determine miRNAs whose expression was different between two sample groups without a direct comparison of the expression levels, we performed quantile rank analysis (Smalheiser et al., 2012). For each miRNA, we computed its mean expression for each sample group separately. In each group, the miRNAs were sorted by their mean group expression in a decreasing order to determine their ranks. Next, we computed for each pair of groups the rank differences of the miRNAs.
These values were approximately normally distributed with mean = 0 for all three group pairs and standard deviation of 195.5, 211.5 and 224.8 for spermatozoa/seminal plasma, spermatozoa/testicular tissue and seminal plasma/testicular tissue respectively (Figure 3e). Finally, the probability to achieve a more extreme rank change by chance was determined for all rank differences. Table S5 contains the 18 miRNAs with at least one significant rank difference comparing two groups.
These findings suggest that there is a certain similarity between spermatozoa and seminal plasma regarding the relative miRNA expression changes as compared to testicular tissue. However, there are also rank differences between spermatozoa and seminal plasma for six other miRNAs, which seem to be specific for these two groups.

| Most expressed and most stable miRNAs
We determined the most expressed and stable miRNAs. The most expressed miRNAs were defined by a high mean expression value and the most stable miRNA by a low standard deviation. This analysis was done for each tissue/cell type separately. Among the most expressed miRNAs, three miRNAs namely miR-4281, miR-2861 and miR-638 were shared across all sample groups and 7 miRNAs were shared between spermatozoa and seminal plasma as indicated in Table 2. As for the most stable miRNAs, we did not find a stable miRNA that was shared among the three tested groups (Table 2).

| miRNA target genes prediction
The Gene Ontology (GO) was considered to elucidate the biological function of the identified miRNAs for each sample group separately. To this end, we extracted miRNA targets using the miRWalk 2.0 (Dweep & Gretz, 2015) and we used the list of unique genes for the enrichment analysis. The most enriched GO terms for the spermatozoa, seminal plasma and testicular tissue are listed in Table S6.

| Validation of microarray results by RT-qPCR
Of the miRNAs present in each group, we validated 11 miRNAs from spermatozoa, 5 miRNAs from seminal plasma and 8 miRNAs from testicular tissue by RT-qPCR. These miRNAs were chosen based on their abundance level changes and based on their known associations with spermatogenesis. The selected miRNA expression level was largely concordant with the microarray expression data as expected (Figure 4a-c).

| D ISCUSS I ON
In this study, we investigated the miRNA expression profile in spermatozoa and seminal plasma of men with proven fertility and   -3202, miR-3190, miR-3137, miR-3682, miR-1305, miR-3156, miR-3127, miR-3125, miR-595, miR-1183, miR-1299 and miR-3945 have not been reported to be involved in any biological and/or cellular function associated with spermatogenesis or related processes. It remains to be seen whether the newly identified miR-NAs are specific for germ cells and/or if their aberration is associated with impaired spermatogenesis.
As for the normal testicular tissue, several studies reported that testicular-expressed miRNA changes depending on the stage of spermatogenesis (Hayashi et al., 2008;Maatouk et al., 2008) and the presence of various cell types in the testis (i.e., germ and somatic cells) increase the heterogeneity of the testis. According to the studies, each of these cell populations has a specific miRNA profile and/or a stage-specific miRNA profile. In our study, however, the testicular tissue samples were taken from 8 different men diagnosed with OA and with morphologically and/or histologically normal spermatogenesis, and each of these testicular tissues might contain different stages of normal cell division and development.
There is rather limited information on the cellular function of the expressed and stable miRNAs for the process of spermatogenesis.
The Let-7 family, a highly conserved miRNA family in the testis that contributes significantly to the control of male germ cell lineage, is specifically deregulated in testicular cancer (Johnson et al., 2007;Jung, Gupta, Shin, Uhm, & Lee, 2010). In addition, Let-7 family was observed to be abundantly detected in testicular tissue, spermatozoa and embryonic spent culture media of IVF correlated with successful pregnancy (Abu-Halima, Abu-Halima et al., 2013.
The expression level of miR-574-5p is reported to be significantly higher in the semen sample of subfertile men as compared to controls (Liu, Cheng, Gao, Wang, & Liu, 2012). The expression level of miR-574-5p is also reported to be lower in spermatozoa samples from oligoasthenozoospermic subfertile men as compared to those from normozoospermic fertile men (Abu-Halima et al., 2013). In men with Sertoli cell-only syndrome, the miR-574-5p expression level was significantly lower compared to controls (Abu-Halima, , indicating that sustained lower expression of miR-574-5p is pivotal for the progression of infertility. MiR-638 plays a role in pig spermatogenesis by regulating immature Sertoli cell growth and apoptosis by targeting the SPAG1 gene (Hu et al., 2017). The miR-320 family is predicted to target the protocadherin alpha proteins that are expressed in the mammalian testes. Protocadherin proteins play a crucial function in maintaining germ cells adhesion to Sertoli cells in the interstitial compartment of the testes (Dai et al., 2011;Marcon, Babak, Chua, Hughes, & Moens, 2008).
The results demonstrate higher samples correlation and thus homogeneity of miRNA pattern within the spermatozoa and seminal plasma groups as compared to testicular tissue group. The observed heterogeneity within the latter group could be explained by the fact that testis contains various cell types, and its composition may vary between the samples. Quantile rank analysis identifies miRNAs whose expression changes significantly relative to the other observed miRNAs. These analyses identified 18 miRNAs with significant rank changes of the mean expression values. Five of these miRNAs namely miR-143, miR-200c, miR-510, miR-3651 and miR-214 were shared between spermatozoa/testicular tissue and seminal plasma/testicular tissue. This may indicate a higher similarity between spermatozoa and seminal plasma as compared to the testicular tissue. That is also supported by the fact that from the most expressed miRNAs for spermatozoa and seminal plasma, seven miRNAs namely miR-4281, miR-2861, miR-638, miR-1225-5p, miR-320c, miR-1207-5p and miR-1202 were shared, while for testicular tissue three miRNAs miR-4281, miR-2861 and miR-638 were shared among all groups ( Table 2). The expression patterns of most miRNAs in testicular tissue samples have also be identified as abundant in testis as reviewed in Kotaja (2014).
The correlation analysis was performed to see whether specific miRNAs are commonly regulated. As for miRNAs clustered in the same genomic region, a correlation of expression might indicate a common promoter sequence orchestrating the expression of miR-NAs. Correlated miRNAs located on different chromosomes are potentially regulated by the same transcription factors (TFs) as reviewed extensively by Chen and Rajewsky (2007). However, it remains to be elucidated which TFs regulate miRNA genes in spermatogenesis. Identification of co-expressed miRNAs in spermatozoa, seminal plasma and testicular tissue is an important first step towards the identification of potential miRNA-regulating TFs.
Gene Ontology (GO) analysis revealed many genes involved in reproduction, development and differentiation. Specifically, in spermatozoa and testicular tissue, many of the enriched GO terms were directly or indirectly associated with spermatogenesis, sperm motility and flagellum, male germ cell, gamete generation and fertilisation.
As for the seminal plasma, most of the enriched terms were directly or indirectly associated with regulation of protein secretion, focal adhesion, peptidase activity, cytoplasmic membrane-bounded vesicle, positive regulation of male gonad development and Sertoli cell differentiation.

| CON CLUS ION
Our results suggest that there is a certain similarity between spermatozoa and seminal plasma for the relative miRNA expression changes with respect to testicular tissue. Taken together, several miRNAs of unknown function were identified or demonstrated different levels in each sample type, suggesting that aberration of these miRNAs may adversely affect the rate of spermatogenesis.

AUTH O R CO NTR I B UTI O N
MA carried out the experimental work and wrote the article. VG, CB and AK performed bioinformatics analysis. MH collected and diagnosed the semen samples. EM designed the experimental work and edited the article.