Mettl3-mediated m⁶A regulates spermatogonial differentiation and meiosis initiation

Mettl3 is essential for male fertility and spermatogonial differentiation during spermatogenesis

To explore the function of Mettl3 in mouse spermatogenesis, we first used immunostaining to examine the expression of METTL3 in the mouse testis at 6 and 12 days post-partum (P6 and P12). METTL3 was expressed in both germ cells and somatic cells during testis development (Supplementary information, Figure S1A). Furthermore, METTL3 expression was relatively higher in undifferentiated spermatogonia that expressed PLZF (promyelocytic leukemia zinc-finger protein) at P6 (Supplementary information, Figure S1B). To investigate the function of Mettl3-mediated m⁶A modification in gametogenesis in mice, we specifically knocked out the Mettl3 gene in the germ cells. Using CRISPR-Cas9 system-assisted homologous recombination, two loxp sites were inserted into the intron 1 and intron 4 of the Mettl3 gene to generate mice carrying the floxed allele (Mettl3^flox/+) (Supplementary information, Figure S1C). To disrupt the Mettl3 gene in germ cells, the Mettl3^flox/+ mice were further mated with the Vasa-Cre mice that specifically expressed the Cre recombinase in germ cells driven by a Vasa promoter as early as embryonic day 15.5 (E15.5)²⁹ (Supplementary information, Figure S1C). Six genotypes of Mettl3 allele, including Mettl3^+/+, Mettl3^flox/+, Mettl3^flox/flox, Vasa-Cre, Mettl3^flox/+Vasa-Cre and Mettl3^flox/−Vasa-Cre, were obtained and verified by PCR and Sanger sequencing (Supplementary information, Figure S1D-S1H). The Mettl3^flox/−Vasa-Cre mice showed specific deletion of exons 2∼4 and loss of METTL3 expression in PLZF-positive spermatogonia, confirming the conditional knockout of Mettl3 (referred as Mettl3^cKO) (Figure 1A and Supplementary information, Figure S1E, S1H). The Mettl3^+/+, Mettl3^flox/+, Mettl3^flox/flox, Vasa-Cre and Mettl3^flox/+Vasa-Cre mice were healthy and phenotypically normal, and were used as control in the following experiments (referred as Mettl3^Ctrl).

Next, we analyzed the phenotypes of the Mettl3^cKO mice. The Mettl3^cKO mice were normal in growth, but were completely infertile and had much smaller testes with 80% reduction in testis weight at 8 weeks (Figure 1B and 1C). Rare differentiated cell types of spermatocytes were observed in the Mettl3^cKO seminiferous tubules at P8, P10 and P12, and no spermatids were observed at 8 weeks by hematoxylin and eosin (H&E) staining, indicating a severe defect in spermatogenesis in the Mettl3^cKO mice (Figure 1D). Co-staining of the germ cell marker DDX4 and the undifferentiated spermatogonia marker PLZF showed the presence of undifferentiated spermatogonia that were normal, but the number of germ cells was significantly reduced in the Mettl3^cKO testes at P8, P10 and P12 (Figure 1E, 1F and Supplementary information, Figure S1I-S1L). We next traced advanced stages of spermatogonial differentiation using the differentiated spermatogonial cell marker KIT at P8, the initial stage of spermatogonial differentiation, and P12. KIT-positive cells were significantly reduced in the Mettl3^cKO testes at both P8 and P12 (Figure 1G-1J). Collectively, these results demonstrated that Mettl3 is essential for male fertility in mice and its deletion leads to defects in the initiation of spermatogonial differentiation.

Mettl3 deletion causes a severe defect in meiosis during spermatogenesis

