Comprehensive Analysis of the Transcriptome-Wide m6A Methylation in Mouse Pachytene Spermatocytes and Round Spermatids

Spermatogenesis, an efficient and complex system in male germline development, requires a series of elaborately regulated genetic events in which diploid spermatogonia differentiate into haploid spermatozoa. N6-methyladenosine (m6A) is an important epigenetic RNA modification that occurs during spermatogenesis. ALKBH5 is an m6A eraser and knocking out Alkbh5 increases the level of total m6A methylation and causes male infertility. In this study, comprehensive analyses of MeRIP-seq and RNA-seq data revealed differences between wild-type (WT) and Alkbh5 knockout (KO) mice. In pachytene spermatocytes (PA), 8,151 m6A peaks associated with 9,959 genes were tested from WT and 10,856 m6A peaks associated with 10,016 genes were tested from KO mice. In the round spermatids (RO), 10,271 m6A peaks associated with 10,109 genes were tested from WT mice and 9,559 m6A peaks associated with 10,138 genes were tested from KO mice. The peaks were mainly concentrated in the coding region and the stop codon of the GGAC motif. In addition, enrichment analysis showed significant m6A methylation genes in related pathways in spermatogenesis. Furthermore, we conducted joint analyses of the m6A methylome and RNA transcription, suggesting an m6A regulatory mechanism of gene expression. Finally, seven differentially expressed mRNAs from RNA-seq data in both PA and RO were verified using qPCR. Overall, our study provides new information on m6A modification changes between WT and KO in PA and RO, and may provide new insights into the molecular mechanisms of m6A modification in germ cell development and spermatogenesis.

Introduction

In mammals, spermatogenesis consists of a series of tightly controlled processes, including mitosis, meiosis, and spermiogenesis (Li et al., 2014; Griswold, 2016). The appropriate function of spermatogonial stem cells (SSCs) determines the sustained production of sperm (Oatley and Brinster, 2012). SSCs proliferate through mitosis, and the resulting daughter cells are divided into two types: those that maintain their characteristics for mitosis, and those that divide to form A paired (Apr) and A aligned (Aal) spermatogonial cells (Oatley and Brinster, 2008). The Aal spermatogonia undergo a series of mitotic divisions, followed by differentiation into type B spermatogonia. The type B spermatogonia generate primary spermatocytes which enter meiotic division which includes the pachytene spermatocyte (PA) stage, and homologous chromosome recombination which leads to the formation of haploid round spermatids (RO). Finally, RO transform into mature spermatozoa (Ahmed and de Rooij, 2009; Handel and Schimenti, 2010; O'Donnell et al., 2011; Naro et al., 2017).

N6-methyladenosine (m6A) modification is implicated in the regulation of expression at the post-transcriptional level, including pre-mRNA splicing, mRNA nuclear export, stability, decay, and translation (Wei et al., 1975; Dominissini et al., 2012; Meyer et al., 2012; Wang et al., 2014; Alarcon et al., 2015; Wang et al., 2015). In mammalian cells, m6A is a dynamic and reversible modification process involving writers, erasers, and readers (Jia et al., 2013; Fu et al., 2014; Yue et al., 2015; A. Alemu et al., 2016; Zhao et al., 2017). METTL3, an m6A writer, modulates m6A modification that acts primarily as a major catalytic component, plays a functional role in mice spermatogenesis, as METTL3-deficient mice are infertile (Lin et al., 2017). ALKBH5, an m6A eraser, was recognized as the second m6A demethylase (Ensfelder et al., 2018). Upon deletion of Alkbh5, the fertility of male mice was impaired due to apoptosis of pachytene and metaphase-stage spermatocytes and aberrant spermiogenesis. Significantly reduced numbers of PA and RO were observed in Alkbh5 knockout (KO) mice (Zheng et al., 2013). Another study showed that Alkbh5 deficiency resulted in aberrant splicing and increased levels of m6A modification in spermatocytes and RO, generated transcripts with a longer 3ʹ-untranslated region (3ʹ-UTR), and shorter transcripts, suggesting that transcripts with longer 3ʹ-UTR have much higher levels of m6A than those with shorter 3ʹ-UTR near the stop codon (Tang et al., 2018). Furthermore, the study showed that spermatogenesis in Alkbh5 knockout (KO) mice stagnated at the round spermatozoa phase, and occasionally, some sperm were found in the caudal epididymis of KO mice. These KO sperm were malformed and displayed various structural abnormalities. However, currently, the function of m6A methylation in PA and RO is unclear.

In the present study, we aimed to investigate the role of m6A modification, to facilitate future functional studies of mammalian m6A, and to explore the pathogenic mechanism of spermatogenic disorder. To do this, we examined the m6A peaks in the RNA of PA and RO purified from Alkbh5 KO and wild-type (WT) adult mouse testes using the STA-PUT method. RNA-Seq and MeRIP-Seq were used to uncover the impact of the Alkbh5 gene on transcriptome-wide m6A methylation levels in mouse spermatogenesis and to find out whether the functional role of m6A methylome varies between different spermatogenic stages. We also identified the molecular functions within the cell type-specific differential methylation peaks in spermatogenesis, which provides a useful reference for further study into underlying mechanisms in male germline development.

Materials and Methods

Mice

The animal study was reviewed and approved by the Institutional Review Board of Nantong University (Nantong, China). The Alkbh5^+/− mice were a gift from Dr. Yamei Niu, Institute of Basic Medical Sciences, Chinese Academy of Medical Science, and were purchased from The Jackson Laboratory (Bar Harbor, ME). Homozygous Alkbh5 mice were produced by intercrossing of heterozygous Alkbh5 mice. Mice were housed in a temperature-and humidity-controlled cages with 12-h light/dark cycles. All WT and Alkbh5 KO mice used in this study were C57BL/6J (B6) mice. Male Alkbh5 KO mice have been previously reported (Zheng et al., 2013).

Gross Morphology and Sperm Analysis

Testes and whole epididymis were dissected from 12-week-old WT and Alkbh5 KO mice, and washed and immersed in phosphate buffered saline (PBS). After weighing the testis and recording its morphology, the cauda epididymitis was punctured several times with a needle, and the sperm was released into a PBS solution at 37°C for 30 min. Sperm motility was determined using computer-assisted sperm analysis (CASA) and the measurements were repeated thrice. Sperm concentration was measured in triplicate using a hemacytometer.

Isolation of Spermatogenic Cells at Various Stages of Development

PA and RO were obtained from five adult WT and ten adult Alkbh5 KO mice testes using the STA-PUT method as described previously (Song et al., 2009). Cell purity and morphology were examined using a microscope, fractions containing similar cell types were mixed and centrifuged to collect purified PA and RO. The purity of the PA and RO cell fraction was higher than 90% as determined by the cell morphology analysis, as previously described (Bryant et al., 2013). The above experiments were repeated three times.

Histology

Testes and epididymis from WT and KO mice were fixed in 4% paraformaldehyde solution, dehydrated through gradient mixtures of ethyl alcohol and water, cleared by passing through xylene, embedded in paraffin, and sectioned to 5 μm-thickness. The sections were stained with haematoxylin and eosin (HE) solution. Finally, the micrographs were obtained using a light microscope (Leica, DM 2000, Germany).

Western Blot Analysis

Total protein lysates were extracted from whole testes of WT and KO mice using RIPA lysis buffer. The Pierce BCA Protein Assay Kit (P0012, Beyotime) was used to determine the protein concentration. Protein lysates (30 μg) were resuspended in sodium dodecyl sulphate (SDS) sample buffer, subjected to electrophoresis on a 10% SDS-polyacrylamide gel, and then transferred onto nitrocellulose membranes. The membranes were blocked with skimmed milk at a 5% concentration and incubated at 20°C for 1 h. Next, the membrane was incubated with anti-Alkbh5 (1:500, HPA007196; Sigma) and anti-β-tubulin (1:3000, ab108342; Abcam) antibodies overnight at 4°C, washed in PBS-T the next day, incubated with goat anti-rabbit HRP (1:5000, ab205718; Abcam) at 20°C for 1 h. Proteins were finally visualized using an Odyssey infrared imaging system (Li-Cor Biosciences).

Quantitative Real-Time PCR Assay

Total RNA was extracted from the testicular tissues of WT and KO mice using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). A total of 1 μg of total RNA was reverse transcribed using the HiScript Ⅱfirst Strand cDNA Synthesis Kit (Vazyme Biotech Co., Ltd). Real-time PCR was performed using TB Green™ Premix Ex Taq™ II (TakaRa Biotechnology, Kyoto, Japan) according to the manufacturer’s instructions. All primers for genes were synthesized by Genewiz and are shown in Supplementary Table S1. The relative abundance of mRNA was calculated using the 2^−ΔΔCt method and normalized to β-actin mRNA levels.

