Men who inject opioids exhibit altered tRNA-Gly-GCC isoforms in semen

Abstract

In addition to their role in protein translation, tRNAs can be cleaved into shorter, biologically active fragments called tRNA fragments (tRFs). Specific tRFs from spermatocytes can propagate metabolic disorders in second generations of mice. Thus, tRFs in germline cells are a mechanism of epigenetic inheritance. It has also been shown that stress and toxins can cause alterations in tRF patterns. We were therefore interested in whether injecting illicit drugs, a major stressor, impacts tRFs in germline cells. We sequenced RNA from spermatocytes and from semen-derived exosomes from people who inject illicit drugs (PWID) and from non-drug using controls, both groups of unknown fertility status. All PWID injected opioids daily, but most also used other illicit drugs. The tRF cleavage products from Gly-GCC tRNA were markedly different between spermatocytes from PWID compared to controls. Over 90% of reads in controls mapped to shorter Gly-GCC tRFs, while in PWID only 45% did. In contrast, only 4.1% of reads in controls mapped to a longer tRFs versus 45.6% in PWID. The long/short tRF ratio was significantly higher in PWID than controls (0.23 versus 0.16, P = 0.0128). We also report differential expression of a group of small nucleolar RNAs (snoRNAs) in semen-derived exosomes, including, among others, ACA14a, U19, and U3-3. Thus, PWID exhibited an altered cleavage pattern of tRNA-Gly-GCC in spermatocytes and an altered cargo of snoRNAs in semen-derived exosomes. Participants were not exclusively using opioids and were not matched with controls in terms of diet, chronic disease, or other stressors, so our finding are not conclusively linked to opioid use. However, all individuals in the PWID group did inject heroin daily. Our study indicates a potential for opioid injection and/or its associated multi-drug use habits and lifestyle changes to influence epigenetic inheritance.

Introduction

The use and abuse of opioids has become a major public health threat in the last decade. There are familial patterns of drug abuse demonstrating an inherited propensity toward drug dependence, some of which could be inherited by specific genetic sequences (Donaldson et al., 2017). In animal models, consequences of drug use across multiple generations without variation in DNA allelic sequence have been observed (Yohn et al., 2015). In these models, there was also no offspring exposure to drugs, therefore implying that epigenetic changes due to parental drug abuse can be passed to future generations. Sperm cells express opioid receptors at both mRNA and protein levels, suggesting that they can be directly influenced by opioid use (Agirregoitia et al., 2006). Many studies of other outcomes, such as metabolic phenotypes and behavioral changes in both mice and in humans, also support the inheritance of epigenetic information, specifically via sperm (Chen et al., 2016).

Epigenetic information includes DNA methylation status, chromatin structure, and small regulatory RNA expression (Zhang et al., 2020b). Understanding how changes in these parameters affect cell function is relatively straightforward, but how epigenetic changes in sperm can influence future generations remains a puzzle.

In addition to their role in protein translation, tRNAs can be cleaved into shorter, biologically active regulatory RNA fragments called tRNA fragments (tRFs). The generation of tRFs is highly regulated and alterations in tRF patterns are associated with many biological processes, including viral infection, cancer, metabolic disorders, and exposure to stress or toxicants (Jiang and Yan, 2019). There is evidence that tRFs can impair protein translation non-specifically, acting as regulatory microRNAs (miRNAs), or repress retroelements in germ cells and embryos (Sharma et al., 2016). A recent study showed that a single i.m. injection of heroin in monkeys changed the expression pattern of genes involved in tRNA biosynthesis in the hippocampus, but the authors did not examine genetic changes in sperm (Choi et al., 2020).

Spermatocytes contain a substantial amount of tRFs, accounting for a much larger fraction of small RNAs (sRNAs) than miRNAs (Sharma et al., 2016). Recent studies in mice fed high-fat diets demonstrate an essential role for sperm tRFs in conferring heritable metabolic disorders on offspring (Chen et al., 2016; Zhang et al., 2018). This suggests that tRFs could be a major mechanism of epigenetic inheritance. As yet, the impact of opioid use on tRFs in germline cells has not been explored.

In this study, we conducted an unbiased small RNA-sequencing (RNA-Seq) profiling using semen from people who inject illicit drugs (PWID), who are assigned male at birth and compared it with semen from non-users. Our results provide evidence for a heritable impact of opioid use and lay the groundwork for further mechanistic studies on its transgenerational health implications.

Materials and methods

Biorepository of semen samples

The protocol for obtaining semen and blood samples from patients was approved by the Institutional Review Boards of the University of Washington in Seattle, WA, with informed consent signed by each donor (IRB number: 00000835). Subjects were heroin users recruited at two institutions in Seattle, WA: Evergreen Treatments Services (‘ETS’) and Harborview Medical Center (‘HMC’) Addiction Clinic.

The inclusion criteria were: assigned male at birth, being 18–55 years old, HIV-uninfected and active users of illicit (unprescribed) opioids. Subjects intending to seek medication-assisted therapy for heroin addiction with either methadone or buprenorphine/naloxone were approached by intake staff. Subjects self-reported drug and alcohol use through a questionnaire. Blood was also withdrawn for HIV and hepatitis B/C testing. Urine was obtained for qualitative testing of morphine, oxycodone, methadone, buprenorphine, cocaine, methamphetamine, and cannabis. Subjects with a negative result for heroin in the urine test, with clinical signs of infections (e.g. fever), with HIV infection or with active hepatitis B/C infection were excluded. Cryopreserved samples from non-drug-using men were collected from the University of Washington, Department of Urology, Andrology Clinic (directed by Dr John Amory), following institutional ethical recommendations. Samples were produced by masturbation after 2–3 days of sexual abstinence, collected in sterile urine collection cups containing cold Roswell Park Memorial Institute (RPMI) cell culture media and transported on ice to the laboratory within 3 h, as described previously (Vojtech et al., 2014). For both PWID and control cohorts, sperm banking (the process of collecting, freezing, and storing semen samples) was carried out using the same protocols and reagents, as described in the sections below.

Purification of mature spermatozoa

Semen samples were allowed to liquefy for 20 min at room temperature. Seminal plasma, containing exosomes, was separated from the cell fraction by centrifugation at 1000 ×g for 10 min at 20°C. Cell debris was removed by subsequent centrifugation at 2400×g for 30 min at 20°C followed by 0.45 and 0.22 μm syringe filtration (Millex HA, Darmstadt, Germany). The cell pellet was cryopreserved (both in PWID and control semen samples) in 1 ml of freezing medium containing 90% heat-inactivated fetal bovine serum (Nucleus Biologics, San Diego, CA, USA) and 10% dimethyl sulfoxide. Cryovials were placed in a –80°C freezer in a Mr. Frosty™ freezing device containing isopropanol; the approximate freezing rate was –1°C/min. For long-term storage, cryovials were transferred to vapor-phase nitrogen.

