Evidence of transposon-derived primary and secondary piRNAs in fetal testis

We next looked for evidence of piRNAs across developmental timepoints spanning GW13-22. Small RNA sequencing was performed on bulk fetal testis tissue from eight samples at a depth of 35-40 million reads and processed according to our custom pipeline (Fig. S2). Trimmed 18-40 nt reads were first aligned to a curated list of human TE consensus sequences obtained from the GIRI Repbase database (Bao et al., 2015), and the remainder were aligned to the human genome (hg38). An overall alignment rate of >97% was obtained across all samples (Table S1). At all timepoints, the majority of reads (~75%) corresponded to miRNAs, with smaller fractions derived from tRNAs, snoRNAs and protein-coding transcripts (Fig. 2A, left panel; Fig. S3A; Table S1). Strikingly, fewer than 1% of reads for all samples mapped to Repbase, and these resolved into two size classes: 27 nt and 22 nt (Fig. S3B). Reads comprising the 22 nt peak showed extremely low sequence diversity, with the vast majority mapping to the L2C retrotransposon (data not shown). The 27 nt peak was highest at GW20, and likely corresponds to putative TE-derived piRNAs, which range from 26-32 nt in other species. As our small RNA populations were derived from bulk fetal testis tissue without cell type or biochemical enrichment, we focused on the 26-32 nt small reads for the remainder of the analysis (Fig. S2). This increased the fraction of Repbase-mapping reads to ~6% in the 26-32 nt RNAs but not in 18-22 nt small RNAs, which were mostly miRNAs (Fig. 2A, middle and right panels). The number of Repbase-mapping reads was dynamic over development, with a maximum at GW20-21 (Fig. 2B, Table S1), corresponding to the window of HIWI2 nuclear translocation. Hallmarks of pre-pachytene piRNAs were identified in these 26-32 nt Repbase-mapping reads, including a bias toward 1st position uridine (Fig. 2C, Fig. S3C) – a key signature of primary piRNA biogenesis (Iwasaki et al., 2015). Additionally, an adenine bias at the 10th position of mapped 26-32 nt reads at all timepoints with the exception of the earliest, GW13A (Fig. 2C, Fig. S3D), suggests active secondary piRNA biogenesis, also known as ‘ping-pong’ amplification (Czech and Hannon, 2016). Corroborating the specificity of Repbase-mapping reads, nucleotide distribution analysis on miRNA-annotated reads demonstrated a reduced, yet still apparent, 5′ U bias (Fig. S3E), similar to miRNAs in fetal ovary and adult testis (Williams et al., 2015), and no adenine bias at position 10 (Fig. S3F). Another signature of ping-pong biogenesis, ten nucleotide overlaps between reverse complementary piRNAs, was detected in all samples to varying extents, with the exception of GW13A, and most strongly at GW20 (Fig. 2D). Together, these analyses identify small RNAs in human fetal testes with defining characteristics of pre-pachytene piRNAs, although at much lower abundance than found in mice or human fetal ovaries (Roovers et al., 2015; Vasiliauskaitė et al., 2017; Williams et al., 2015).

Open in a separate window

Fig. 2.

Identification of transposon derived piRNAs in the human fetal testis. (A) Relative percentages of different RNA subtypes present in gestational week (GW) 20 sample after size filtering for 18-40 nt, 26-32 nt and 18-22 nt small RNAs. Values in bold italic highlight Repbase-mapped reads, which include putative piRNAs. (B) Normalized number of reads (26-32 nt) mapped to human transposon consensus sequences in Repbase database across all samples. (C) Relative nucleotide distribution at a given read position of pooled Repbase-mapping reads (26-32 nt) at GW20. (D) Ping-pong overlap signatures calculated on pooled Repbase-mapping reads (26-32 nt). (E) Top transposable elements mapped to by Repbase-mapping reads (26-32 nt) at GW20. (F) Distribution of read sizes mapping to the L1HS transposon across all samples. Positive and negative values represent reads that map to the positive and negative strand, respectively.

