Transcriptome-Wide m6A Analysis Provides Novel Insights Into Testicular Development and Spermatogenesis in Xia-Nan Cattle

Introduction

Xia-Nan (XN) cattle, the first specialized beef breed produced by the crossbreeding of French Charolais (male) and Nanyang cattle (female) in China, has important features including fast growth rate and high meat production performance. Enhancing the reproductive performance of XN cattle is an important breeding objective to increase the efficiency and sustainability of this breed as a major beef producer in the Chinese meat market. Testis, as a basic male reproductive organ, plays a critical role in spermatogenesis and steroidogenesis. Bull’s fertility has been always considered a key issue for both cattle breeders and scientists. Different methods and common practices are traditionally used to check male reproductive performance, including a physical assessment of the bull (e.g., volume of the testicles) and semen evaluation (e.g., sperm morphology, concentration, and motility). Although these procedures offered a considerable background on testicular and epididymal function and quantitative production of sperms, infertility due to major reasons was not well recognized. Nowadays, molecular genetics techniques can efficiently target semen quality providing a powerful approach to evaluate the fertility potential of mammalian males (Puglisi et al., 2016; Vendelbo et al., 2021). Establishing an association between spermatogenesis and gene expression profile may enhance a better understanding of the causes of infertility. Improving the testicular function for acceptable production yield of high-quality semen (e.g., healthy spermatozoa) is still the basis for increasing the opportunity for the success rate of obtaining healthy offspring. Spermatogenesis can be divided into three major functional stages including the proliferative stage, meiosis, and maturation stage (Griswold, 2016). Notably, spermatogenesis is strictly regulated by the expression of stage-specific genes at both transcription and post-transcription levels (Ding et al., 2020). Thus, the testes at different developmental stages in XN cattle were taken as the experimental object in the present study, taking into consideration that identifying key regulators and signaling pathways related to testis development and spermatogenesis will provide valuable insights into improving semen quality.

N⁶-methyladenosine (m6A) is the most abundant internal RNA modification that has been recently suggested as a critical post-transcriptional mRNA regulator in most organisms, and is strongly associated with production traits (Shi et al., 2019; Zhang et al., 2020b; Xiong et al., 2021). In mammals, the RNA m6A modification is installed by a methyltransferase complex that mainly includes methyltransferase-like 3 (METTL3) and METTL14, which are responsible for catalyzing m6A modification (Liu et al., 2014). Wilms’ tumor 1-associated protein (WTAP) is another essential member of the core component that interacts with METTL3 and is required for its localization in the nucleus (Schwartz et al., 2014). RNA m6A modification can be removed by two known demethylases: AlkB homolog 5 (ALKBH5) (Zheng et al., 2013) and fat mass and obesity-associated factor (FTO) (Jia et al., 2011). RNA m6A modification affects almost all stages of RNA metabolism, like alternative splicing, RNA degradation, nuclear RNA export, and translation (Niu et al., 2013), and these influences can be recognized by a category of proteins primarily composed of the YTH domain family (YTHDF1-3) and IGF2BPs (IGF2BP1-3) (Dominissini et al., 2012; Fu et al., 2014; Wang et al., 2014). Also, m6A modification can influence the spermatogenic function, for example, the combined deletion of METTL3 and METTL14 leads to impaired murine spermiogenesis (Lin et al., 2017). Moreover, m6A demethylase ALKBH5 deficiency causes compromised spermatogenesis and apoptosis in mouse testis (Zheng et al., 2013). The YTHDC2 is an m6A-binding protein gene and its knockout is related to infertile in mice (Hsu et al., 2017). Considering the various functions of m6A modification mentioned above in different species, it would seem logical to assume that m6A modification may also affect bovine testis development. During the last few years, continuously incremented works were carried out to investigate the mechanisms of spermatogenesis of cattle by comparative analysis of the gene expression associated with reproductive traits on the molecular level. However, little is known until now about the potential impacts of m6A modification on spermatogenesis in ruminants, in general, and there are no reports that document the association of m6A modification and bovine testis development in XN cattle, in specific.

Therefore, to identify the functional m6A and explore the spermatogenesis mechanism of the testes caused by the m6A RNA modification, we detected the m6A methylomes of bovine testes from birth to adulthood. This enabled us to obtain the transcriptome-wide m6A profiles in bovine testes at different stages and to provide reasonable insights into the roles of m6A modification in the mechanisms underlying spermatogenesis during post-natal testicular development.

Materials and Methods

Experiment Design and Sample Collection

The main steps and bioinformatics used for data analysis in the present study is shown in Supplementary Figure S1. A total of sixty-nine clinically healthy XN cattle bulls from Kerchin Cattle Industry (Nanyang) Co., China was randomly selected. These animals were grouped into three groups based on their sexual maturity. The first group indicates the bulls are prepuberty (at birth, n = 23), the second group represents the bulls are puberty (about 1 year old showing heat for the first time, n = 23), and the last group represents the bulls are postpuberty (about 2 years of age, n = 23). In addition, the semen within 1 month of the bulls in puberty and postpuberty was taken by using artificial vagina. Sperm concentration was assessed using a haemocytometer (Bane, 1952). Sperm motility was evaluated based on the methods described (Björndahl et al., 2003). Further, the t-test was used to determine the significant difference levels of the semen quality parameters (semen volume, sperm motility and sperm concentration) between the puberty and postpuberty groups. For the puberty and postpuberty groups, the representative individuals of each group were selected in two ways: 1) the semen parameters existed significance differences between them; 2) the semen quality parameters closed to mean value for each group. In the prepuberty, the samples were randomly chosen because the bulls have no semen parameters. Finally, a total of nine bulls (triplicate in each group) were selected based on the above-mentioned criterion. The selected animals were given general anesthesia (Zoletil 50, Virbac Co., France), a combination of zolazepam (5–9 mg/kg, i.m) and tiletamine, and xylazine hydrochloride (1.5–2 mg/kg, i.m) before sampling. Samples of the left testes were collected posterior to castration that performed by professional veterinarians. The length, width, and weight of the left testes were measured. Each testis was divided into three pieces and immediately subjected to snap-freezing in liquid nitrogen and stored at −80°C until RNA extraction.