Mature spermatozoa were purified using an adapted Percoll gradient protocol (Claassens et al., 1998). Isotonic Percoll solution was prepared by mixing 9 ml of Percoll (Sigma-Aldrich, St. Louis, MO, USA) with 1 ml of 10× concentrated Ham’s F10 solution (MP Biomedicals, Solon, OH, USA). The bottom fraction of the Percoll gradient (90%) contained 1 ml of 90% isotonic Percoll mixed with 10% X-VIVO 15 (Lonza, Basel, Switzerland) or human tubal fluid ‘HTF’ culture media (in-house). The upper fraction (45%) contained 1 ml of 45% isotonic Percoll mixed with 55% X-VIVO 15 (Lonza, Basel, Switzerland) or HTF culture media (in-house). The gradient was prepared in 15 ml Falcon conical tubes. Semen samples from both groups (PWID and controls) were thawed quickly in a 37°C water bath and 1 ml of X-VIVO 15 was added dropwise to each cryovial. Thawed semen samples were washed in 25 ml of X-VIVO 15, centrifuged at 400×g for 10 min at 20°C, resuspended in ∼500 µl of X-VIVO15 and layered on top of the gradient and centrifuged again at 400×g for 20 min (brake off). After centrifugation, ∼1.5 ml of the supernatant was aspirated off and discarded. The pellet and remaining ∼1 ml of Percoll were washed with 14 ml of X-VIVO 15 or PBS and spun down (400×g for 10 min at 20°C). The final pellet was resuspended in 500 µl of X-VIVO 15 and used to prepare smears on coverslips for Wright and hematoxylin and eosin (H&E) staining and/or aliquoted and frozen at –80°C for later RNA purification. Both stainings were performed at the Histology and Imaging Core (HIC) of the University of Washington, South Lake Union campus (Seattle, WA, USA).

Isolation of seminal extracellular vesicles

Exosomes were purified from the entire cell- and debris-free seminal plasma (ranging from 0.9 to 4.0 ml in our donors) by ultracentrifugation over a sucrose cushion, as previously described (Vojtech et al., 2014). Up to 2.5 ml of supernatant was added to ultracentrifuge tubes and under-layered with 300 μl of a 20 mM Tris/30% sucrose/deuterium oxide (D₂O) cushion (pH 7.4) (Sigma, St. Louis, MO, USA). Samples were ultra-centrifuged at 100 000 ×g for 90 min at 4°C in an SW 50 swinging bucket rotor (Beckman, Brea, CA, USA). The upper layer was collected and ultracentrifuged again at 100 000 ×g for 14 h at 4°C over a 20 mM Tris/25% sucrose/D₂O cushion (pH 7.4). The upper layer of exosome-depleted seminal plasma was stored at −80°C. The 30% and 25% sucrose cushions containing the exosome fraction were pooled and brought to 15 ml with Dulbecco’s PBS (DPBS, GIBCO Thermo Fisher Scientific, Bothell, WA, USA). The exosomes were washed by centrifuging at 2400 ×g at 20°C through an Amicon Ultracel 100 kDa cellulose centrifugal filter (Millipore Sigma, St. Louis, MO, USA) with 10 ml of PBS and concentrated to a final volume of 425 μl–3.2 ml. Exosomes were stored at −80°C.

Exosome quantification and size determination

Concentration and size distribution of the isolated exosomes were measured by nanoparticle tracking analysis using a Nanosight NS300 instrument (Nanosight, Salisbury, UK), according to the manufacturer’s instructions. In brief, seminal extracellular vesicle (SEV) samples were vortexed and serially diluted to a final dilution of 1:8000–1:80 000 in filtered molecular grade H₂O. Blank-filtered H₂O was run as a negative control. Polystyrene latex standards (97 nm ± 3 nm) of known concentration (Life Technologies, Carlsbad, CA, USA) were analyzed along with the diluted exosome solution to validate the instrument. Each sample analysis was conducted for 60 s using Nanosight automatic analysis settings. Samples were evaluated in triplicate and concentration values were averaged.

RNA isolation

Total RNA for SEV and seminal plasma sequencing studies was isolated using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) with the following modification: samples were mixed with 5× volumes of Trizol LS (Invitrogen Thermo Fisher Scientific, Bothell, WA, USA), vortexed and incubated at room temperature for 5 min. Then, 0.2 volume of chloroform was added to each sample, and the miRNeasy protocol was then followed according to the manufacturer’s instructions.

For RNA studies on sperm cells, small and large RNA was isolated in one fraction (total RNA) using the NucleoSpin™ kit (Macherey-Nagel Kit, Duren, Germany, catalog: 740971), following the manufacturer’s instructions. The DNA was digested on the column by an RNAse-free recombinant DNase provided in the kit.

All RNA samples were quantified using fluorometry (Qubit™ RNA HS Assay Kit, Q32852 and Qubit™ Fluorometer 3.0, Life Technologies, Carlsbad, CA, USA). RNA size distributions were assessed by Agilent 4200 Tapestation at the Genomics Services core, Fred Hutchinson Cancer Research Center (FHCRC) in Seattle, WA, USA.

Quantification of sperm-specific and other markers by RT-qPCR assay