Among the Repbase-mapping reads, the L1 family was the most abundantly represented class of repeat elements at GW20 (Fig. 2E). This result is consistent with the derivation of piRNAs from actively transcribed TEs, given that L1 elements are known to be active in humans (Hancks and Kazazian, 2012). The evolutionarily youngest L1 family (Smit et al., 1995), L1HS, was the most highly represented element in all samples (with the exception of GW13A), followed by ‘L1’ elements, which here represent the broad primate-specific L1PA clade (Kojima, 2018) (Table S2). Reads were mapped to both the forward and reverse strand of L1HS, particularly at later stages (Fig. 2F), consistent with a progressive ping-pong amplification process. Analysis of other highly represented classes of TEs revealed different strand-specific biases, yet consistently showed GW20 as the sample with the largest number of relative mapping reads (Fig. S3G-I). We conclude that pre-pachytene piRNAs are produced and amplified in the human fetal testes, with a large fraction derived from TEs, which is similar to the developing mouse testis (Aravin et al., 2007, 2008).

L1 transposons are expressed in AGCs coordinately with the transposon repression network at the single cell level

The prevalence of L1-derived piRNAs and the possibility of piRNA-mediated nuclear silencing prompted a closer examination of the relationship between transposons and the transposon-repression network in an unbiased, population-wide manner. We took advantage of a recently published single cell RNA-seq dataset from human fetal testis tissue (Li et al., 2017) to assess the heterogeneity of transposon-derived and gene-coding transcripts. By mapping to repeat elements (based on RepeatMasker annotations) as well as the genome, we obtained profiles for expression of TEs as well as genes for each cell (Fig. S4). We focused on GW19 in light of the high level of piRNA expression and dynamic HIWI2 localization. t-distributed stochastic neighbor embedding (t-SNE) analysis based on gene expression profiles identified three distinct clusters corresponding to PGCs, AGCs and somatic support cells based on known lineage- and stage-specific markers (Fig. 3A,B, Fig. S5A-C). As anticipated, PGC and AGC clusters were distinguished by high expression of OCT4 (POU5F1) and VASA (DDX4), respectively (Fig. 3B). The machinery of primary piRNA biogenesis [MAEL and HILI (PIWIL2)] was expressed in both PGCs and AGCs, though higher in the latter, whereas HIWI2 mRNA was largely limited to AGCs (Fig. 3B,C). This developmental delay in transcription of HIWI2 as compared with HILI recapitulates the sequential expression of the proteins observed in the fetal testis tissues (Fig. 1A,B). With a few exceptions, TE expression was dramatically elevated in AGCs compared with PGCs, and absent in the soma (Fig. 3B,C). The most highly expressed TEs in AGCs belong to the L1 clade, specifically L1HS, L1PA2 and L1PA3. This corresponds well with the abundance of piRNAs derived from the L1 clade (Fig. 2E). Other highly expressed TEs include members of the Alu family, especially the evolutionarily young AluYa5 and AluYb8 subfamilies, which are known to be active in humans (Batzer and Deininger, 2002). Paired correlation analyses between L1HS, and either HILI or HIWI2 both yielded positive correlations (Fig. 3D,E). Together, this analysis shows that at the single cell level, the expression of genes of the transposon repression network is upregulated during germ cell development in concert with the expression of transposons (Fig. 3C).

Open in a separate window

Fig. 3.

Single cell sequencing reveals dramatic upregulation of L1 family retrotransposons in advanced germ cells. (A) tSNE clustering on transcriptome reference was used to generate three distinct cell populations in the GW19 dataset. (B) Violin plots showing expression of germ cell markers (top), transposon repression genes (middle) and L1 family retrotransposons (bottom) in AGCs, PGCs and soma. (C) Heatmap displaying the most differentially expressed transposons sorted by myAUC score (top), and expression of germ cell markers (middle) and transposon repression genes (bottom). (D-F) Pairwise correlation analysis between L1HS and HILI (D) and HIWI2 (E) in GCs; and HSP90α (F) in AGCs (i), PGCs (ii) and combined GCs (ii). Pearson's correlation coefficient scores and P-values are listed above graph.