Meiosis, as the key feature of spermatogenesis, initiates around P10 and spermatocytes enter meiotic prophase around P12. To investigate the function of Mettl3 in meiosis, we examined the expression of meiotic marker genes at P10 and P12 by immunostaining. STRA8, a marker of meiosis initiation, showed a much-reduced expression level in Mettl3^cKO testes at P10 and P12. Both the ratio of STRA8-positive seminiferous tubules and the number of STRA8-expressing cells per tubule were significantly reduced in Mettl3^cKO testes compared to Mettl3^Ctrl testes at P10 and P12 (Figure 2A). The number of cells expressing the synaptonemal complex component SYCP3 in Mettl3^cKO testes was also much less than that in Mettl3^Ctrl testes (Figure 2B). We next traced the progression of meiotic prophase that can be divided into leptotene, zygotene, pachytene stages according to the distribution patterns of SYCP3 and γH2AX in spermatocyte nuclear spreads. Pachytene-stage spermatocytes were only detectable in Mettl3^Ctrl mice and not in Mettl3^cKO mice, while leptotene- and zygotene/zygotene-like-stage spermatocytes were observed in both Mettl3^Ctrl and Mettl3^cKO mice (Figure 2C-2E and Supplementary information, Figure S2A). Mettl3^Ctrl testes had around 29%, 31% and 40% of spermatocytes at the leptotene, zygotene and pachytene stages, respectively, whereas 84% of spermatocytes in Mettl3^cKO testes were in the zygotene/zygotene-like stage but none were in pachytene (Figure 2F).

To confirm the blockage of meiosis at the zygotene/zygotene-like stage, we examined testes for the presence of additional marker proteins. We found that the synapsis marker SYCP1 was also not detected in Mettl3^cKO spermatocytes (Supplementary information, Figure S2B). In contrast to the much greater number of foci of the DSB repair protein RAD51 at the zygotene/zygotene-like stage in Mettl3^Ctrl spermatocytes, very few RAD51 foci were detected in Mettl3^cKO spermatocytes, indicating defective DSB repair during meiosis upon Mettl3 knockout (Figure 2C-2E). γH2AX expression is strong at the leptotene stage (referred to as positive), gradually decreases around the zygotene stage (referred to as partial positive), and the protein is enriched in the sex-body at the pachytene stage (referred as only sex-body positive). While all three expression patterns could be seen in nuclear spreads of Mettl3^Ctrl spermatocytes, the only sex-body-positive pattern was not observed in Mettl3^cKO nuclear spreads (Figure 2G and Supplementary information, Figure S2C). Spermatocytes blocked in meiosis will undergo apoptosis; this was detectable by TUNEL assays showing a significantly increased apoptotic signal in Mettl3^cKO testes at P12 (Figure 2H and 2I). These results demonstrate that Mettl3-deficient spermatocytes are unable to reach the pachytene stage of meiotic prophase and that Mettl3 is required for normal meiosis.

Mettl3 deletion alters expression pattern of spermatogenesis-related genes

