Whole-genome methylation analysis of testicular germ cells from cryptozoospermic men points to recurrent and functionally relevant DNA methylation changes

Abstract

Background: Several studies have reported an association between male infertility and aberrant sperm DNA methylation patterns, in particular in imprinted genes. In a recent investigation based on whole methylome and deep bisulfite sequencing, we have not found any evidence for such an association, but have demonstrated that somatic DNA contamination and genetic variation confound methylation studies in sperm of severely oligozoospermic men. To find out whether testicular germ cells (TGCs) of such patients might carry aberrant DNA methylation, we compared the TGC methylomes of four men with cryptozoospermia (CZ) and four men with obstructive azoospermia, who had normal spermatogenesis and served as controls (CTR).

Results: There was no difference in DNA methylation at the whole genome level or at imprinted regions between CZ and CTR samples. However, using stringent filters to identify group-specific methylation differences, we detected 271 differentially methylated regions (DMRs), 238 of which were hypermethylated in CZ (binominal test, p < 2.2 × 10^-16). The DMRs were enriched for distal regulatory elements (p = 1.0 × 10^-6) and associated with 132 genes, 61 of which are differentially expressed at various stages of spermatogenesis. Almost all of the 67 DMRs associated with the 61 genes (94%) are hypermethylated in CZ (63/67, p = 1.107 × 10^-14). As judged by single-cell RNA sequencing, 13 DMR-associated genes, which are mainly expressed during meiosis and spermiogenesis, show a significantly different pattern of expression in CZ patients. In four of these genes, the promoter is hypermethylated in CZ men, which correlates with a lower expression level in these patients. In the other nine genes, eight of which downregulated in CZ, germ cell-specific enhancers may be affected.

Conclusions: We found that impaired spermatogenesis is associated with DNA methylation changes in testicular germ cells at functionally relevant regions of the genome. We hypothesize that the described DNA methylation changes may reflect or contribute to premature abortion of spermatogenesis and therefore not appear in the mature, motile sperm.

Keywords: Cryptozoospermia; DNA methylation; Differentially methylated regions; Male infertility; Single cell RNA sequencing; Spermatogenesis; Testicular germ cells.

Fig. 1

Testicular germ cell samples selection for whole genome bisulfite sequencing. a Schematic representation of the experimental design. b Box plot showing the methylation values of H19, MEST, DDX4 and XIST measured using deep bisulfite sequencing (DBS) of the supernatant fraction at day 3–4 of culture for the normal controls (CTR, teal, n = 4) and the cryptozoospermic (CZ, purple, n = 4) samples. No significant difference was found in the methylation values of the four genes between the two groups. c Box plot showing the results of the ploidy analysis of the day 0 single-cell suspension of the normal controls (CTR, teal, n = 4) and the cryptozoospermic (CZ, purple, n = 4) samples used for whole genome bisulfite sequencing (WGBS). d Stacked bar plots showing the percentages of tubules containing germ cells (most advanced germ cell type shown), only Sertoli cells, or tubular shadows in each biopsy from which the samples for WGBS were prepared. e Box plot showing the results of the ploidy analysis of the supernatant fraction at day 3–4 of culture for the normal controls (CTR, teal, n = 4) and the cryptozoospermic (CZ, purple, n = 4) samples. No significant difference was found in the cellular composition of the supernatant fraction between the two groups

Fig. 2

CTR-CZ DMRs are associated with 132 genes. a Global comparison of methylomes from control (CTR) and cryptozoospermic (CZ) testicular germ cells. PCA generated for ca. 20 million CpG loci where all samples show methylation values. Only loci with minimum coverage of five in all samples and minimum mapping quality of 10 are considered. CTR testicular germ cell samples in teal, CZ samples in purple. b Cluster analysis of the methylation values for the 271 CTR-CZ DMRs in the eight TGC samples. c Enrichment/depletion of DMRs for functional genomic regions. LMR, low-methylated region; UMR, unmethylated regions; CGI, CpG islands. d DMRs are associated with 132 genes by overlapping genes, promoters and/or “double-elite” enhancers (GeneHancer)

Fig. 3

Characterization of the expression of the DMR associated genes in the scRNA-seq datasets. a Right panel: uniform manifold approximation and projection (UMAP) plot showing the germ cell subset of the scRNAseq data published in [26]. Left panel: UMAP plots showing the cells obtained from three normal controls (CTR, n = 14,098) and from the three cryptozoospermic patient samples (CZ, n = 5939). The cells are color coded according to their identity defined by the assignment published in [26]. b UMAP plot showing the integrated CTR-CZ germ cell dataset aligned along the latent time. The cells are color-coded according to their progression along the latent time. The undifferentiated spermatogonia cluster was set as starting point of the differentiation process. c Heatmap showing the normalized expression of the 55 DMR-associated genes with more than 500 counts in the CTR dataset. The cells are plotted along the latent time with the undifferentiated spermatogonia as starting point on the left side. The clustering analysis identified 5 clusters and left 12 genes unclustered. d Line plots showing the normalized expression along the latent time of the 55 DMR-associated genes with high expression grouped according to the belonging cluster in the CTR (teal) dataset

Fig. 4

DNA methylation levels and functional characterization of the DMRs associated with genes showing spermatogenesis-regulated expression. Box plots show the distribution of the methylation values for CTR (teal, n = 4) and CZ (purple, n = 4) (Additional file 1: Table S8). The overlaps of each DMR with specific genomic features are shown: exons (yellow), introns (light-blue), lncRNAs (grey), promoters (orange), UMRs (unmethylated regions, green), LMRs (low-methylated regions, dark-blue) and GeneHancer “double-elite” regulatory regions (pink). The 13 genes in bold were shown to be differentially expressed between CTR and CZ (see Fig. 5)

Fig. 5

Characterization of the 13 differentially expressed DMR-associated genes between CTR and CZ patients. a Line plots showing the expression along the latent time of the 11 DMR associated differentially expressed genes in CTR (teal) and CZ (purple) dataset determined by tradeSeq analysis. The dashed lines mark the knots dividing the three knot groups. The grey areas identify the knot groups in which statistical significance is reached for each gene. Values can be found in Additional file 1: Table S11. b Box plots showing the average expression values of the DMR-associated genes that resulted to be differentially expressed using MAST while comparing the CTR (teal, n = 3) and CZ scRNA-seq datasets (purple, n = 3). Values can be found in Additional file 1: Table S12