MeRIP Sequencing and RNA Sequencing

Tissue samples were isolated using grinding beads and a high-speed reciprocating shaking device (Tissuelyser-24) in a low-temperature grinding environment, and total RNA from three sets of PA and RO was isolated using TRIzol™ Reagent (cat. No 15596018; Invitrogen). The isolated RNA was treated with DNase I (M0303S; NEB) to remove genomic DNA contamination. Extracted RNA was precipitated in glycogen (25 μg/ml final) and isopropanol at −30°C for 2 h. The pellet was rinsed twice with 70% ethanol, dried, and dissolved in ultra-pure water, then the RNA was fragmented. Briefly, 2 μl 10× RNA fragmentation buffer (100 mM Tris-HCl, 100 mM ZnCl2 in nuclease-free water) was added, and 2 μl 0.5 M EDTA was added immediately after 6 min incubation at 70°C to terminate the reaction (Zeng et al., 2018). The fragmented RNA was purified and recovered using the Zymo RNA Clean and Concentrator-5 Kit (cat. no R1013; Zymo Research). A portion of the purified fragmented RNA was saved as an input control. The remainder of purified fragmented RNA was incubated with 2 μg anti-m6A antibody (ABE572; Sigma-Aldrich) in an IP reaction system containing Dynabeads Protein A (10002D; Invitrogen) and Dynabeads Protein G (10004D; Invitrogen) overnight at 4°C with gentle shaking. After magnetic separation of the beads, the supernatant was removed and ×5 precipitation buffer and RNA enzyme inhibitor were added. The reaction proceeded at 4°C for 1–3 h, followed by 2–3 washes with a low salt precipitation buffer and 2–3 washes with a high salt buffer solution. Thereafter, Dynabeads were placed in an elution buffer (200 μl) at 50°C for 30 min to separate out the m6A-antibody immunoprecipitated RNA fragments. The mRNA was purified via phenol-chloroform extraction and ethanol precipitation. The products were removed by ribosomal RNA, first-strand cDNA was synthesized according to the SMART protocol, and fragments were enriched through PCR amplification for cDNA library construction. The magnetic bead library fragments were purified by DNA to extract ultramicro-RNA methylated m6A to detect the library. A Bioptic QseP70 Analyzer was used to inspect the library and determine whether the size distribution of the library conforms to the theoretical size. Lastly, the PE150 sequencing mode was performed on the NovaSeq high-throughput sequencing platform.

Data Analysis