RNA Isolation, Library Construction and Sequencing

Total RNA for each testis sample was isolated using the Trizol reagent (Invitrogen, CA, United States) following the manufacturer’s protocols. The quality and quantity of total RNA were determined by Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., CA, United States) and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, United States) with RIN number >7.0. Poly(A) RNA from total RNA was isolated with Arraystar Seq-Star™ poly(A) mRNA Isolation Kit (Arraystar, MD, United States). The RNA was further fragmented into fragments with an average length of 100 nt using RNA Fragmentation Reagents (Sigma, MO, United States). The fragmented RNA segments were divided into two groups. One group was used to perform m6A RNA immunoprecipitation (IP) using the GenSeq™ m6A RNA IP Kit (GenSeq Inc., China). The other group was utilized to construct the input samples without immunoprecipitation. The IP and input libraries were both constructed with NEBNext^® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., MA, United States). The quality of all libraries was measured by the Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., CA, United States). The sequencing of nine cDNA libraries was performed on the Illumina HiseqTM 4000 by Gene Denovo Biotechnology Co., Ltd (Guangzhou, China). The raw data were deposited in the NCBI SRA database (BioProject ID: PRJNA776655).

Bioinformatics Analysis of m6A-Seq and RNA-Seq Data

TrimGalore v0.6.6 (Krueger, 2021) software was used to eliminate the reads containing adaptor contaminants, low-quality bases, and undetermined bases. Meanwhile, the sequence quality of IP and input of all samples were validated by the FASTP v0.20.1 (Chen et al., 2018) software. Subsequently, the high-quality clean reads were mapped to the Bos taurus reference genome (ARS-UCD1.2) by HISAT2 ver.2.1.0 (Kim et al., 2015) software with default parameters. Peak calling for mapped reads of IP and input libraries were performed using the exomePeak2 (Meng et al., 2014) package in R. The m6A intensity was visualized using the IGV software (http://www.igv.org/). The identified m6A peaks were carried out to conduct the motif enrichment analysis by MEME (Bailey et al., 2009) and HOMER (Sven et al., 2010) software with default parameters. The ChIPseeker v1.0 (Guangchuang et al., 2015) software was used to annotate the called peaks by intersection with gene architecture. The differential m6A peaks [fold changes (FC) ≥ 2 and adjusted p-value < 0.05] between the pairwise comparison groups were analyzed using the exomePeak2 package in R. These differential peaks were annotated using the Ensembl database (Bos taurus/ARS-UCD1.2). Moreover, the expression level for all mRNAs from input libraries was calculated using StringTie ver. 1.3.5 (Pertea et al., 2016) software. The expression level of each transcript was normalized by the Trimmed Mean of M-values (TMM) implemented in the edgeR R-package. The differential expression analysis for pairwise contrasts was performed using the DESeq2 (Love et al., 2014) package in R. The adjusted p-value ≤ 0.05 and FoldChange >1.5 were defined as the cutoff criteria for the differentially expressed mRNAs (DEGs). Gene enrichment analysis was performed by the Gene Ontology (GO) functional analysis and the Kyoto Encyclopedia of Genes and Genomics (KEGG) pathway enrichment using the KEGG Orthology-Based Annotation System (KOBAS) 3.0 with cutoff criteria of p ≤ 0.05, aiming to identify their biological significance. The plot results were visualized using the ggplot2 (Wickham, 2016) package in R.

Quantitative Real-Time PCR Confirmation

Four differentially m6A methylated (DMGs) genes were selected and analyzed by qRT-PCR. Primers were designed using Primer 5.0 software (Supplementary Table S1) and synthesized by Sangon Biotech (Shanghai) Co. Ltd. RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States) was used to reverse transcribe the total RNAs into cDNA following the manufacturer’s protocols. Then qPCR was conducted using QuantiNova SYBR Green PCR Kit (QIAGEN, Shanghai, China). The GAPDH gene was used for normalizing the relative abundance of genes. The 2^−ΔΔCt method (Livak and Schmittgen, 2001) was used to analyze the data for all samples in triplicate technical replicates.

Statistical Analysis

The data are expressed as mean ± standard error of the mean (SEM). The student’s t-test and one-way analysis of variance (ANOVA) was performed to determine the significance of the differences between the contrasting groups by using Graphpad Prism 8 software. Differences between means were considered statistically significant when adjusted pairwise comparison between means reached p-value ≤ 0.05 (Bonferroni).

Results

Testis Source Description