To systematically illuminate gene expression changes resulting from the Mettl3 conditional knockout in spermatogenesis, we compared the transcriptome between Mettl3^Ctrl and Mettl3^cKO testes at P6 and P12 through high-throughput RNA sequencing. The correlation analysis showed that the expression of the majority of genes in testes was not influenced by Mettl3 deletion (Supplementary information, Figure S3A, S3B and Table S4). We found that a total of 258 genes were differentially expressed in the Mettl3^cKO testes at P6. These included 157 down-regulated genes that were significantly enriched in genes required for the regulation of spermatogenesis including the meiotic cell cycle and synaptonemal complex assembly (Supplementary information, Figure S3C). We found that 941 genes were differentially expressed in Mettl3^cKO testes at P12. These included 699 down-regulated genes that were mainly enriched in genes required for spermatogenesis and regulation of meiosis (Figure 3A and Supplementary information, Table S4). We next examined expression patterns of 50 genes specifically expressed in germ cells according to a previous study³⁰. Unsupervised cluster analysis showed that the Mettl3^cKO expression pattern at P12 was much more similar to Mettl3^Ctrl at P6 than Mettl3^Ctrl at P12, consistent with the severely impaired spermatogonial differentiation observed in Mettl3^cKO testes (Figures 3B, 3C and Supplementary information, Table S4). Genes with reported functions in spermatogenesis^31,32,33 were mapped into an interaction network to illustrate regulationary relationships by Pathway Studio (Figure 3D). The genes significantly down-regulated in Mettl3^cKO testes included 9 genes involved in spermatogonial stem cell maintenance (Dazl, Ddx4, Plzf, Pax7, Nanos2, Id4, Pou3f1, Taf4b, Bcl6); 6 genes in spermatogonial differentiation (Sohlh1, Sohlh2, Neurog3, Kit, Dmrt1, Sox3); 8 genes in cell cycle (Ccna1, Ccna2, Ccnb1, Ccnb2, Ccnb3, Ccne2, Cdk1, Rad51); and 9 genes in meiosis (Stra8, Sycp1, Sycp3, Spo11, Rad18, Dmc1, Rec8, Mlh1, Smc1b) (Figure 3E, 3F and Supplementary information, Figure S3D). Taken together, this indicate that the Mettl3 deficiency globally influenced the expression pattern of spermatogenesis-related genes regulating multiple biological processes including spermatogonial stem cell maintenance, differentiation and meiosis.

Mettl3-mediated m⁶A regulates alternative splicing of genes functioning in spermatogenesis

Previous studies have shown that m⁶A could influence biological process through regulating alternative splicing events^11,19. We therefore employed rMATS software to analyze the differential alternative splicing events between Mettl3^Ctrl and Mettl3^cKO transcriptomes³⁴. Among the differential splicing events, we mainly detected exon skipping to be enriched, consistent with previous studies^11,19 (Figure 4A and Supplementary information, Table S5). Gene ontology analysis revealed that the differentially spliced genes between Mettl3^Ctrl and Mettl3^cKO samples at P6 were enriched in those for DNA replication and cell division processes whereas those at P12 were enriched for spermatogenesis (Supplementary information, Figure S4A).

We then detected the m⁶A sites and explored their distribution in P12 testes of wild-type mice through m⁶A individual-nucleotide-resolution cross-linking and immunoprecipitation (miCLIP) and sequencing (miCLIP-seq). A total of 12 516 putative m⁶A residues within 5 092 genes were identified in testes (Supplementary information, Table S6). These sites were enriched in the m⁶A consensus motif (Figure 4B) and were predominantly distributed in the CDS and 3′ UTR regions, especially near stop codon regions (Figure 4C and 4D). Notably, the average m⁶A site number was higher in genes with skipped exons (Figure 4E). An alternative exon inclusion level measurement also revealed that genes containing m⁶A had a lower inclusion level in Mettl3^cKO samples compared to that in Mettl3^Ctrl samples whereas genes without m⁶A modification did not show this tendency (Figure 4F and Supplementary information, Figure S4B). These results are consistent with previous reports showing the role of m⁶A modification in promoting exon inclusion^11,19.

Next, we examined differential splicing events between Mettl3^Ctrl and Mettl3^cKO samples in spermatogenesis-related genes. Our findings showed that 9 genes had significantly decreased exon inclusion levels in Mettl3^cKO testes (Figure 4G, 4H and Supplementary information, Figure S4C-S4I). The m⁶A-RIP-qPCR analysis was used to validate the presence of m⁶A residues in alternative exons or neighboring regions, and the results supported occurrence of alternative splicing events regulated by m⁶A (Supplementary information, Figure S4J-S4Q). For example, Dazl, essential for spermatogenesis, had a shorter transcript with loss of 17 amino acids next to its RNA recognition motif in Mettl3^cKO testes (Supplementary information, Figure S4C). Similar splicing events were also found in Sohlh1 regulating spermatogonial stem cell differentiation^35,36,37, and Nasp- and Cdk11b-regulating cell division^38,39,40,41 (Figure 4I and Supplementary information, Figure S4D, S4E). Taken together, these data revealed that Mettl3-mediated m⁶A modification regulates the alternative splicing of genes that function in spermatogenesis.