Dynamic relationships between L1 and the repression pathway at transcript and protein level

To confirm the correlation between L1 transcripts and the transposon repression network, we performed immunostaining. Using a previously validated antibody against the L1-ORF1p RNA-binding protein (Rodić et al., 2014), we detected expression primarily in VASA⁺ AGCs starting at GW16, whereas earlier it was found very sparsely and at low levels (Fig. S6). L1-ORF1p was restricted to the cytoplasm and concentrated in granules, some of which colocalized with HILI foci (Fig. 4A). However, L1-ORF1p granules were distinct from those containing HIWI2 (Fig. 4B). As suggested by the heatmap in Fig. 3C, a small subset of PGCs express high levels of L1 transcript but not the transposon repression network. We similarly observed rare instances of robust L1-ORF1p expression in OCT4⁺ PGCs (Fig. 4C, top panel) and in cells lacking VASA and harboring low levels of OCT4 (Fig. 4C, bottom panel). It is not clear whether these cells are transitioning between PGCs and AGCs or whether they will successfully turn on the repression machinery, but this observation raises the issue of whether transposon protein expression precedes that of HIWI2 and piRNA biogenesis.

Open in a separate window

Fig. 4.

L1 expression in advanced germ cells. (A) L1 and HILI immunofluorescence in a GW16 human testis. Arrowheads indicate puncta of L1/HILI colocalization. (B) Expression of L1 and HIWI2 at GW17. Scale bars: 20 μm. Arrows indicate regions of mutually exclusive staining between L1 and HIWI2 puncta. (C) L1 expression at GW21 (top panel) and GW18 (bottom panel), co-stained with OCT4 or a pan-GC marker (OCT4 and VASA), respectively. (Top) Arrowheads indicate an L1⁺, OCT4⁺ germ cell; arrows highlight a L1⁺, OCT4⁻ cell. (Bottom) Arrows indicate L1⁺ and pan-GC negative cell. (D) L1⁺/HIWI2^low germ cells (arrows), a HIWI2^high/L1^low germ cell (arrowheads) and double-positive AGCs (stars) in human fetal testis at GW19 and GW22. (E) xy scatter plots of HIWI2 versus L1 fluorescence intensities (left panels) and HIWI2 versus VASA (right panels), with Pearson's correlation coefficients and corresponding P values. For GW19 (top), n=526 VASA⁺ objects, counted from three stitched xy sections. For GW22 (bottom), n=472 VASA⁺ objects, from two stitched xy sections. (F) L1^high/HSP90α^low cells (arrows), HSP90α^high/L1^low cells (arrowheads) and double-positive cells (stars) in human fetal testis at GW16-22. (G) xy scatter plots of HSP90α versus L1 fluorescence intensities (left) and HSP90α versus VASA fluorescence intensities (right). Pearson's correlation coefficients and corresponding P values are indicated. For GW16 (top), n=768 VASA⁺ objects, counted from three stitched xy sections. For GW19 (middle), n=1809 VASA⁺ objects, counted from four stitched xy sections. For GW22 (bottom), n=2098 VASA⁺ objects, counted from three stitched xy sections. Scale bars: 20 μm in A,B,D,F; 10 μm in C. Dashed red boxes in E and G indicate the ‘arms’ of single-positive cell populations.

Examination of the relative expression of L1-ORF1p, HIWI2 and VASA in AGCs revealed a more dynamic relationship during fetal development than previously identified in the single cell RNA-seq analysis. The majority of L1-ORF1p expressing AGCs at GW22 co-expressed HIWI2 with significant correlation that approached the r value between HIWI2 and VASA (Fig. 4D,E). Earlier, however, at GW19, we observed VASA⁺ cells with exclusive expression of either L1-ORF1p or HIWI2 and within the AGC population, the correlation between the two proteins was lacking at this timepoint (Fig. 4D,E). Compared with the high correlation between L1 and HIWI2 transcripts in single cell RNA-seq at GW19 (Fig. 3E, Fig. S5D,E), these results indicate a delay in the coordinated expression of HIWI2 and L1-ORF1p proteins during fetal development. It is likely that the VASA single-positive cells at GW22 (Fig. 4E, right) represent newly differentiating PGCs that have not yet turned on HIWI2.