Semen quality parameters of the bulls in puberty and postpuberty, including the semen volume, sperm motility and sperm concentration were analyzed, and their results was listed in Table 1. Our data showed the semen volume and sperm motility in postpuberty were markedly higher (p < 0.05), while sperm concentration was not different as compared to bulls in puberty. The Supplementary Figure S2 further showed the detailed information on the semen volume, sperm motility and sperm concentration of bulls in puberty and postpuberty groups. Meanwhile, we selected six animals (three bulls for each group) to subsequently testes collection according to the principle that semen quality parameters closed to mean value. In addition, 3 bulls in prepuberty were randomly chosen.

TABLE 1

TABLE 1. Evaluations of semen quality and testicular phenotypic parameters of XN bull testes at different ages (mean ± SEM).

For the above-selected 9 bulls, we measured their weight, length, and width of left testes, and their estimates are listed in Table 1. The results showed that all the tested testicular phenotypic parameters varied significantly (p < 0.05) among the groups studied. According to their semen quality and phenotypic difference, 9 testes samples of XN males at three developmental stages (prepuberty = TY0, in puberty = TY1, and postpuberty = TY2) were employed for further m6A-seq analysis. Each stage was in triplicate.

Transcriptome-Wide m6A-Seq Revealed m6A Modification Patterns During Testis Development

In the present study, the bull testes in prepuberty, puberty and postpuberty were used for m6A-Seq and RNA-Seq assays, with three replicates for each group. For the m6A-Seq, we obtained approximately 4.10 million raw reads for each sample, and about 3.39 million clean reads were mapped to the bovine reference genome for each animal (Supplementary Table S2). Regarding RNA-Seq, approximately 4.20 million raw reads for each individual were generated, and about 4.06 million valid reads were mapped to the reference genome for each individual (Supplementary Table S2). The proportions of mapped reads ranged from 86.39 to 97.12%, correspondingly (Supplementary Table S2). These results demonstrated that the high-quality sequence data obtained are proper to be used in subsequent analysis.

Using the exomePeak2 analysis, a total of 12,947, 20,016, and 16,200 m6A peaks were detected in the TY0, TY1, and TY2 groups, respectively (Figure 1A). Among them, a total of 2,351, 4,259, and 1,701 specific peaks were observed in the TY0, TY1, and TY2, respectively, reflecting the significant difference among the studied groups in total m6A modification trends. Likewise, a total of 7,716, 10,170, and 9,063 genes were annotated in the TY0, TY1, and TY2 groups, respectively (Figure 1B). In addition, 8,774 peaks were consistently observed in three groups, and 6,206 genes (54.3% of total genes) within the studied groups were modified by m6A.

FIGURE 1

FIGURE 1. Transcriptome-wide m6A analysis in bovine testes. (A) The number of common and specific m6A peaks in three groups. (B) The Venn diagram shows the m6A-related transcripts in three groups. (C) Distribution of m6A peaks across chromosomes in the three groups. (D) The top motifs enriched from m6A peaks were identified among the studied groups.

Considering genome coverage pattern, the m6A distribution analysis revealed that the m6A peaks for each group differentially distributed on bovine chromosomes (Figure 1C). Most of the m6A peaks enriched in chromosome 19 in all studied groups. The motif analysis results indicated that the three studied groups had the classic m6A RRACH consensus sequences (Figure 1D). This enabled us to obtain the high credibility of the m6A peaks and revealed the presence of a prevailing methylated modification mechanism.

Analysis of m6A Modification Distribution in Testes Transcriptome

To investigate the preferential locations of m6A in transcripts, we explored the profiles of m6A peaks in the mRNA transcriptome by coordinating the bovine reference genome. The transcript was divided into the following regions: the 5′untranslated regions (5′UTRs), near the start codon, CDS, near the stop codon, and the 3′untranslated regions (3′UTRs). The analysis revealed that enrichment of m6A modified peaks was highest in the CDS followed by the 3′UTRs, 5′UTRs, stop codon region, and start codon region in the studied groups (Figure 2A). The distribution of m6A peak density in each group exhibited similar trends (Figure 2B). In addition, we observed that the number of m6A modified peaks enriched in the CDS region was higher in the both puberty and postpuberty groups than in the prepuberty group, but the opposite trend was observed for enrichment near the stop codon. Thereafter, the enrichment degree of m6A peaks significantly varied among all the pairwise groups (p-value < 2.21E-16; Figure 2C). The distribution of m6A modified peaks with each gene was explored and it showed that almost >50% of affected genes hold only one m6A peak and the majority of genes harbored one to three m6A peaks (Figure 2D).

FIGURE 2

FIGURE 2. Overview of m6A methylation profiles in bovine testes. (A) Pie charts demonstrating m6A peak distribution in the gene structures of mRNAs. (B) Metagene plots displaying the regions of m6A peaks identified across the transcripts in TY0, TY1, and TY2 groups. (C) Violin plot displays the distribution of enrichment degree of m6A peaks in each group. (D) The number of m6A peaks per gene in the TY0, TY1, and TY2 groups.

Differentially Methylated Gene Analysis

To explore the potential function of the m6A modification in bovine testes, a pairwise comparison was applied to scan the DMGs. Compared to the TY0 group, 2,495 significantly differential m6A peaks within 2,036 mRNAs were found in the TY1 group, while 2,047 differential peaks within 1,719 mRNAs were detected in the TY2 group (Figure 3A; Supplementary Table S3). In addition, we found 22 differential m6A peaks within 11 mRNAs in the TY2 group compared to the TY1 group. Several randomly selected mRNAs, including SLU7, GOLGA7, METTL14, and TGFB1, showed significantly hypermethylated peaks as presented in Figure 3B.

FIGURE 3