The bioinformatics data analysis were performed by Guangzhou Epibiotek Co., Ltd. Cutadapt (v2.5) was used to trim adapters and filter for sequences. The remaining sequence reads were aligned to the mouse Ensembl genome GRCm38. We used Hisat2 aligner (v2.1.0) under the parameters: “--rna-strandness RF” (Kim et al., 2015). m6A peak calling was performed using the exomePeak R package (v2.13.2) under parameters: “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = NA, FRAGMENT_LENGTH = 200” (Meng et al., 2014). Differential m6A peaks were identified using the exomePeak R package under the indicated parameters. GO and REACTOME analyses of the differentially expressed genes (DEGs) were performed using the clusterprofile R package (v3.6.0). The resulting data were visualized using the Guitar R package (v1.16.0) under the parameter “-len 6-RNA”. Homer (V4.10.4) was used to select the m6A peak with a p-value < 0.05 for de novo motif analysis.

Results

Alkbh5 Knockout Mice Display Impaired Spermatogenesis

Western blot analysis revealed that the ALKBH5 protein was indeed deleted in KO mice (Figure 1A). To investigate the cause of infertility in Alkbh5 KO mice, we first analysed testicular development from gross morphology and histological perspectives. The shape and mass of the testes in the KO mice were smaller and lower respectively than those in the WT group (Figure 1B). Histological analysis of testes sections revealed that the seminiferous tubules from KO mice contained multiple vacuoles of varying sizes and fewer spermatogenic cells compared to those from WT mice (Figure 1C). Furthermore, caput and cauda epididymis from Alkbh5-KO mice showed only a few degenerating immature spermatids or spermatocytes (Figure 1D). To further analyse this, we isolated epididymal spermatozoa from adult WT and KO mice, even though sperm counts in the epididymis of KO mice were low. Using HE staining, we also found that the spermatozoa from KO mice showed various structural abnormalities (Supplementary Figures S1). As expected, KO mice exhibited reduced sperm counts and impaired sperm motility and morphology (Supplementary Figures S2).

FIGURE 1

FIGURE 1. Western blot analysis of the ALKBH5 protein. β-Tubulin was used as the loading control (A). Representative right testes and right epididymis from 12-week-old WT and Alkbh5 KO mice (left image). Average testis weight from 12-week-old WT and Alkbh5 KO mice (right image) Bars indicate means ± SEM (n = 10, ****p < 0.0001) (B). Seminiferous tubule architecture of WT and Alkbh5 KO mice were stained with hematoxylin-eosin (H&E) (Scale bars, 50 μm) (C). Cauda and caput epididymis from WT and Alkbh5 KO mice stained with H&E (Scale bars, 50 μm) (D).

Overview of the m6A Methylation Map in Pachytene Spermatocytes (PA) and Round Spermatids (RO)

We performed MeRIP-seq analysis of spermatogenic cells in PA and RO at various stages from Alkbh5-KO and WT mice (Meng et al., 2013). Alkbh5-KO PA had 1,929 significantly upregulated methylation peaks, and 326 downregulated peaks (|log2FC| > 1 and FDR <0.05). The Alkbh5-KO RO had 1,589 significantly upregulated methylation peaks, while 1,602 were downregulated compared with that in the WT RO (Figures 2A,B), suggesting that knocking out Alkbh5 may have a considerable effect on the upregulation of methylation in PA, while the number of downregulated methylation peaks was much smaller than that of upregulated peaks in PA. However, RO results differ from those of PA where the number of upregulated methylation peaks was found to be similar to that of downregulated methylation peaks. The top 10 altered methylation peaks in PA and RO are presented in Supplementary Table S2. Next, the distribution of methylation peaks in PA from both Alkbh5-KO and WT were similar in the 3ʹ-UTR region of mRNA and most parts of the long non-coding (lnc) RNA (Figure 2C). Specifically, methylation peaks in the Alkbh5-KO PA indicate a significantly different pattern from peaks in the WT PA, with a relative increase in the number of methylation peaks in the stop codon (21.3 vs. 19.7%), a start codon (4.7% vs. 4.3%), 3ʹ-UTR (20% vs. 17.8%), 5ʹ-UTR (1.6% vs. 1.4%), and TSS (0.6% vs. 0.3%), and a notable decrease in the coding DNA sequences (CDS) (56.5% vs. 51.8%) (Figure 2E). Compared to that in PA, the distribution of different regions in RO differed considerably. The knocking out of Alkbh5 contributed to a larger change in the distribution of lncRNA in RO than in PA (Figure 2D). The distribution of different regions in RO is conversely different from that in PA. We found that methylation peaks in the Alkbh5-KO RO indicate a significantly different pattern from peaks in the WT RO with a relative increase in the number of methylation peaks in the CDS (51.6% vs. 49.7%), 5ʹ-UTR (1.4% vs. 1.2%), and a relative decrease in the 3ʹ-UTR (21.4% vs 22%), TSS (0.4% vs 0.7%), stop codon (21.1% vs 21.4%), and start codon (4.3% vs. 5.1%) (Figure 2F). In addition, the GGAC motif was highly enriched within m6A sites in PA and RO (Figure 2G).

FIGURE 2

FIGURE 2. Volcano plot showing differential methylation peaks between WT PA and Alkbh5-KO PA (A). Volcano plot showing differential methylation peaks between WT RO and Alkbh5-KO RO (B). m6A peaks were most enriched in the coding sequence near the stop codon on mRNA and the distribution on lncRNA in PA (C). m6A peaks were most enriched in the coding sequence near the stop codon on mRNA and the distribution on lncRNA in RO (D). Distribution of different regions in Alkbh5-KO and WT PA (E). Distribution of different regions in Alkbh5-KO and WT RO (F). Top five consensus sequence motifs enriched in m6A peaks identified from Alkbh5-KO PA (G).

Differential Methylation Peaks are Enriched in Pathways Related to Spermatogenesis