Mouse

The mice used in this study were B6D2F1 (C57BL/6 × DBA2), DBA2 and C57BL6 strains. Specific pathogen-free-grade mice were purchased from Beijing Vital River laboratory animal center and housed in the animal facilities of the Institute of Zoology, Chinese Academy of Sciences. All animal experiments were carried out in accordance with the guidelines for the Use of Animals in Research issued by the Institute of Zoology, Chinese Academy of Sciences.

Mouse breeding

Mettl3^flox/+ mice were generated by the CRISPR-Cas9 system-assisted homologous recombination. Vasa-Cre transgenic mice were a gift from professor Jia-Hao Sha, Nanjing Medical University. Mettl3^flox/+Vasa-Cre mice were obtained by mating Mettl3^flox/+ and Vasa-Cre mice and Mettl3^flox/flox mice were obtained by mating Mettl3^flox/+ and Mettl3^flox/+ mice. Mettl3^flox/+Vasa-Cre and Mettl3^flox/flox were mated to generate Mettl3^flox/−Vasa-Cre mice (Mettl3^cKO).

Construction of gene targeting vector

The T7-Cas9 vector was from our lab. The T7-sgRNA vector was digested with BsaI enzyme (NEB) and the linearized vector was gel purified. A pair of oligonucleotides for each targeting site (Mettl3 intron 1 and intron 4) was annealed and ligated to linearized vector. All oligonucleotides were listed in the Supplementary information, Table S1.

In vitro preparation of sgRNAs, Cas9 mRNA and oligonucleotide DNAs for microinjection

All RNAs prepared for microinjection were in vitro transcribed. SgRNAs and Cas9 mRNA were transcribed by the HiScribe T7 In Vitro Transcription Kit (NEB). The Cas9 mRNA was capped with the m7G (5′) ppp (5′) G RNA Cap Structure Analog (NEB) by T7 transcriptase. All RNAs were dissolved in DEPC water and stored in −80 °C.

Oligonucleotides used for CRISPR-Cas9 system-assisted homologous recombination through embryo injection were synthetized in Integrated DNA Technologies. All oligonucleotides were listed in the Supplementary information, Table S1.

Intracytoplasmic RNA microinjection

The intracytoplasmic RNA microinjection was performed according to a previous report. Briefly, one-cell-stage embryos were collected at 0.5 day post coitum (dpc). Each embryo was microinjected with 25 ng/μl L-sgRNA, 25 ng/μl R-sgRNA, 100 ng/μl Cas9 mRNA, 50 ng/μl L-loxp oligonucleotide and 50 ng/μl R-loxp oligonucleotide into the cytoplasm. After microinjection, surviving embryos were implanted on the same day into the oviduct of pseudopregnant B6D2F1 female mice. Full-term pups were obtained by natural labor at 19.5 dpc.

Immunohistochemical analysis

For immunostaining analysis: testes were fixed in 4% paraformaldehyde overnight at 4 °C. Following dehydration, the samples were embedded in paraffin and sectioned (5 μm). After dewaxing and hydration, the sections were boiled in an antigen retrieval buffer (0.01 M citric acid/sodium citrate, pH 6.0) for 10 min using a microwave oven, washed in PBS three times and permeated in 0.2% Triton X-100 for 15 min. The sections were blocked with 5% BSA in PBS for 2 h at room temperature before adding primary antibodies. The samples were washed three times in PBS after incubating at 4 °C overnight and then incubated with secondary antibodies for 2 h at room temperature. After washing in PBS three times, the sections were incubated in 2 μg/ml of DAPI (Molecular Probes, D1306) for 10 min at room temperature and washed in PBS three times. Sections were mounted on slices with Dako Fluorescence Mounting Medium (Dako Canada, ON, Canada). The immunohistochemical analysis of STRA8 were imaged with a Leica Aperio VERSA 8 microscope (Leica Biosystems, Germany). The immunofluorescence staining was imaged with a laser scanning confocal microscope LSM780 (Carl Zeiss, Germany). All the antibodies were listed in the Supplementary information, Table S2.

