Sperm mRNAs and microRNAs as candidate markers for the impact of toxicants on human spermatogenesis: an application to tobacco smoking

Abstract

Spermatozoa contain a complex population of RNAs including messenger RNAs (mRNAs) and small RNAs such as microRNAs (miRNA). It has been reported that these RNAs can be used to understand the mechanisms by which toxicological exposure affects spermatogenesis. The aim of our study was to compare mRNA and miRNA profiles in spermatozoa from eight smokers and eight non-smokers, and search for potential relationships between mRNA and miRNA variation. All men were selected based on their answers to a standard toxic exposure questionnaire, and sperm parameters. Using mRNA and miRNA microarrays, we showed that mRNAs from 15 genes were differentially represented between smokers and non-smokers (p < 0.01): five had higher levels and 10 lower levels in the smokers. For the microRNAs, 23 were differentially represented: 16 were higher and seven lower in the smokers (0.004 ≤ p < 0.01). Quantitative RT-PCR confirmed the lower levels in smokers compared to non-smokers for hsa-miR-296-5p, hsa-miR-3940, and hsa-miR-520d-3p. Moreover, we observed an inverse relationship between the levels of microRNAs and six potential target mRNAs (B3GAT3, HNRNPL, OASL, ODZ3, CNGB1, and PKD2). Our results indicate that alterations in the level of a small number of microRNAs in response to smoking may contribute to changes in mRNA expression in smokers. We conclude that large-scale analysis of spermatozoa RNAs can be used to help understand the mechanisms by which human spermatogenesis responds to toxic substances including those in tobacco smoke.

• Human
• miRNA
• mRNA
• spermatogenesis
• spermatozoa
• tobacco

Introduction

It is now well established that spermatozoa contain messenger RNAs (mRNAs), which reflect germ cell transcriptional activity during different stages of spermatogenesis [Miller et al. 1999; Ostermeier et al. 2002]. In spite of the heterogeneity of the sperm RNA population, some transcripts seem to be strictly regulated [Lalancette et al. 2009]. Thus, the analysis of transcripts present in men’s spermatozoa represents a non-invasive method to study the processes involved during their spermatogenesis. More recently, it has also been shown that human spermatozoa bear small RNAs such as microRNAs (miRNAs) [Krawetz et al. 2011; Ostermeier et al. 2005]. MicroRNAs are short non-coding molecules that bind specifically, in conjunction with a protein complex, to several target messenger RNAs by complementary binding to the 3′-untranslated region [Eriksson et al. 2003] to trigger their degradation or transcriptional repression [Guo et al. 2010; Thomas et al. 2010]. Several studies have shown that many miRNAs are involved in the regulation of gene expression through the targeting of mRNAs during several processes such as cell proliferation, apoptosis, development, or cancer [Huang et al. 2011], including testis cancer [Voorhoeve et al. 2006]. One study has reported that spermatozoa mRNA profiles can be used to study the interaction between the environment and the germinal cells [Linschooten et al. 2009]. Moreover, microRNA expression levels have been associated with environmental exposures, especially cigarette smoke, making them potential biomarkers for better understanding the mechanisms of environmental diseases [Hou et al. 2011; Izzotti et al. 2009; Izzotti and Pulliero 2014].

Among the environmental agents to which humans can be exposed, cigarette smoke is one of the most common. Cigarette smoking releases more than 4,000 components into the human body. Dozens of them are classified as carcinogenic, mutagenic, and/or reproductive toxicants under European regulations [INRS 2008], and are thus likely to impact the quality of the gametes. Several studies have shown deleterious effects of tobacco smoking on semen quality [Rubes et al. 1998; Shrivastava et al. 2014]. Moreover, a link between paternal tobacco addiction and an increased risk of development of some childhood cancers has been reported [Ji et al. 1997]. The mechanisms by which cigarette smoking can impact spermatogenesis are multiple. Benzo(a)pyrene is formed during the combustion of tobacco, and one of its strongly reactive metabolites, the benzo(a)pyrene-diol-epoxide (BPDE), has the capacity to bind DNA covalently, leading to the formation of stable, potentially mutagenic, adducts. Indeed, the presence of such adducts can modify the three-dimensional conformation of the DNA. DNA adducts are recognized and repaired by the nucleotide excision-repair system, but repairs are error prone. Moreover, the large size of adducts and the resulting DNA distortion can cause DNA double-strand breaks, whose repair is likely to involve a loss of genetic material. These DNA adducts have been detected in the spermatozoa of smokers [Perrin et al. 2011a,b; Zenzes et al. 1999a]. In addition, these adducts were also found in human embryos conceived by couples where the man but not the woman smoke, [Zenzes et al. 1999b], suggesting that the transmission of DNA adducts to the human embryos occurs via the spermatozoon. Moreover, a recent study of paternal exposure to benzo(a)pyrene showed that it might affect the expression of several microRNAs in mouse embryo. This suggests that smoking induced miRNA variations in human spermatozoa may also affect human embryo development and may have a potential trans generational effect [Brevik et al. 2012]. It must however be noted that the specific nature of spermatogenesis impairment in smokers has never been investigated. The metabolic pathways affected are probably multiple, and these may interact with each other. Large-scale molecular techniques are now available to explore these mechanisms, which will allow the development of assays that can be used to define the molecular basis of spermatogenic dysfunction [Anton and Krawetz 2012].

