- Email: [email protected]
Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids
GENE-41615; No. of pages: 9; 4C: Gene xxx (2016) xxx–xxx
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Research paper
Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids Hui Zhang a, Chuanchao Zhang a, Jinting Yan a, Zhongshuai Sun a, Shuhui Song b, Yazhou Sun a, Caixia Guo a,⁎ a b
CAS Key Laboratory of Genomic and Precision Medicine, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China Core Genomic Facility, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China
a r t i c l e
i n f o
Article history: Received 27 July 2016 Received in revised form 12 September 2016 Accepted 3 October 2016 Available online xxxx Keywords: Transcriptome RNA-seq Spermatogenesis MMS DNA damage response Differentially expressed genes
a b s t r a c t Spermatogenesis is threatened by DNA alkylating agents, one major category of DNA damaging agents. Currently, little is known about the alterations in transcriptome profiling of the mouse spermatogenic cells in response to DNA alkylation at distinct stages of spermatogenesis. In this study, RNA sequencing (RNA-seq) was performed in pachytene spermatocytes (PS) and round spermatids (RS) at 0 or 30 min following Methyl Methanesulfonate (MMS) treatment and with untreated controls. A large number of differentially expressed genes (DEGs) were identified by comparison of the three groups in PS and RS, respectively. Functional analyses of all DEGs highlighted the protein ubiquitination pathway and DNA damage response (DDR) network being the two main biological processes in common in the two cell types. Further analyses of the DEGs with 2-fold or more changes between 30 min repair and control group indicated that several cytokine signaling pathways were the most strongly affected in PS and DDR related pathways in RS, respectively. Gene ontology (GO) analyses directly showed differential biological process (BP) affected between PS and RS, with “regulation of transcription” being most overrepresented in PS and “cellular response to stress” in RS, respectively. Moreover, 374 DDR-related genes in PS and 158 in RS among all DEGs were filtered and clustered, which showed dynamic expression patterns in PS and RS. Our analyses provide a transcriptional landscape for male germ cells in response to MMS during spermatogenesis. © 2016 Elsevier B.V. All rights reserved.
1. Introduction In adult mammals, spermatogenesis is a continuous and well-orchestrated process, which can be divided into three principal stages: a mitotic stage in spermatogonia, a meiotic stage in spermatocytes and a post-meiotic stage in spermatids. During the lengthy meiotic prophase, the primary spermatocyte progresses through leptotene, zygotene, pachytene and diplotene stages to generate two secondary spermatocytes and then is rapidly converted into four haploid round spermatids. Subsequently, during the post-meiotic stage, round spermatids differentiate into elongated spermatid and mature spermatozoa. Maintenance of genomic integrity during male germ cell development is crucial for production of healthy offspring, which requires Abbreviations: RNA-seq, RNA sequencing; PS, pachytene spermatocytes; RS, round spermatids; MMS, methyl methanesulfonate; DEGs, differentially expressed genes; DDR, DNA damage response; GO, Gene ontology; BP, biological process; BER, base excision repair; NER, nucleotide excision repair; TM, Tail Moment; RPKM, Reads Per Kilobase of exon model per Million mapped reads; DAVID, Database for Annotation, Visualization and Integrated Discovery; IPA, Ingenuity Pathway Analysis. ⁎ Corresponding author. E-mail address: [email protected] (C. Guo).
efficient surveillance and repair mechanisms for the fidelity of genetic material transmissions. Spermatogenesis is constantly challenged by various endogenous and exogenous factors which induce DNA damage; male germ cells combat DNA damage through DDR involving a range of DNA repair mechanisms and alteration of gene expression (Ciccia and Elledge, 2010). Alkylating agents represent an important class of DNA damaging agents. Among them, S-adenosylmethionine (SAM) is a ubiquitous methyl donor, constantly contributing to endogenous DNA alkylation. Tobacco-specific N-nitrosamines and cancer chemotherapeutics such as cyclophosphamide, melphalan, and ifosfamide are familiar exogenous alkylating agents which probably threaten spermatogenesis of human beings (Meistrich, 2013). In the laboratory, methyl methanesulfonate (MMS) has served as one of the most extensively used alkylating agents in DNA repair studies. MMS predominantly produce N-methylation adducts including N7-methylguanine (7meG) and N3-methyladenine (3meA). In somatic cells, adducts induced by MMS are mainly repaired by direct DNA repair, by base excision repair (BER), and by nucleotide excision repair (NER) (Fu et al., 2012). The effects of MMS on spermatogenesis have been initially detected in several models animals including mouse, rat and drosophila. Due to the unavailability of efficient scientific tools, most of the studies have
http://dx.doi.org/10.1016/j.gene.2016.10.006 0378-1119/© 2016 Elsevier B.V. All rights reserved.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
2
H. Zhang et al. / Gene xxx (2016) xxx–xxx
mainly focused on the mutagenic effects of MMS on male germ cells, and emphasized expression changes of a small number of genes. Studies by Robaire et al. have described the alterations of stress responsive genes to cyclophosphamide using microarray in rat male germ cells (Aguilar-Mahecha et al., 2001). Several other large scale studies have elucidated the stage specificity and complexity of testicular transcription machinery with microarray experiments or RNA-seq during normal spermatogenesis (Johnston et al., 2008). However, the complete landscape of the gene expression patterns and DDR molecular events following MMS treatment in specific spermatogenic cells remains poorly understood. In this study, comprehensive transcriptome analyses in response to MMS treatment in PS and RS were conducted using RNA-seq and functional annotation of DEGs and highlight important biological processes, including DDR related pathways. 2. Materials and methods 2.1. Isolation of spermatogenic cells We used male C57BL/6Jx129 mice for the preparations of different testis cell populations. All animal experiments were reviewed and approved by the institute of Zoology, Institutional Animal Care and Use Committee and were conducted according to the committee's guidelines. The spermatogenic cell populations were isolated from the testes of 3-month-old mice using a StaPut 2–4% bovine serum albumin gradient as described previously (Bellve et al., 1977; Romrell et al., 1976). The PS and RS fractions were identified both by their size as well as by typical cellular morphological by phase contrast microscopy (Leica). The purity of PS exceeded 80% and that of RS exceeded 90% (Supplemental Fig. S1). Cell viability was N 95% as determined with the method of trypan blue exclusion. 2.2. Cell culture and MMS treatment Primary PS and RS cells (4 × 106 cells per sample) were treated with 0.12 mM MMS or left untreated as control for 1 h at 32 °C in DMEM/F12 culture medium. The cells were then placed on ice for reaction termination and washed two times with medium. For repair incubation, cells were maintained in DMEM/F12 medium supplemented with 10% FCS, pyruvate (0.1 mg/ml) and antibiotics (100 IU/ml penicillin; 0.1 mg/ml streptomycin) in a humidified atmosphere at 32 °C, 5% CO2 for up to 2 h. The cells were collected at different time points for comet assay or RNA extraction. 2.3. Comet assay Alkaline comet assays were performed according to the manufacturer's protocol (Trevigen). Samples were stained with SYBR Gold and images were taken by a fluorescent microscope (Leica) following alkaline electrophoresis. The average of Tail Moment (TM) was established by measuring at least 100 cells per slide using Comet Assay Software Project Casp-1.2.2 (University of Wroclaw). The reported TMs represented mean ± standard error (SE). The experiments were repeated three times.
software allowing up to four mismatches. We counted the number of reads mapped to exons, introns, splice junctions, intron-exon adjacent regions of the annotated genes, and intergenic regions using in-house built pipeline wapRNA (Zhao et al., 2011). Gene expression levels were quantified by RPKM (Reads Per Kilobase of exon model per Million mapped reads), and then DEGs were defined using the R-package DEGseq with the method MARS (MA-plot-based method with Random Sampling) model (Wang et al., 2010), and with the BenjaminiHochberg false-discovery rate (FDR) adjustment for multiple testing. Hierarchical clustering of all expressed genes was conducted using R software. 2.5. Function analysis of differentially expressed genes GO analyses of DEGs were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 (http://david. abcc.ncifcrf.gov/). GO BP was selected as the functional annotation category for this analysis. Enrichment map of DEGs was generated by Cytoscape 2.8.2 software installed with Enrichment Map plugin. Each node represents a pathway, and the size of the node reflects the number of genes in the pathway. The canonical pathways and interaction network were determined using Ingenuity Pathway Analysis (IPA) v9.0-3211 (Ingenuity Systems, Inc.) Ensemble accession was used as the identifier and Ingenuity knowledge gene database was used as a reference for the pathway and network analysis. 2.6. Quantitative real-time RT-PCR (qRT-PCR) validation qRT-PCR was performed to validate a subset gene expression changes predicted by RNA-seq, primers used in qRT-PCR were listed in Supplemental Table 1. One microgram of total RNA from each sample was reverse-transcribed in 20 μl of reaction medium using the First Strand cDNA Synthesis kit (Promega) following the manufacturer's instructions. PCR was detected using SYBR Green Real time PCR Master Mix (Toyobo) on a Real-time thermal cycler CFX96 (Bio-Rad). The target gene expression levels were calculated relative to the expression of mouse Gapdh, employing an optimized comparative Ct (ΔΔCt) value method. The experiments were performed in triplicate. Data are expressed as mean ± standard deviation (SD). 3. Results 3.1. Cellular repair of MMS-induced DNA lesions To investigate the cellular response to DNA damage reagent MMS treatment in spermatogenic cells of adult mouse, the purified PS and RS (Supplemental Fig. 1) were exposed to 0.12 mM MMS at 32 °C for 1 h and then were harvested at different time points. Efficient and
2.4. RNA sequencing Total RNA was extracted from MMS-treated PS and RS cells at 0 min and 30 min after repair incubation and the control group using Trizol reagent (Invitrogen), respectively. The total six samples were indicated as PScon, PS0, PS30, RScon, RS0 and RS30. Sequencing libraries were prepared according to the manufacturer's instructions (Illumina). After trimming low-quality reads, the high-quality reads were aligned against the reference genome of Mus musculus (NCBIM 37; ENSEMBL release 62; http://asia.ensembl.org/info/data/ftp/index.html) using BWA 0.5.8
Fig. 1. Repair of MMS-induced DNA lesions in PS and RS. PS and RS were exposed to 0.12 mM MMS for 1 h, respectively, and incubated for up to 2 h to allow repair. DNA damage (tail moment; mean ± SE) was measured with alkaline comet assay at indicated time points. At least 100 cells were scored for each sample.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
3
Table 1 Summary of sequence data and read alignment statistics for individual samples. Sample
Total reads
Uniquely mapped
Counted on exons
Counted on introns
Expressed genes
PScon PS0 PS30 RScon RS0 RS30
45,889,931 45,698,664 46,550,854 44,207,333 44,846,475 45,911,712
36,382,567 (79.28%) 34,855,055 (76.28%) 36,749,237 (78.94%) 35,206,820 (79.64%) 34,839,289 (77.69%) 34,994,342 (76.22%)
29,928,689 29,300,496 30,079,158 30,241,642 29,304,027 29,965,255
5,431,178 4,905,678 5,687,934 4,328,639 4,816,948 4,639,136
21,430 21,411 21,500 20,994 21,214 21,280
rapid repair of damaged DNA was observed by both PS and RS using the comet assay, and the repair kinetics were similar between the two cell types (Fig. 1). 3.2. Gene expression profiling of PS and RS in response to MMS treatment To gain a complete understanding of gene expression profiles of PS and RS with MMS-treatment, RNA-seq was employed. In each cell type, 3 samples were selected and designated as the control, 0 or 30 min group following treatment with MMS, representing the untreated, damaged and repaired status, respectively (a total of 6 samples were labeled as PScon, PS0, PS30; RScon, RS0, RS30). The total number of reads generated from each sample varied from 44,207,333 to 46,550,854, which were aligned to the reference genome in a unique manner, with mapping rates between 76.22% and 79.61% (Table 1). In total, 27,930 genes were expressed at least in one stage. Hierarchical clustering of these 27,930 genes showed a clear difference in gene expression pattern of sample groups from PS and RS and also revealed a farther correlation between 30 min repair and the control group in both PS and RS (Fig. 2A). DEG analysis showed the amounts of genes changed significantly (P b 0.001), including 3778 for PS0/PScon, 7680 for PS30/PS0 and 7482 for PS30/PScon in PS (Fig. 2B); 5897 for RS0/RScon, 3012 for RS30/RS0 and 4596 for RS30/RScon in RS (Fig. 2B). The total number of overlapping DEGs in PS and RS was 9705 and 7071, respectively
(Fig. 2C and D). Moreover, the number of DEGs with 2-fold or greater change and the range of fold change (FC) in each comparison of PS and RS show the highest difference in the number of DEGs between the 30 min repair and control groups in both PS (868) and RS (295) cells, consistent with the results of hierarchical clustering (Table 2). Special attention was thus paid to the comparisons between the 30 min repair and control groups in both PS and RS in the following analyses. The lists of top 30 up- and down-regulated DEGs comparing PS or RS are presented (Supplemental Tables 2 and 3, respectively).
3.3. Functional analysis of differentially expressed genes Canonical pathways and molecular networks related to DEGs were identified using the “Core Analysis” function in IPA for each comparison. The protein ubiquitination pathway is among the top canonical pathways in all comparisons of both PS and RS (Fig. 3). Other top canonical pathways affected are oxidative phosphorylation, mitochondrial dysfunction and EIF2 signaling in PS cells (Fig. 3A) and AMPK signaling, mitotic roles of polo-like kinase, cyclins and cell cycle regulation in RS cells, respectively (Fig. 3B). Network analyses revealed several interactome networks in common for PS and RS cells, such as developmental and hereditary disorder, RNA post-transcriptional modification, cancer, and processes related to the DDR (i.e., DNA replication, recombination and repair, cell cycle) (Supplemental Tables 4 and 5).
Fig. 2. Dynamic changes in gene expression and overlapping DEGs. (A) Hierarchical clustering for whole expressed genes in PS and RS sample groups. The default colors in R were used. Upregulated genes are indicated in white and downregulated in red. (B) Bar graph shows total number of DEGs and up- and down regulated gene numbers among comparisons in PS and RS. (C) and (D) The number of overlapping DEGs between comparisons in PS and RS, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
4
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Table 2 The number of DEGs with 2-fold change or greater. Comparisons
Total
Up regulated
Down regulated
Fold change range
PS0/PScon PS30/PS0 PS30/PScon RS0/RScon RS30/RS0 RS30/RScon
139 692 868 293 60 295
132 282 643 85 20 118
7 410 225 208 40 177
14.45 to −26.16 48.90 to −45.00 26.24 to −25.76 13.79 to −31.33 6.67 to −25.33 71.64 to −26.13
DEGs with 2-fold or greater changes in expression between the 30 min repair and control groups were further subjected to the same analyses using IPA. The most significant top 10 canonical pathways of PS and RS cells are presented in Fig. 4A and B. Notably, inflammatory responses including IL-17A, IL-17F, IL-8 and TNFR2 signaling are among the top canonical pathways in PS while cell cycle control and checkpoint regulation are among that in RS. The top one network in PS was embryonic development, organismal development, and developmental disorder, hereditary disorder, and immunological disease in RS, respectively
(Fig. 4C and D). Meanwhile, GO analyses based on the DAVID tool revealed that a large proportion of these genes in PS was enriched in BP, such as regulation of transcription, regulation of RNA metabolic process and embryonic organ development (Fig. 4E). In RS, the enriched genes belong to several different functional categories, including regulation of transcription, cellular response to stress, chromosome organization and immune effecter processes (Fig. 4F). 3.4. Comparisons of spermatocytes and spermatids We next focused on the difference in DEGs in response to MMS treatment that distinguishes PS from RS. Of total DEGs between the 30 min repair and control groups, 59% (4407/7482) was unique to PS and 33% (1521/4596) to RS (Fig. 5A). When considering those with 2fold or greater changes, the proportions specific to PS and RS were sharply increased to 94% (814/868) and 82% (241/295), respectively (Fig. 5B). Further GO BP enrichment analysis of 814 genes unique to PS revealed that “regulation of transcription (such as C/EBPbeta, EP300, Junb and Pou3f3)”, “regulation of RNA metabolic process (such as Atf7
Fig. 3. Top canonical pathways associated with DEGs through comparisons among the three groups by IPA. (A) DEGs from PS. (B) DEGs from RS.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
5
Fig. 4. Functional analysis of the expressed gene with 2-fold change or more between the 30 min repair and control groups in PS and RS. (A) and (B) The top 10 canonical pathways in which genes were enriched by IPA. (C) and (D) Interaction networks of genes identified by IPA. The genes in green were down-regulated and the genes in red were up-regulated with respect to the control group. (E) and (F) Function analysis of genes based on the DAVID GO BP analysis result. The cutoff parameters for enrichment analysis with Cytoscape software are: P b 0.005, FDR q b 0.05, overlap cutoff N 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
and Fos)”, and “embryonic development (such as Nkx2-6 and Lmo4)” were overrepresented (Fig. 5C red strip). Among the 241 DEGs specific to RS, the processes “cellular response to stress (such as Dmc1, Wrnip1 and Udg)”, “regulation of cell cycle (such as Ccnd1and Cdt1)” and “DNA integrity checkpoint (such as Msh2 and Mapk14)” were enriched (Fig. 5C blue strip). 54 DEGs common in PS and RS were enriched in “regulation of transcription (such as Abca2 and Atf3)”, “regulation of RNA metabolic process (such as Fosb and Bhlhe40)” and “response to protein stimulus (such as Cyr61 and Egr1)” (Fig. 5C grey strip). These results suggest that distinct biological processes and gene expression patterns are required for dealing with MMS-induced DNA damage repair during and after meiosis.
3.5. DNA damage and repair gene expression patterns in spermatocytes and spermatids following MMS treatment Top networks of all DEGs involve processes related to DDR (Supplemental Tables 4 and 5), we sought to gain as complete a picture as possible of the dynamic regulation of DNA damage and repair processes in PS and RS in the following analyses. We briefly summarized the genes involved in DDR events by integrating GO annotation with a literature search, in which 1251 genes were annotated by GO:0007049 (cellcycle control); 573 were annotated by GO:0006974 (cellular response to DNA damage stimulus); 550 cell cycle genes were expressed during spermatogenesis from the literature (Roy Choudhury et al., 2010) and
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
6
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Fig. 5. Comparisons of PS and RS in gene expressions and functions. (A) and (B) Venn diagram showing the distribution of all and 2-fold or greater change DEGs between 30 min and control groups in PS and RS, respectively. The numbers of DEGs in PS are indicated in the red circle, while the numbers of DEGs in RS are highlighted in the blue circle. (C) GO BP categories enriched in DEGs from Venn analysis in 3B. Red represents DEGs unique to PS while blue to RS, and grey representing common DEGs in PS and RS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
174 mouse homologs of human DNA repair genes from the Wood laboratory website (http://sciencepark.mdanderson.org/labs/wood/dna repair genes.html) (Fig. 6A). Based on these criteria, 1021 and 736 DDR related genes were differentially expressed in PS and RS following MMS treatment, respectively. These DDR related genes were then sorted into two groups: 1) genes that are “similarly expressed” in PS or RS (no change of N = 1.5-fold) and 2) genes that were “differentially expressed” in PS or RS (one or more changes of N = 1.5-fold) (Fig. 6A). We chose the fold change of 1.5 fold as the criterion in this analysis because we wanted to identify the alteration of DDR related genes in spermatogenic cells following MMS treatment in a more comprehensive manner, which was also used in a previous study (Roy Choudhury et al., 2010). Our results show that, 374 and 158 genes were sorted as “differentially expressed” in PS and RS, respectively (Fig. 6A). We further clustered these “differentially expressed” genes in PS and RS with K-means and identified six distinct expression patterns of DDRrelated genes in PS or RS, including three patterns showing gene expression changes in one direction (cluster 1, 4, 6 in PS; cluster 2, 4, 5 in RS) and the other three patterns in two opposite directions (cluster 2, 3, 5 in PS; cluster 1, 3, 6 in RS) (Fig. 6B and Supplemental Tables 6 and 7), which aids our understanding the activation or suppression of DDR pathways. Several noteworthy observations from these analyses follow. In PS, although clustered in distinct patterns, several important genes
required for DNA double-strand breaks (DSBs) by homologous recombination (HR) all showed significantly decreased expression when comparing the 30 min repair and control groups, including Rad50, Exo1, and Rev1 in cluster 1 and Brca2, Dna2, Atm in cluster 5. DNA mismatch repair gene Mlh3 was also down-regulated as shown in cluster1. As expected, Aptx in cluster 2 and Apex2 in cluster 6, both of which participate in BER, are up-regulated. The gene of highest alteration is Irf1 in cluster 4, a transcription factor that regulates the transcription of genes involved in DNA damage/repair pathways (Frontini et al., 2009). In RS, the genes showing lasting increased expression patterns in cluster 4 included H2afx, Ercc2, and Gtf2h4. H2afx (H2ax) is variant gene coding for histone H2A, whose phosphorylated form is an early marker of DSBs and present in the whole population of spermatids throughout chromatin remodeling. Both Ercc2 and Gtf2h4 are components of the core-TFIIH basal transcription factor and involved in transcription-coupled NER. Interestingly, the expression of Jun and Junb are up-regulated both in PS and RS, which is consistent with a previous study showing increased mRNA levels of Jun and Junb following chemotherapeutic treatment (Delbès et al., 2009). In addition, an increase in the expression of ALKBH family members including Alkbh7 in cluster 2 of PS and Alkbh4 in cluster 4 of RS was observed. Some of the ALKBH proteins like ALKBH2 and ALKBH3 have been suggested to catalyze direct reversal of certain N-alkyl lesions, while other members might function beyond
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
7
Fig. 6. Expression patterns of DDR-related genes in response to MMS treatment in PS and RS. (A) Pipeline for the analysis of DDR-related gene. GO terms and literature was used to identify all DDR-related genes that were differentially expressed among sample groups in PS or RS. Next, those genes were sorted by fold change of 1.5 or greater. (B) and (C) K-means clustering was used to group the 374 DDR-related DEGs in PS and 158 in RS, respectively, into six different expression patterns.
DNA repair. ALKBH7 lacks DNA repair activity and has been reported to facilitate cell death following MMS treatment in 293T cells (Fu et al., 2013). 3.6. Validation of RNA-seq results for DDR related genes We further confirm the validity of the expression patterns of DDR related genes sorted by our RNA-seq based studies of PS and RS cells
treated with MMS. First, the expression change of thirteen wellknown DDR genes from our RNA-seq dataset matched well with the analysis of mRNA expression between PS and RS cells in a microarray dataset (GSE4193) from NCBI GEO database (Supplemental Fig. 2). Then, validation of DEGs between the 30 min repair group and the control was accomplished by qRT-PCR. Seven genes from PS and six genes from RS were amplified, which were important in BER (Apex2, Aptx, Ung, Xrcc1 and Lig3), NER and/or single-strand break repair (Aptx and
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
8
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Fig. 7. Validation of RNA-seq data. To validate differential expression between 30 min repair and control groups revealed by RNA-seq analysis, (A) seven DDR-related genes were amplified from PS and (B) six genes from RS by qRT-PCR. To verify the expression patterns of DDR-related genes in chronological order observed by RNA-seq analysis, (C) Mdc1 and Exo1 were amplified from PS and Gadd45b from RS by qRT-PCR.
Parp1), recombinational repair (Rad50 and Mlh3) and other DNA repair pathway (Fance, Fancg). The results indicated good agreement between the RNA-seq and qRT-PCR data (Fig. 7A and B). Finally, three more genes, Mdc1 and Exo1 from PS and Gadd45b from RS, were selected to verify the expression pattern in chronological order by qRT-PCR. The results showed that these genes exhibited a similar expression pattern as compared to RNA-seq, implying that Mdc1 and Exo1 were inhibited in PS cells as well as Gadd45b was activated in RS cells in response to MMS treatment (Fig. 7C). 4. Discussion Mammalian germ cells encounter various forms of DNA damage. In order to ensure genome integrity and the fidelity of transmission of hereditary information during spermatogenesis, damage should be completely repaired within a short time. However, little is known about the alterations to the transcription profile of male germ cells and the differences between the meiotic and post-meiotic stages in response to damage. In this RNA-seq study, we analyzed the dynamic transcriptional alterations in PS and RS affected by MMS. Comparative analyses of DEGs derived from PS and RS were conducted, in which the most obvious differences in gene expression profiles and major biological processes were observed between the 30 min repair and untreated control groups in both cell types. Although a large number of DEGs were identified with DEGseq, the DEGs manifesting higher fold change are relatively few, while those with 2-fold or more change account for only 13.4% of total genes in PS and 7.7% in RS, respectively. Canonical pathway analysis revealed that the ubiquitination pathway is the main pathway affected in common by the DEGs in every comparison of both PS and RS, which is not surprising when considered the importance of ubiquitin-proteasome system in spermatogenesis and DNA repair process (Richburg et al., 2014). It was reported that seventy E3 ubiquitin ligases were expressed in the mouse testis and that deficiency in several of these enzymes leads to infertility in male mice. Deficiency in CUL4A, for example, resulted in persistent DSBs in PS (Kopanja et al., 2011). In RS, E3 ligases have been shown to be fundamental during the spermiogenesis in its morphological differentiation through condensation of the spermatid DNA, supporting our analyses (Nakamura, 2013). When considering expression of DEGs
that are 2-fold or greater, the pathways affected are somewhat different, especially in PS in which several inflammatory signaling pathways are overrepresented with IL-17A being the most significantly overrepresented. Some inflammatory processes have been suggested to regulate normal spermatogenesis (O'Bryan and Hedger, 2008), where the role of IL-17A signaling, to our knowledge, has never been reported. The role of IL-17A during normal spermatogenesis or under alkylating stimuli requires further investigation. In addition, Network analysis indicated DDR is a common interactome network both in PS and RS following MMS treatment. Three major repair mechanisms specific for alkylation damage include direct DNA repair by the AlkB dioxygenase and O6-methylguanine-DNA methyltransferase (MGMT); and by well-known mechanisms of BER and NER. In K-means clustering analyses, continued activation of ALKBH7 in PS and ALKBH4 in RS, the mammalian homologs to bacterial AlkB dioxygenase, probably uncovers novel functions of ALKBH family members in germ cells following MMS treatment. Based on our RNA-seq data, activation of BER in conjunction with inhibition of the HR pathway occurred in response to MMS treatment in PS, indicating precise coordination between different repair pathways. The NER pathway may mainly contribute to the repair of DNA lesions in RS induced by MMS. Olsen et al. (2001) reported that BER was proficient and dominates in repairing MMS-induced DNA lesions in both human and rat in primary spermatocytes and round spermatids. However, the authors considered that it cannot be excluded based on a comet assay alone that other repair pathways besides BER participate in repair. Our results indicated the functionality of NER in mouse RS when facing MMS, which needs further elucidation. A previous report indicated that PS and RS exhibit stage-specific transcriptome differences during normal spermatogenesis. Our results showed dramatic differences in the patterns of DEGs that were 2-fold or greater between the 30 min repair and control groups when comparing PS and RS. Functional categorization of these unique DEGs revealed regulation of transcription-related terms that were specifically enriched in PS and DDR-related terms specifically enriched in RS. Previous studies with microarray have demonstrated that acute exposure of male rats to cyclophosphamide altered the expression of stress responsive genes during spermatogenesis in a germ cell-specific manner. RS were highly responsive to cyclophosphamides, while PS was not (Aguilar-Mahecha
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
et al., 2001). They explained that the transcriptional response to damage may be more rapid in PS than in RS, which was partially supported by our data with MMS in mice. In somatic cells, the cellular transcriptional response to MMS treatment has been conducted initially in yeast by microarray experiments. The authors showed 325 transcripts induced and 76 transcripts repressed by MMS (Jelinsky and Samson, 1999). A genome-wide screen on haploid gene deletion mutants identified a subset of genes that display a specific MMS response (Chang et al., 2002; Svensson et al., 2011). More recently, systematic examination of the proteins induced by MMS has been performed in human HeLa cells using mass spectrometry-based proteomics, suggesting that the ubiquitin-proteasome system plays an important role in the chromatin response to MMS treatment (Aslanian et al., 2014). Our RNA-seq data analysis indicates that a stage-specific response is characteristic for mouse spermatogenic cells following MMS treatment. In conclusion, we have depicted the landscape of the transcriptome in response to MMS “during” and “post” meiosis in mice, and identified several novel candidate genes and biological processes that are affected in PS and RS, which provides new evidence for transcriptional regulation during spermatogenesis during alkylation insult. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2016.10.006. Conflict of interest statement None declared. Acknowledgments The authors thank Dr. Chunsheng Han (Institute of Zoology, CAS) for instruction on germ cells isolation, Dr. Tie-Shan Tang (Institute of Zoology, CAS) for helpful discussion, M.S. Lijia Yu and Dr. Baofa Sun (Beijing Institute of Genomics, CAS) for assistance on data analysis, Dr. Paula L. Fischhaber (California State University Northridge, USA) for proofreading the manuscript. This work was supported by the Chinese National 973 Project [2013CB945000]; the National Natural Science Foundation of China [81630078, 31471331, and 31670822]; the Strategic Priority Research Program of the CAS [XDB14030300]; and the State Key Laboratory of Membrane Biology.
9
Aslanian, A., Yates 3rd, J.R., Hunter, T., 2014. Mass spectrometry-based quantification of the cellular response to methyl methanesulfonate treatment in human cells. DNA Repair (Amst) 15, 29–38. Bellve, A.R., Millette, C.F., Bhatnagar, Y.M., O'Brien, D.A., 1977. Dissociation of the mouse testis and characterization of isolated spermatogenic cells. J. Histochem. Cytochem. 25, 480–494. Chang, M., Bellaoui, M., Boone, C., Brown, G.W., 2002. A genome-wide screen for methyl methanesulfonate-sensitive mutants reveals genes required for S phase progression in the presence of DNA damage. Proc. Natl. Acad. Sci. U. S. A. 99, 16934–16939. Ciccia, A., Elledge, S.J., 2010. The DNA damage response: making it safe to play with knives. Mol. Cell 40, 179–204. Delbès, G., Chan, D., Pakarinen, P., Trasler, J.M., Hales, B.F., Robaire, B., 2009. Impact of the chemotherapy cocktail used to treat testicular cancer on the gene expression profile of germ cells from male Brown-Norway rats. Biol. Reprod. 80, 320–327. Frontini, M., Vijayakumar, M., Garvin, A., Clarke, N., 2009. A ChIP–chip approach reveals a novel role for transcription factor IRF1 in the DNA damage response. Nucleic Acids Res. 37, 1073–1085. Fu, D., Calvo, J.A., Samson, L.D., 2012. Balancing repair and tolerance of DNA damage caused by alkylating agents. Nat. Rev. Cancer 12, 104–120. Fu, D., Jordan, J.J., Samson, L.D., 2013. Human ALKBH7 is required for alkylation and oxidation-induced programmed necrosis. Genes Dev. 27, 1089–1100. Jelinsky, S.A., Samson, L.D., 1999. Global response of Saccharomyces cerevisiae to an alkylating agent. Proc. Natl. Acad. Sci. U. S. A. 96, 1486–1491. Johnston, D.S., Wright, W.W., Dicandeloro, P., Wilson, E., Kopf, G.S., Jelinsky, S.A., 2008. Stage-specific gene expression is a fundamental characteristic of rat spermatogenic cells and Sertoli cells. Proc. Natl. Acad. Sci. U. S. A. 105, 8315–8320. Kopanja, D., Roy, N., Stoyanova, T., Hess, R.A., Bagchi, S., Raychaudhuri, P., 2011. Cul4A is essential for spermatogenesis and male fertility. Dev. Biol. 352, 278–287. Meistrich, M.L., 2013. Effects of chemotherapy and radiotherapy on spermatogenesis in humans. Fertil. Steril. 100, 1180–1186. Nakamura, N., 2013. Ubiquitination regulates the morphogenesis and function of sperm organelles. Cells 2, 732–750. O'Bryan, M.K., Hedger, M.P., 2008. Inflammatory networks in the control of spermatogenesis: chronic inflammation in an immunologically privileged tissue? Adv. Exp. Med. Biol. 636, 92–114. Olsen, A.K., Bjortuft, H., Wiger, R., Holme, J., Seeberg, E., Bjoras, M., Brunborg, G., 2001. Highly efficient base excision repair (BER) in human and rat male germ cells. Nucleic Acids Res. 29, 1781–1790. Richburg, J.H., Myers, J.L., Bratton, S.B., 2014. The role of E3 ligases in the ubiquitin-dependent regulation of spermatogenesis. Semin. Cell Dev. Biol. 30, 27–35. Romrell, L.J., Bellve, A.R., Fawcett, D.W., 1976. Separation of mouse spermatogenic cells by sedimentation velocity. A morphological characterization. Dev. Biol. 49, 119–131. Roy Choudhury, D., Small, C., Wang, Y., Mueller, P.R., Rebel, V.I., Griswold, M.D., McCarrey, J.R., 2010. Microarray-based analysis of cell-cycle gene expression during spermatogenesis in the mouse. Biol. Reprod. 83, 663–675. Svensson, J.P., Pesudo, L.Q., Fry, R.C., Adeleye, Y.A., Carmichael, P., Samson, L.D., 2011. Genomic phenotyping of the essential and non-essential yeast genome detects novel pathways for alkylation resistance. BMC Syst. Biol. 5, 157. Wang, L., Feng, Z., Wang, X., Wang, X., Zhang, X., 2010. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. Zhao, W., Liu, W., Tian, D., Tang, B., Wang, Y., Yu, C., Li, R., Ling, Y., Wu, J., Song, S., Hu, S., 2011. wapRNA: a web-based application for the processing of RNA sequences. Bioinformatics 27, 3076–3077.
References Aguilar-Mahecha, A., Hales, B.F., Robaire, B., 2001. Acute cyclophosphamide exposure has germ cell specific effects on the expression of stress response genes during rat spermatogenesis. Mol. Reprod. Dev. 60, 302–311.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Research paper
Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids Hui Zhang a, Chuanchao Zhang a, Jinting Yan a, Zhongshuai Sun a, Shuhui Song b, Yazhou Sun a, Caixia Guo a,⁎ a b
CAS Key Laboratory of Genomic and Precision Medicine, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China Core Genomic Facility, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China
a r t i c l e
i n f o
Article history: Received 27 July 2016 Received in revised form 12 September 2016 Accepted 3 October 2016 Available online xxxx Keywords: Transcriptome RNA-seq Spermatogenesis MMS DNA damage response Differentially expressed genes
a b s t r a c t Spermatogenesis is threatened by DNA alkylating agents, one major category of DNA damaging agents. Currently, little is known about the alterations in transcriptome profiling of the mouse spermatogenic cells in response to DNA alkylation at distinct stages of spermatogenesis. In this study, RNA sequencing (RNA-seq) was performed in pachytene spermatocytes (PS) and round spermatids (RS) at 0 or 30 min following Methyl Methanesulfonate (MMS) treatment and with untreated controls. A large number of differentially expressed genes (DEGs) were identified by comparison of the three groups in PS and RS, respectively. Functional analyses of all DEGs highlighted the protein ubiquitination pathway and DNA damage response (DDR) network being the two main biological processes in common in the two cell types. Further analyses of the DEGs with 2-fold or more changes between 30 min repair and control group indicated that several cytokine signaling pathways were the most strongly affected in PS and DDR related pathways in RS, respectively. Gene ontology (GO) analyses directly showed differential biological process (BP) affected between PS and RS, with “regulation of transcription” being most overrepresented in PS and “cellular response to stress” in RS, respectively. Moreover, 374 DDR-related genes in PS and 158 in RS among all DEGs were filtered and clustered, which showed dynamic expression patterns in PS and RS. Our analyses provide a transcriptional landscape for male germ cells in response to MMS during spermatogenesis. © 2016 Elsevier B.V. All rights reserved.
1. Introduction In adult mammals, spermatogenesis is a continuous and well-orchestrated process, which can be divided into three principal stages: a mitotic stage in spermatogonia, a meiotic stage in spermatocytes and a post-meiotic stage in spermatids. During the lengthy meiotic prophase, the primary spermatocyte progresses through leptotene, zygotene, pachytene and diplotene stages to generate two secondary spermatocytes and then is rapidly converted into four haploid round spermatids. Subsequently, during the post-meiotic stage, round spermatids differentiate into elongated spermatid and mature spermatozoa. Maintenance of genomic integrity during male germ cell development is crucial for production of healthy offspring, which requires Abbreviations: RNA-seq, RNA sequencing; PS, pachytene spermatocytes; RS, round spermatids; MMS, methyl methanesulfonate; DEGs, differentially expressed genes; DDR, DNA damage response; GO, Gene ontology; BP, biological process; BER, base excision repair; NER, nucleotide excision repair; TM, Tail Moment; RPKM, Reads Per Kilobase of exon model per Million mapped reads; DAVID, Database for Annotation, Visualization and Integrated Discovery; IPA, Ingenuity Pathway Analysis. ⁎ Corresponding author. E-mail address: [email protected] (C. Guo).
efficient surveillance and repair mechanisms for the fidelity of genetic material transmissions. Spermatogenesis is constantly challenged by various endogenous and exogenous factors which induce DNA damage; male germ cells combat DNA damage through DDR involving a range of DNA repair mechanisms and alteration of gene expression (Ciccia and Elledge, 2010). Alkylating agents represent an important class of DNA damaging agents. Among them, S-adenosylmethionine (SAM) is a ubiquitous methyl donor, constantly contributing to endogenous DNA alkylation. Tobacco-specific N-nitrosamines and cancer chemotherapeutics such as cyclophosphamide, melphalan, and ifosfamide are familiar exogenous alkylating agents which probably threaten spermatogenesis of human beings (Meistrich, 2013). In the laboratory, methyl methanesulfonate (MMS) has served as one of the most extensively used alkylating agents in DNA repair studies. MMS predominantly produce N-methylation adducts including N7-methylguanine (7meG) and N3-methyladenine (3meA). In somatic cells, adducts induced by MMS are mainly repaired by direct DNA repair, by base excision repair (BER), and by nucleotide excision repair (NER) (Fu et al., 2012). The effects of MMS on spermatogenesis have been initially detected in several models animals including mouse, rat and drosophila. Due to the unavailability of efficient scientific tools, most of the studies have
http://dx.doi.org/10.1016/j.gene.2016.10.006 0378-1119/© 2016 Elsevier B.V. All rights reserved.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
2
H. Zhang et al. / Gene xxx (2016) xxx–xxx
mainly focused on the mutagenic effects of MMS on male germ cells, and emphasized expression changes of a small number of genes. Studies by Robaire et al. have described the alterations of stress responsive genes to cyclophosphamide using microarray in rat male germ cells (Aguilar-Mahecha et al., 2001). Several other large scale studies have elucidated the stage specificity and complexity of testicular transcription machinery with microarray experiments or RNA-seq during normal spermatogenesis (Johnston et al., 2008). However, the complete landscape of the gene expression patterns and DDR molecular events following MMS treatment in specific spermatogenic cells remains poorly understood. In this study, comprehensive transcriptome analyses in response to MMS treatment in PS and RS were conducted using RNA-seq and functional annotation of DEGs and highlight important biological processes, including DDR related pathways. 2. Materials and methods 2.1. Isolation of spermatogenic cells We used male C57BL/6Jx129 mice for the preparations of different testis cell populations. All animal experiments were reviewed and approved by the institute of Zoology, Institutional Animal Care and Use Committee and were conducted according to the committee's guidelines. The spermatogenic cell populations were isolated from the testes of 3-month-old mice using a StaPut 2–4% bovine serum albumin gradient as described previously (Bellve et al., 1977; Romrell et al., 1976). The PS and RS fractions were identified both by their size as well as by typical cellular morphological by phase contrast microscopy (Leica). The purity of PS exceeded 80% and that of RS exceeded 90% (Supplemental Fig. S1). Cell viability was N 95% as determined with the method of trypan blue exclusion. 2.2. Cell culture and MMS treatment Primary PS and RS cells (4 × 106 cells per sample) were treated with 0.12 mM MMS or left untreated as control for 1 h at 32 °C in DMEM/F12 culture medium. The cells were then placed on ice for reaction termination and washed two times with medium. For repair incubation, cells were maintained in DMEM/F12 medium supplemented with 10% FCS, pyruvate (0.1 mg/ml) and antibiotics (100 IU/ml penicillin; 0.1 mg/ml streptomycin) in a humidified atmosphere at 32 °C, 5% CO2 for up to 2 h. The cells were collected at different time points for comet assay or RNA extraction. 2.3. Comet assay Alkaline comet assays were performed according to the manufacturer's protocol (Trevigen). Samples were stained with SYBR Gold and images were taken by a fluorescent microscope (Leica) following alkaline electrophoresis. The average of Tail Moment (TM) was established by measuring at least 100 cells per slide using Comet Assay Software Project Casp-1.2.2 (University of Wroclaw). The reported TMs represented mean ± standard error (SE). The experiments were repeated three times.
software allowing up to four mismatches. We counted the number of reads mapped to exons, introns, splice junctions, intron-exon adjacent regions of the annotated genes, and intergenic regions using in-house built pipeline wapRNA (Zhao et al., 2011). Gene expression levels were quantified by RPKM (Reads Per Kilobase of exon model per Million mapped reads), and then DEGs were defined using the R-package DEGseq with the method MARS (MA-plot-based method with Random Sampling) model (Wang et al., 2010), and with the BenjaminiHochberg false-discovery rate (FDR) adjustment for multiple testing. Hierarchical clustering of all expressed genes was conducted using R software. 2.5. Function analysis of differentially expressed genes GO analyses of DEGs were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 (http://david. abcc.ncifcrf.gov/). GO BP was selected as the functional annotation category for this analysis. Enrichment map of DEGs was generated by Cytoscape 2.8.2 software installed with Enrichment Map plugin. Each node represents a pathway, and the size of the node reflects the number of genes in the pathway. The canonical pathways and interaction network were determined using Ingenuity Pathway Analysis (IPA) v9.0-3211 (Ingenuity Systems, Inc.) Ensemble accession was used as the identifier and Ingenuity knowledge gene database was used as a reference for the pathway and network analysis. 2.6. Quantitative real-time RT-PCR (qRT-PCR) validation qRT-PCR was performed to validate a subset gene expression changes predicted by RNA-seq, primers used in qRT-PCR were listed in Supplemental Table 1. One microgram of total RNA from each sample was reverse-transcribed in 20 μl of reaction medium using the First Strand cDNA Synthesis kit (Promega) following the manufacturer's instructions. PCR was detected using SYBR Green Real time PCR Master Mix (Toyobo) on a Real-time thermal cycler CFX96 (Bio-Rad). The target gene expression levels were calculated relative to the expression of mouse Gapdh, employing an optimized comparative Ct (ΔΔCt) value method. The experiments were performed in triplicate. Data are expressed as mean ± standard deviation (SD). 3. Results 3.1. Cellular repair of MMS-induced DNA lesions To investigate the cellular response to DNA damage reagent MMS treatment in spermatogenic cells of adult mouse, the purified PS and RS (Supplemental Fig. 1) were exposed to 0.12 mM MMS at 32 °C for 1 h and then were harvested at different time points. Efficient and
2.4. RNA sequencing Total RNA was extracted from MMS-treated PS and RS cells at 0 min and 30 min after repair incubation and the control group using Trizol reagent (Invitrogen), respectively. The total six samples were indicated as PScon, PS0, PS30, RScon, RS0 and RS30. Sequencing libraries were prepared according to the manufacturer's instructions (Illumina). After trimming low-quality reads, the high-quality reads were aligned against the reference genome of Mus musculus (NCBIM 37; ENSEMBL release 62; http://asia.ensembl.org/info/data/ftp/index.html) using BWA 0.5.8
Fig. 1. Repair of MMS-induced DNA lesions in PS and RS. PS and RS were exposed to 0.12 mM MMS for 1 h, respectively, and incubated for up to 2 h to allow repair. DNA damage (tail moment; mean ± SE) was measured with alkaline comet assay at indicated time points. At least 100 cells were scored for each sample.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
3
Table 1 Summary of sequence data and read alignment statistics for individual samples. Sample
Total reads
Uniquely mapped
Counted on exons
Counted on introns
Expressed genes
PScon PS0 PS30 RScon RS0 RS30
45,889,931 45,698,664 46,550,854 44,207,333 44,846,475 45,911,712
36,382,567 (79.28%) 34,855,055 (76.28%) 36,749,237 (78.94%) 35,206,820 (79.64%) 34,839,289 (77.69%) 34,994,342 (76.22%)
29,928,689 29,300,496 30,079,158 30,241,642 29,304,027 29,965,255
5,431,178 4,905,678 5,687,934 4,328,639 4,816,948 4,639,136
21,430 21,411 21,500 20,994 21,214 21,280
rapid repair of damaged DNA was observed by both PS and RS using the comet assay, and the repair kinetics were similar between the two cell types (Fig. 1). 3.2. Gene expression profiling of PS and RS in response to MMS treatment To gain a complete understanding of gene expression profiles of PS and RS with MMS-treatment, RNA-seq was employed. In each cell type, 3 samples were selected and designated as the control, 0 or 30 min group following treatment with MMS, representing the untreated, damaged and repaired status, respectively (a total of 6 samples were labeled as PScon, PS0, PS30; RScon, RS0, RS30). The total number of reads generated from each sample varied from 44,207,333 to 46,550,854, which were aligned to the reference genome in a unique manner, with mapping rates between 76.22% and 79.61% (Table 1). In total, 27,930 genes were expressed at least in one stage. Hierarchical clustering of these 27,930 genes showed a clear difference in gene expression pattern of sample groups from PS and RS and also revealed a farther correlation between 30 min repair and the control group in both PS and RS (Fig. 2A). DEG analysis showed the amounts of genes changed significantly (P b 0.001), including 3778 for PS0/PScon, 7680 for PS30/PS0 and 7482 for PS30/PScon in PS (Fig. 2B); 5897 for RS0/RScon, 3012 for RS30/RS0 and 4596 for RS30/RScon in RS (Fig. 2B). The total number of overlapping DEGs in PS and RS was 9705 and 7071, respectively
(Fig. 2C and D). Moreover, the number of DEGs with 2-fold or greater change and the range of fold change (FC) in each comparison of PS and RS show the highest difference in the number of DEGs between the 30 min repair and control groups in both PS (868) and RS (295) cells, consistent with the results of hierarchical clustering (Table 2). Special attention was thus paid to the comparisons between the 30 min repair and control groups in both PS and RS in the following analyses. The lists of top 30 up- and down-regulated DEGs comparing PS or RS are presented (Supplemental Tables 2 and 3, respectively).
3.3. Functional analysis of differentially expressed genes Canonical pathways and molecular networks related to DEGs were identified using the “Core Analysis” function in IPA for each comparison. The protein ubiquitination pathway is among the top canonical pathways in all comparisons of both PS and RS (Fig. 3). Other top canonical pathways affected are oxidative phosphorylation, mitochondrial dysfunction and EIF2 signaling in PS cells (Fig. 3A) and AMPK signaling, mitotic roles of polo-like kinase, cyclins and cell cycle regulation in RS cells, respectively (Fig. 3B). Network analyses revealed several interactome networks in common for PS and RS cells, such as developmental and hereditary disorder, RNA post-transcriptional modification, cancer, and processes related to the DDR (i.e., DNA replication, recombination and repair, cell cycle) (Supplemental Tables 4 and 5).
Fig. 2. Dynamic changes in gene expression and overlapping DEGs. (A) Hierarchical clustering for whole expressed genes in PS and RS sample groups. The default colors in R were used. Upregulated genes are indicated in white and downregulated in red. (B) Bar graph shows total number of DEGs and up- and down regulated gene numbers among comparisons in PS and RS. (C) and (D) The number of overlapping DEGs between comparisons in PS and RS, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
4
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Table 2 The number of DEGs with 2-fold change or greater. Comparisons
Total
Up regulated
Down regulated
Fold change range
PS0/PScon PS30/PS0 PS30/PScon RS0/RScon RS30/RS0 RS30/RScon
139 692 868 293 60 295
132 282 643 85 20 118
7 410 225 208 40 177
14.45 to −26.16 48.90 to −45.00 26.24 to −25.76 13.79 to −31.33 6.67 to −25.33 71.64 to −26.13
DEGs with 2-fold or greater changes in expression between the 30 min repair and control groups were further subjected to the same analyses using IPA. The most significant top 10 canonical pathways of PS and RS cells are presented in Fig. 4A and B. Notably, inflammatory responses including IL-17A, IL-17F, IL-8 and TNFR2 signaling are among the top canonical pathways in PS while cell cycle control and checkpoint regulation are among that in RS. The top one network in PS was embryonic development, organismal development, and developmental disorder, hereditary disorder, and immunological disease in RS, respectively
(Fig. 4C and D). Meanwhile, GO analyses based on the DAVID tool revealed that a large proportion of these genes in PS was enriched in BP, such as regulation of transcription, regulation of RNA metabolic process and embryonic organ development (Fig. 4E). In RS, the enriched genes belong to several different functional categories, including regulation of transcription, cellular response to stress, chromosome organization and immune effecter processes (Fig. 4F). 3.4. Comparisons of spermatocytes and spermatids We next focused on the difference in DEGs in response to MMS treatment that distinguishes PS from RS. Of total DEGs between the 30 min repair and control groups, 59% (4407/7482) was unique to PS and 33% (1521/4596) to RS (Fig. 5A). When considering those with 2fold or greater changes, the proportions specific to PS and RS were sharply increased to 94% (814/868) and 82% (241/295), respectively (Fig. 5B). Further GO BP enrichment analysis of 814 genes unique to PS revealed that “regulation of transcription (such as C/EBPbeta, EP300, Junb and Pou3f3)”, “regulation of RNA metabolic process (such as Atf7
Fig. 3. Top canonical pathways associated with DEGs through comparisons among the three groups by IPA. (A) DEGs from PS. (B) DEGs from RS.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
5
Fig. 4. Functional analysis of the expressed gene with 2-fold change or more between the 30 min repair and control groups in PS and RS. (A) and (B) The top 10 canonical pathways in which genes were enriched by IPA. (C) and (D) Interaction networks of genes identified by IPA. The genes in green were down-regulated and the genes in red were up-regulated with respect to the control group. (E) and (F) Function analysis of genes based on the DAVID GO BP analysis result. The cutoff parameters for enrichment analysis with Cytoscape software are: P b 0.005, FDR q b 0.05, overlap cutoff N 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
and Fos)”, and “embryonic development (such as Nkx2-6 and Lmo4)” were overrepresented (Fig. 5C red strip). Among the 241 DEGs specific to RS, the processes “cellular response to stress (such as Dmc1, Wrnip1 and Udg)”, “regulation of cell cycle (such as Ccnd1and Cdt1)” and “DNA integrity checkpoint (such as Msh2 and Mapk14)” were enriched (Fig. 5C blue strip). 54 DEGs common in PS and RS were enriched in “regulation of transcription (such as Abca2 and Atf3)”, “regulation of RNA metabolic process (such as Fosb and Bhlhe40)” and “response to protein stimulus (such as Cyr61 and Egr1)” (Fig. 5C grey strip). These results suggest that distinct biological processes and gene expression patterns are required for dealing with MMS-induced DNA damage repair during and after meiosis.
3.5. DNA damage and repair gene expression patterns in spermatocytes and spermatids following MMS treatment Top networks of all DEGs involve processes related to DDR (Supplemental Tables 4 and 5), we sought to gain as complete a picture as possible of the dynamic regulation of DNA damage and repair processes in PS and RS in the following analyses. We briefly summarized the genes involved in DDR events by integrating GO annotation with a literature search, in which 1251 genes were annotated by GO:0007049 (cellcycle control); 573 were annotated by GO:0006974 (cellular response to DNA damage stimulus); 550 cell cycle genes were expressed during spermatogenesis from the literature (Roy Choudhury et al., 2010) and
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
6
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Fig. 5. Comparisons of PS and RS in gene expressions and functions. (A) and (B) Venn diagram showing the distribution of all and 2-fold or greater change DEGs between 30 min and control groups in PS and RS, respectively. The numbers of DEGs in PS are indicated in the red circle, while the numbers of DEGs in RS are highlighted in the blue circle. (C) GO BP categories enriched in DEGs from Venn analysis in 3B. Red represents DEGs unique to PS while blue to RS, and grey representing common DEGs in PS and RS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
174 mouse homologs of human DNA repair genes from the Wood laboratory website (http://sciencepark.mdanderson.org/labs/wood/dna repair genes.html) (Fig. 6A). Based on these criteria, 1021 and 736 DDR related genes were differentially expressed in PS and RS following MMS treatment, respectively. These DDR related genes were then sorted into two groups: 1) genes that are “similarly expressed” in PS or RS (no change of N = 1.5-fold) and 2) genes that were “differentially expressed” in PS or RS (one or more changes of N = 1.5-fold) (Fig. 6A). We chose the fold change of 1.5 fold as the criterion in this analysis because we wanted to identify the alteration of DDR related genes in spermatogenic cells following MMS treatment in a more comprehensive manner, which was also used in a previous study (Roy Choudhury et al., 2010). Our results show that, 374 and 158 genes were sorted as “differentially expressed” in PS and RS, respectively (Fig. 6A). We further clustered these “differentially expressed” genes in PS and RS with K-means and identified six distinct expression patterns of DDRrelated genes in PS or RS, including three patterns showing gene expression changes in one direction (cluster 1, 4, 6 in PS; cluster 2, 4, 5 in RS) and the other three patterns in two opposite directions (cluster 2, 3, 5 in PS; cluster 1, 3, 6 in RS) (Fig. 6B and Supplemental Tables 6 and 7), which aids our understanding the activation or suppression of DDR pathways. Several noteworthy observations from these analyses follow. In PS, although clustered in distinct patterns, several important genes
required for DNA double-strand breaks (DSBs) by homologous recombination (HR) all showed significantly decreased expression when comparing the 30 min repair and control groups, including Rad50, Exo1, and Rev1 in cluster 1 and Brca2, Dna2, Atm in cluster 5. DNA mismatch repair gene Mlh3 was also down-regulated as shown in cluster1. As expected, Aptx in cluster 2 and Apex2 in cluster 6, both of which participate in BER, are up-regulated. The gene of highest alteration is Irf1 in cluster 4, a transcription factor that regulates the transcription of genes involved in DNA damage/repair pathways (Frontini et al., 2009). In RS, the genes showing lasting increased expression patterns in cluster 4 included H2afx, Ercc2, and Gtf2h4. H2afx (H2ax) is variant gene coding for histone H2A, whose phosphorylated form is an early marker of DSBs and present in the whole population of spermatids throughout chromatin remodeling. Both Ercc2 and Gtf2h4 are components of the core-TFIIH basal transcription factor and involved in transcription-coupled NER. Interestingly, the expression of Jun and Junb are up-regulated both in PS and RS, which is consistent with a previous study showing increased mRNA levels of Jun and Junb following chemotherapeutic treatment (Delbès et al., 2009). In addition, an increase in the expression of ALKBH family members including Alkbh7 in cluster 2 of PS and Alkbh4 in cluster 4 of RS was observed. Some of the ALKBH proteins like ALKBH2 and ALKBH3 have been suggested to catalyze direct reversal of certain N-alkyl lesions, while other members might function beyond
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
7
Fig. 6. Expression patterns of DDR-related genes in response to MMS treatment in PS and RS. (A) Pipeline for the analysis of DDR-related gene. GO terms and literature was used to identify all DDR-related genes that were differentially expressed among sample groups in PS or RS. Next, those genes were sorted by fold change of 1.5 or greater. (B) and (C) K-means clustering was used to group the 374 DDR-related DEGs in PS and 158 in RS, respectively, into six different expression patterns.
DNA repair. ALKBH7 lacks DNA repair activity and has been reported to facilitate cell death following MMS treatment in 293T cells (Fu et al., 2013). 3.6. Validation of RNA-seq results for DDR related genes We further confirm the validity of the expression patterns of DDR related genes sorted by our RNA-seq based studies of PS and RS cells
treated with MMS. First, the expression change of thirteen wellknown DDR genes from our RNA-seq dataset matched well with the analysis of mRNA expression between PS and RS cells in a microarray dataset (GSE4193) from NCBI GEO database (Supplemental Fig. 2). Then, validation of DEGs between the 30 min repair group and the control was accomplished by qRT-PCR. Seven genes from PS and six genes from RS were amplified, which were important in BER (Apex2, Aptx, Ung, Xrcc1 and Lig3), NER and/or single-strand break repair (Aptx and
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
8
H. Zhang et al. / Gene xxx (2016) xxx–xxx
Fig. 7. Validation of RNA-seq data. To validate differential expression between 30 min repair and control groups revealed by RNA-seq analysis, (A) seven DDR-related genes were amplified from PS and (B) six genes from RS by qRT-PCR. To verify the expression patterns of DDR-related genes in chronological order observed by RNA-seq analysis, (C) Mdc1 and Exo1 were amplified from PS and Gadd45b from RS by qRT-PCR.
Parp1), recombinational repair (Rad50 and Mlh3) and other DNA repair pathway (Fance, Fancg). The results indicated good agreement between the RNA-seq and qRT-PCR data (Fig. 7A and B). Finally, three more genes, Mdc1 and Exo1 from PS and Gadd45b from RS, were selected to verify the expression pattern in chronological order by qRT-PCR. The results showed that these genes exhibited a similar expression pattern as compared to RNA-seq, implying that Mdc1 and Exo1 were inhibited in PS cells as well as Gadd45b was activated in RS cells in response to MMS treatment (Fig. 7C). 4. Discussion Mammalian germ cells encounter various forms of DNA damage. In order to ensure genome integrity and the fidelity of transmission of hereditary information during spermatogenesis, damage should be completely repaired within a short time. However, little is known about the alterations to the transcription profile of male germ cells and the differences between the meiotic and post-meiotic stages in response to damage. In this RNA-seq study, we analyzed the dynamic transcriptional alterations in PS and RS affected by MMS. Comparative analyses of DEGs derived from PS and RS were conducted, in which the most obvious differences in gene expression profiles and major biological processes were observed between the 30 min repair and untreated control groups in both cell types. Although a large number of DEGs were identified with DEGseq, the DEGs manifesting higher fold change are relatively few, while those with 2-fold or more change account for only 13.4% of total genes in PS and 7.7% in RS, respectively. Canonical pathway analysis revealed that the ubiquitination pathway is the main pathway affected in common by the DEGs in every comparison of both PS and RS, which is not surprising when considered the importance of ubiquitin-proteasome system in spermatogenesis and DNA repair process (Richburg et al., 2014). It was reported that seventy E3 ubiquitin ligases were expressed in the mouse testis and that deficiency in several of these enzymes leads to infertility in male mice. Deficiency in CUL4A, for example, resulted in persistent DSBs in PS (Kopanja et al., 2011). In RS, E3 ligases have been shown to be fundamental during the spermiogenesis in its morphological differentiation through condensation of the spermatid DNA, supporting our analyses (Nakamura, 2013). When considering expression of DEGs
that are 2-fold or greater, the pathways affected are somewhat different, especially in PS in which several inflammatory signaling pathways are overrepresented with IL-17A being the most significantly overrepresented. Some inflammatory processes have been suggested to regulate normal spermatogenesis (O'Bryan and Hedger, 2008), where the role of IL-17A signaling, to our knowledge, has never been reported. The role of IL-17A during normal spermatogenesis or under alkylating stimuli requires further investigation. In addition, Network analysis indicated DDR is a common interactome network both in PS and RS following MMS treatment. Three major repair mechanisms specific for alkylation damage include direct DNA repair by the AlkB dioxygenase and O6-methylguanine-DNA methyltransferase (MGMT); and by well-known mechanisms of BER and NER. In K-means clustering analyses, continued activation of ALKBH7 in PS and ALKBH4 in RS, the mammalian homologs to bacterial AlkB dioxygenase, probably uncovers novel functions of ALKBH family members in germ cells following MMS treatment. Based on our RNA-seq data, activation of BER in conjunction with inhibition of the HR pathway occurred in response to MMS treatment in PS, indicating precise coordination between different repair pathways. The NER pathway may mainly contribute to the repair of DNA lesions in RS induced by MMS. Olsen et al. (2001) reported that BER was proficient and dominates in repairing MMS-induced DNA lesions in both human and rat in primary spermatocytes and round spermatids. However, the authors considered that it cannot be excluded based on a comet assay alone that other repair pathways besides BER participate in repair. Our results indicated the functionality of NER in mouse RS when facing MMS, which needs further elucidation. A previous report indicated that PS and RS exhibit stage-specific transcriptome differences during normal spermatogenesis. Our results showed dramatic differences in the patterns of DEGs that were 2-fold or greater between the 30 min repair and control groups when comparing PS and RS. Functional categorization of these unique DEGs revealed regulation of transcription-related terms that were specifically enriched in PS and DDR-related terms specifically enriched in RS. Previous studies with microarray have demonstrated that acute exposure of male rats to cyclophosphamide altered the expression of stress responsive genes during spermatogenesis in a germ cell-specific manner. RS were highly responsive to cyclophosphamides, while PS was not (Aguilar-Mahecha
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006
H. Zhang et al. / Gene xxx (2016) xxx–xxx
et al., 2001). They explained that the transcriptional response to damage may be more rapid in PS than in RS, which was partially supported by our data with MMS in mice. In somatic cells, the cellular transcriptional response to MMS treatment has been conducted initially in yeast by microarray experiments. The authors showed 325 transcripts induced and 76 transcripts repressed by MMS (Jelinsky and Samson, 1999). A genome-wide screen on haploid gene deletion mutants identified a subset of genes that display a specific MMS response (Chang et al., 2002; Svensson et al., 2011). More recently, systematic examination of the proteins induced by MMS has been performed in human HeLa cells using mass spectrometry-based proteomics, suggesting that the ubiquitin-proteasome system plays an important role in the chromatin response to MMS treatment (Aslanian et al., 2014). Our RNA-seq data analysis indicates that a stage-specific response is characteristic for mouse spermatogenic cells following MMS treatment. In conclusion, we have depicted the landscape of the transcriptome in response to MMS “during” and “post” meiosis in mice, and identified several novel candidate genes and biological processes that are affected in PS and RS, which provides new evidence for transcriptional regulation during spermatogenesis during alkylation insult. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2016.10.006. Conflict of interest statement None declared. Acknowledgments The authors thank Dr. Chunsheng Han (Institute of Zoology, CAS) for instruction on germ cells isolation, Dr. Tie-Shan Tang (Institute of Zoology, CAS) for helpful discussion, M.S. Lijia Yu and Dr. Baofa Sun (Beijing Institute of Genomics, CAS) for assistance on data analysis, Dr. Paula L. Fischhaber (California State University Northridge, USA) for proofreading the manuscript. This work was supported by the Chinese National 973 Project [2013CB945000]; the National Natural Science Foundation of China [81630078, 31471331, and 31670822]; the Strategic Priority Research Program of the CAS [XDB14030300]; and the State Key Laboratory of Membrane Biology.
9
Aslanian, A., Yates 3rd, J.R., Hunter, T., 2014. Mass spectrometry-based quantification of the cellular response to methyl methanesulfonate treatment in human cells. DNA Repair (Amst) 15, 29–38. Bellve, A.R., Millette, C.F., Bhatnagar, Y.M., O'Brien, D.A., 1977. Dissociation of the mouse testis and characterization of isolated spermatogenic cells. J. Histochem. Cytochem. 25, 480–494. Chang, M., Bellaoui, M., Boone, C., Brown, G.W., 2002. A genome-wide screen for methyl methanesulfonate-sensitive mutants reveals genes required for S phase progression in the presence of DNA damage. Proc. Natl. Acad. Sci. U. S. A. 99, 16934–16939. Ciccia, A., Elledge, S.J., 2010. The DNA damage response: making it safe to play with knives. Mol. Cell 40, 179–204. Delbès, G., Chan, D., Pakarinen, P., Trasler, J.M., Hales, B.F., Robaire, B., 2009. Impact of the chemotherapy cocktail used to treat testicular cancer on the gene expression profile of germ cells from male Brown-Norway rats. Biol. Reprod. 80, 320–327. Frontini, M., Vijayakumar, M., Garvin, A., Clarke, N., 2009. A ChIP–chip approach reveals a novel role for transcription factor IRF1 in the DNA damage response. Nucleic Acids Res. 37, 1073–1085. Fu, D., Calvo, J.A., Samson, L.D., 2012. Balancing repair and tolerance of DNA damage caused by alkylating agents. Nat. Rev. Cancer 12, 104–120. Fu, D., Jordan, J.J., Samson, L.D., 2013. Human ALKBH7 is required for alkylation and oxidation-induced programmed necrosis. Genes Dev. 27, 1089–1100. Jelinsky, S.A., Samson, L.D., 1999. Global response of Saccharomyces cerevisiae to an alkylating agent. Proc. Natl. Acad. Sci. U. S. A. 96, 1486–1491. Johnston, D.S., Wright, W.W., Dicandeloro, P., Wilson, E., Kopf, G.S., Jelinsky, S.A., 2008. Stage-specific gene expression is a fundamental characteristic of rat spermatogenic cells and Sertoli cells. Proc. Natl. Acad. Sci. U. S. A. 105, 8315–8320. Kopanja, D., Roy, N., Stoyanova, T., Hess, R.A., Bagchi, S., Raychaudhuri, P., 2011. Cul4A is essential for spermatogenesis and male fertility. Dev. Biol. 352, 278–287. Meistrich, M.L., 2013. Effects of chemotherapy and radiotherapy on spermatogenesis in humans. Fertil. Steril. 100, 1180–1186. Nakamura, N., 2013. Ubiquitination regulates the morphogenesis and function of sperm organelles. Cells 2, 732–750. O'Bryan, M.K., Hedger, M.P., 2008. Inflammatory networks in the control of spermatogenesis: chronic inflammation in an immunologically privileged tissue? Adv. Exp. Med. Biol. 636, 92–114. Olsen, A.K., Bjortuft, H., Wiger, R., Holme, J., Seeberg, E., Bjoras, M., Brunborg, G., 2001. Highly efficient base excision repair (BER) in human and rat male germ cells. Nucleic Acids Res. 29, 1781–1790. Richburg, J.H., Myers, J.L., Bratton, S.B., 2014. The role of E3 ligases in the ubiquitin-dependent regulation of spermatogenesis. Semin. Cell Dev. Biol. 30, 27–35. Romrell, L.J., Bellve, A.R., Fawcett, D.W., 1976. Separation of mouse spermatogenic cells by sedimentation velocity. A morphological characterization. Dev. Biol. 49, 119–131. Roy Choudhury, D., Small, C., Wang, Y., Mueller, P.R., Rebel, V.I., Griswold, M.D., McCarrey, J.R., 2010. Microarray-based analysis of cell-cycle gene expression during spermatogenesis in the mouse. Biol. Reprod. 83, 663–675. Svensson, J.P., Pesudo, L.Q., Fry, R.C., Adeleye, Y.A., Carmichael, P., Samson, L.D., 2011. Genomic phenotyping of the essential and non-essential yeast genome detects novel pathways for alkylation resistance. BMC Syst. Biol. 5, 157. Wang, L., Feng, Z., Wang, X., Wang, X., Zhang, X., 2010. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. Zhao, W., Liu, W., Tian, D., Tang, B., Wang, Y., Yu, C., Li, R., Ling, Y., Wu, J., Song, S., Hu, S., 2011. wapRNA: a web-based application for the processing of RNA sequences. Bioinformatics 27, 3076–3077.
References Aguilar-Mahecha, A., Hales, B.F., Robaire, B., 2001. Acute cyclophosphamide exposure has germ cell specific effects on the expression of stress response genes during rat spermatogenesis. Mol. Reprod. Dev. 60, 302–311.
Please cite this article as: Zhang, H., et al., Transcriptome analysis of the responses to methyl methanesulfonate treatment in mouse pachytene spermatocytes and round spermatids, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.10.006