FIGURE 3. Distribution of significantly differential m6A peaks between the pairwise comparison groups. (A) Volcano plots showing the differential peaks between the studied groups. (B) Data visualization analysis of differential m6A peaks in the selected mRNAs (SLU7, GOLGA7, METTL14, and TGFB1) among the studied groups. (C) GO analysis for differentially methylated genes. (D) KEGG enrichment analysis for differentially methylated genes.

Moreover, the GO and KEGG enrichment analyses for all DMGs were performed to demonstrate the important function of m6A modification in bovine testes. GO analysis revealed that all the DMGs were mainly annotated into the nucleus, cytoplasm (ontology: cellular component), metal ion binding, RNA binding and ATP binding (ontology: molecular function), and protein ubiquitination, ubiquitin-dependent protein catabolic process (ontology: biological process) (Figure 3C). KEGG enrichment analysis demonstrated that the DMGs were significantly implicated with the MAPK signaling pathway, Herpes simplex virus 1 infection, Axon guidance, and FoxO signaling pathway, besides other methylated genes (Figure 3D).

RNA-Seq Identification of Differentially Expressed Genes

Using the RNA-Seq technique, a total of 8,280, 8,589, and 146 DEGs were detected between TY0 vs. TY1, TY0 vs. TY2, and TY1 vs. TY2, respectively (Figure 4A). Correspondingly, a set of 3,785, 4,001, and 77 up-regulated as well as 4,495, 4,588, and 69 down-regulated were found, respectively (Figure 4B). The hierarchical cluster expression pattern of the DEGs was shown in Figure 4C. Furthermore, all the DEGs were mainly enriched to 17 GO terms and 20 KEGG pathways (Figure 4D). Notably, most DEGs were highly annotated into the reproduction, reproductive process, multicellular organism reproduction, developmental process involved in reproduction, sexual reproduction, spermatogenesis, male gamete generation (ontology: biological process) and motile cilium (ontology: cellular component) and calcium ion binding (ontology: molecular function) of GO biological process. In addition, the PI3K-Akt signaling pathway is the most significant enrichment pathway for all DEGs identified, followed by MAPK signaling pathway and cAMP signaling pathway etc.

FIGURE 4

FIGURE 4. Differential expression analysis of RNA-Seq data among the studied groups. (A) Volcano plots showing the differential expressed genes between the studied groups. (B) Barplot showing the number of up-and down-regulated DEGs. (C) Heatmap plot of all DEGs among the studied groups. (D) GO and KEGG enrichment analysis of all DEGs.

Conjoint Analysis of m6A-Seq and RNA-Seq Data

To explore the potential relationship between m6A modification and gene expression, we performed a cross-analysis of the m6A-seq and RNA-seq data. As shown in Figure 5A, a positive correlation between differentially methylated peaks and gene expression levels is obvious (p = 0.0001, Spearman r = 0.1320). Based on the principle that absolute value of both X and Y axes was greater than 2, the results of the four-quadrant diagram analysis revealed that there were 502 differentially methylated genes, of which 124 genes belonged to the hypo-up, 2 genes to the hyper-up, 368 genes to the hypo-down, and 8 genes to the hyper-down quadrants (Figure 5B). Further, all of these genes were utilized for GO and KEGG pathway analyses (Figure 5C). Most of them were significantly enriched for 37 GO terms and 3 KEGG pathways. Interestingly, 84 genes were annotated into the nucleus of the GO term. The STRING analysis revealed that these genes were clustered into 8 PPIs networks (Figure 5D). Of them, 30 candidate genes were located in the largest PPIs network (Table 2). Further, 4 hub DMGs (PLK4, PTEN, EGR1, and PSME4) were selected and analyzed by RT-qPCR. The results of qPCR showed that the expression level of the 4 hub DMGs displayed a similar tendency with that of the RNA-Seq (Figures 5E,F). The above results suggested that these candidate genes may have crucial roles in spermatogenesis during bovine testis development.

FIGURE 5

FIGURE 5. Conjoint analysis of m6A-Seq and RNA-Seq data. (A) Dot plot of Log2 FC (mRNA expression) versus Log2 FC (differential m6A methylation) revealing a positive association between total m6A methylation and level of mRNA expression. (B) Four quadrant plots showing differentially expressed genes with differentially methylated m6A peaks. (C) GO and KEGG pathway enrichment analysis of the genes with a significant change in both m6A and mRNA levels. (D) Protein-protein interactions (PPIs) of the genes enriched in the nucleus of GO terms. The red circle represents the node with a high degree, while the blue circle indicates the node with a low degree. (E) The relative mRNA levels were determined by the qPCR of four hub genes in three groups. (F) The genes change levels based on RNA-Seq data.

TABLE 2

TABLE 2. List of 30 genes with significant changes in m6A and mRNA transcript abundance in bovine testes.

Discussion