To investigate the biological significance of m6A differences between Alkbh5-KO and WT spermatogenic cells, GO and REACTOME pathway analyses for differentially methylation peaks were performed (Ashburner et al., 2000; Draghici et al., 2007). GO analysis can be divided into two parts: biological processes (BP) and molecular functions (MF). First, we examined the differences between Alkbh5-KO PA and WT PA. We found that both upregulated and downregulated methylation peaks were closely related to covalent chromatin modification and histone modification. Moreover, we found that upregulated methylation peaks are linked with activation of MAPK activity and that downregulated methylation peaks are linked with negative regulation of the MAPK cascade (Figures 3A,B). Next, we analysed the differences between RO from Alkbh5-KO and WT. We found that the mechanism of Alkbh5 in RO may be little different from that in PA, as some upregulated and downregulated methylation peaks are also related to chromatin modification and histone modification (Figures 3C,D). We performed REACTOME-pathway analysis to further study the potential reason for the reproduction dysfunction of Alkbh5-KO mice. In the comparison between Alkbh5-KO and WT PA, we found that both the upregulated and downregulated methylation peaks were related to chromatin modifying enzymes and chromatin organization. Moreover, we found that the upregulated methylation peaks were related to the M phase. Comparing Alkbh5-KO and WT RO, we found that the upregulated methylation peaks are more related to mitotic prometaphase and that the downregulated methylation peaks are more related to RHO GTPase cycle (Supplementary Figures S3).

FIGURE 3

FIGURE 3. Differentially methylated mRNAs analysis. GO enrichment analysis of upregulated (A) and downregulated (B) methylated mRNAs between Alkbh5-KO PA and WT PA. GO enrichment analysis of upregulated (C) and downregulated (D) methylated mRNAs between Alkbh5-KO RO and WT RO.

Overview of Transcriptome Profiles Revealed Important Mechanisms in a Spermatogenetic Malfunction

Using the datasets from RNA-seq results (MeRIP-seq input library), we compared the transcriptome profiles of Alkbh5-KO with WT spermatogenic cells. We found no remarkable difference in the density and distribution of expression levels between these samples, which could comprehensively indicate that the sequencing result is reliable (Figure 4A). We found 785 upregulated and 440 downregulated genes in PA (|log2FC| > 1 and FDR <0.05) using the R package DESeq2. We also found 564 upregulated genes and 539 downregulated genes in the RO. The top 10 differentially expressed genes in PA and RO are presented in Supplementary Table S3. In brief, the number of upregulated genes in Alkbh5-KO PA was much greater than that in WT PA. However, the number of upregulated genes in Alkbh5-KO RO is similar to that in WT RO (Figures 4B,C). Thereafter, we performed GO analysis of the DEGs between the Alkbh5-KO and WT spermatogenic cells. GO analysis can be divided into two parts: BP and MF. First, we considered the DEGs between the Alkbh5-KO PA and WT PA. From the BP results, we found that these DEGs are closely related to reproduction processes, such as spermatid differentiation, spermatid development and germ cell development. In addition, the upregulated genes are more related to spermatogenesis and the downregulated genes are more related to the meiotic cycle. From this, it may be speculated that knocking out Alkbh5 may cause damage to the reproductive functions of mice. From the MF results, we found that the upregulated genes were linked with protein disulphide isomerase activity and that the downregulated genes were linked with dynein light intermediate chain binding (Figures 5A,B). Next, we compared the DEGs closely related to reproduction processes between Alkbh5-KO RO and WT RO. However, we found that the upregulated genes were more related to protein binding and that the downregulated genes were more related to ATPase activity (Figures 5C,D). Additionally, we screened out DEGs and performed REACTOME enrichment analysis in PA and RO. In the PA and RO, we found that both the downregulated genes were related to fertilization and reproduction. Moreover, we found that the upregulated genes were related to the M phase. Comparing Alkbh5-KO and WT RO, we found that the upregulated genes are more related to stress responses (Supplementary Figure S4).

FIGURE 4

FIGURE 4. Differential gene analysis FPKM density distribution (A). Volcano plot showing differential genes between Alkbh5-KO PA and WT PA (B), Alkbh5-KO RO and WT RO (C).

FIGURE 5

FIGURE 5. GO enrichment analysis of upregulated (A) and downregulated (B) genes between Alkbh5-KO PA and WT PA. GO enrichment analysis of upregulated (C) and downregulated (D) genes between Alkbh5-KO RO and WT RO.

The Expression of Genes Related to Spermatogenesis is Markedly Altered in Alkbh5 KO

To further investigate the influence of Alkbh5 on spermatogenesis, we performed an association analysis of the genes with significant changes in methylation and mRNA expression levels. The results were divided into four parts: hypermethylated-down, hypermethylated-up, hypomethylated-down, and hypomethylated-up. Figure 6A shows the top 10 differentially methylated genes. To discover the relationship between differential peaks and genes, we counted the number of upregulated and downregulated mRNAs and peaks in KO in PA and RO and calculated their intersection. We found a total of 54 upregulated overlapped genes and 10 downregulated overlapped genes in PA. However, in RO, 37 genes were upregulated, and 55 were downregulated. This revealed the relatively large differences between PA and RO (Figure 6B). Next, we selected the top 20 differential protein-coding genes in Alkbh5-KO PA and WT PA (Table 1), as well as in Alkbh5-RO and WT RO (Table 2). Thereafter, we intersected the 20 genes in each table. We found that among the top 20 DEGs in each table, there were approximately 11 genes that appeared in both comparisons. We found that these genes were related to the development of spermatogenesis. For this reason, we used qPCR to discover the comprehensive changes in these genes, and whether they exist in the whole testis. We found that Mroh4, Prm1, Lrrc7, Dbil5, H1fnt, Catsperd, and Odf4 were greatly downregulated in the testes in the Alkbh5 KO (Figure 6C).

FIGURE 6

FIGURE 6. Association analysis of the genes with significant changes in methylation and mRNA expression level in PA and RO (A). Venn diagrams showed the intersection of upregulated and downregulated methylation peaks and mRNAs of PA and RO respectively (B). mRNA expression level of seven genes in Alkbh5-KO and WT mice Testicular tissue detected by qPCR (****p < 0.0001, ***p = 0.0001, **p < 0.01.) (C).