Genotyping of mice

First, all mice were genotyped with the tail DNA. Two pair of primers were used to detect the loxp insertion into the Mettl3 intron 1 (L-loxp-forward and L-loxp-reverse) and intron 4 (R-loxp-forward and R-loxp-reverse). The product sizes were 222 bp and 335 bp with the loxp sequence insertion into Mettl3 intron 1 and intron 4, respectively; whereas the product sizes from WT were 182 bp and 295 bp, respectively. Cre recombinase was detected by the Vasa-Cre primers and the product of PCR was 240 bp. The PCR reaction conditions were as follows: 95 °C, 2 min; 95 °C, 30 s, 61 °C, 30 s, 72 °C, 30 s, 35 cycles; 72 °C, 5 min. Then the germ cells were used to confirm the deletion of Mettl3 with the primers of L-loxp-forward and R-loxp-reverse and the product of Mettl3 deletion was 318 bp whereas the WT product was 2 554 bp. The PCR reaction conditions were as follows: 95 °C, 5 min; 95 °C, 30 s, 61 °C, 30 s, 72 °C, 2.5 min, 35 cycles; 72 °C, 10 min. All primers are listed in the Supplementary information, Table S3.

Histological analysis

For histological analysis: testes from Mettl3^Ctrl and Mettl3^cKO male mouse were fixed in Modified Davidson's Fluid (Formaldehyde:Glacial acetic acid:Ethanol:H₂O = 6:3:1:10) overnight at room temperature. The testes were dehydrated stepwise through an ethanol series (70%, 80%, 90%, 100% ethanol) and processed for paraffin embedding. 5 μm sections were cut with a Leica slicing machine (Leica RM2235, Leica Biosystems, Germany) and mounted on poly-D-lysine coated glass slices (Zhong Shan Golding Bridge biotechnology, Beijing, China). After dewaxing and hydration, the sections were stained with H&E using standard methods and imaged with a Leica Aperio VERSA 8 microscope (Leica Biosystems, Germany).

TUNEL analysis

The TUNEL assay was performed using the TUNEL BrightRed Apoptosis Detection Kit (Vazyme) according to the manufacturer's instructions. Briefly, sections were permeabilized by protein K and labeled with rTdT reaction mix for 1 h at 37 °C, reaction was stopped by 1× PBS. After washing in PBS, the sections were incubated in 2 μg/ml of DAPI (Molecular Probes, D1306) for 10 min at room temperature and washed in PBS. Sections were mounted on slices with Dako Fluorescence Mounting Medium (Dako Canada, ON, Canada). Images were obtained using a laser scanning confocal microscope LSM780 (Carl Zeiss).

RNA extraction and qRT-PCR

Total RNA was extracted with TRIzol reagent (Invitrogen, 15596-018) from whole testes. After removing the genomic DNA with the RNase-Free DNase I Kit (Promega, M6101), 1 μg total RNA reverse-transcribed into cDNAs using the reverse Transcription System (Promega, A3500). qRT-PCR was performed using a SYBR Premix Ex Taq kit (TaKaRa, RR420A) on Agilent Stratagene Mx3005P. Relative gene expression was analyzed based on the 2^−ΔΔCt method with Actin as internal control. At least three independent experiments were analyzed. All primers were listed in the Supplementary information, Table S3.

m⁶A-miCLIP-seq