HSP90α is a chaperone protein, with an inducible isoform enriched in germ cells, which has been implicated in transposon repression and piRNA biogenesis in flies and mice (Gangaraju et al., 2011; Ichiyanagi et al., 2014; Izumi et al., 2013; Specchia et al., 2010; Xiol et al., 2012). Levels of HSP90α transcript were high but variable in AGCs and PGCs (Fig. 3B,C), with a negative correlation observed between HSP90α and L1HS in AGCs but not PGCs (Fig. 3F). Consistent with this, immunofluorescence revealed that a significant number of AGCs at GW16 and GW19 expressed exclusively L1-ORF1p or HSP90α proteins but not both (Fig. 4F, top and middle panels); accordingly, the XY scatter plots for L1-ORF1p and HSP90α intensities displayed ‘arms’ at the extreme end of each axis of AGCs with mutually exclusive expression (dashed red boxes in Fig. 4G). This relationship changed by GW22, when most AGCs expressing L1-ORF1p also expressed HSP90α (Fig. 4F, bottom panel), and a highly significant correlation was found (Fig. 4G). Together, these observations suggest that three AGC subpopulations co-exist at GW16 – L1-high, HSP90α-high and those low for both – and by GW22 a subpopulation of HSP90α single positive AGCs persist, but the majority of L1⁺ cells also express moderate to high levels of HSP90α.

Animals and xenograft model

Wild-type embryos were generated by mating CD1 females with homozygous Oct4-ΔPE-GFP males (MGI:3057158, multiple-copy transgene insertion; Anderson et al., 1999). Embryos were dissected from timed matings and staged using anatomical landmarks from E16.5 to E18.5.

Immunodeficient CD1 nud/nud homozygous (Charles River) mice were used as host for the xenograft. Prior to surgery, human fetal testes were cut into four fragments and held in ice-cold RPMI 1640 medium until xenografting, with one fragment set aside for fixation. Hosts were anesthetized, a 2 cm back midline incision was made to exteriorize each kidney and a 2 mm opening was made in the renal capsule. One fetal testis fragment was xenografted into the left kidney capsule of each host. The peritoneum was closed with c-9 silk surgical sutures (Henry Schein), and the incision was stapled closed with stainless steel wound clips (MikRon). Following surgery, animals were monitored daily for food and water intake, and any obvious signs of stress. All mouse work was carried out under the University of California, San Francisco, Institutional Animal Care and Use Committee guidelines in an approved facility of the Association for Assessment and Accreditation of Laboratory Animal Care International.

Immunofluorescence

We performed immunofluorescent staining on sectioned slides using our general staining protocol. Slides were washed three times with 1× PBS for 5 min at room temperature to dissolve OCT. Samples were permeabilized with PBT (1× PBS+0.5% Triton X -100) for 10 min at room temperature, followed by three short rinses with 1× PBS, and then blocked for 1 h at room temperature with fresh blocking buffer (10% heat inactivated donkey serum in 1× PBS+0.1% Triton X-100). Subsequent primary antibody incubations were performed in humidified chambers at 4°C overnight, with primary antibodies (see Table S3) diluted appropriately in blocking buffer. The next day, slides were washed three times in 1× PBS for 5 min at room temperature, and then treated with secondary antibodies (1:200) and DAPI nuclear counterstain (1:1000) for 1 h at room temperature. Slides were washed three times in 1× PBS, and then mounted with Vectashield Antifade Mounting Medium (Vector Laboratories, H-1000). Some antibodies required antigen retrieval with either sodium citrate buffer (pH 6.0) or Tris-EDTA buffer (pH 8.8) (see Table S3). Antigen retrieval conditions consisted of boiling slides in appropriate solution for 10 min, followed by washes in 1× PBS and processing according to the general staining protocol. 5meC antibody staining required Tris-EDTA (pH 8.8) antigen retrieval followed by PBT permeabilization and incubation in 3.5 N HCl for 1 h prior to the general staining protocol. Xenograft samples incubated with a primary antibody from a mouse host required antigen retrieval and blocking with Mouse on Mouse blocking reagent (Vector Laboratories, MKB-2213). Following our staining protocols, images were acquired from a Keyence fluorescence microscope and/or Leica Sp5 upright confocal microscope w/AOBS.