TABLE 1

TABLE 1. The top 20 differential protein-coding genes in Alkbh5-KO PA and WT PA.

TABLE 2

TABLE 2. The top 20 differential protein-coding genes in Alkbh5-KO RO and WT RO.

Discussion

Modification of m6A demethylases plays a critical role in mammalian spermatogenesis. ALKBH5 was highly expressed in mouse testis as an m6A eraser, and loss of ALKBH5 significantly increased the level of m6A in testicular cells. Alkbh5-KO mice show abnormal spermatogenesis and subfertility (Zheng et al., 2013). An independent study suggested that correct splicing of longer 3ʹ-UTR transcripts may require ALKBH5-mediated m6A demethylation (Tang et al., 2018). However, previous studies on Alkbh5 only focused on the overall changes in spermatogenesis but not on the specific changes in RNA and m6A methylation changes at each stage of spermatogenesis. To address this question, we performed high-throughput sequencing to reveal the m6A landscape in the spermatogenic stages of PA and RO in Alkbh5-KO and WT mice. Using MeRIP-seq, we identified that Alkbh5-KO PA had 1,929 noticeably upregulated methylation peaks relative to that in the WT PA, while 326 were downregulated. The Alkbh5-KO RO had 1,589 noticeably upregulated methylation peaks relative to that in the WT RO, while 1,602 were downregulated. We then used GO and REACTOME pathway enrichment analyses to explore the potential functions of differentially methylated transcripts. Finally, we found a difference in the genes in methylation peaks and synchronous differential expression genes enriched in the spermatogenesis-related pathway.

In this study, since the analytical results of MeRIP-seq include the results of RNA-seq, our analysis reveals the impact of Alkbh5 on spermatogenesis from the perspective of differences in methylation and gene expression. We found that the effects of Alkbh5 on PA and RO differed considerably. First, we found that ALKBH5 is a demethylase has different effects on the methylation of transcripts in PA and RO. In PA, we found that the number of significantly upregulated methylation peaks was far larger than that of significantly downregulated methylation peaks in the Alkbh5 KO, which is consistent with our previous knowledge about Alkbh5 that the degree of RNA methylation significantly increases in the testes of Alkbh5 KOs (Tang et al., 2018). However, in RO, we found that in the Alkbh5 KO, the number of significantly upregulated methylation peaks was very close to the number of significantly downregulated methylation peaks. We found that the distribution of methylation peaks in PA and RO differed considerably. In particular, in the Alkbh5 knockout, we found that the proportion of methylation peaks distributed in the CDS region of PA increased significantly, while the proportion of methylation peaks distributed in the CDS region of RO reduced. Moreover, lncRNAs, as epigenetic regulators, play an essential role in spermatogenesis (Wichman et al., 2017), such as the definition of sperm parameters (An et al., 2017) and sperm capacitation (Song et al., 2011). Our observations regarding the distribution of methylation peaks of lncRNA, PA, and RO were also affected by knocking out of Alkbh5. The expression of lncRNA was characterized by a sharp increase in the number of expression sites and the overall expression level with the development of spermatogenesis, and reached a peak in the round spermatogenesis stage (Davis et al., 2017).

Here, we used GO enrichment analysis to further explore differentially methylation peaks. We found that after Alkbh5 knockout, both the upregulated and downregulated methylation peaks in PA and RO were closely related to covalent chromatin modification and histone modification. Histone modification influences chromatin structure and function. Previous studies have indicated various histone modifications are involved in the regulation of spermatogenesis (Okada et al., 2007; Liu et al., 2010; Okada et al., 2010; Iwamori et al., 2011), but the role of covalent chromatin modification in spermatogenesis remains unclear. However, recent studies demonstrated that m6A modification activates gene transcription by altering histone modifications and chromatin structure (Wang et al., 2018). In mouse embryonic stem cells, knocking out METTL3 reduced m6A modification levels and stabilized opening of chromatin that facilitates and sustains gene transcription (Liu et al., 2020). Based on our results, we propose the following conjecture: m6A modification level and related binding proteins on the mRNA of spermatogenic cells at PA and RO can together regulate open state of chromatin and the transcription level of downstream genes, which play an important role in the smooth progress of spermatogenesis, and an abnormal regulation can lead to male infertility. Moreover, genes with upregulated m6A-modified sites in PA were related to activation of MAPK activity and downregulated m6A-modified sites were related to the negative regulation of the MAPK cascade. The MAPK cascade plays important roles in male reproductive processes, such as spermatogenesis, sperm maturation, and the acrosome reaction (Li et al., 2009). TCP1, which is involved in the regulation of the MAPK cascade, is involved in the assembly of actin and tubulin filaments (Dun et al., 2012). Actin filaments exist in mammalian reproductive cells and participate in some of the changes that occur during spermatogenesis, such as sperm maturation and capacitation (Howes et al., 2001). Our results revealed that changes in the m6A modifications of the MAPK signalling pathway may affect sperm maturation and capacitation. Taken together, these factors may contribute to fertility in KO mice.