Single-base resolution high-throughput sequencing of the mouse testes methylome was carried out according to previously reported methods⁶⁰ with some modifications. Briefly, mRNA was purified using a Dynabeads^® mRNA Purification Kit (Life Technologies, 61006) and fragmented into around 100 nt fragments using the RNA fragmentation reagent (Life Technologies, AM8740). 5 μg mRNA was incubated with 10 μg of anti-m⁶A antibody (Abcam, ab151230) in 450 μl immunoprecipitation buffer (50 mM Tris, pH 7.4, 100 mM NaCl, 0.05% NP-40) and incubated at 4 °C for 2 h with gentle rotation. The mixture was transferred to a clear flat-bottom 96-well plate (Corning) on ice and irradiated three times with 0.15 J/cm² at 254 nm in a CL-1000 Ultraviolet Crosslinker (UVP). After irradiation, the mixture was collected together and then immunoprecipitated by incubation with dynabeads protein A (Life Technologies, 10008D) at 4 °C by rotating for another 2 h. After extensive washing, end-repair and linker ligation was performed on-bead. Finally, the bound RNA fragments were eluted from the beads with 2 μg/μl proteinase K in PK buffer (100 mM Tris-HCl, pH 7.4, 50 mM NaCl, 10 mM EDTA) by digestion at 55 °C for 20 min. RNA was isolated from the eluate by phenol-chloroform (Life Technologies, AM9730) extraction and ethanol precipitation. Purified RNA was reverse transcribed with Superscript III reverse transcriptase (Life Technologies, 18080093) according to the manufacturer's protocol as described previously^18,61. cDNA was size-selected on a 6% TBE-Urea gel (Life Technologies, EC68655BOX) and regions corresponding to 70-100 nt were used for further analysis. Thus, purified cDNA was circularized with Circligase II (Epicentre, CL9021K) and re-linearized with BamHI (NEB, R0136M). Libraries were PCR amplified with Accuprime Supermix 1 enzyme for 18 cycles and size-selected on an 8% TBE gel (Life Technologies, EC62155BOX). Sequencing was carried out on an Illumina HiSeq 3000 platform according to the manufacturer's instructions.

m⁶A-RIP-qPCR

The procedure of m⁶A immunoprecipitation (m⁶A-RIP) was modified from previously reported methods²³. In brief, purified mRNAs from P6/P12 mouse testes were digested by DNase I and then fragmented into ∼300 nt fragments by 30 s incubation at 94 °C in RNA fragmentation reagent (Life technologies, AM8740). This reaction was then stopped with stop solution (Life technologies, AM8740), followed by standard ethanol precipitation and collection. 12 μg anti-m⁶A polyclonal antibody (Synaptic Systems) was pre-incubated with 50 μl Dynabeads Protein A (Life technologies, 10013D) in IPP buffer (150 mM NaCl, 0.1% NP-40, 10 mM Tris-HCl, pH 7.4) for 1 h at room temperature. After pre-heating at 75 °C for 5 min and chilling on ice immediately, 6 μg mRNAs were added to the above prepared antibody-beads mixture and incubated for 4 h at 4 °C by rotating. After extensive washing, bound RNAs were extracted by proteinase K digestion, phenol-chloroform (Life Technologies, AM9730) extraction and ethanol precipitation. The enrichment of m⁶A was quantified through qPCR as reported²³. m⁶A-RIP-qPCR primers are listed in the Supplementary information, Table S3.

RNA-seq library preparation and data processing

Total RNA was extracted from cultured cells by TRIzol reagent (Invitrogen, 15596-018). For RNA-seq library construction, the PolyA+ tailed RNA purification was performed for each sample using a Dynabeads mRNA purification kit (Ambion, 61006). The cDNA library was generated with a KAPA Stranded mRNA-Seq Kit (KAPA, KR0960) protocol. Sequencing was performed on an Illumina HiSeq 3000 platform with 101 bp paired-end-sequencing reactions. The RNA-seq reads of each sample were mapped to the mouse mm9 genome assembly independently by the HISAT2 software using the annotated gene structures as templates. Default parameters of HISAT2 were used except with the options “--dta-cufflinks” and “--rna-strandness RF” opening⁶². Reads with unique genome location were retained for gene expression calculation using Cufflinks (version 2.0.2) with the option “--GTF”⁶³. Genes with no less than 0.5 FPKM in at least one sample were used for the differentially expressed genes analysis.