Immunofluorescence quantification

Small RNA-sequencing of human fetal testis

All samples were processed and sequenced by Quick Biology. Total RNA extractions were performed on bulk fetal testis tissue using Trizol, and quality control was performed on total RNA using the Agilent Total RNA Nano Bioanalyzer kit. Resulting RNA integrity number (RIN) scores, are provided in Table 2. Libraries were prepared using the QIAseq miRNA Library Kit, and library quality was assessed using the Agilent Bioanalyzer High Sensitivity DNA assay. Paired-end sequencing was performed on an Illumina HiSeq 4000.

Small RNA-sequencing data analysis

All samples were sequenced to a depth of approximately 35-40 million reads in order to meet the minimum read requirement for small RNA-seq data of 30 million as outlined by the ENCODE consortium (www.encodeproject.org/rna-seq/small-rnas/#restrictions). Raw reads in fastq files were trimmed to remove the adapter sequence (AACTGTAGGCACCATCAAT) using cutadapt version 1.8.1 (Martin, 2011), and reads were subsequently size filtered to 18-40 nt using cutadapt. FastQC (version 0.11.5) was utilized to broadly assess read quality. Trimmed reads were aligned to Repbase (version 22.11), a collection of consensus sequences of known human transposons, curated by GIRI: the Genetic Information Research Institute (Bao et al., 2015). Alignments were performed with bowtie version 1.2.2 (Langmead et al., 2009) using the following options: -v 1, -M 1, –best, –strata, signifying a maximum of one mismatch between the read and reference sequence, with multimappers being randomly assigned to only one location.

All reads that did not align to Repbase were subsequently aligned (with bowtie version 1.2.2 and same options) to the human genome (GRCh38) in order to determine the contribution of other RNA types to the total small RNA pool. Aligned reads were intersected with genomic features and counted using the htseq-count function of HTSeq version 0.9.1 (Anders et al., 2015), with the following parameters: -f sam -t exon -s yes -i gene_type -m intersection-nonempty –nonunique all -a 0. The reference file containing genomic feature coordinates was a custom GTF file generated by merging the GENCODE 28 comprehensive gene annotation list with GENCODE 28 predicted tRNA genes. Feature counts for miRNAs, tRNAs, snoRNAs and protein-coding transcripts were combined with the number of Repbase-mapping reads to generate a breakdown of RNA subtypes. For a comprehensive explanation of RNA subtypes and relative percentages, see Table S1.

Repbase-mapped reads were further analyzed for known piRNA signatures, including a 5′ uridine bias and 10th position adenine bias, by calculating the nucleotide frequency distribution per base position along the read for all pooled Repbase-mapped reads of a given sample. Read overlap signatures, which are characteristic of the piRNA ping-pong biogenesis cycle, were calculated on the MISSISSIPPI Galaxy server (mississippi.snv.jussieu.fr) using the ‘Small RNA Signatures’ tool (Galaxy Version 3.1.0), developed by Antoniewski (2014). Transposable element annotations and transposable element-specific read size distributions were calculated directly from alignment files post Repbase mapping.

Single cell RNA-sequencing reanalysis

We are grateful to Brandon Sigamony, Daniel Nguyen, Suping Peng and other members of the Laird lab for their feedback and support, as well as to Drs Gerald Cunha and Miguel Ramalho-Santos for comments on the manuscript.