A mammalian-like piRNA pathway in Axolotl reveals the origins of piRNA-directed DNA methylation

Abstract

The piRNA pathway protects animal germlines from active transposons. Mammals employ a cytoplasmic pathway to destroy transposon transcripts during germline reprogramming. This post-transcriptional mechanism is ancient and found throughout the animal kingdom. A nuclear piRNA pathway mediates transposon DNA re-methylation, which is believed to be bespoke to mammals. However, when exactly piRNA-directed DNA methylation evolved remains unknown. We found that a mammalian-like piRNA pathway evolved early in tetrapod evolution and is found and expressed in its current configuration in the axolotl salamander. Analysis of axolotl testes and oocytes revealed diverse repertoires of piRNAs and pervasive post-transcriptional targeting of young transposons. We identified high levels of genome methylation in axolotl spermatozoa, with full-length transposons being heavily methylated. Our findings reveal that the mammalian nuclear piRNA pathway has ancient vertebrate origins, and it has likely been safeguarding the germline throughout most of tetrapod evolution. Thus, the emergence of piRNA-directed DNA methylation is a pivotal epigenetic evolutionary event that may have laid the foundation for germline reprogramming and genomic imprinting.

Keywords: DNA Methylation; Genomic Imprinting; Germline; Germline Reprogramming; Transposon; piRNA.

Figure EV1. Identification and distribution of potentially active transposons in the axolotl genome.

(A) Distribution of substitution rate and sequence length for every single copy in each transposon family. Defined potentially active copies are highlighted in red. (B) Composition of potentially active copies. (C) Scatter plot of percentage of potentially active copy number and total copy number for transposon classes. Transposon classes without potentially active copies are omitted.

Figure EV2. piRNA repertoires and transposon targeting in axolotl germlines.

(A) The overlap of unique piRNA sequences among the four samples. (B) Mean piRNA signal level over testis or oocyte for each transposon family. ****, active copy number >1000; ***, active copy number >100; **, active copy number >10; *, active copy number >0. (C) Tracks for piRNAs targeting both sense and antisense strands of representative transposon consensus sequence. For all panels, testis, n = 2; oocyte, n = 2.

Figure EV3. Features of ping-pong cycle in axolotl germlines.

(A) Relative frequency of the nucleotide distance between 5′ ends of complementary piRNA pairs over representative transposon families. Mean and S.E.M. are presented. (B) Nucleotide composition of the first and tenth position in piRNAs over representative transposon families. A, Adenine; U, Uracil; C, Cytosine; G, Guanine. For all panels, testis, n = 2; oocyte, n = 2.

Figure EV4. DNA methylation landscape of minor transposon classes and CpG islands in axolotl spermatozoa.

(A) Percentages of CpG methylation levels over all copies of minor transposon classes in the axolotl. For boxplots, the middle line represents the median; boxes represent the 25th (bottom) and 75th (top) percentiles; whiskers represent median ± 1.5× interquartile range; and outside values are not shown. Rhombus, mean level. (B) Length distribution of all CGI regions (n = 248,988). (C) CGI enrichment over TSS and adjacent regions. For statistical tests, ***P value < 0.001, n.s. not significant; Chi-square test. (D) Distribution of CpG methylation level and substitution rate for all transposon copies or copies overlapped with CGI of minor transposon classes. (E) Percentages of transposons overlapped with shuffled region. For all panels, spermatozoa, n = 3; pooled for analysis.

Figure 1. Modern mammalian piRNA pathway might originate from common tetrapod ancestors with the axolotl.

(A) Schematic of piRNA pathway in mammals. Major piRNA biogenesis factors and nucleus factors are listed in box with dashed lines. (B) The piRNA pathway and germline specification mechanism for vertebrates. Left, phylogenetic tree of representative vertebrate species. Middle, genome size, transposon percentage, and the types of germline specification mechanism. Germline specification mechanisms for each species are from Extavour and Akam, 2003 and Hansen and Pelegri, 2021 (see “Methods”). Right, the existence of piRNA factors in vertebrate genomes. For lungfishes, South American lungfish (SA), African lungfish (AF), Australian lungfish (AU). Source data are available online for this figure.

Figure 2. The axolotl has a mammalian-like piRNA pathway.

(A) Heatmap for expression of piRNA pathway factors in axolotl somatic and gonadal tissues. Source data are available online for this figure.