Improvement of our knowledge regarding the testis functions developments is vital for gaining our understating regard to the spermatogenesis process. The role of m6A modification influences the mechanisms underlying testis development and spermatogenesis. To our best knowledge, this work is the first comprehensive high-throughput study of RNA methylation in the testes of XN young calves and mature bulls. To explore the potential function of the m6A modification affecting spermatogenesis, we specified three age points (TY0, TY1, and TY2) with significantly different physiological statuses in bovine testis development to analyze the transcriptome-wide m6A profile. The generated data in the present study displayed that a diverse pattern of mRNA methylation have occurred over the testis development. These mRNA m6A sites were mainly concentrated around the CDS and 3′UTRs region (95.28%), in agreement with the distributional characteristics of the mammalian transcriptomes (Wen et al., 2020). Our data showed that although the m6A peaks for each group markedly varied on different chromosomes, they had a highest expressed distribution in chromosome 19 within groups. The great differences in the genome coverage of reads among chromosomes observed in the present study may be attributed to the variation in the gene density of the bovine chromosomes. In addition, these m6A sites tended to occur in the conserved motif sequence “RRACH”, similar to that reported in other animals, such as the goat (Wang et al., 2020), sheep (Lu et al., 2019), and yak (Zhang et al., 2020a). This finding supported our hypothesis on the presence of a predominant methylated modification mechanism based on the m6A peaks identified.

In the study, a total of 2,278 unique DMGs were identified among all the pairwise groups. Differential m6A methylation has proved to be responsible for tissue or organ differentiation and development (Wen et al., 2020). Regarding the function and pathway of DMGs, we performed GO and KEGG enrichment analysis for DMGs. Importantly, GO analysis revealed that all the DMGs were mainly annotated into the nucleus, cytoplasm, metal ion binding, RNA binding, ATP binding, protein ubiquitination, and ubiquitin-dependent protein catabolic process. Earlier reports showed that metal ion binding was involved in porcine spermatogenesis (Luo et al., 2015) and horse testis development (Han et al., 2020). A growing body of evidence supports the fact that protein ubiquitination also has a critical role in the regulation of spermatogenesis in testes of the mouse (Guo et al., 2021) and buffalo (Huang et al., 2020; Zhang et al., 2021). Similarly, RNA binding was also essential for spermatogenesis (Venables and Eperon, 1999; Paronetto and Sette, 2010; Jan et al., 2017). Moreover, KEGG enrichment analysis revealed that the DMGs were significantly implicated in the MAPK signaling pathway, Herpes simplex virus 1 infection, Axon guidance, and FoxO signaling pathway. Ni et al. (2019) found that the MAPK signaling pathway regulates dynamics of tight junctions and adherents junctions, proliferation, and meiosis of germ cells, proliferation and lactate production of Sertoli cells. Accumulating evidence has demonstrated that the FoxO signaling pathway was involved in the regulation of spermatogenesis process (Huang et al., 2016b; Ge et al., 2019; Hu et al., 2021). Therefore, above GO terms and pathways results suggested that the DMGss identified were related to bovine testis development and spermatogenesis process.

Likewise, our results showed that there are 8,952 unique DEGs were observed among all the pairwise groups. We observed that most DEGs were highly annotated into the reproduction, reproductive process, multicellular organism reproduction, developmental process involved in reproduction, sexual reproduction, spermatogenesis, male gamete generation (ontology: biological process) of the GO biological process, indicating that the identified DEGs in the present study were reliable. Notably, the PI3K-Akt signaling pathway was found to be the most significant enrichment pathway for all DEGs. Many studies have demonstrated that the PI3K-Akt signaling pathway was involved in the regulation of spermatogenesis (Duan et al., 2016; Deng et al., 2020; Ni et al., 2020). In brief, our findings suggested that the DEGs were strongly associated with the testis development and spermatogenesis.

A surprising finding was that the enrichment degree of m6A peaks (Figure 2C), intensity of the DMGs (Figure 4A) and DEGs (Figure 5A) exhibited marked differences between before (TY0) and after puberty (TY1 and TY2 groups), suggesting that m6A modification and their gene expression mainly involved in the bovine testicular development and spermatogenesis. To further explore the potential genes underlying bovine testis development and spermatogenesis, the integrated analysis of m6A-Seq and mRNA-Seq data were performed. Results revealed that a positive correlation existed in differentially methylated peaks and gene expression levels, and 502 genes with concomitant changes in the mRNA expression and m6A modification. Interestingly, 84 genes were annotated into the nucleus of GO term, of which 30 candidate genes were located in the largest PPIs network. Importantly, four key node genes (PLK4, PTEN, EGR1, and PSME4) were observed on the largest PPIs network. Notably, the results of qPCR also supported that the expression level of the 4 selected DMGs exhibited marked differences between before and after puberty, which displayed a similar tendency with that of the RNA-Seq. These results indicated the putative role of these four genes in bovine testis development and spermatogenesis. It was noted that the four genes have already been reported to regulate the mammalian testis development and spermatogenesis. For example, the PLK4 gene was highly expressed in testes during both pre- and post-natal stages and had a role in the initiation of spermatogenesis (Harris et al., 2011). Moreover, Miyamoto et al. (2016) reported that a mutation in PLK4 causing azoospermia in a man with Sertoli cell-only syndrome. Dorostghoal et al. (2020) found that sperm miR-26a-5p and its target PTEN transcript content may contribute to the etiology of male infertility in unexplained infertile patients. More specifically, Neirijnck et al. (2019) reported that ablation of the PTEN appears dispensable for Sertoli cell proliferation and spermatogenesis, where inactivation of PTEN gene in the absence of Insr and Igf1r rescued the Sertoli cell proliferation rate during late fetal development, testis size, and sperm production. Besides, the EGR1 gene has a role in maintaining spermatogonia stem cells’ self-renewal and is a target to better understand the molecular basis of spermatogenesis (Sisakhtnezhad and Heshmati, 2018). Also, the EGR1 gene can act as an activator of the sex-determining region Y box 18 promoter (Petrovic et al., 2010). Additionally, the loss of PSME4 gene led to a marked reduction in male fertility, due to defects in spermatogenesis observed in meiotic spermatocytes and also during the maturation of postmeiotic haploid spermatids (Khor et al., 2006). Huang et al. (2016a) reported that the double knockout of PSME3 and PSME4 genes in mice resulted in completely infertile males. These findings suggested that m6A modifications play an essential role during bovine testis development and spermatogenesis.