Alternative splicing analysis

Reads with a unique genome location in SAM format were used for the splicing events identification. The AS events between the two sample groups were identified using rMATS version 3.2.5³⁴. Five types of AS events from RNA-seq data with two replicates were identified. In each rMATS run, the first group was compared to the second group to identify differentially spliced events with FDR less than 0.05. The alternative exon accumulation levels in each sample were calculated using methods previously described. RT-PCR validation of alternative splicing events were performed with cDNA reverse-transcribed from mouse testes RNA. Actin was used as an internal control. All primers are listed in the Supplementary information, Table S3.

Analysis of miCLIP-sequencing data

The method of pre-processing for reads was performed essentially as previously reported⁶⁴. The adaptor was trimmed by fastx_clipper fastx_toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Low quality bases were filtered by fastq_filter.pl, a custom perl script from CLIP Tool Kit (CTK)⁶⁵ and reads shorter than 24 nt discarded. After trimming, reads without a barcode marker sequence or within large fragment (length > 60 nt) are discarded to eliminate unspecific RNA tags. Based on the same criteria of a previously reported approach to process paired-end data, the forward reads were demultiplexed based on 5′ barcodes for individual replicates by fastq2collapse to remove PCR amplified reads and the reverse reads were reverse-complemented and processed like their forward mates. Finally, the random barcodes of remaining reads were stripped by stripBarcode.pl and moved to read headers for downstream processing by the CIMS pipeline.

Remaining reads were mapped to the mouse genome (version mm9) with BWA (v0.7.10)⁶⁶ and allowed ≤ 0.06 error rate (substitutions, insertions or deletions) per read as online CTK Documentation (https://zhanglab.c2b2.columbia.edu/index.php/CTK_Documentation). To identify the m⁶A locus, the mode of mutation calling was performed as previously reported with minor modification⁶⁰. For each mutation position, the coverage of a unique tag (k) and mutations (m) was determined by the CIMS.pl program⁶⁷. The positions within k ≥ 15 and m/k ≤ 50% were kept and only mutation positions within the RRACH motif were identified as m⁶A for downstream analysis to exclude the potential m⁶Am modification⁶⁰. The logo plot for bases adjacent to m⁶A was generated with WebLogo⁶⁸.

Other tools

The unsupervised hierarchical cluster analysis performed by the hcluster function of the amap package in R. Point plots were produced by the smoothScatter functions of R. And the Pearson correlation coefficients shown in the figures. The heatmap was produced by the heatmap.2 function of R. Gene ontology analysis of differentially expressed genes was analyzed in DAVID^69,70 and biological processes were selected based on P-values smaller than 0.05 and the figures were produced by ggplot2. The interaction relationship of spermatogenesis related genes was constructed by Pathway Studio⁷¹ based on the regulation reported in public reference and the network was showed by the Cytoscape⁷². The cumulative frequency curve and violin plot were also produced by ggplot2. The RNA-seq reads distribution were exhibited in IGV browser⁷³.

Quantification and statistical analysis

Statistical parameters including statistical analysis, statistical significance and n value are reported in the Figure legends. Statistical analyses were performed using Prism Software (GraphPad). For statistical comparison, Student's t-test was employed. A value of P < 0.05 was considered significant.

Data and software availability

The accession number for the sequencing data reported in this paper is NCBI GEO: GSE99773. This parent directory includes the following datasets: GEO: GSE99771 (RNA-seq) and GEO: GSE99772 (m⁶A-miCLIP-seq). The public RNA-seq datasets (GSE53249, GSE56010, GSE61997 and GSE37005) for alternative splicing events identification by rMATS were downloaded from GEO database.