Total RNA (15–35 ng) was reverse transcribed with a High-Capacity cDNA Reverse Transcription Kit, using random primers (Applied Biosystems, Waltham, MA, USA). To quantify the levels of sperm-specific marker (sex determining region Y: SRY), we used a pre-designed quantitative PCR (qPCR) gene expression assay (Integrated DNA Technologies, Coralville, IA, USA, Ref Seq # NM_003140, exon location: 1-1). A GAPDH pre-designed assay (Ref Seq # NM_002046, exons: 2-3) was included to normalize the data. In both cases, the dye/quencher modifications were 6-FAM/ZEN/3IABkFQ and the primer:probe ratio was 2:1. To detect potential leukocyte contamination in sperm samples, we quantified the protein tyrosine phosphatase, receptor type, C (PTPRC marker, also known as CD45) (Ref Seq # NM_001267798, exon location: 3a-4). Reactions were set up in duplicates in 96-well plates with 1× PrimeTime Gene Expression Master Mix (Integrated DNA Technologies, Coralville, IA, USA) and run on a QuantStudio 5 Flex System (0.2 ml 96-well block, Applied Biosystems, Waltham, MA, USA), with the following conditions: denaturation/enzyme activation for 3 min at 95°C, 40 cycles of 15 s denaturation at 95°C and 60 s annealing/extension at 60°C followed by a hold at 4°C.

Transglutaminase 4 (TGM4) and protamine 2 (PRM2) assays were run using SYBR™ Green chemistry. The primers used were: TGM4-sense, 5′-TCCCTCAAGCCCACAGATA-3′; TGM4, antisense: 5′-ACTGCCTGTCCACTTGTATG-3′; PRM2-sense, 5′-GAAAGTCACCTGCCCAAGAA-3′ and PRM2-antisense, 5′-CTCAGATCTTGTGGGCTTCTC-3′. The PCR reactions were run using 1× PowerUp™ SYBR™ Green Master Mix (Applied Biosystems, Waltham, MA, USA) and 500 nM of each primer. The cycling conditions were: UDP activation for 2 min 50°C, Dual-Lock™ DNA polymerase activation for 2 min at 95°C, 40 cycles of 15 s denaturation at 95°C and 60 s annealing/extension at 60°C followed by a dissociation step (melt curve analysis). PCR reactions were run in duplicates. No-template control PCR reactions were run with each experiment and we observed no amplification. For both Taqman and SYBR Green assays, results were analyzed using Quant Studio v 1.4.1 (Applied Biosystems, Waltham, MA, USA).

tRF and tiRNA PCR array

To quantify the level of 185 prevalent tRFs or tRNA halves (tiRNAs) in sperm cells from PWID and controls, we used the nrStart™ tRF&tiRNA PCR array (catalog: AS-NR-002-1, ArrayStar, Rockville, MD, USA). The array covers all the nuclear anticodon isoacceptors catalogued in the genomic tRNA database (GtRNAdb), providing insights on both the isoacceptor and the specific isodecoder expression. The total RNA purified from sperm cells as described above was submitted to ArrayStar PCR service. An rtStar™ tRF&tiRNA Pretreatment Kit (Cat# AS-FS-005, Arraystar, Rockville, MD, USA) was used for RNA sample pretreatment; this process removed chemical modifications such as methylation and acylation from tRF and tiRNA. After phenol–chloroform extraction and ethanol precipitation, rtStar™ first-strand cDNA synthesis kit (3′ and 5′ adaptor) (Cat# AS-FS-003, Arraystar, Rockville, MD, USA) was used as to synthesize cDNA. Equivalent amounts of cDNA from each donor (in a 1:40 dilution) were mixed with Arraystar SYBR^® Green qPCR Master Mix (ROX+) (AS-MR-006-5, Arraystar, Rockville, MD, USA) and loaded onto the 384-well PCR array. As controls, we included: genomic DNA contamination control (GDC), spike-in, positive (PPC) and blank controls. The cycling conditions were: 95°C for 10 min, 40 cycles of 95°C for 10 s and 60°C for 1 min, followed by melting curve analysis. The ΔΔCt method was used to analyze the results. Raw cycle threshold (CT) values for each RNA fragment were normalized to the average of the CTs of the housekeeping genes SNORD43 and SNORD45. CT values over 35 or undetectable were set to 35. ΔΔCt was defined as: ΔCt (PWID) – ΔCt (controls). Fold-difference for each gene from controls to PWID was calculated as 2^–ΔΔCt. The P values were calculated based on two-tailed heteroscedastic Student’s t-test of the replicate 2^–ΔCt values for each gene in the control and PWID groups.

Phospho-RNA-seq library preparation and sequencing

To sequence sperm RNA fragments with alternate 5′ and 3′ phosphorylation states such as those found in tRFs, we applied phosphor-RNA-seq (Giraldez et al., 2019). To perform phospho-RNA-seq, sperm total RNA samples were pretreated with T4 polynucleotide kinase (Anza™ T4 PNK Kit, catalog: IVGN2304, ThermoFisher Scientific, Bothell, WA, USA) using an RNA input of 8.5 μl in a final reaction volume of 10 μl and incubated at 37°C for 30 min following the manufacturer’s instructions. After the enzymatic treatment, RNA samples were purified with sequential washes in silica columns according to the manufacturer’s instructions (Zymo, Irvine, CA, USA). Libraries for SEV RNA, seminal plasma RNA and PNK-treated sperm RNA were then prepared using the TruSeq small RNA kit according to the manufacturer’s instructions, with an input of 5 μl of RNA. Size selection was performed using pre-cast 5% acrylamide gels (Bio-Rad, Hercules, CA, USA). For SEV and plasma RNA, all products from 145 to 160 bp were selected. For PNK-treated sperm RNA, all products from 140 to 170 bp were included. Libraries were multiplexed and sequenced using the Illumina MiSeq (San Diego, CA, USA), specifying 50 bp single-end run. SEV extracellular RNA (exRNA) and seminal plasma exRNA samples were additionally sequenced on the Illumina NextSeq platform available at the University of Michigan Advanced Genomics Core (Ann Arbor, MI, USA).

Computational methods/small RNA-seq analysis

Reads were processed, aligned and quantified using the sRNAnalyzer sRNA-seq analysis pipeline (Wu et al., 2017). The sequences databases were downloaded from the sRNAnalyzer resources. The database and pipeline configuration files used are provided in the Supplementary Information. Read counts from the sRNAnalyzer pipeline output were loaded into R for differential expression analysis, filtering to require ≤1 mismatch and alignment in the sense (+) orientation. Differential expression was performed using DESeq2 and EdgeR using a blocking design to test for differences between opioid users and non-users while controlling for possible batch effects. Different RNA families (e.g. miRNA, piRNA, small nucleolar RNA (snoRNA), tRNA, etc.) were analyzed separately, since we could not assume that the overall composition of RNA classes remained the same across conditions. Expression differences with a false discovery rate of ≤ 0.05 were considered differentially expressed.

tRNA processing analysis

Reads processed by the sRNAnalyzer pipeline were re-aligned using SALMON (Patro et al., 2017) to generate BAM alignment files. Custom scripts were used to convert sRNAnalyzer database sequences into SALMON indices and annotation files. Analysis of tRNA read alignments were performed in R. Code for differential expression, and alignment analysis is available upon request.

tRF-Gly crystal digital PCR assays

Crystal digital PCR (dPCR) assays were used to validate sequencing data on the Naica^® system (Stilla Technologies, Villejuif, France). Total RNA was purified using a NucleoSpinTM kit (Macherey-Nagel Kit, Duren, Germany, catalog: 740971), following the manufacturer’s instructions. Genomic DNA was eliminated using the rDNase included in the kit.

cDNA synthesis was performed in a small scale using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Waltham, MA, USA) and adapting a protocol previously published (Bar et al., 2008). Two microliters of an aqueous solution containing total RNA (up to 2.4 ng) were reversed transcribed in a 5 µl reaction, containing each stem–loop specific primer at 0.75 μM (Table I), 1× RT buffer, 1× RNAse inhibitor, 1 mM of each dNTP and 16.5 IU of MultiScribe reverse transcriptase. The thermal cycler conditions were: 65°C for 1 min, 4°C for 1 min, 16°C for 30 min, 42°C for 30 min, 85°C for 5 min and hold at 4°C. Control reactions without reverse transcriptase were included to account for potential amplification from genomic DNA.

The sequences for the stem–loop RT primers as well as those for primers and probes (all from Integrated DNA Technologies, Coralville, IA, USA) are listed in Table I. In tRF-Gly short and tRF-Gly long Taqman assays, the dye/quencher modifications were 6-FAM/ZEN/3IABkFQ and the primer:probe ratio was 2:1. In each stem RT primer sequence, letters denoted in bold and italics indicate the nucleotides that form the seed in each RT reaction.

The reaction mixture in each tube contained 1× Perfecta Multiplex ToughMix (Quantabio, Beverly, MA, USA), 100 nM fluorescein, 1.25 μl of a 20× multiplex mix of primers (final concentration of 0.5 μM each) and Taqman probes (final concentration: 0.25 μM each) and 5 μl of diluted cDNA, for a final volume of 25 µl. To detect tRF-Gly long molecules, 1:15–1:30 dilutions of cDNA were used and 1:150–1:300 were used to detect short tRF-Gly copies. No-template (water) controls were included as contamination controls. The Sapphire chip (Stilla Technologies, Villejuif, France) was loaded into the Naica Geode thermocycler to compartmentalize the droplets and to perform the PCR reaction. The thermal conditions were: 95°C for 3 min, followed by 45 cycles of 95°C for 10 s and 60°C for 30 s. Data were acquired with Crystal Reader software and analyzed using Crystal Miner v2.3.5 (Stilla Technologies, Villejuif, France). In each reaction, at least 20 000 droplets were analyzed.

miR-100-5p crystal dPCR assay

Two microliters of an aqueous solution containing total RNA (up to 2.4 ng) were reverse-transcribed using the Taqman Advanced miRNA Assay (catalog: A25576, Applied Biosystems, Foster City, CA, USA), following the manufacturer’s protocol. This kit uses 3′ poly-A tailing and 5′ ligation of an adaptor sequence to extend the mature miRNAs present in the sample on each end prior to reverse transcription. Universal RT primers recognize the universal sequences present on both the 5′ and 3′ extended ends of the mature miRNAs.

The reaction mixture in each tube contained 1× Perfecta Multiplex ToughMix, 100 nM fluorescein, 1.25 μl of a 20× 5′-FAM hsa-miR-100-5p Taqman assay (assay: 478224_miR, catalog number: A25576, miRBase Accession #: MIMAT0000098, Applied Biosystems) and 5 μl of 1:30–1:300 diluted miR-Amp reaction, for a final volume of 25 µl. Crystal dPCR and data analysis were performed as described above.

Cross amplification qPCR control with synthetic long tRF-Gly motif

To evaluate cross amplification, a synthetic 53-ultramer RNA oligonucleotide encompassing the longest tRF Gly sequence was purchased from Integrated DNA Technologies (IDT, San Diego, CA, USA). The sequence of this synthetic construct was: 5′ GCAUGGGUGGUUCAGUGGUAGAAUUCUCGCCUGCCACGCGGGAGGCCCGGGUU 3′. Serial dilutions ranging from 100 μM to 1 pM were prepared in RNA-free DNA-free water. cDNA synthesis was carried out on a small scale, as explained above, using 2 μl of diluted RNA oligomer and either the short or long tRF Gly stem–loop primer.

We ran 6.75 μl of 1:30 diluted RT reactions using a qPCR system and the same primers and probes for long and short tRF-Gly molecules as described before. The master mix contained the following: 7.5 μl of 2× TaqMan Fast Advanced Master mix (Applied Biosystems, Foster City, CA, USA) and 0.75 μl of each 20× Taqman assay mix. We ran two technical replicates for each sample. Cycling conditions were as follows: initial denaturation and activation of Taq DNA polymerase at 95°C for 20 s and 40 cycles of 1 s of denaturation at 95°C and 20 s of annealing/extension at 60°C. Thermocycling was performed on a QuantStudio 5 (Applied Biosystems, Foster City, CA, USA).

Statistical analysis

Statistical analysis was carried out using the GraphPad Prism suite (Prism 9 for Windows 64-bit, v9.1.0(221), 15 March 2021, RRID: SCR_002798). Statistical tests include one- and two-tailed Mann–Whitney U-tests and two-tailed Student’s t-test. Differences were considered statistically different when P < 0.05.

Results

Standardization of sperm cell purification

We collected semen and blood samples from 13 PWID who were assigned male at birth, as well as from a group of 10 individuals not using drugs as controls. In our cohort of illicit drug users, almost no individuals used only a single drug (Supplementary Table SI). By both urinalysis and by self-reports, users of illicit opioids also used other drugs, including marijuana, cocaine, alcohol, methamphetamine, and benzodiazepines. However, all participants identified illicit opioids as their primary drug.

Using these samples, we investigated the hypothesis that opioid use changes the sRNA cargo carried by sperm. The experimental strategy is outlined in Fig. 1; the methods used for fractioning were similar to the ones we applied previously (Vojtech et al., 2014). For each donor semen sample, we separated and cryopreserved different fractions: the cellular portion/constituent of ejaculates, the seminal plasma and the SEVs (Fig. 1, Step 1).

An important potential confounding factor is the presence of cellular DNA from somatic cells in semen. Semen contains a complex mix of cells, including those derived from the prostate and the epididymis, but regarding transgenerational inheritance, only changes in mature, fertile sperm should be considered. We developed and validated a method to purify mature sperm from additional control (N = 6) frozen semen cell samples, using Percoll density gradient centrifugation (Fig. 1, Step 2). Pure, mature sperm cells were isolated by centrifugation over a Percoll gradient (90% and 45%), which separates them from somatic cells, defective and dead sperm, and other debris. We used Wright and H&E staining to monitor the efficiency of purification; a typical staining of the initial unfractionated semen sample and Percoll-purified fraction is shown in Fig. 2a and b. The final samples lacked mucoid cells or bacteria, and the number of contaminating epithelial cells was strongly reduced.

After standardizing the purification step for sperm cells, we isolated total RNA using commercially available silica columns (Fig. 1, Step 3). Using electrophoretic analysis (Fig. 2c), we found that this method of total RNA purification enriched for RNA fragments that ranged from 25 nucleotides to 5–6 kb. Mature-sperm RNA also lacked defined 18S and 28S rRNA peaks (Fig. 2c), in line with other reports (Ostermeier et al., 2005; Johnson et al., 2011). We next purified total RNA from all the semen samples of the biorepository following the Percoll gradient and silica-based method and applied a previous experimental pipeline for seminal plasma and SEV RNA (Vojtech et al., 2014) (Fig. 1, Step 4).

We performed gene expression analysis by qPCR to characterize the purity of the purified sperm cells. As expected after the Percoll gradient, the sperm-purified fraction was enriched in the sperm-specific transcript SRY, lacked the leukocyte pan-marker CD45 and showed a decrease in the prostate marker TGM4 (Fig. 2d). To conclusively test how efficient Percoll centrifugation decreased the proportion of contaminating cells, we set up additional experiments in which semen samples were spiked with increasing numbers of cervical epithelial cells and peripheral blood mononuclear cells, as surrogate markers for epithelial cells and leukocytes, respectively (Supplementary Fig. S1). Gene expression analysis revealed that upon Percoll purification, SRY and PRM2 mRNAs (sperm protamine 2 gene, highly expressed in sperm) (Steger et al., 2000) were selectively enriched, thus confirming the effectiveness of the purification method.

We recovered significantly higher amounts of RNA from control sperm (median: 54.0 ng, n = 10) than from PWID sperm (median: 16.4 ng, n = 13, P = 0.0167) (Fig. 2e); only in one donor from the first group was the concentration of RNA below the lower limit of detection, whereas four samples from PWID had undetectable amounts of RNA. This difference in yield agreed with the fact that the semen cellular pellets in PWID were smaller than those in non-drug using controls.

Altered cleavage site in tRF-Gly among semen samples from PWIDs

We applied a method called ‘phospho‐RNA‐seq’ (Akat et al., 2019; Giraldez et al., 2019) to profile global sRNAs in semen samples from nine PWID and nine controls (Fig. 1, Step 5); samples with undetectable RNA concentrations were not subjected to this analysis. Typical RNAseq library preparation with adapter ligation relies on 5′-monophosphate and 3′ hydroxyl groups. Because cleavage products of RNAs generated by many ribonucleases have different 5′ and 3′ ends, we utilized phospho-RNA-seq. In phospho-RNA-seq, RNA is treated with T4 PNK to allow adapter ligation to detect RNAs that are missed by other RNAseq library preparation methods (Giraldez et al., 2019). We size-selected to sequence sRNAs ∼20–50 nucleotides long. After applying a high-stringency bioinformatic pipeline (described in the Materials and methods section), we found that among the sRNAs present in sperm cells, the most abundant biotypes were from protein-coding genes, tRNAs and from retained introns (Supplementary Fig. S2a).

We focused our analysis on tRFs because they are the most abundantly expressed sRNAs in semen (Hua et al., 2019) and they have already been shown to be involved in transgenerational inheritance. The tRNA isotypes that we detected with the highest frequency were: tRNA-Val, tRNA-Glu and tRNA-Gly (Supplementary Fig. S2b).

Next, we analyzed different patterns of tRF size distributions. For tRF-Gly-GCC, we found differential cleavage patterns in men using opioids compared to non-drug users. We did not see size-distribution differences for other tRFs in our dataset. Prior studies have implicated tRF-Gly-GCC as mechanistically important for epigenetic inheritance (Chen et al., 2016; Sharma et al., 2016; Zhang et al., 2018). Figure 3a depicts the aggregate data derived from all the samples in each group. We observed a bimodal distribution in both the 3′ end position and the overall length of fragments. The data revealed that these tRF fragments are predominantly truncated at position ∼32 in control and at position ∼51 in PWID samples.

When reads were mapped to the tRNA structure, we found that over 90% of reads in non-drug using controls map to shorter tRFs cleaved in the anticodon loop (Fig. 3b), while in PWID only 45% of reads mapped to this short tRF. In contrast, only 4.1% of reads in controls mapped to a longer tRF species cleaved in the proximity of the TψC loop of the tRNA, while in PWID men, 45.6% mapped to this long tRF (Fig. 3b).

Validation of phospho-RNA-seq results by crystal dPCR

As phospho-RNA-seq revealed specific differences in tRNA-Gly cleavage patterns and because tRNA-Gly-GCC fragments have been involved in regulation of gene expression (Sharma et al., 2016; Boskovic et al., 2020), we used a PCR approach to confirm our next-generation sequencing findings. We measured both the absolute concentration and the cleavage patterns of these tRF-Gly-GCC in our samples using crystal dPCR (Fig. 1, Step 6).

For the tRF Gly, we designed specific RT primers with complementary regions for either a long fragment (53 nucleotides) and or a short tRF fragment (32 nucleotides) (from references Schuster et al., 2016b; Nätt et al., 2019), plus a standard stem–loop sequence to extend the 3′ end of the complementary DNA for quantitative dPCR amplification. We designed and tested Taqman assays to quantify both the short and long tRF fragments; agarose gel electrophoresis confirmed the presence of 72 bp and 93 bp amplification products that matched the expected sizes (not shown). Owing to the high sensitivity of dPCR, we were able to detect both tRF forms in samples that had concentrations of RNA below the lower limit of detection.

Since the ‘short-RT’ primer contained a seven-nucleotide seed sequence that was common to both the 32-nucleotide tRF as well as to longer tRFs that extended beyond the anticodon loop, we reasoned that part of the amplification obtained using the short Taqman assay would originate from cross-hybridization of the RT-primer with longer tRFs. To correct for this cross-priming artifact, we designed a synthetic 51-ribonucleotide long tRF, which spanned the most frequent tRF found among PWIDs (Fig. 3a, first sequence). We did serial dilutions of this template and performed RT reactions using either the 53-nucleotide (‘long-RT’) primer or the 32-nucleotide (‘short-RT’) primer. We ran both long and short qPCR Taqman assays on corresponding long-RT and short-RT reactions. Supplementary Fig. S3a and b shows the results of these experiments. As expected, the data show that the efficiency of amplification of the synthetic 51-nucleotide long tRFs was much higher for the long RT-primer/long Taqman assay than for the short RT-primer/short Taqman assay; for a given amount of template, the qPCR curve had a lower Ct for the first pair (Supplementary Fig. S3a) than for the second pair (Supplementary Fig. S3b). The quantitative data are shown in Supplementary Fig. S3c; the difference in Ct ranged from 4 to 5 cycles, which corresponds to ∼16- and 32-fold differences in amplification efficiency. Based on these data, we determined that the qPCR assay to measure short tRF molecules can detect also up to ∼6% of the total number of longer tRFs; we used this correction factor to adjust dPCR data.

As shown in Fig. 4a, we observed that semen samples from PWID had a significantly lower concentration of long tRF Gly fragments compared to controls (median PWID: 1626 absolute copy number, n = 13 versus median controls: 37 111 absolute copy number n = 10; P = 0.0128). Likewise, there was a significantly lower concentration of short tRF Gly fragments in PWID than controls (median PWID: 8313, n = 13 versus median controls: 277 273, n = 10; P = 0.0001). For both groups (PWID and controls), the cDNA reactions were carried out with equivalent amounts of RNA. When we calculated the ratio of long/short tRF in each donor, the dPCR results confirmed the RNA-seq data: in the sperm of people who inject heroin, there is a relative enrichment of tRF-Gly long fragments (median PWID = 0.2305, n = 13 versus median control: 0.1602, n = 10, P = 0.0128, Fig. 4c). In summary, we found that both long and short tRF-Gly were markedly decreased in PWIDs compared to controls, with a significant relative enrichment of the long isoforms in PWID over controls.

We also determined the expression of miR100-5p, a marker that has been used to normalize qRT-PCR datasets from semen samples (Corral-Vazquez et al., 2017), and measured it at much lower concentrations in PWIDs than controls (median 900 molecules per µl cDNA in 12 PWID versus 7183 in 10 controls, P = 0.0249; Fig. 4d); this relative downregulation of miR100-5p in PWID made it unsuitable as a normalizer. Finally, we found that the dPCR results correlated very well with those obtained by phospho-RNA-seq. We illustrate this concordance by showing the results from three semen donors in Supplementary Fig. S4a–c. For example, in Supplementary Fig. S4a, 168 of the reads had sequences that extend beyond the anticodon loop (larger than 32 nucleotides), whereas only 22 were shorter; the dPCR long:short ratio for this PWID was 2.72. Conversely, in the control shown in Supplementary Fig. S4c, only 1 read had 53 nucleotides, 41 were <32 nucleotides long; the dPCR long:short ratio for this control individual was 0.14.

Differentially expressed tRFs in PWID

tRNA-derived non-coding RNA fragments fall largely into two main categories: tiRNAs (or tRNA halves, 30–50 nucleotides long) and tRFs (tRNA-derived fragments, 16–28 nucleotides long); each group is generated through precise biosynthetic pathways, contains a certain nucleotide composition and elicits specialized functions (Anderson and Ivanov, 2014). Using phospho-RNA-seq, we found that the tRF cleavage products from Gly-GCC tRNA were markedly different between spermatocytes from PWID compared to non-drug using controls. Because tRFs are heavily post-transcriptionally modified (Lyons et al., 2018) and these modifications may hinder the sensitivity of detection of sRNAs by RNA-seq methods (Giraldez et al., 2019), we used a complementary approach that relies on qPCR. For this purpose, we identified tRF fragments with differential expression in PWID compared to controls using the human tRF and tiRNA PCR array, which includes a demethylation step to ensure continuing RT during cDNA synthesis. This methodology has been successfully used to delineate differentially expressed tRFs in cancer (Li et al., 2018; Kong et al., 2020; Zhang et al., 2020a). We quantified the expression of 185 tRFs and tiRNAs, using the levels of SNORD43 and SNORD95 as endogenous normalization control genes.

The results of the quantification are shown in Supplementary Table SII and Fig. 5. The fraction of fragments that were undetectable (i.e. assays with cycle thresholds above 35) was comparable between PWIDs and controls (average of tRFs with absent calls: 69% versus 65%, respectively). Overall, the majority of genes were either down regulated in opioid-users or expressed at similar levels in both groups. We found three fragments to be significantly down-regulated in PWID (> 2-fold decrease with P < 0.05): tiRNA-5031-PheGAA, 3′ tiR_056_ValTAC (mt) and the transcript 3026/27/28A (Fig. 5a). The 3026/27/28A transcript is classified as a tRF-3 fragment (Supplementary Table SIII), the potential precursors are tRF-Gly-GCC, Gly-CCC and Gly-TCC; it ends at the 3′ end of mature tRNAs and starts by cleavage within the T-loop (Kumar et al., 2015). We found four more fragments whose expression was lower in PWID (> 2-fold), but the difference did not reach statistical significance (Fig. 5b): 3'tiR_007_GluTTC (n), transcript 3020/21A, 3'tiR_082_ThrTGT (mt) and RNU6. Only one gene trended toward upregulation in PWID compared to controls (tiRNA-5034-ValCAC-3, 2.19-fold up regulation), though the difference was not significant (P = 0.2168; Supplementary Table SII).

Chronic heroin use changes the exRNA cargo of semen

In our previous work, we reported that human semen contains trillions of SEV carrying exRNA molecules known to target immune-related mRNA (Vojtech et al., 2014). Here, we investigated whether chronic heroin use changed the composition of exRNA. For this purpose, we purified and sequenced sRNAs from SEVs and from seminal plasma from PWID and controls. We applied the sRNAnalyzer pipeline (Wu et al., 2017) to map the reads to sRNA databases, and used two independent algorithms (EdgeR (Robinson et al., 2010) and DESeq2 (Love et al., 2014)) to conduct differential expression analysis.

Figure 6 shows the results of differentially expressed markers detected by both algorithms, with a false discovery rate < 0.05. The results obtained from independent runs were consistent. Both algorithms agreed for 53.8% of differentially expressed genes in SEV and 59.1% in seminal plasma. With a few exceptions (E3, miR181d-5p, miR-429-3p), we observed a trend in both sample types for genes to be expressed at lower levels in PWIDs than in controls. Four miRNAs (miR103a-1-3p, miR-152-3p, miR30-a-3p and miR-361-5p) were detected at significantly lower levels in PWIDs in both sample types.

We ran a similar analysis, asking what genes were identified as being differentially expressed by only one of the two algorithms. The results are shown in Supplementary Fig. S5. In SEV, we found differential expression of a group of snoRNAs that includes, among others, ACA14a, U19, and U3-3. Again here, the tendency was to find genes expressed less in PWIDs than in controls.

Finally, we found two tRNAs differentially expressed in seminal plasma in PWID compared to controls: tRNA-Arg-CGG by both algorithms (Fig. 6), and tRNA-Val-GTG by only one algorithm (Supplementary Fig. S5). The direction again showed that there were fewer tRNAs in PWIDs. In summary, there was not a dramatic shift in exRNA cargo in semen from PWID compared to non-drug using controls, though we identified a few differentially expressed exRNAs.

Discussion

The opioid epidemic has broad public health reverberations across society, and it is important to consider how it will impact future generations. In this study, we aimed to investigate how opioid use altered a key epigenetic mechanism: the composition of sRNA in human sperm cells, seminal plasma, and SEVs, making use of our existing biorepository of semen samples from PWID and non-drug-using controls. We made three main observations. First, we found that a tRF from isoacceptor Gly-GCC tRNA had differential cleavage patterns in PWID compared to controls. Second, in SEVs, we report differential expression of a group of snoRNAs that includes, among others, ACA14a, U19, and U3-3. Third, we described a common set of four miRNAs (miR103a-1-3p, miR-152-3p, miR30-a-3p, and miR-361-5p) that were expressed at lower levels in both SEVs and seminal plasma of PWIDs. Altogether, we described an altered cargo of sRNAs in human sperm associated with chronic opioid use.

Importantly, however, a limitation of our study was that poly-drug use among PWID is extremely common; in fact, we could not enroll individuals who used illicit opioids exclusively. Although we assessed urine for other drugs, we were not able to adjust for their use because it was too varied, including changing combinations of benzodiazepines, marijuana, cocaine, and methamphetamines. Furthermore, we could not match PWID and non-drug users in terms of diet, environmental exposures, chronic diseases or other stressors, all factors that could influence epigenetic markers. Thus, we were not able to conclusively link our findings specifically to opioid use. However, all individuals in the PWID group did inject heroin daily, and this fact, and/or its associated factors, thus appeared to cause the epigenetic changes we observed.

Long-term opioid exposure has been shown to induce epigenetic modifications in the brain such as changes in histone acetylation, histone methylation, DNA methylation and non-coding RNA expression (reviewed in Browne et al. (2020)). Nevertheless, few studies have examined the chronic effect of opioids on the epigenetic status of human gametes. Endogenous opioid peptides play an important role in male reproductive function, and sperm express functional delta, kappa and mu opioid receptors (Agirregoitia et al., 2006; Subirán et al., 2011). Addition of synthetic opioids to sperm ex-vivo causes reduced motility (Xu et al., 2013), implying that sperm are directly affected by the presence of opioids. In rats, chronic paternal morphine exposure adversely affected fertility (Cicero et al., 2002). In humans, opioid use is associated with an increase in global DNA methylation in lymphocytes, and increased methylation of CpG sites in the mu opioid receptor in sperm and lymphocytes (Nielsen et al., 2009; Chorbov et al., 2011; Ebrahimi et al., 2018), but the global effects of opioid use in sperm, particularly on the sRNA transcriptome, have not yet been investigated. Furthermore, no study has closely examined whether heroin exposure or consumption of other synthetic illicit opioids, such as fentanyl, could lead to differential biogenesis/metabolism of tRF fragments.

Many reports indicate that small non-coding RNAs in sperm can sense environmental cues (Rodgers et al., 2013; Gapp et al., 2014; Rodgers et al., 2015; de Castro Barbosa et al., 2016; Short et al., 2016, 2017) and, in recent years, tRFs have emerged as new epigenetic markers (Chen et al., 2016; Sharma, 2019; Shen et al., 2019). The first insights into the biology of tRFs came from animal studies. It has been well documented that sperm-derived RNAs play a role in early embryo development and can mediate transgenerational epigenetic inheritance (Peng et al., 2012). Small-RNA sequencing studies have revealed that sperm are remarkably enriched for 5′ tRFs cleaved at very specific sites, mapping primarily to only a few tRNA isoacceptor types (Peng et al., 2012; Sharma et al., 2016; Nätt et al., 2019). The demonstrated high degree of regulation in generation of tRFs suggests that specifically cleaved versions might play important biological roles. These 5′ tRFs are abundant in mouse sperm heads, indicating they can be delivered to oocytes (Peng et al., 2012). In a few recent landmark studies, 5′ tRFs have been implicated in transgenerational effects linked to diet. Injecting oocytes with tRFs isolated from the sperm of male mice fed high-fat diets caused metabolic disorders and altered metabolic gene expression in the offspring (Chen et al., 2016). This has also been shown with males exposed to high-fat diets in utero (via feeding pregnant mice): these males had increased expression of 5′ tRFs in sperm. Injection of these 5′ tRFs into a third generation of healthy oocytes led to obesogenic traits in the pups (Sarker et al., 2019). The mechanism of this effect is still under investigation, but evidence suggests that tRFs can act as post-transcriptional gene regulators and may alter metabolic gene expression and cause a transcriptional cascade change in embryos that reprograms gene expression in key metabolic tissues (Chen et al., 2016; Magee and Rigoutsos, 2020; Xie et al., 2020). In addition, tRFs have been linked to control of endogenous retrotransposons (Sharma et al., 2016; Schorn et al., 2017). Control of retrotransposons during development is particularly important, since they are normally controlled by methylation, which is erased during reprogramming to a pluripotent state.

In humans, a subset of tRFs is also highly enriched in sperm (Hua et al., 2019). We published an investigation of the small regulatory RNA cargo carried by extracellular vesicles in semen and found that they are also highly enriched for a specific subset of tRFs (Vojtech et al., 2014). These vesicles can fuse with sperm to deliver cargo, so they are also important to consider for sperm function (Ronquist, 2012). So far, there are only a few studies investigating differences in tRF patterns of sperm from men with different environmental exposures. One study found that 10 tRFs were differentially expressed between sperm samples from men in IVF clinics who had high rates of good quality embryos compared to men with low rates of good quality embryos (Hua et al., 2019), demonstrating their utility as sperm quality markers for IVF. A recent analysis of 15 men in a dietary intervention study found significant upregulation of tRFs in sperm following just 1 week of a high-sugar diet, compared to sperm tRF profiles during the healthy diet phase of the study (Nätt et al., 2019). These changes are similar to those seen in mouse studies and correlated with observed changes in sperm motility. The responsiveness of the sperm tRF repertoire to other exposures has not yet been investigated. Though we were not able to assess the diets of individuals in our cohorts, our study tentatively suggests that drug use, like diet, may be another environmental factor with the potential to shift tRF patterns and potentially mediate transgenerational inheritance. Arguably, injection drug use is a major stressor that also tends to be prolonged and increases in severity over time.

Here, we report that the sRNA reads from mature sperm cells were mostly enriched in tRNAs. Among them, the tRNAs that we mapped with higher frequency from the phospho-RNA-seq libraries were tRNA-Val, tRNA-Glu, tRNA-Gly, and tRNA-Lys, which were also found in other reports (Schuster et al., 2016b; Nätt et al., 2019). Next, we established a method to quantitatively distinguish between differently cleaved tRFs from the same parent tRF, using sequence-specific RT primers and dPCR. We found that the most dramatic difference in cleavage pattern was observed in the tRF from the isoacceptor Gly-GCC tRNA, and dPCR experiments confirmed the results obtained from phospho-small-RNA-seq. The longest read had 53 nucleotides and its 3′ end maps to the TΨC loop, which makes it unclassifiable according to the current nomenclature (Magee and Rigoutsos, 2020). The shift toward longer tRFs has also been found for other stimuli. For example, exposure to environmental toxins has been shown to change the average size of tRFs in rats (Stermer et al., 2019). Notably, a shift toward higher expression of longer tRF species with chronic alcohol exposure was also observed in mouse sperm (Rompala et al., 2018). When cells are exposed to certain stressors (e.g. lack of amino acids, oxidative stress, UV, viral infections), several enzymes, such as Dicer, Angiogenin and RNAase Z, respond and result in the biogenesis of tRFs (Magee and Rigoutsos, 2020; Xie et al., 2020). It is still unknown which enzymatic complexes generate the longer tRF-Gly-GCC fragments we found in our datasets. Intriguingly, our longest reads, 50–53 nucleotides long, had 3′ ends close to or within the T-loop. Similar T loop cleavage tRNA-Gly-GCC products of ∼55–60 nucleotides were reported in epididymis but not in testis samples (Sharma et al., 2016), and also in mouse brain (gene ID: chr1.trna706-GlyGCC (Flores et al., 2017)). These results also agree with those by Nätt et al., who found that sperm tRFs that changed in response to a high-sugar diet had a TTCGA sequence in their T-loop (Nätt et al., 2019). This motif is a target of pseudouridine synthase 7 (PUS7), which catalyzes pseudouridylation of uridines (Guzzi et al., 2018; Guzzi and Bellodi, 2020). Notably, the isoacceptor Gly-GCC has been shown to experimentally bind to PUS7 (Guzzi et al., 2018). This suggests that different stimulations, such as opioid exposure or nutritional flux, may drive a common mechanism that results in the selective cleavage of the T-loop, likely accompanied by differential pseudouridylation in defined sets of tRFs.

The specific consequences of this imbalance in tRF Gly-GCC isoforms in the fertilized embryo is unknown. Based on the tRF activity targeting mRNAs at their 5′-untranslated regions (Schuster et al., 2016a,b), it is conceivable that they may affect gene expression networks. They can also modulate other molecular processes such as cap-independent translation, localization of mRNAs into stress granules and others (Guzzi and Bellodi, 2020; Magee and Rigoutsos, 2020; Xie et al., 2020). Mechanistically, tRF-Gly specifically has been reported as a master regulator of endogenous retroviral transposons and other non-coding RNA elements (snoRNAs, small Cajal body-specific RNAs and small nuclear RNAs (Sharma et al., 2016; Boskovic et al., 2020). One of these downstream effectors, the U7 noncoding RNA, plays a role in processing histone pre-mRNAs. By controlling U7 levels, tRF-Gly may indirectly influence histone mRNA and protein levels, thereby affecting gene activity and chromatin compaction (Boskovic et al., 2020). Therefore, tRF-Gly-GCC may have an impact on the epigenome of embryonic stem cells and the preimplantation embryo.

In addition to profiling tRFs in mature sperm cells by phospho-RNA-seq and confirming the findings by crystal dPCR, we applied a PCR array to quantify differentially expressed tRF fragments in seven PWIDs and eight controls. This platform interrogates 185 recently identified tiRNAs (5′ and 3′) and fragments (tRF-1, tRF-3, and tRF-5), which have been compiled in recent databases and publications (Wang et al., 2013; Selitsky et al., 2015; Olvedy et al., 2016). We reasoned that this methodology could provide additional information because it combines the high sensitivity of PCR and it employs a pretreatment step that improves the conversion of RNA into cDNA. In agreement with the results derived from phospho-RNA-seq (Supplementary Fig. 2b) and crystal dPCR (Fig. 4), we found that the relative mRNA expression levels of most of the sRNA fragments screened was lower in PWID compared to controls. We report that three fragments had statistically significantly lower expression (> 2-fold, P < 0.05) in PWIDs: tiRNA-5031-PheGAA, 3'tiR_056_ValTAC (mt) and transcript 3026/27/28a. Phospho-RNA-seq also identified tRNA-Val as one of the most abundant and differentially expressed tRNA (Supplementary Fig. S2b). The transcript 3026/27/28a is classified as a tRF-3 fragment since it ends at the 3′ end of mature tRNAs and starts by cleavage within the T-loop; the potential precursors are tRF-Gly-GCC, Gly-CCC and Gly-TCC. Taken together with the results of dPCR shown in Fig. 4b, the data suggest that the level of expression of short tRF-Gly-GCC fragments (both the 5′ half and the one derived from the 3′ end) is lower in PWID. Only by using complementary approaches that include dPCR were we able to reveal the enrichment in long tRFs in this population.

From the qPCR array results, we also found that the transcript RNU6 was expressed almost to ∼1/3 level in sperm cells from PWID compared to controls, but the difference did not reach statistical significance (Supplementary Table SII). RNU 6 (U6 snRNA) is a non-coding small nuclear RNA that has been shown to play a role in splicing (Didychuk et al., 2018). Mechanistically, sRNAs, such as U7 snRNA, have been implicated in the regulation of histone proteins by tRF fragments (Boskovic et al., 2020). It can be posited that U6 and other snRNAs in sperm cells of PWIDs could be mediators of epigenomic effects of tRFs.

In the last part of this study, we investigated whether opioid use in PWID changes the composition of exRNA in SEVs and in seminal plasma. Although seminal exosomes (SEVs) are generated by all cell types present in the male genital tract, the prostate epithelium is the source of most SEVs (Renneberg et al., 1997), and opioids are present in the prostate gland (Hosseini et al., 2012). Many reports describe a pronounced change in cellular miRNA following exposure to opioids (reviewed in Rodríguez (2012)). In line with the hypothesis that tRF can influence the transcription of small non-coding RNAs (Boskovic et al., 2020), we found that the snoRNAs ACA14a and U19 were expressed at lower levels in SEVs from PWIDs than from non-drug-using controls. Interestingly, ACA14a belongs to the class of snoRNAs containing a H/ACA box that participates in the conversion of uridine to pseudouridine (Ψ) (Kiss et al., 2004). U19 (SNORA74B) belongs to the H/ACA class of snoRNAs that form ribonucleoprotein complexes, also with pseudouridylation activity (Courtes et al., 2014). It is tempting to speculate that the dysregulation of these snoRNAs in epididymal extracellular vesicles from PWIDs can cause changes in the pseudouridylation of tRFs in developing sperm cells from the testis, affecting stability and function (Charette and Gray, 2000). Similar intergenerational transmission via extracellular vesicles has been reported before (Chan et al., 2020).

Another limitation of this study is that we do not know the fertility status of any of the men in our repository, neither men who inject illicit opioids nor the non-drug using controls. It is well known that opioid use causes low testosterone, low sperm counts, and hypogonadism (Drobnis and Nangia, 2017). We have spent considerable effort developing a well-validated method to purify live mature sperm from semen samples to ensure we are comparing, as much as possible, those sperm cells with the greatest chance of actually fertilizing an oocyte. In addition, for our samples, we do not have access to other potentially confounding factors (smoking status, diet, obesity, days of abstinence, history of epididymitis or other scrotal infections, use of antibiotics, sex hormonal supplements or chronic medications, etc.). Owing to the difficulty of recruiting men using opioids with demonstrated fertility or infertility, we consider our cohort and our method of sperm purification to be a good achievable comparison for pilot studies to understand basic epigenetic changes mediated by opioid use. The unequivocal effects of specific illicit drugs, such as heroin or fentanyl, on the repertoire of tRFs (and potentially other epigenetic parameters) will have to be confirmed using in vivo models such as mice or non-human primates.

The findings in this and follow-up studies could help answer several questions of importance for individual and societal family planning:

• Can therapeutic use of opioid pain medications (e.g. hydrocodone) result in epigenetic changes in sperm and thereby negatively affect later generations?

• Are opioid-induced epimutations reversible when the drug is no longer used and, if so, what is the timeframe for partial and/or complete reversion?

• Can opioid substitution therapies (OST), such as methadone or buprenorphine, reverse epimutations that are induced by illicit opioids?

Our findings emphasize the significance of further mechanistic studies into how these tRF fragments effect the transcriptome/epigenome of human embryonic cells. Future experiments treating epididymal cells with illicit opioids/OSTs may reveal their effects on the sRNA cargo of the cells’ secreted vesicles and how uptake of these vesicles by maturing sperm cells may, in turn, lead to changes in the fertilized embryo. Examining chromatin or histone modifications and changes in DNA methylation status in sperm from PWID may reveal additional mechanisms of potential epigenetic changes resulting from drug use. Together, such mechanistic studies are important to provide guidance for choosing the best therapeutic options to treat PWID, with the specific goal of limiting the transgenerational impact of opioid use.

Supplementary data

Supplementary data are available at Molecular Human Reproduction online.

Data availability

The data underlying this article are available in the article and in its online supplementary material. Raw data will be shared upon reasonable request to the corresponding author.

Acknowledgements

The authors want to acknowledge Dr John Amory from the University of Washington, Andrology Clinic, for providing some semen samples from non-drug using men. We want to also thank all semen donors who contributed to this study.

Authors’ roles

G.G., M.T, F.H., and L.V. conceived the original idea and designed the experiments. G.G, R.M.S., E.S., Y.K., C.L., and L.V. conducted the experiments. G.G, R.M.S., E.S., M.T., F.H., and L.V. analyzed and interpreted the data. G.G, L.V., and F.H. wrote the manuscript. R.M.S., E.S., C.L., and M.T. provided critical discussion, reviewed, and revised the manuscripts. All authors provided final approval off the manuscript prior to submission.

Funding

G.G. received funds from the University of Washington Alcohol & Drug Abuse Institute (ADAI) to cover the expenses of this project. This study was also partially funded by grant R01 DA040386-05 (National Institute of Health, National Institute on Drug Abuse, NIDA) awarded to F.H. and by NIH/NHLBI (National Institute of Health, National Heart, Lung and Blood Institute) grant U01-HL126499 to M.T. This publication was supported by the National Center For Advancing Translational Sciences of the National Institutes of Health under Award Number KL2 TR002317. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Conflict of interest

All authors declare that they have no conflicts of interest.

References