Conclusion

We identified 8,774 m6A peaks and 6,206 m6A genes among the studied groups by using MeRIP-Seq analysis. The cross-analysis of m6A and mRNA expression revealed 502 genes with concomitant changes in the mRNA expression and m6A modification. Further, 30 candidate genes were found to be located in the largest network of protein-protein interactions. Notably, four key node genes (PLK4, PTEN, EGR1, and PSME4) have been implicated in the regulation of mammalian testis development and spermatogenesis. The present study is pioneer to present a map of RNA m6A modification in bovine testes at different ages, and provides novel insights into m6A topology and associated molecular mechanisms underlying bovine spermatogenesis, and provides a basis for future research on mammalian testis development and spermatogenesis.

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: NCBI SRA; PRJNA776655.

Ethics Statement

The animal study was reviewed and approved by the Henan Agricultural University Animal Care and Use Committee. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

S-HL, T-TY, Z-CW and Y-YG collected the phenotype data and testes tissue samples. S-HL, X-YM and T-TY isolated RNA. T-XD, S-HL, H-ER and X-YM, created and carried out the analysis, interpreted the data. S-HL, H-ER and T-XD wrote the manuscript. T-XD, S-HL, and X-LH developed the study and participated in its design and coordination. J-CL, K-LQ, TF, ML, FL, T-YG and L-GY reviewed the paper. All authors read and approved the manuscript.

Funding

This work was supported by the China Agriculture Research System of MOF and MARA and China Postdoctoral Science Foundation (2020M672233) and Henan Key Research and Development Program (192102110067).

Conflict of Interest

T-tY was employed by Henan Dairy Herd Improvement Co., Ltd., Y-yG was employed by Henan Dingyuan Cattle Breeding Co., Ltd.

The remaining 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

The authors thank Kerchin Cattle Industry (Nanyang) Co., Ltd. (China) for providing experimental animals and biological samples.

Supplementary Material

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

References

Bailey, T. L., Boden, M., Buske, F. A., Frith, M., Grant, C. E., Clementi, L., et al. (2009). MEME Suite: Tools for Motif Discovery and Searching. Nucleic Acids Res. 37 (Web Server issue), W202–W208. doi:10.1093/nar/gkp335

PubMed Abstract | CrossRef Full Text | Google Scholar

Bane, A. (1952). A Study on the Technique of Hemocytometric Determination of Sperm Motility and Sperm Concentration in Bull Semen. Cornell Vet. 42 (4), 518–531. doi:10.1016/0035-9203(51)90016-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Björndahl, L., Söderlund, I., and Kvist, U. (2003). Evaluation of the One-step Eosin-Nigrosin Staining Technique for Human Sperm Vitality Assessment. Hum. Reprod. (4), 813–816. doi:10.1093/humrep/deg199

CrossRef Full Text | Google Scholar

Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). Fastp: an Ultra-fast All-In-One FASTQ Preprocessor. Bioinformatics 34 (17), i884–i890. doi:10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, C.-Y., Lv, M., Luo, B.-H., Zhao, S.-Z., Mo, Z.-C., and Xie, Y.-J. (2020). The Role of the PI3K/AKT/mTOR Signalling Pathway in Male Reproduction. Curr. Mol. Med. 21 (7), 539–548. doi:10.2174/1566524020666201203164910

CrossRef Full Text | Google Scholar

Ding, H., Liu, M., Zhou, C., You, X., Su, T., Yang, Y., et al. (2020). Integrated Analysis of miRNA and mRNA Expression Profiles in Testes of Duroc and Meishan Boars. BMC Genomics 21 (1), 686. doi:10.1186/s12864-020-07096-7

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

Dorostghoal, M., Galehdari, H., Moramezi, F., and Danyari, R. (2020). Sperm miR‐26a‐5p and its Target PTEN Transcripts Content in Men with Unexplained Infertility. Andrologia 8 (5), 1167–1173. doi:10.1111/andr.12801

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, P., Quan, C., Huang, W. T., and Yang, K. D. (2016). PI3K-Akt/LKB1-AMPK-mTOR-p70S6K/4EBP1 Signaling Pathways Participate in the Regulation of Testis Development and Spermatogenesis: An Update. Zhonghua Nan Ke Xue 22 (11), 1016–1020. doi:10.13263/j.cnki.nja.2016.11.011

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

Ge, P., Zhang, J., Zhou, L., Lv, M. Q., Li, Y. X., Wang, J., et al. (2019). CircRNA Expression Profile and Functional Analysis in Testicular Tissue of Patients with Non-obstructive Azoospermia. Reprod. Biol. Endocrinol. 17 (1), 100. doi:10.1186/s12958-019-0541-4

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

Guo, Y., Zhang, H., Yao, L., Li, Y., Situ, C., Sha, J., et al. (2021). Systematic Analysis of the Ubiquitome in Mouse Testis. Proteomics.

Google Scholar