Figure 3. The landscape of potentially active axolotl transposons.

(A, B) Genomic composition of transposon classes in axolotl. Both the whole genome (A) and each chromosome (B) are displayed. (C) Distribution of substitution rate and sequence length for every single copy in the representative family. Defined potentially active copies are highlighted in red. (D) Lollipop plots display the copy number for each transposon family. Total copy number (left) and potentially active copy number (right) are shown. The red dashed line indicates transposon families with copy numbers greater than 1000. Bar plot indicates potentially active copy number percentage within transposon family. Source data are available online for this figure.

Figure 4. The axolotl genome encodes a diverse repertoire of transposon-targeting piRNAs.

(A) Nucleotide (nt) length distribution from small RNA libraries. Mean and S.E.M. are presented. (B) Classification of piRNAs mapped to the axolotl genome. (C) Mean piRNA signal level over four testis and oocyte samples for each transposon family. ****, active copy number >1000; ***, active copy number >100; **, active copy number >10; *, active copy number >0. (D) Tracks for piRNAs targeting both sense and antisense strands of representative transposon consensus sequences. (E) Number of annotated piRNA clusters in the genome. (F) piRNA clusters ranked by the cumulative fraction of piRNA RPMs (orange). The RPKM of piRNAs from each cluster is displayed (gray). The red dashed line indicates the number of top-ranked clusters that collectively account for 90% of cluster-derived piRNAs (90th percentile). (G) Length distribution of annotated piRNA clusters in the genome. (H) The number of piRNA clusters that overlap between the two genders. For all panels, testis, n = 2; oocyte, n = 2. Source data are available online for this figure.

Figure 5. The axolotl piRNA clusters exhibit sex-specific expression patterns.

(A) Genome-wise distribution of piRNA signal and piRNA clusters (inner circle), with zoom-in view for chr3p and chr13p (outer circle). (B) Tracks on both sense and antisense strands for piRNA clusters are highlighted in (A). Genic regions and transposon regions are shown. For all panels, testis, n = 2; oocyte, n = 2. Source data are available online for this figure.

Figure 6. piRNA pathway post-transcriptionally silences transposons in both male and female germlines.

(A) Schematic of the ping-pong cycle and transposon silencing. Two features of piRNAs generated from the ping-pong cycle for transposon transcript cleavage are displayed. One is the 10 nt overlap between 5′ of complementary piRNA pairs, the other is the 1U/10A preference. (B, D) Relative frequency of the nucleotide distance between 5′ of complementary piRNA pairs. 5′ nucleotide overlap frequency over all families from major transposon classes (B) and representative transposon families (D) are shown. Mean and S.E.M. are presented in (D). (C, E) Nucleotide composition of the first and tenth position in piRNAs. Nucleotide proportion over all families from major transposon classes (C) and representative transposon families (E) are shown. A adenine, U uracil, C cytosine, G guanine. For all panels, testis, n = 2; oocyte, n = 2. Source data are available online for this figure.

Figure 7. The genome of axolotl spermatozoa is highly methylated.

(A) Percentages of CpG methylation levels over the whole genome or certain genomic elements. (B) Metaplots of mean CpG methylation over gene bodies and adjacent 2 kb. (C) Percentages of CpG methylation levels over all genomic transposons or certain transposon classes. Methylation levels over all copies or only active copies are shown. (D) Metaplots of mean CpG methylation for selected transposon families over all transposon/active transposon bodies and adjacent 2 kb. (E) Distribution of CpG methylation levels across the whole genome or CGI regions. (F) Genomic annotation for CGIs or shuffled regions. (G) CpG methylation levels over all transposon/TSS regions or regions overlapped with CGIs. (H) Distribution of CpG methylation level and substitution rate for all transposon copies or copies overlapped with CGIs. (I) Percentages of transposons overlapped with CGIs. For boxplots, the middle line represents the median; boxes represent the 25th (bottom) and 75th (top) percentiles; whiskers represent median ± 1.5× interquartile range; and outside values are not shown. Rhombus, mean level. For statistical tests, ***P value <0.001, n.s. not significant; (E, F) Chi-square test, (C, G) unpaired two-sided t test. The exact P values for (C, E–G) are provided in Source Data. For all panels, spermatozoa, n = 3; pooled for analysis. Source data are available online for this figure.