Our aim was to conduct a preliminary innovative study on spermatozoa RNAs, to compare mRNA and miRNA profiles in spermatozoa of smokers and non-smokers, and search for potential relationships between mRNA and miRNA variation. Our results show that this is a promising approach for beginning to understand how spermatogenesis is impacted in men exposed to potential reprotoxic molecules.

Results

Study group, semen parameters, and RNA extracts

The two study groups included eight smokers and eight non-smokers of similar ages (respective mean age 40.25 vs. 37.87, p = 0.393). The experimental design and the techniques used on each sperm sample are illustrated in Figure 1 (see details in the Materials and Methods section). Semen parameters for all men are available in Table 1. The quantity of spermatozoa available for this study varied from 25 to 145 million. RNA recovery ranged from 6.48 to 81.71 ng/10⁶ spermatozoa (Table 2), a yield variability already described for RNA extraction from spermatozoa [Lalancette et al. 2009]. Comparison of semen parameters between smokers and non-smokers showed no statistical difference in sperm concentration (74.6+/−25.58 vs. 155.9+/−104.56, p = 0.065), or percentage of progressive sperm (motility a+b) (35.62+/−11.48 vs. 42.50+/−11.95, p = 0.260), respectively, but the percentage of morphologically normal sperm was lower in smokers (25.75+/−7.83 and 43.50+/−17.28 (p = 0.024), respectively.

Figure 1. Illustration of the experimental design of the study. (A) Methodology of participant selection. (B) Diagram illustrating the techniques used on each participant’s sperm sample. The eight smokers and the eight non-smokers are identified using letters A to H as indicated in Table 1. S: smoker; NS: non-smoker; miRNA: microRNAs; qPCR: quantitative PCR.

Table 1. Sperm parameters and profiles of smoker and non-smoker patients used for transcriptome study, miRNA chips, and miRNA quantitative PCR.

Group	Germetheque Number	Sperm concentration (M/ml)	Round cells (M/ml)	Motility (a/b/c/d)*	% NF**	Age/ Years	CIG/day (starting year)	Other drugs	Transcriptome microarray	miRNA microarray	miRNA qPCR
S A	R 010900069	79.6	9.6	0/30/30/40	26	48	20 (1980)	0	Yes	–	–
S B	R 010900090	102	3.8	0/35/30/35	33	38	15 (1991)	0	–	Yes	Yes
S C	R 010900093	52	8.8	5/25/30/40	18	52	15 (1979)	0	Yes	Yes	–
S D	R 010900171	89	1.8	20/30/20/30	26	38	10 (1987)	0	–	Yes	–
S E	R 010900209	85.6	15.6	10/10/20/60	15	36	15 (1989)	Alcohol occ	Yes	Yes	Yes
S F	R 010900217	67.6	3.2	15/30/10/45	38	36	10 (1999)	0	Yes	–	–
S G	R 010900229	96	3.6	0/25/25/50	30	40	10 (1989)	Alcohol occ	Yes	–	–
S H	R 010900455	25	6	15/35/10/40	20	34	15 (1994)	0	Yes	–	–
NS A	R 010900143	39.2	4.4	20/20/30/30	17	37	0	Alcohol occ	Yes	–	–
NS B	R 010900164	109	17.6	10/40/10/40	52	37	0	0	–	Yes	–
NS C	R 010900175	364.4	15.6	5/35/30/30	51	43	0	0	Yes	–	Yes
NS D	R 010900252	138.8	4.6	0/25/35/40	54	39	0	0	–	Yes	–
NS E	R 010900399	113	9.4	20/35/25/20	67	43	0	Alcohol occ	–	Yes	–
NS F	R 010900409	259.8	30	10/20/10/60	52	30	0	Alcohol occ	Yes	Yes	Yes
NS G	R 010900443	125	14	30/10/10/50	25	37	0	Alcohol occ	Yes	–	–
NS H	R 010900468	98	7.6	30/30/10/30	30	37	0	Alcohol occ	Yes	–	–

miRNA: microRNAs; S: smokers; NS: non-smokers. The 8 smokers and the 8 non-smokers are identified using letters A to H, respectively. *WHO (2009) criteria; **normal forms (David criteria); CIG: cigarettes; Alcohol occ: alcohol less than 3 glasses/week.

Table 2. RNA quantity obtained from each sperm sample.

Group	Germetheque Number	Total spz available (x10^6)	Total RNA (ng)	Yield (ng per 10^6 spz)
S	R 010900069	39.8	2397.55	60.24
S	R 010900090	51	2072.45	40.64
S	R 010900093	52	3718.00	71.50
S	R 010900171	45.5	1979.25	43.50
S	R 010900209	42.8	3497.37	81.71
S	R 010900217	44.6	2096.20	47.00
S	R 010900229	76.8	2344.07	30.52
S	R 010900455	25	1138.00	45.52
NS	R 010900143	39.2	1854.33	47.30
NS	R 010900164	54.5	2524.44	46.32
NS	R 010900175	145.7	1821.25	12.50
NS	R 010900252	69	2725.50	39.50
NS	R 010900399	56.5	2070.00	36.64
NS	R 010900409	90	3810	42.33
NS	R 010900443	125	810.00	6.48
NS	R 010900468	98	2752.91	28.09

S: smokers; NS: non-smokers; Spz: spermatozoa.

Population selection and mRNA levels in spermatozoa in smokers and non-smokers

Samples from six smokers and five non-smokers were used on mRNA microarrays (Table 1). After microarray data normalization, we selected 18,836 genes that were detected in more than five out of the 11 samples. We further performed a t test to look for gene expression differences between smokers and non-smokers. We selected 15 genes having p < 0.01 (Table 3), and we performed unsupervised hierarchical clustering of samples based on the expression of the selected genes. As shown in Figure 2, this procedure successfully classified six out of six smoker samples, and five out of five non-smoker samples. Moreover, we performed a two-way ANOVA to take into account a moderate level of alcohol as a factor (see Materials and Methods). Gene expression differences between smokers and non-smokers were not affected by alcohol for BAX, S100A13, B3GAT3, and PKD2 (p < 0.01).

Figure 2. Unsupervised hierarchical clustering of 15 differentially represented transcripts in spermatozoa after p value filtering of t test (p < 0.01) on six smokers and five non-smokers. Smoker samples are identified by the yellow bar and are located on the left of the cluster. Heatmap representation of over represented (red) and under represented (green) genes. The color scale bar indicates Log2 ratio of intensities.

Table 3. List of transcripts with differential level of detection in spermatozoa between smokers and non-smokers after t test (p < 0.01).

Gene	Description	Smoker median	Non-smoker median	Fold change (S/NS)	p Value	Proposed biological function (Gene Ontology)
BAX	ref|Homo sapiens BCL2-associated X protein (BAX), transcript variant epsilon, mRNA [NM_138764]	2.83 3.40	3.69 3.81	−1.82 −1.33	1.11E-03 5.53E-03	Apoptosis Regulation of apoptosis (GO:0035904)
THC2670523	Unknown	2.90	3.60	−1.62	1.08E-03	–
S100A13	ref|Homo sapiens S100 calcium binding protein A13 (S100A13), transcript variant 1, mRNA [NM_001024210]	3.59	4.23	−1.56	8.18E-04	Positive regulation of cell proliferation (GO:0008284) Regulation of cell shape (GO:0008360)
LOC340947	gb|PREDICTED: Homo sapiens similar to eukaryotic translation initiation factor 3, subunit 6 interacting protein (LOC340947), mRNA [XM_001131609]	2.72	3.32	−1.52	9.55E-04	–
OASL	ref|Homo sapiens 2′-5′-oligoadenylate synthetase-like (OASL), transcript variant 1, mRNA [NM_003733]	3.22	3.74	−1.43	1.00E-03	Defense response to virus (GO:0051607) Double-stranded RNA binding (GO:0003725)
HNRNPL	ref|Homo sapiens heterogeneous nuclear ribonucleoprotein L (HNRNPL), transcript variant 1, mRNA [NM_001533]	3.80	4.27	−1.39	3.18E-03	mRNA processing (GO:0006397) RNA binding (GO:0003723) Nucleotide binding (GO:0000166)
CNGB1	ref|Homo sapiens cyclic nucleotide gated channel beta 1 (CNGB1), transcript variant 1, mRNA [NM_001297]	3.99	4.37	−1.30	5.42E-04	Photoreceptor cell maintenance (GO:0045494) Protein localization to organelle (GO:0033365) Membrane (GO:0016020)
A_24_P814246	Unknown	4.67	5.04	−1.29	3.70E-03	–
ODZ3	ref|Homo sapiens odz, odd Oz/ten-m homolog 3 (Drosophila) (ODZ3), mRNA [NM_001080477]	4.16	4.52	−1.28	3.20E-03	Homophilic cell adhesion (GO:0007156) Signal transduction (GO:0007165) Membrane cellular component (GO:0016020) Protein heterodimerization activity (GO:0042803)
B3GAT3	ref|Homo sapiens beta-1,3-glucuronyltransferase 3 (glucuronosyltransferase I) (B3GAT3), mRNA [NM_012200]	3.68	4.04	−1.28	9.64E-03	Membrane (GO:0016020) 3-beta-glucuronosyltransferase activity (GO:0015018)
A_24_P187340	Unknown	4.43	4.10	−1.26	2.57E-03	–
THC2661208	thc|ALU7_HUMAN (P39194) Alu subfamily SQ sequence contamination warning entry, partial (9%) [THC2661208]	3.36	2.74	1.54	4.48E-03	–
FBXL8	ref|Homo sapiens F-box and leucine-rich repeat protein 8 (FBXL8), mRNA [NM_018378]	3.26	2.58	1.60	8.30E-03	Protein binding (GO:0005515)
PKD2	ref|Homo sapiens polycystic kidney disease 2 (autosomal dominant) (PKD2), mRNA [NM_000297]	3.21	2.40	1.75	7.94E-03	Branching involved in ureteric bud morphogenesis (GO:0001658) Calcium ion transport (GO:0070588) Cellular response to osmotic stress (GO:0071470) Cellular response to reactive oxygen species (GO:0034614) Centrosome duplication (GO:0051298) Detection of mechanical stimulus (GO:0050982) Determination of left/right symmetry (GO:0007368) Embryonic placenta development (GO:0001892)
AL137292	gb|Homo sapiens mRNA; cDNA DKFZp434L1016 (from clone DKFZp434L1016) [AL137292]	3.53	2.70	1.78	2.22E-03	–

miRNA levels in spermatozoa from smokers and non-smokers

We performed miRNA microarray analysis with spermatozoa RNA from four smokers and four non-smokers (Table 1), three of which were included in the mRNA microarray analysis. We first identified the miRNAs that showed the highest level of expression in the eight samples. To this end, we selected the miRNAs with an average level of detection corresponding to the 10% most represented miRNAs, and a variance corresponding to the 20% least variable. The 19 most represented miRNAs in spermatozoa of the eight patients are listed in Table 4.

Table 4. List of miRNAs present in human spermatozoa with an average level of expression included in the 10% of most represented miRNAs.

	Average	Standard deviation
hsa-miR-638	6,290	0,000
hsa-miR-2861	6,120	0,083
hsa-miR-1915	5,888	0,076
hsa-miR-1280	2,752	0,060
hsa-miR-634	2,319	0,085
hsa-miR-149	2,318	0,058
hsa-miR-636	2,313	0,092
hsa-let-7f-1*	2,222	0,086
hsa-miR-23c	2,199	0,084
hsa-miR-1539	2,171	0,053
hsa-miR-18b*	2,161	0,058
hsa-miR-135a	2,154	0,092
hsa-miR-4307	2,139	0,068
hsa-miR-135b	2,133	0,091
hsa-miR-302a	2,113	0,048
hsa-miR-30d	2,098	0,092
hsa-miR-519e	2,096	0,092
hsa-miR-136	2,091	0,092
hsa-miR-520e	2,083	0,086

miRNA: microRNAs.

Thus, we selected 522 microRNAs with an expression above background in four out of eight samples (the background being defined by the average expression of all microRNAs), and identified 23 microRNAs differentially represented between smokers and non-smokers with a p value between 0.004 and 0.094 (median 0.0042) (Table 5), seven of them being down regulated, and 16 up regulated, in smokers. Among these 23 miRNAS, 12 had a p value < 0.05. In spite of the relatively low level of global expression in spermatozoa, the unsupervised hierarchical clustering of samples based on the 23 selected microRNAs allowed us to discriminate samples by smoking status (Figure 3). As for mRNA analysis,we performed a two-way ANOVA to take into account moderate alcohol consumption. Differences between smokers and non-smokers were not affected by amoderate level of alcohol for hsa-miR-126*, hsa-miR-92a-2*, hsa-miR-26a-1, hsa-miR-517c, and hsa-miR-340* (p < 0.01).

Figure 3. Unsupervised hierarchical clustering of 23 differentially represented miRNAs in spermatozoa after p value filtering of t test (p < 0.01) on four smokers and four non-smokers. Smoker samples are identified by the yellow bar and are located on the left of the cluster. Heatmap representation of over represented (red) and under represented (green) genes. The color scale bar indicates Log2 ratio of intensities.

Table 5. List of miRNAs with differential level of detection in spermatozoa between smokers and non-smokers (p < 0.1).

Probe name	Systematic name	Smoker median	Non-smoker median	Fold change (S/NS)	P value
A_25_P00015143	hsa-miR-1246	2.68	2.20	1.39	0.030
A_25_P00015182	hsa-miR-548n	2.13	1.87	1.20	0.089
A_25_P00016168	hsa-miR-3647-3p	2.20	1.98	1.16	0.017
A_25_P00010600	hsa-miR-126*	2.09	1.90	1.14	0.004
A_25_P00015462	hsa-miR-762	2.43	2.24	1.14	0.090
A_25_P00013230	hsa-miR-92a-2*	2.06	1.89	1.13	0.010
A_25_P00013373	hsa-miR-23b*	2.02	1.87	1.11	0.043
A_25_P00013650	hsa-miR-33b*	2.30	2.16	1.10	0.037
A_25_P00012611	hsa-miR-522	2.13	1.99	1.10	0.066
A_25_P00012514	hsa-miR-181d	2.08	1.94	1.10	0.094
A_25_P00010635	hsa-miR-584	2.02	1.89	1.09	0.079
A_25_P00013345	hsa-miR-221*	1.98	1.86	1.09	0.010
A_25_P00013277	hsa-miR-129*	2.25	2.15	1.07	0.074
A_25_P00012446	hsa-miR-451	2.06	1.97	1.06	0.055
A_25_P00012242	hsa-miR-186	1.99	1.90	1.06	0.028
A_25_P00015174	hsa-miR-548g	2.01	1.92	1.06	0.072
A_25_P00013546	hsa-miR-340*	1.99	2.06	−1.05	0.011
A_25_P00013195	hsa-miR-26a-1*	1.90	2.04	−1.10	0.062
A_25_P00013978	hsa-miR-296-5p	2.08	2.25	−1.13	0.072
A_25_P00015353	hsa-miR-517c	1.87	2.07	−1.15	0.011
A_25_P00012588	hsa-miR-520d-3p	2.03	2.24	−1.16	0.048
A_25_P00016003	hsa-miR-3940	1.74	2.12	−1.30	0.075
A_25_P00016234	hsa-miR-3679-5p	3.82	4.32	−1.41	0.043

miRNA: microRNAs; S: smokers; NS: non-smokers.

Quantitative real-time RT-PCR (RT-qPCR) analysis of smoker (R01090090, R010900209) and non-smoker samples (R010900175, R0109409) was performed for miRNAs differentially represented on the microarray. We choose two miRNAs up regulated in smokers (hsa-miR-129*, hsa-miR-1246), and three miRNAs down regulated in smokers (hsa-miR-296-5p, hsa-miR-3940, hsa-miR-520d-3p). Quantitative PCR (qPCR) confirmed that hsa-miR-296-5p, hsa-miR-3940 (miRBase mature sequence hsa-miR-3940-5p), hsa-miR-520d-3p were down regulated in smokers (Figure 4). The qPCR did not confirm the up regulation for hsa-miR-129* and hsa-miR-1246 in smokers. The three down regulated miRNAs are at a level that is two to four times less than that present in smokers.

Figure 4. Expression levels of hsa-miR-296-5p, hsa-miR-3940, and hsa-miR-520 d-3p in the spermatozoa of smoking and non-smoking men analysed using quantitative real-time PCR (RT-qPCR). Values represent mean +/- SEM. t test was performed between smokers and non-smokers (*p < 0.01 and **p < 0.001).

Putative miRNA target gene expression in human spermatozoa

Putative targets for the 23 miRNAs differentially represented on the microarray according to smoker status, and for hsa-miR-638 (highly represented) used as control were determined using TargetScan (version 6.2, June 2012) [Lewis et al. 2005]. To determine whether miRNAs might be modulators of mRNAs that were differentially detected in smokers, we examined whether our selected mRNAs identified as affected by smoking included predicted targets of each of the miRNAs differentially represented between smokers and non-smokers. We found that nine out of the 15 are predicted to be targets of 11 miRNAs out of the 23 differentially detected in smokers (Figure 5). This suggests that alterations in the levels of a relatively small number of miRNAs in response to smoking could potentially contribute to smoking-associated changes in mRNA expression. Interestingly, the microarray data show that the direction of variation of mRNAs compared to microRNAs is in agreement with a potential interaction during human spermatogenesis for five down-regulated genes (B3GAT3, HNRNPL, OASL, ODZ3, and CNGB1) and one up-regulated gene, PKD2. Indeed, PKD2 is a putative target of hsa-miR-520d-3p that we found to be down regulated in smokers by microarray and qPCR analysis.

Figure 5. Unsupervised hierarchical clustering of the differentially represented miRNAs (Y axis) and mRNAs (X axis) during spermatogenesis of studied patients, according to the smoking status. Only the miRNAs having target gene, in the TargetScan database v 6.2, corresponding to our 15 differentially represented mRNAs are represented. Black boxes correspond to miRNA target genes; red boxes: gene/miRNA over represented in Smokers; green boxes: gene/miRNA under represented in Smokers.

Discussion

In the present study, we examined the possibility that sperm mRNAs and miRNAs could be used as candidate markers to study the impact of toxicants such as cigarette smoke on human spermatogenesis. To the best of our knowledge, this is the first study of mRNA and miRNA profiles during spermatogenesis in men exposed to environmental toxicants, and let us highlight two interesting points.

First, we show that by using total RNA extracts from spermatozoa of smokers and non-smokers, we were able to select transcripts and miRNAs differentially represented in human spermatozoa based on smoking status. Among the 15 transcripts shown to be differential, PKD2, BAX, and HNRNPL are known to be involved in spermatogenesis regulation. Evidence from online resources indicate that the proteins they encode are present in human adult germ cells (www.proteinatlas.org, www.humanproteomemap.org [Fagerberg et al. 2014; Kim et al. 2014]). Polycystic Kidney Disease 2 (PKD2) overrepresented in smokers, encodes, the polycystin-2 (PC-2) involved in mesonephric tubule morphogenesis. PKD isoforms located in cilia have been demonstrated to play a role in the transduction of mechanosensation via calcium signalling [Nauli et al. 2003], and have been shown to operate in sperm as a determinant of directional movement inside the drosophila female reproductive tract [Gao et al. 2003; Yang et al. 2011]. The calcium channel PKD2 mediates sensory perception on cilia and flagella, and these defects can contribute to ciliopathic diseases. In human, asthenozoospermia and necrozoospermia have been described in patients bearing autosomal dominant polycystic kidney diseases [Vora et al. 2008]. If PKD2 may be involved in the modulation of flagellar properties, and given that there is no difference in motility between smokers and non-smokers in our study, one could hypothesize that an increased level of PKD2 transcripts in smokers may be a part in a compensatory mechanism that maintains normal motility despite a negative impact of smoking. The two other transcripts BAX and HNRNPL appear under represented in smokers. BAX encodes a well-known pro-apoptotic factor described as involved in male germ cell apoptosis [Carreau et al. 2011; Saygin et al. 2011; Sinha Hikim et al. 2003]. HNRNPL encodes an mRNA processing factor shown to be expressed in spermatids and involved in post-transcriptional expression during testicular spermatozoa maturation [Chapman et al. 2013; Li et al. 2012]. HNRNPL is supposed to be a key factor involved in the bioprocess of apoptosis during spermatogenesis. Interestingly, a recent study of spermatozoal RNA profile in smokers reported that apoptosis appeared to be the most influenced process, which means that damaged cells may not be eliminated as efficiently during the germ cell cycle in smokers [Linschooten et al. 2009]. Nevertheless, although differential transcripts between smokers and non-smokers in our study are not the same as those reported in this previous study, both support the hypothesis that tobacco smoking impacts the germ cell apoptotic process. Indeed, it has been shown that although discordances may be found between data sets for the same biological process, common patterns of expression and outcome predictions can be identified in data sets generated by independent labs [Hu et al. 2006]. However, it should be stressed that further studies with independent samples are required to clarify how smoking affects spermatozoa mRNA profiles.

Regarding miRNAs, we were able to distinguish spermatozoal miRNA profiles according to the smoking status despite the small amount of miRNAs that are present in our total human spermatozoa extracts. We show that 23 miRNAs were differentially represented between smokers and non-smokers. Among them, hsa-miR-1246 was previously shown to be up regulated in spermatozoa of smokers [Marczylo et al. 2012]. It is interesting to note that this previous study and ours find that miRNAs in spermatozoa of smokers seem to be more frequently up regulated than down regulated. This pattern is different than the one described in cancer or in airway epithelium of smokers where miRNAs are preferentially down regulated [Schembri et al. 2009]. This suggests that cigarette smoking could impact human spermatogenesis via mechanisms distinct from those acting during tumorigenesis. Using RTqPCR, three miRNAs were confirmed to be down regulated in smokers. Their roles in human spermatogenesis are still unclear, but two of them are known to be involved in human spermatogenesis or in tobacco smoke induced cell modifications. Hsa-miR-296-5p was shown to be involved in tumour angiogenesis and was present in embryonic stem cells, suggesting it may have a role in regulation of early embryonic development. Hsa-miR-520d-3p was shown to be present in human testis and was down regulated in the testis of patients with non-obstructive azoospermia. Altogether there are presumptions for a potential role of those miRNAs in the spermatogenetic modifications induced by cigarette smoking.

Second, we show that it is possible to find a relationship between the dysregulation of mRNAs and miRNAs in human spermatozoa according to smoking status. Among miRNAs that were down regulated in smokers, it is interesting to note that some of their theoretical targets are present among transcripts that appear up regulated. Indeed, hsa-miR-520d-3p potentially regulates PKD2, found to be up regulated in the same group and reported to be involved in regulating spermatogenesis [Gao et al. 2003; Yang et al. 2011]. Among miRNAs that were up regulated in smokers, hsa-miR-23b is predicted to target ODZ3, B3GAT3, and CNGB1, which we found to be down regulated. Hsa-miR-23b was reported to be up regulated in spermatozoa of asthenozoospermic patients [Abu-Halima et al. 2013]. ODZ3 encodes a protein expressed in neurons that may function as a cellular signal transducer [Beckmann et al. 2011; Nunes et al. 2005]. B3GAT3 encodes a member of the glucuronyltransferase gene family [Gulberti et al. 2003; Ouzzine et al. 2000], and CNGB1 a cyclic nucleotide-gated cation channel [Becirovic et al. 2010]. Although a biological relationship could not be demonstrated, our results suggest that alterations in the levels of a small number of miRNAs in response to smoking could potentially contribute to smoking–associated changes in mRNA expression, that potentially impact cellular transducers and possible consequences for apoptosis or spermatozoa motility.

Our study of mRNAs and miRNAs from human spermatozoa that contain very low levels of RNAs has nevertheless some limitations. In particular, the technique we used for RNAs extraction was not modified to fully recover the miRNAs and thus may explain our small yield of miRNAs and differences compared to miRNAs previously reported in spermatozoa [Krawetz et al. 2011; Hammoud et al. 2014]. Indeed, it has been reported that the concentration of total RNA in the starting material is critical in determining the yield of small RNA [Kim et al. 2012]. As a very low level of RNA is present in human spermatozoa, a technical adaptation may improve the significance of the results.

In conclusion, this pilot study shows that large-scale analysis of sperm RNAs can be used to help understand the mechanisms that mediate the response to tobacco smoke exposure during human spermatogenesis. However, the identification of the complete molecular pathways affected is a real challenge because of the low RNA levels present in human spermatozoa. While it has been shown that smoking-induced differences in spermatozoal mRNA can be detected with few sperm samples [Linschooten et al. 2009], the extension to a larger number of samples will improve the significance of our results. Moreover, new tools such as RNA-seq may soon provide the means to more fully characterize the RNAs present in human spermatozoa and thus to elucidate the fine regulatory mechanisms that mediate responses to environmental agent exposure during human spermatogenesis.

Materials and Methods

The Germetheque Biobank provided all samples to the authors. Germetheque obtained informed consent from each patient for inclusion of samples in the biobank and for their use for studies on human fertility in accordance with the Helsinki Declaration of 1975 on human experimentation. The Germetheque Scientific Committee approved the present study design. Consequently, no IRB approval was sought.

Participants and tobacco exposure

Smoker and non-smoker men consulted for primary infertility in the infertility clinic of the University Hospital of La Conception (Marseilles, France). A standardized questionnaire was used to evaluate cigarette smoking (with details about smoking levels, duration, and contribution of passive smoke), alcohol consumption, and occupational exposure [De Fleurian et al. 2009].

None of the participants had any known heritable condition that could favour an abnormal semen profile, no current medication known to impair spermatogenesis, no drug consumption (except cigarette), which can influence spermatogenesis. They had not taken medication for at least six months and their Body Mass Index was consistent with a normal body mass. Moreover, they were not occupationally exposed to chemical or physical agents likely to have effects on the spermatogenesis. Alcohol consumption of more than three drinks/w was a factor of exclusion. Patients smoking at least 10 cigarettes/d for at least six months defined the smokers group. Patients who had never smoked or had stopped cigarette smoking for at least 10 y were considered the non-smokers.

Sperm analysis and preparation

Semen was collected via masturbation after a period of sexual abstinence of 2-6 d. All patients gave an informed consent for the conservation of the remnant sperm after semen analysis in the Germetheque biobank. After 30 min of liquefaction, semen analysis was performed according to the WHO [2009] criteria and according to the French David classification for the morphology analysis [Auger and Eustache 2000]. Within 30 min, a sample of one to two ml of sperm was then centrifuged at 3,000 rpm for 15 min. Seminal plasma was removed; the pellet was snap-frozen and stored in the Germetheque biobank in liquid nitrogen until use.

Experimental design

Sperm samples were selected from those stored in the Germetheque Biobank (n = 158), and were selected according to the result of the occupational exposures questionnaire, the number of spermatozoa available, and the sperm parameters. Therefore, eight smokers and eight non-smokers were selected, with a standard profile for morphological abnormalities. Whole mRNA microarray was performed on RNA extracts from six smokers and five non-smokers, miRNA microarray on RNA from four smokers and four non-smokers, and miRNA RT-qPCR analysis on RNA from two smokers and two non-smokers.

RNA extraction

Subsequent to storage, samples were thawed and washed twice in 2 ml of PBS, then suspended in round cell lysis buffer (0.1% SDS, 0.5% Triton X-100 in RNase free H2O). RNA was extracted using the RNeasy mini kit (Qiagen GmbH, Germany) with minor modifications. Sperm cells were suspended in 600 µl of lysis buffer for 10⁷ cells. This cell suspension was lysed using a Politron as previously described [Saade et al. 2007], to maximize the recovery of spermatozoa RNA. The lysates were then processed according to the manufacturer’s instructions using the animal cell protocol. After RNase-free DNase treatment (Qiagen GmbH), the RNA was eluted with 30 µl of RNase free H₂O. RNA concentration was determined using a nanodrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA quality was assessed on an Agilent 2100 Bioanalyzer with the Agilent RNA 6000 nano kit (Agilent Technologies, Santa Clara, CA, USA), showing that the RNA is highly fragmented and that the 18S and 28S ribosomal RNAs are absent (see supplementary Figure S2).

Quality control of RNA extracts

One aliquot of 400 -500 ng of each sample was subjected to reverse transcriptase PCR using oligo-dT 25 primers and SuperScript II (Invitrogen, Carlsbad, CA, USA). The DNA was amplified by polymerase chain reactions using selected oligonucleotides in a final volume of 20 µl. To assure the purity of RNA extracts and be sure that there were no round cell extracts, we choose primers such as CD45, c-kit, or E-cadherin which are positive markers for leukocytes, testicular germ cells, and epithelial cells, respectively [Lambard et al. 2004]. Genomic DNA contamination was tested under conditions described above using a set of primers specific to PRM2 that span an intron [Miller et al. 1999] (see supplementary Figure S1). The positive controls for leukocytes and immature germ cell contamination were as follow: total RNA extracted from blood sample and total RNA extracted from human testis (636533 Clontech, purchased from Ozyme, Saint Quentin en Yvelines, France). Samples were also submitted to PCR.

Whole human genome Agilent microarray gene expression

Cy3 labelled cRNA was generated from 500 ng total RNA of patient’s spermatozoa, using the Low RNA Input Fluorescent Linear Amplification Kit (Agilent Technologies) following manufacturer’s instructions. Equal amounts of purified Cy3 labelled cRNA were hybridized for 17 h at 65°C to the Agilent whole human genome 4x44k microarray (G4112F). Each Agilent slide was hybridized with two smoker samples and two non-smoker samples according to the manufacturer’s instructions, taking care to change the order between smokers and non-smokers for each slide. The microarray was washed using the Gene Expression Wash Buffer Kit (Agilent Technologies) and scanned as for the standard Agilent protocol. The data were processed using Feature Extraction software. Data have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE44133.

Microarray data normalization

Whole genome microarray quantification files derived from the Agilent Feature Extraction Software were analyzed using the Agilent R package (ftp://tagc.univ-mrs.fr/public/AgiND/) for quality control and normalization. The quantile method was used for normalization. Among probes obtained from Feature Extraction files, those for which at least 60% of samples are above background using the “gIsWellAboveBg” feature were conserved for statistical analysis (n = 18,836).

To reduce the effect of noise in the miRNA microarray data, we aggregated probe signals to have a more robust measure of expression, by computing the median of the intensity values “gProcessedSignal” of the set of probes corresponding to each miRNA. The samples were then normalized by the quantile method allowing a comparison of experimental samples between each other.

Statistical analysis

Statistical analyses were performed using the TIGR MeV (MultiExperiment Viewer) v4.5 and the GeneANOVA programs [Didier et al. 2002]. A student t test was applied to identify gene and miRNA whose expression differed between the smoker and the non-smoker groups. In addition, a two-way ANOVA was used to take simultaneously into account the effect of smoking and alcohol consumption on gene expression, each factor being binary: smokers/non-smokers and alcohol/non-alcohol.

MicroRNAs microarray preparation and hybridization

Total spermatozoa RNA extracts (100 ng) were labeled using the Agilent miRNA complete labeling and hybridization kit (5190-0456). Hybridization was performed on Agilent human miRNA microarray 8x60K (G-4872A-031181) sourced from the Sanger miRBASE. All procedures were performed according to the manufacturer’s instructions. Microarrays with labeled samples were hybridized 20 h at 55°C and washed using the Gene Expression Wash Buffer Kit (5188-5327). Arrays were scanned using a High-Resolution Microarray Scanner and images were analyzed using Feature extraction software. Data were corrected as described above using AgiND R library with the same parameters and filters. Data are accessible through GEO Series accession number GSE44134.

Real time miRNA qPCR

To validate miRNA microarray results, a subset of miRNAs identified as more concentrated or less concentrated in smokers by microarray analysis were analyzed for smoker-related differential expression. RT-qPCR was performed using a standard TaqMan MicroRNA assay according to the manufacturer’s standard protocol (Applied Biosystems, Carlsbad, CA, USA) [Chen et al. 2005]. Twenty nanograms of the isolated total RNA was reverse transcribed into cDNA using the TaqMan Multiscribe Reverse Transcription kit (4366596 Applied Biosystems). The MiRNA sequences and assay IDs are listed in Table 6. Each reaction was performed in triplicate in two smokers and two non-smokers. Quantification was performed using validated TaqMan probes, the TaqMan MicroRNA assay and the TaqMan Gene Expression Master Mix (Applied Biosystems). The fold change in miRNA level was calculated using has-miR-548q, as endogenous control.

Table 6. MicroRNAs (miRNA) sequences and assay IDs used for miRNA quantitative real-time PCR.

miRNA	miR sequence	Assay ID
hsa-mir-129*	49-AAGCCCUUACCCCAAAAAGUAU-70	002298
hsa-mir-1246	11-AAUGGAUUUUUGGAGCAGG-29	462575_mat
hsa-mir-296-5p	14-AGGGCCCCCCCUCAAUCCUGU-34	000527
hsa-mir-3940	23-GUGGGUUGGGGCGGGCUCUG-42	462391_mat
hsa-mir-520d-3p	54-AAAGUGCUUCUCUUUGGUGGGU-75	002743
hsa-miR-548q	GCUGGUGCAAAAGUAAUGGCGG	241103_mat

Supplementary material available online Supplementary Figures

Acknowledgments

We are grateful to the patients who gave their informed consent for the use of their samples for research. We thank C. Metton and M.J. Fays-Bernardin for technical assistance and Germetheque support. We also thank M. Mitchell for the text revision, and M. Guichaoua for her precious advise.

Declaration of interest

This work was supported by grants from INSERM: Institut National de la Santé Et de la Recherche Médicale and APHM: Assistance Publique – Hôpitaux de Marseille. The Germetheque biobank was supported by grants from the ANR (Agence Nationale de la Recherche), the Agence de la biomedicine, the Centre Hospitalier Universitaire de Toulouse, and APHM (Assistance Publique Hôpitaux de Marseille). All authors declare no declarations of interest and no competing financial interest.

Author contributions

Played a role in the sample collection, the biological experiments, the study design and the manuscript writing: CMG. Had a major role in the bioinformatics data analysis: CL, AB, PR. In charge of microarray and quantitative PCR experiments: GV, MY, NB. Played a role in the patient recruitment: JP, ISM. Direct responsibility for the manuscript: CN, PR. All authors contributed to the writing and revision of the manuscript.