In the DEG analysis, we detected several the DEGs related to many important biological pathways. REACTOME pathway analysis revealed that transcripts related to fertilization and reproduction genes were downregulated in PA and RO during spermatogenesis, showing that the expression of genes related to spermatogenesis in PA and RO both are affected by the increase in methylation level after knockout of the methylase. This further illustrates that m6A methylation is necessary during spermatogenesis. Using MeRIP-seq results, we analysed the differences from the perspective of RNA. The RNA differential analysis results were similar to the methylation results. We found that the number of upregulated genes in PA was significantly greater than that in the downregulated genes in the Alkbh5 knockout. However, the number of upregulated and downregulated genes in RO was close. In the enrichment analysis, we found that the differential genes of PA and RO are closely related to spermatogenesis. Therefore, this result reveals that knocking out Alkbh5 affects gene expression in PA and RO. Further molecular experiments should be performed to verify the differentially methylated genes in mouse PA and RO. Next, we found that the top 20 DEGs in Alkbh5-KO PA, WT PA, Alkbh5-RO, and WT RO were closely related to spermatogenesis. Moreover, after intersecting these genes, we found that more than half of these genes appeared in both comparisons. In the qPCR experiment, we found that transcripts related to spermatogenesis were significantly downregulated because of Alkbh5 KO. One of these was, Prm1, which is a protamine substitute for histones in the chromatin of sperm during the haploid phase of spermatogenesis, compacts sperm DNA into highly condensed, stable, and inactive complexes (Fagerberg et al., 2014). In addition, we found Odf4, which encodes a protein that is localized in the outer dense fibres of the tails of mature sperm. This protein is thought to play an important role in the sperm tail (Nakamura et al., 2002). Furthermore, H1FNT, which was also downregulated in Alkbh5 KO, is a germ cell-specific linker histone variant expressed during spermiogenesis, specifically in round and elongating spermatids (Tanaka et al., 2006). In conclusion, Alkbh5 greatly contributes to the development of spermatogenesis at both the epigenetic and mRNA levels.

While the causes of approximately 50% of human male infertility are unknown, genetic causes, in particular, often remain undiagnosed (Zhang et al., 2011). Our current mouse study provides basis to a genetic cause for male infertility. Alkbh5 is also widely expressed in the human testis. In humans, Alkbh5 expression is higher in PA compared with other stages, implying that it affects semen quality and a perturbation in its function increases the risk of male infertility (Landfors et al., 2016). However, throughout evolution, humans have also developed features in spermatogenesis that are distinguishable from those of mice. It is also unknown whether the m6A methylation level and mRNA expression level in PA and RO will be the same in humans as in mice. This study sheds light on possible targets for future research into developing treatments for certain types of human infertility.

In summary, we analysed the differential m6A methylome in PA and RO in Alkbh5-KO mice relative to the corresponding WT mice, and demonstrated a strong association between m6A modification and spermatogenesis. Furthermore, in the Alkbh5 KO, m6A levels were increased, affecting many spermatogenesis-related genes, very likely leading to spermatogenesis disorders. In addition, the data show that the m6A methylation level and mRNA expression level in mouse PA and RO were significantly different; further research is needed to uncover the underlying mechanisms that affect male fertility.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: BioProject: PRJNA761579, SRA accession: SRP336482, GEO: GSE186217.

Ethics Statement

The animal study was reviewed and approved by the Institutional Review Board of Nantong University (Nantong, China).

Author Contributions

SH and XS participated in the study design, statistical analyses, and writing of the paper. HT assisted with data collection. JC reviewed and edited manuscript. FS revised the manuscript and were responsible for supervision of study.

Funding

This work was supported by grants from the National Key Research and Development Program of China (No. 2021YFC2700200 to FS) and Postgraduate Research and Practice Innovation Program of Jiangsu Province (KYCX20_2841).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Dr. Yamei Niu from the Institute of Basic Medical Sciences, Chinese Academy of Medical Science, for her generous gift of the heterozygous Alkbh5 mice. We would like to thank Home for Researchers (www.home-for-researchers.com) for paper writing. We would also like to thank Editage (www.editage.cn) for English language editing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.832677/full#supplementary-material

References

A. Alemu, E., He, C., and Klungland, A. (2016). ALKBHs-Facilitated RNA Modifications and De-modifications. DNA Repair 44, 87–91. doi:10.1016/j.dnarep.2016.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahmed, E. A., and de Rooij, D. G. (2009). Staging of Mouse Seminiferous Tubule Cross-Sections. Methods Mol. Biol. 558, 263–277. doi:10.1007/978-1-60761-103-5_16

PubMed Abstract | CrossRef Full Text | Google Scholar

Alarcón, C. R., Goodarzi, H., Lee, H., Liu, X., Tavazoie, S., and Tavazoie, S. F. (2015). HNRNPA2B1 Is a Mediator of m6A-dependent Nuclear RNA Processing Events. Cell 162 (6), 1299–1308. doi:10.1016/j.cell.2015.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

An, T., Fan, H., Liu, Y. F., Pan, Y. Y., Liu, Y. K., Mo, F. F., et al. (2017). The Difference in Expression of Long Noncoding RNAs in Rat Semen Induced by High-Fat Diet Was Associated with Metabolic Pathways. PeerJ 5, e3518. doi:10.7717/peerj.3518

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25 (1), 25–29. doi:10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Bryant, J. M., Meyer-Ficca, M. L., Dang, V. M., Berger, S. L., and Meyer, R. G. (2013). Separation of Spermatogenic Cell Types Using STA-PUT Velocity Sedimentation. JoVE 80, 1. doi:10.3791/50648

CrossRef Full Text | Google Scholar

Davis, M. P., Carrieri, C., Saini, H. K., Dongen, S., Leonardi, T., Bussotti, G., et al. (2017). Transposon‐driven Transcription Is a Conserved Feature of Vertebrate Spermatogenesis and Transcript Evolution. EMBO Rep. 18 (7), 1231–1247. doi:10.15252/embr.201744059

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the Human and Mouse m6A RNA Methylomes Revealed by m6A-Seq. Nature 485 (7397), 201–206. doi:10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

Draghici, S., Khatri, P., Tarca, A. L., Amin, K., Done, A., Voichita, C., et al. (2007). A Systems Biology Approach for Pathway Level Analysis. Genome Res. 17 (10), 1537–1545. doi:10.1101/gr.6202607

PubMed Abstract | CrossRef Full Text | Google Scholar