Han, H., Chen, Q., Gao, Y., Li, J., Li, W., Dang, R., et al. (2020). Comparative Transcriptomics Analysis of Testicular miRNA from Cryptorchid and Normal Horses. Animals (Basel) 10 (2). doi:10.3390/ani10020338

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, R. M., Weiss, J., and Jameson, J. L. (2011). Male Hypogonadism and Germ Cell Loss Caused by a Mutation in Polo-Like Kinase 4. Endocrinology 152 (10), 3975–3985. doi:10.1210/en.2011-1106

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. 2010. Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol. Cel 38(4):576–589.doi:10.1016/j.molcel.2010.05.004

CrossRef Full Text | Google Scholar

Hsu, P. J., Zhu, Y., Ma, H., Guo, Y., Shi, X., Liu, Y., et al. (2017). Ythdc2 Is an N6-Methyladenosine Binding Protein that Regulates Mammalian Spermatogenesis. Cell Res 27 (9), 1115–1127. doi:10.1038/cr.2017.99

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, T., Luo, S., Xi, Y., Tu, X., Yang, X., Zhang, H., et al. (2021). Integrative Bioinformatics Approaches for Identifying Potential Biomarkers and Pathways Involved in Non-obstructive Azoospermia. Transl Androl. Urol. 10 (1), 243–257. doi:10.21037/tau-20-1029

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, L., Haratake, K., Miyahara, H., and Chiba, T. (2016a). Proteasome Activators, PA28γ and PA200, Play Indispensable Roles in Male Fertility. Sci. Rep. 6, 23171. doi:10.1038/srep23171

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, P., Zhou, Z., Shi, F., Shao, G., Wang, R., Wang, J., et al. (2016b). Effects of the IGF-1/PTEN/Akt/FoxO Signaling Pathway on Male Reproduction in Rats Subjected to Water Immersion and Restraint Stress. Mol. Med. Rep. 14 (6), 5116–5124. doi:10.3892/mmr.2016.5880

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y. L., Zhang, P. F., Hou, Z., Fu, Q., Li, M. X., Huang, D. L., et al. (2020). Ubiquitome Analysis Reveals the Involvement of Lysine Ubiquitination in the Spermatogenesis Process of Adult buffalo (Bubalus Bubalis) Testis. Biosci. Rep. 40. doi:10.1042/BSR20193537

PubMed Abstract | CrossRef Full Text | Google Scholar

Jan, S. Z., Vormer, T. L., Jongejan, A., Röling, M. D., Silber, S. J., de Rooij, D. G., et al. (2017). Unraveling Transcriptome Dynamics in Human Spermatogenesis. Development 144 (20), 3659–3673. doi:10.1242/dev.152413

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, G., Fu, Y., Zhao, X., Dai, Q., Zheng, G., Yang, Y., et al. (2011). N6-Methyladenosine in Nuclear RNA Is a Major Substrate of the Obesity-Associated FTO. Nat. Chem. Biol. 7 (12), 885–887. doi:10.1038/nchembio.687

PubMed Abstract | CrossRef Full Text | Google Scholar

Khor, B., Bredemeyer, A. L., Huang, C.-Y., Turnbull, I. R., Evans, R., Maggi, L. B., et al. (2006). Proteasome Activator PA200 Is Required for normal Spermatogenesis. Mol. Cel. Biol. 26 (8), 2999–3007. doi:10.1128/mcb.26.8.2999-3007.2006

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

Krueger, F. (2021). Trim Galore! in. Available at: https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/.

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. Cel Res 27 (10), 1216–1230. doi:10.1038/cr.2017.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Yue, Y., Han, D., Wang, X., Fu, Y., Zhang, L., et al. (2014). A METTL3-METTL14 Complex Mediates Mammalian Nuclear RNA N6-Adenosine Methylation. Nat. Chem. Biol. 10 (2), 93–95. doi:10.1038/nchembio.1432

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 25 (4), 402–408. doi:10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Z., Ma, Y., Li, Q., Liu, E., Jin, M., Zhang, L., et al. (2019). The Role of N6-Methyladenosine RNA Methylation in the Heat Stress Response of Sheep (Ovis aries). Cell Stress and Chaperones 24 (2), 333–342. doi:10.1007/s12192-018-00965-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Z., Liu, Y., Chen, L., Ellis, M., Li, M., Wang, J., et al. (2015). microRNA Profiling in Three Main Stages during Porcine Spermatogenesis. J. Assist. Reprod. Genet. 32 (3), 451–460. doi:10.1007/s10815-014-0406-x

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

Miyamoto, T., Bando, Y., Koh, E., Tsujimura, A., Miyagawa, Y., Iijima, M., et al. (2016). APLK4mutation Causing Azoospermia in a Man with Sertoli Cell-Only Syndrome. Andrology 4 (1), 75–81. doi:10.1111/andr.12113

PubMed Abstract | CrossRef Full Text | Google Scholar

Neirijnck, Y., Kühne, F., Mayère, C., Pavlova, E., Sararols, P., Foti, M., et al. (2019). Tumor Suppressor PTEN Regulates Negatively Sertoli Cell Proliferation, Testis Size, and Sperm Production In Vivo. Endocrinology 160 (2), 387–398. doi:10.1210/en.2018-00892

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, F.-D., Hao, S.-L., and Yang, W.-X. 2020. Molecular Insights into Hormone Regulation via Signaling Pathways in Sertoli Cells: With Discussion on Infertility and Testicular Tumor. Gene 753: 144812. doi:10.1016/j.gene.2020.144812

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, F. D., Hao, S. L., and Yang, W. X. (2019). Multiple Signaling Pathways in Sertoli Cells: Recent Findings in Spermatogenesis. Cell Death Dis 10, 541. doi:10.1038/s41419-019-1782-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, Y., Zhao, X., Wu, Y.-S., Li, M.-M., Wang, X.-J., and Yang, Y.-G. (2013). N6-methyl-adenosine (m6A) in RNA: An Old Modification with A Novel Epigenetic Function. Genomics, Proteomics & Bioinformatics 11 (1), 8–17. doi:10.1016/j.gpb.2012.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Paronetto, M. P., and Sette, C. (2010). Role of RNA-Binding Proteins in Mammalian Spermatogenesis. Int. J. Androl. 33 (1), 2–12. doi:10.1111/j.1365-2605.2009.00959.x

CrossRef Full Text | Google Scholar

Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level Expression Analysis of RNA-Seq Experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11 (9), 1650–1667. doi:10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Petrovic, I., Kovacevic-Grujicic, N., and Stevanovic, M. (2010). Early Growth Response Protein 1 Acts as an Activator of SOX18 Promoter. Exp. Mol. Med. 42 (2), 132–142. doi:10.3858/emm.2010.42.2.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Puglisi, R., Gaspa, G., Balduzzi, D., Severgnini, A., Vanni, R., Macciotta, N., et al. (2016). Genomewide Analysis of Bull Sperm Quality and Fertility Traits. Reprod. Dom Anim. 51 (5), 840–843. doi:10.1111/rda.12747

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, S., Mumbach, M. R., Jovanovic, M., Wang, T., Maciag, K., Bushkin, G. G., et al. (2014). Perturbation of m6A Writers Reveals Two Distinct Classes of mRNA Methylation at Internal and 5′ Sites. Cel Rep. 8 (1), 284–296. doi:10.1016/j.celrep.2014.05.048

CrossRef Full Text | Google Scholar

Shi, H., Wei, J., and He, C. (2019). Where, when, and How: Context-dependent Functions of RNA Methylation Writers, Readers, and Erasers. Mol. Cel 74 (4), 640–650. doi:10.1016/j.molcel.2019.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Sisakhtnezhad, S., and Heshmati, P. (2018). Comparative Analysis of Single‐cell RNA Sequencing Data from Mouse Spermatogonial and Mesenchymal Stem Cells to Identify Differentially Expressed Genes and Transcriptional Regulators of Germline Cells. J. Cel Physiol 233 (7), 5231–5242. doi:10.1002/jcp.26303

CrossRef Full Text | Google Scholar

Venables, J., and Eperon, I. (1999). The Roles of RNA-Binding Proteins in Spermatogenesis and Male Infertility. Curr. Opin. Genet. Dev. 9 (3), 346–354. doi:10.1016/s0959-437x(99)80052-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Vendelbo, N. M., Mahmood, K., Sarup, P., Kristensen, P. S., Orabi, J., and Jahoor, A. (2021). Genomic Scan of Male Fertility Restoration Genes in a 'Gülzow' Type Hybrid Breeding System of Rye (Secale Cereale L.). Int. J. Mol. Sci. 22 (17). doi:10.3390/ijms22179277

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, Y., Zheng, Y., Guo, D., Zhang, X., Guo, S., Hui, T., et al. (2020). m6A Methylation Analysis of Differentially Expressed Genes in Skin Tissues of Coarse and Fine Type Liaoning Cashmere Goats. Front. Genet. 10, 1318. doi:10.3389/fgene.2019.01318

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, K., Zhang, Y., Li, Y., Wang, Q., and Sun, J. (2020b). Comprehensive Analysis of Transcriptome-wide m6A Methylome in the Anterior Capsule of the Lens of High Myopia Patients. Epigenetics, 1–14. doi:10.1080/15592294.2020.1834917

CrossRef Full Text | Google Scholar

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer.

Google Scholar

Xiong, X., Hou, L., Park, Y. P., Molinie, B., Ardlie, K. G., Aguet, F., et al. (2021). Genetic Drivers of m6A Methylation in Human Brain, Lung, Heart and Muscle. Nat. Genet. 53 (8), 1156–1165. doi:10.1038/s41588-021-00890-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., and He, Q. Y. (2015). ChIPseeker: an R/Bioconductor Package for ChIP Peak Annotation, Comparison and Visualization. Bioinformatics 31 (14), 2382–2383. doi:10.1093/bioinformatics/btv145

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P. f., Huang, Y. l., Fu, Q., He, W. t., Xiao, K., and Zhang, M. (2021). Integrated Analysis of Phosphoproteome and Ubiquitylome in Epididymal Sperm of buffalo ( Bubalus Bubalis ). Mol. Reprod. Dev. 88 (1), 15–33. doi:10.1002/mrd.23432

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Liang, C., Wu, X., Pei, J., and Yan, P. (2020a). Integrated Study of Transcriptome-wide m6A Methylome Confirms a Potential Mechanism for Transcriptional Regulation during Yak Adipocytes Differentiation.

Google Scholar

Zhang, Z., Luo, K., Zou, Z., Qiu, M., Tian, J., Sieh, L., et al. (2020b). Genetic Analyses Support the Contribution of mRNA N6-Methyladenosine (m6A) Modification to Human Disease Heritability. Nat. Genet. 52, 939–949. doi:10.1038/s41588-020-0644-z

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. Cel 49 (1), 18–29. doi:10.1016/j.molcel.2012.10.015

CrossRef Full Text | Google Scholar