Dun, M. D., Aitken, R. J., and Nixon, B. (2012). The Role of Molecular Chaperones in Spermatogenesis and the post-testicular Maturation of Mammalian Spermatozoa. Hum. Reprod. Update 18 (4), 420–435. doi:10.1093/humupd/dms009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ensfelder, T. T., Kurz, M. Q., Iwan, K., Geiger, S., Matheisl, S., Müller, M., et al. (2018). ALKBH5-induced Demethylation of Mono- and Dimethylated Adenosine. Chem. Commun. 54 (62), 8591–8593. doi:10.1039/c8cc03980a

CrossRef Full Text | Google Scholar

Fagerberg, L., Hallström, B. M., Oksvold, P., Kampf, C., Djureinovic, D., Odeberg, J., et al. (2014). Analysis of the Human Tissue-specific Expression by Genome-wide Integration of Transcriptomics and Antibody-Based Proteomics. Mol. Cell Proteomics 13 (2), 397–406. doi:10.1074/mcp.M113.035600

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Dominissini, D., Rechavi, G., and He, C. (2014). Gene Expression Regulation Mediated through Reversible m6A RNA Methylation. Nat. Rev. Genet. 15 (5), 293–306. doi:10.1038/nrg3724

PubMed Abstract | CrossRef Full Text | Google Scholar

Griswold, M. D. (2016). Spermatogenesis: The Commitment to Meiosis. Physiol. Rev. 96 (1), 1–17. doi:10.1152/physrev.00013.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Handel, M. A., and Schimenti, J. C. (2010). Genetics of Mammalian Meiosis: Regulation, Dynamics and Impact on Fertility. Nat. Rev. Genet. 11 (2), 124–136. doi:10.1038/nrg2723

PubMed Abstract | CrossRef Full Text | Google Scholar

Howes, E. A., Hurst, S. M., and Jones, R. (2001). Actin and Actin-Binding Proteins in Bovine Spermatozoa: Potential Role in Membrane Remodeling and Intracellular Signaling during Epididymal Maturation and the Acrosome Reaction. J. Androl. 22 (1), 62–72.

PubMed Abstract | Google Scholar

Iwamori, N., Zhao, M., Meistrich, M. L., and Matzuk, M. M. (2011). The Testis-Enriched Histone Demethylase, KDM4D, Regulates Methylation of Histone H3 Lysine 9 during Spermatogenesis in the Mouse but Is Dispensable for Fertility1. Biol. Reprod. 84 (6), 1225–1234. doi:10.1095/biolreprod.110.088955

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, G., Fu, Y., and He, C. (2013). Reversible RNA Adenosine Methylation in Biological Regulation. Trends Genet. 29 (2), 108–115. doi:10.1016/j.tig.2012.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12 (4), 357–360. doi:10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Landfors, M., Nakken, S., Fusser, M., Dahl, J.-A., Klungland, A., and Fedorcsak, P. (2016). Sequencing of FTO and ALKBH5 in Men Undergoing Infertility Work-Up Identifies an Infertility-Associated Variant and Two Missense Mutations. Fertil. Sterility 105 (5), 1170–1179. e1175. doi:10.1016/j.fertnstert.2016.01.002

CrossRef Full Text | Google Scholar

Li, M. W. M., Mruk, D. D., and Cheng, C. Y. (2009). Mitogen-activated Protein Kinases in Male Reproductive Function. Trends Mol. Med. 15 (4), 159–168. doi:10.1016/j.molmed.2009.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R.-K., Tan, J.-L., Chen, L.-T., Feng, J.-S., Liang, W.-X., Guo, X.-J., et al. (2014). Iqcg Is Essential for Sperm Flagellum Formation in Mice. PLoS One 9 (5), e98053. doi:10.1371/journal.pone.0098053

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z., Hsu, P. J., Xing, X., Fang, J., Lu, Z., Zou, Q., et al. (2017). Mettl3-/Mettl14-mediated mRNA N6-Methyladenosine Modulates Murine Spermatogenesis. Cell Res 27 (10), 1216–1230. doi:10.1038/cr.2017.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Dou, X., Chen, C., Chen, C., Liu, C., Xu, M. M., et al. (2020). N 6 -methyladenosine of Chromosome-Associated Regulatory RNA Regulates Chromatin State and Transcription. Science 367 (6477), 580–586. doi:10.1126/science.aay6018

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Zhou, S., Liao, L., Chen, X., Meistrich, M., and Xu, J. (2010). Jmjd1a Demethylase-Regulated Histone Modification Is Essential for cAMP-Response Element Modulator-Regulated Gene Expression and Spermatogenesis. J. Biol. Chem. 285 (4), 2758–2770. doi:10.1074/jbc.M109.066845

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, J., Cui, X., Rao, M. K., Chen, Y., and Huang, Y. (2013). Exome-based Analysis for RNA Epigenome Sequencing Data. Bioinformatics 29 (12), 1565–1567. doi:10.1093/bioinformatics/btt171

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, J., Lu, Z., Liu, H., Zhang, L., Zhang, S., Chen, Y., et al. (2014). A Protocol for RNA Methylation Differential Analysis with MeRIP-Seq Data and exomePeak R/Bioconductor Package. Methods 69 (3), 274–281. doi:10.1016/j.ymeth.2014.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, K. D., Saletore, Y., Zumbo, P., Elemento, O., Mason, C. E., and Jaffrey, S. R. (2012). Comprehensive Analysis of mRNA Methylation Reveals Enrichment in 3′ UTRs and Near Stop Codons. Cell 149 (7), 1635–1646. doi:10.1016/j.cell.2012.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, Y., Tanaka, H., Koga, M., Miyagawa, Y., Iguchi, N., de Carvalho, C. E., et al. (2002). Molecular Cloning and Characterization of Oppo 1: a Haploid Germ Cell-specific Complementary DNA Encoding Sperm Tail Protein. Biol. Reprod. 67 (1), 1–7. doi:10.1095/biolreprod67.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Naro, C., Jolly, A., Di Persio, S., Bielli, P., Setterblad, N., Alberdi, A. J., et al. (2017). An Orchestrated Intron Retention Program in Meiosis Controls Timely Usage of Transcripts during Germ Cell Differentiation. Dev. Cell 41 (1), 82–93. e84. doi:10.1016/j.devcel.2017.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Donnell, L., Nicholls, P. K., O’Bryan, M. K., McLachlan, R. I., and Stanton, P. G. (2011). Spermiation. Spermatogenesis 1 (1), 14–35. doi:10.4161/spmg.1.1.14525

PubMed Abstract | CrossRef Full Text | Google Scholar

Oatley, J. M., and Brinster, R. L. (2008). Regulation of Spermatogonial Stem Cell Self-Renewal in Mammals. Annu. Rev. Cell Dev. Biol. 24, 263–286. doi:10.1146/annurev.cellbio.24.110707.175355

PubMed Abstract | CrossRef Full Text | Google Scholar

Oatley, J. M., and Brinster, R. L. (2012). The Germline Stem Cell Niche Unit in Mammalian Testes. Physiol. Rev. 92 (2), 577–595. doi:10.1152/physrev.00025.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Okada, Y., Scott, G., Ray, M. K., Mishina, Y., and Zhang, Y. (2007). Histone Demethylase JHDM2A Is Critical for Tnp1 and Prm1 Transcription and Spermatogenesis. Nature 450 (7166), 119–123. doi:10.1038/nature06236

PubMed Abstract | CrossRef Full Text | Google Scholar

Okada, Y., Tateishi, K., and Zhang, Y. (2010). Histone Demethylase JHDM2A Is Involved in Male Infertility and Obesity. J. Androl. 31 (1), 75–78. doi:10.2164/jandrol.109.008052

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, N., Liu, J., An, S., Nishino, T., Hishikawa, Y., and Koji, T. (2011). Immunohistochemical Analysis of Histone H3 Modifications in Germ Cells during Mouse Spermatogenesis. Acta Histochem. Cytochem. 44 (4), 183–190. doi:10.1267/ahc.11027

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, R., Ro, S., Michaels, J. D., Park, C., McCarrey, J. R., and Yan, W. (2009). Many X-Linked microRNAs Escape Meiotic Sex Chromosome Inactivation. Nat. Genet. 41 (4), 488–493. doi:10.1038/ng.338

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, H., Matsuoka, Y., Onishi, M., Kitamura, K., Miyagawa, Y., Nishimura, H., et al. (2006). Expression Profiles and Single-Nucleotide Polymorphism Analysis of Human HANP1/H1T2 Encoding a Histone H1-like Protein. Int. J. Androl. 29 (2), 353–359. doi:10.1111/j.1365-2605.2005.00600.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, C., Klukovich, R., Peng, H., Wang, Z., Yu, T., Zhang, Y., et al. (2018). ALKBH5-dependent m6A Demethylation Controls Splicing and Stability of Long 3′-UTR mRNAs in Male Germ Cells. Proc. Natl. Acad. Sci. U.S.A. 115 (2), E325–E333. doi:10.1073/pnas.1717794115

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Lu, Z., Gomez, A., Hon, G. C., Yue, Y., Han, D., et al. (2014). N6-methyladenosine-dependent Regulation of Messenger RNA Stability. Nature 505 (7481), 117–120. doi:10.1038/nature12730

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Zhao, B. S., Roundtree, I. A., Lu, Z., Han, D., Ma, H., et al. (2015). N6-methyladenosine Modulates Messenger RNA Translation Efficiency. Cell 161 (6), 1388–1399. doi:10.1016/j.cell.2015.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Li, Y., Yue, M., Wang, J., Kumar, S., Wechsler-Reya, R. J., et al. (2018). N6-methyladenosine RNA Modification Regulates Embryonic Neural Stem Cell Self-Renewal through Histone Modifications. Nat. Neurosci. 21 (2), 195–206. doi:10.1038/s41593-017-0057-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, C.-M., Gershowitz, A., and Moss, B. (1975). Methylated Nucleotides Block 5′ Terminus of HeLa Cell Messenger RNA. Cell 4 (4), 379–386. doi:10.1016/0092-8674(75)90158-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wichman, L., Somasundaram, S., Breindel, C., Valerio, D. M., McCarrey, J. R., Hodges, C. A., et al. (2017). Dynamic Expression of Long Noncoding RNAs Reveals Their Potential Roles in Spermatogenesis and Fertility†. Biol. Reprod. 97 (2), 313–323. doi:10.1093/biolre/iox084

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, Y., Liu, J., and He, C. (2015). RNA N6-Methyladenosine Methylation in post-transcriptional Gene Expression Regulation. Genes Dev. 29 (13), 1343–1355. doi:10.1101/gad.262766.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, Y., Wang, S., Gao, S., Soares, F., Ahmed, M., Guo, H., et al. (2018). Refined RIP-Seq Protocol for Epitranscriptome Analysis with Low Input Materials. Plos Biol. 16 (9), e2006092. doi:10.1371/journal.pbio.2006092

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Zhou, D.-x., Wang, H.-x., and Tian, Z. (2011). An Association Study of SPO11 Gene Single Nucleotide Polymorphisms with Idiopathic Male Infertility in Chinese Han Population. J. Assist. Reprod. Genet. 28 (8), 731–736. doi:10.1007/s10815-011-9571-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B. S., Roundtree, I. A., and He, C. (2017). Post-transcriptional Gene Regulation by mRNA Modifications. Nat. Rev. Mol. Cell Biol 18 (1), 31–42. doi:10.1038/nrm.2016.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, G., Dahl, J. A., Niu, Y., Fedorcsak, P., Huang, C.-M., Li, C. J., et al. (2013). ALKBH5 Is a Mammalian RNA Demethylase that Impacts RNA Metabolism and Mouse Fertility. Mol. Cell 49 (1), 18–29. doi:10.1016/j.molcel.2012.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar