CTCF-mediated 3D chromatin sets up the gene expression program in the male germline

Abstract

Spermatogenesis is a unidirectional differentiation process that generates haploid sperm, but how the gene expression program that directs this process is established is largely unknown. Here we determine the high-resolution three-dimensional (3D) chromatin architecture of mouse male germ cells during spermatogenesis and show that CTCF-mediated 3D chromatin dictates the gene expression program required for spermatogenesis. In undifferentiated spermatogonia, CTCF-mediated chromatin interactions between meiosis-specific super-enhancers (SEs) and their target genes precede activation of these SEs on autosomes. These meiotic SEs recruit the master transcription factor A-MYB (MYBL1) in meiotic spermatocytes, which strengthens their 3D contacts and instructs a burst of meiotic gene expression. We also find that at the mitosis-to-meiosis transition, the germline-specific Polycomb protein SCML2 facilitates the resolution of chromatin loops that are specific to mitotic spermatogonia. Moreover, SCML2 and A-MYB help shape the unique 3D chromatin organization of sex chromosomes during meiotic sex chromosome inactivation. We propose that CTCF-mediated 3D chromatin organization regulates epigenetic priming that directs unidirectional differentiation, thereby determining the cellular identity of the male germline.

Extended Data Fig. 1 |. Germ cell isolation for Hi-C libraries and evaluation of biological replicates.

a, Representative DAPI-stained images of isolated cell fractions for wild-type (WT) pachytene spermatocytes (PS), round spermatids (RS), Scml2-KO PS, Scml2-KO RS, and A-myb-mutant PS. Scale bars: 20 μm. b, Box plots showing the percent purity of cell fractions obtained from sedimentation velocity at unit gravity for each biological replicate (replicates 1 and 2) for the WT PS, WT RS, A-myb-mutant PS, and Scml2-KO PS Hi-C library. Numbers (n) along the top indicate the numbers of fractions used to prepare the corresponding library replicates below. The box indicates the 25th, median and 75th percentiles, and the whiskers indicate the 10th to 90th percentiles. Means and standard deviations for the purities of each cell fraction comprising the WT PS libraries: 80% ± 3.3% (replicate 1) and 86% ± 1.1% (replicate 2); the WT RS libraries: 91% ± 0.4% (replicate 1) and 88% ± 1.6% (replicate 2); A-myb-mutant PS libraries: 95% ± 0.7% (replicate 1) and 90% ± 0.8% (replicate 2); the Scml2-KO PS libraries: 82% ± 3.6% (replicate 1) and 86% ± 1.4% (replicate 2). c, Scatter plots showing the Pearson correlation between the two biological replicates for the Hi-C dataset. Pearson correlation values are computed by HiCExplorer using the two-sided Pearson method. d, Heat map showing the pairwise correlation coefficients among each cell type. Two biological replicates were merged for this analysis. e, Representative images of chromosome spreads after isolation of WT and A-myb-mutant PS. Cells were collected from three or more mice from each genotype and observed individually. Bars: 10 μm.

Extended Data Fig. 2 |. CTCF distribution and chromatin loops in THY1⁺ and PS.

a, Scatter plots showing the Pearson correlation between two biological replicates for CTCF ChIP-seq data in THY1⁺ and PS. Pearson correlation values are computed by deepTools using the two-sided Pearson method. b, Percentage and number of overlaps between anchor sites of CTCF-anchored loops and those of enhancer-promoter loops, and overlaps between anchor sites of Polycomb loops and those of enhancer-promoter loops in THY1⁺ and PS. c, Hi-C pile-up analysis of anchor sites of CTCF-anchored loops, enhancer-promoter loops, Polycomb loops in THY1⁺ and PS. The normalized contact strength in the central pixel is displayed on the top left. d, Ratio of accumulation of CTCF, H3K4me3/H3K27ac, or H3K27me3 at TAD boundaries in THY1⁺ and PS.

Extended Data Fig. 3 |. Hi-C analysis of Scml2-KO PS and RS.

a, Hi-C maps showing chromatin loops (25-kb bins) in KIT⁺, WT PS, and Scml2-KO PS at indicated genomic locations at Chromosome 4 (left) and Chromosome 7 (right). Chromatin loops are indicated by black circles. The dotted circles in the Hi-C map of Scml2-KO cells indicate the loops also detected in KIT⁺. b, Scatter plots showing the Pearson correlation between 2 biological replicates of CTCF ChIP-seq data in Scml2-KO PS. c, Heat map showing the pairwise correlation coefficients for CTCF ChIP-seq in WT THY1⁺, WT PS, and Scml2-KO PS. d, Chromatin loop length (Mb) from each Hi-C data set. The number of loops used in the analysis (WT PS: n = 1,223, Scml2-KO PS: n = 1,428). The box indicates the 25th, median and 75th percentiles, and the dot in the box indicates mean. Statistical analysis is based on Levene’s test and Wilcoxon rank sum test. **** indicates P < 2e⁻¹⁶. e, Hi-C pile-up analysis of anchor sites of chromatin loops in THY1⁺ with each modification in WT PS and Scml2-KO PS. f–h, ChIP-seq data (average tag density plot and heat maps) for indicated modifications at the anchor sites of enhancer-promoter loops (for H3K4me3 in f, for H3k27ac in g) or Polycomb loops in THY1⁺ (h). i, ChIP-seq data (average tag density plot and heat maps) for indicated modifications at the THY1⁺-unique CTCF binding sites. j, Hi-C pile-up analysis showing averaged intersections of Class I or Class II bivalent domains. The genomic coordinates of Class I or Class II bivalent domains were downloaded from. The normalized contact strength in the central pixel is shown on the top left.

Extended Data Fig. 4 |. Hi-C analysis of A-myb mutant PS.

a, Number and overlap of TAD boundaries between Scml2-KO PS and A-myb-mutant PS. b, Number of the pachytene piRNA clusters in A or B compartment. c, Hi-C pile-up analysis of pachytene piRNA clusters in WT PS and A-myb-mutant PS. The genomic coordinates of pachytene piRNA clusters were downloaded from.

Extended Data Fig. 5 |. Hi-C analysis of meiotic SEs.

a, Pie chart showing the percentage of meiotic SEs located around chromatin loops. The anchor site and 0.4 Mb upstream and downstream of the anchor site were defined as the genomic region associated with the PS chromatin loop. X and Y chromosomes were excluded from this analysis. b, Heat maps showing normalized Hi-C interaction frequencies (10-kb bins) in KIT⁺ and PS. Black arrows highlight meiotic SE-to-SE interactions. ChIP-seq data of H3K27ac in KIT⁺ and PS are displayed blue and red, respectively. c, Track view showing meiotic SEs (black bars, top), H3K27ac in PS (red), and Hi-C interaction from the meiotic SEs (blue curves). d, Hi-C pile-up analysis of Hi-C interaction from the meiotic SEs in each cell type. The normalized contact strength in the central pixel is shown on the top left.

Extended Data Fig. 6 |. Meiotic SE interacting genes.

a, A heatmap showing RNA-seq profiles of 587 expressed genes on Hi-C interacting loci of meiotic SEs. b, Gene ontology (GO) analysis of genes interacting with meiotic SEs. c, A heatmap and violin plots showing gene expression of 26 ‘spermatogenesis’ genes that showed Hi-C interaction with meiotic SEs during spermatogenesis. The box in the violin plots indicates the 25th, median, and 75th percentiles, and the dot in the box indicates the mean. All genes that show at least Transcripts per million (TPM) 4 at one developmental stage are shown. d, e, Average tag densities and heatmaps for H3K4me2, H3K4me3, and H3K27ac ChIP-seq enrichment 5,000 bp upstream and downstream of TSSs. Genes adjacent to meiotic SEs (d) and genes interacting with meiotic SEs (e) are analyzed.

Extended Data Fig. 7 |. Meiotic SEs poising with 3D chromatin.

a, Numbers of CTCF binding sites in PS that overlap or do not overlap with meiotic SE and meiotic SE interacting loci. b, Heat maps showing CTCF enrichment in THY1⁺ and PS at CTCF binding sites on meiotic SEs in PS. c, Venn diagram showing the intersection between all the genes at meiotic SE interacting genomic regions (blue) and A-MYB binding genes (green). d, Average tag density and heatmap for A-MYB ChIP-seq enrichment 3,000 bp upstream and downstream of TSSs of meiotic SE interacting genes. e, Principal component analysis plot showing biological replicates of RNA-seq samples. f, Comparison of transcriptomes between Ctrl and A-myb-mutant PS. Differentially expressed genes (DEGs: Log2FoldChange > 2, Padj < 0.05, binominal test with Benjamini-Honchberg correction) in A-myb-mutant PS (n = 3, respectively) after ERCC normalization are colored (red: upregulated in A-myb-mutant PS; blue: downregulated in A-myb-mutant PS). g, Violin plots showing gene expression of 509 expressed genes among 611 meiotic SE-interacting genes. The box indicates the 25th, median and 75th percentiles, and the dot in the box indicates mean. Statistical analysis is based on Levene’s test and Wilcoxon rank sum test. P value is on the graph. h, Venn diagram showing the overlap between CTCF binding sites in WT PS and A-MYB binding sites in whole testes.

Extended Data Fig. 8 |. Characteristics of chromosome interactions in the X chromosome.

a, Numbers of meiotic SEs on autosomes and the X chromosome. b, Numbers of Hi-C interacting loci of meiotic SEs on autosomes and on the X chromosome. c, Hi-C interaction frequency probabilities P of the X chromosome stratified by genomic distance s for each cell type shown (100-kb bins). d, Comparison of Hi-C interaction frequency probabilities P between autosomes and the X chromosome for each cell type (100-kb bins). e, Log₂ ratio comparisons of the Hi-C interaction frequencies (100-kb bins, chromosome 7 and X) for successive cell types. Arrows indicate the distance of chromosome interactions. 10-kb bins normalized Hi-C matrices were used for the zoom-in. f, Hi-C interaction frequency probabilities P stratified by genomic distance s for KIT⁺, WT PS, Scml2-KO PS, WT RS, and Scml2-KO RS shown (100-kb bins). Autosomes and the X chromosome are shown separately. g, Heat maps showing normalized Hi-C interaction frequencies (100-kb bins, chromosome X) in WT RS, Scml2-KO RS. Red and blue Hi-C maps represent a log₂ ratio comparison of Hi-C interaction frequencies between wild-type and Scml2-KO RS. h, Hi-C interaction frequency probabilities P stratified by genomic distance s for KIT⁺, WT PS, and A-myb-mutant PS (100-kb bins). Autosomes and the X chromosome are shown separately.

Extended Data Fig. 9 |. Interchromosomal interactions.

a–d, Heat map showing genome-wide normalized Hi-C interaction frequencies (250-kb bins) in THY1⁺ (a), WT PS (b), Scml2-KO PS (c) and A-myb-mutant PS (d).

Fig. 1 |. 3D chromatin reprogramming and inter-TAD chromatin loop formation in meiosis.

a, Schematic of the stages of mouse spermatogenesis analyzed in this study. b, Hi-C maps showing normalized Hi-C interaction frequencies (100-kb bins) in THY1⁺ spermatogonia, KIT⁺ spermatogonia, PSs, and RSs. For zoomedin views, normalized Hi-C matrices at a 10-kb bin resolution were used. Black arrows indicate chromatin loops. c, Hi-C interaction frequency probabilities (P) stratified by genomic distance (s) for each cell type (100-kb bins). All autosomes were analyzed. Contact probability, P(s). d, Hi-C interaction heat maps (25-kb bins) in THY1⁺ spermatogonia and PSs. Chromatin loops are indicated by black arrows, red lines in THY1⁺ spermatogonia and blue lines in PSs. e, Number of chromatin loops detected from each Hi-C dataset. f, Number of unique chromatin loops, comparing each pairwise developmental stage. g, Chromatin loop length (Mb) from each Hi-C data set. The number of loops used in the analysiswas equal to the number shown in e (THY1⁺, n = 3,562; KIT⁺, n = 3,336; PS, n = 1,223; RS, n = 609). The box limits indicate the 25th and 75th percentiles, the center line indicates the median and the dot in the box indicates the mean. Statistical analysis is based on Levene’s test and Wilcoxon rank-sum test. ****P < 2 × 10⁻¹⁶. h, Hi-C pile-up analysis of chromatin loops in each cell type. Color represents normalized contact strength in the log scale. The normalized contact strength values in the central pixel are shown on the top left.

Fig. 2 |. TAD and chromatin loop reorganization during spermatogenesis.

a, Number of TADs detected in each Hi-C data set (25-kb bins). b, Venn diagram showing the overlap of TAD boundaries in each developmental stage, THY1⁺ and KIT⁺ (left), KIT⁺ and PS (middle) and PS and RS (right). c, Hi-C pile-up analysis of TAD boundaries in each cell type. Using 10-kb bins Hi-C data around the central pixel, normalized contact strength was visualized in a heat map on a log scale, with color intensity reflecting interaction strength. d, Hi-C pile-ups of TADs and their adjacent TADs from 10 kb bin Hi-C data in each cell type. The dotted regions represent interactions between adjacent TADs. e, Modifications at the anchor sites of chromatin loops. Overlap with the anchor sites of chromatin loops in THY1+ spermatogonia and PSs. f, CTCF enrichment at the anchor sites of CTCF-anchored loops in THY1⁺ spermatogonia (detected in e, 2,532 sites). Heat maps for each locus are shown at the bottom. g, Model showing TAD and chromatin-loop reorganization at the mitosis-to-meiosis transition.

Fig. 3 |. SCML2 is required for the resolution of spermatogonia-type 3D chromatin.

a, Heat maps showing normalized Hi-C interaction frequencies (100-kb bins) in WT PS and Scml2-KO PS (left), and WT RS and Scml2-KO RS (right). Red and blue Hi-C maps represent a log₂ ratio comparison of Hi-C interaction frequencies between WT and Scml2-KO. b, The number of TADs detected from each Hi-C data set (25-kb bins) in WT PSs, Scml2-KO PSs, WT RSs and Scml2-KO RSs. c, Venn diagram showing the overlap between all KIT⁺ TAD boundaries and Scml2-KO PS TAD boundaries (left), and the overlap between KIT⁺-specific TAD boundaries and Scml2-KO PS TAD boundaries (right). KIT⁺-specific boundaries are defined by excluding TAD boundaries detected in WT PS. d, Hi-C pile-up analysis of KIT⁺-specific TAD boundaries in WT PS and Scml2-KO PS. e, Number of chromatin loops detected in each Hi-C dataset in WT PS, Scml2-KO PS, WT RS and Scml2-KO RS. f, Hi-C pile-up analysis of chromatin loops in each cell type. The normalized contact strength in the central pixel is shown on the top left. g, Number of specific and common chromatin loops between WT PSs and Scml2-KO PSs; 677 Scml2-KO PS-specific loops overlapped with loops detected in KIT+ spermatogonia. Overlapping loops were detected by Juicer. h, ChIP-seq data (average tag density plot and heat maps) for SCML2 in GS cells at the anchor sites of KIT⁺-specific loops (1,165 loci) or PS-specific loops (623 loci). i, CTCF enrichment in WT PSs and Scml2-KO PSs at the anchor site of CTCF-anchored loops in THY1+ spermatogonia. j, Violin plots of RNA-seq reads converted to log10(transcripts per million (TPM) + 1) value for genes associated with Scml2-KO PS specific loops in KIT⁺, WT PS and Scml2-KO PS. 1,243 genes were identified by extracting genes present in the anchor site of Scml2-KO PS-specific loops. The box indicates the 25^th, median, and 75^th percentiles and the dot in the box indicates mean. Statistical analysis is based on Levene’s test and Wilcoxon rank sum test. P values are on the graph. k, Model of resolution of spermatogonia-type 3D chromatin by SCML2.

Fig. 4 |. A-MYB is required for the formation of meiotic-type 3D chromatin.

a, Left and middle, heat maps showing normalized Hi-C interaction frequencies (100-kb bins) in WT PSs and A-myb-mutant (A-myb-mut) PSs. Right, red and blue Hi-C maps represent a log2 ratio comparison of Hi-C interaction frequencies between WT and A-myb-mutant PSs. b, Number of TADs detected from each Hi-C data set (25-kb bins) in WT PSs and A-myb-mutant PSs. c, Venn diagram showing the overlap between all TAD boundaries in KIT⁺ spermatogonia and those in A-myb-mutant PSs (left), and the overlap between the groups (right). KIT⁺-specific boundaries are defined by excluding boundaries detected in WT PS. d, Hi-C pile-up analysis of KIT⁺-specific TAD boundaries in WT PS and A-myb-mutant PSs. e, Number of chromatin loops detected in each Hi-C data set in WT PSs and A-myb-mutant PSs. Pink area in the graph of A-myb mutant PS indicates that the same loops are detected in WT PS (357 loops). f, Hi-C pile-up analysis of chromatin loops in each cell type. The normalized contact strength in the central pixel is shown on the top left. g, ChIP–seq data (average tag density plot and heat maps) for A-MYB using whole testis at the regions adjacent to TSS overlapping the anchor sites of the KIT⁺-specific chromatin loops (13,605 loci) and PS-specific chromatin loops (6,823 loci). h, Venn diagram showing the intersection of genes located at anchor sites of chromatin loops in PS (blue) and all A-MYB-bound genes (green). The overlap is statistically significant (P = 8.5 × 10⁻⁵⁹), as determined by the hypergeometric test, compared with the proportion of all A-MYB-bound genes to all RefSeq genes (5,929/22,661). i, Model of the establishment of meiotic-type chromatin loops by A-MYB.

Fig. 5 |. Meiotic super-enhancers are poised with 3D chromatin.

a, Track view showing meiotic SEs, H3K27ac, and chromatin loops in PSs (top) and an enlarged view of the boxed area (bottom). b,c, Hi-C pile-up analysis of averaged intersections of mitotic SEs (b) and meiotic SEs (c). d, Track view showing CTCF distribution and Hi-C interactions around meiotic SEs in PS. Pink highlights indicate CTCF binding sites that do not overlap with meiotic SEs; blue highlights indicate CTCF binding sites that overlap with meiotic SEs and their loops. e, CTCF binding and Hi-C maps of THY1⁺ spermatogonia and PSs around meiotic SEs. f, Venn diagram showing the overlap between all CTCF-binding sites in THY1+ spermatogonia and CTCF-binding sites in PSs on meiotic SEs. When considering the ratio of all CTCF-binding sites overlapping between THY1⁺ spermatogonia and PSs to all CTCF-binding sites in PSs (7,094/14,534 loci), the association (406 overlapping loci/617 CTCF-binding sites in PSs on meiotic SEs) is statistically significant (P = 2.67 × 10⁻¹⁸, hypergeometric test). g, Hi-C pile-up analysis of CTCF-binding sites overlapping with meiotic SEs and their interacting genomic regions. h, Track view showing the distributions of A-MYB binding and H3K27ac around meiotic SEs. Hi-C interaction from the meiotic SEs is also shown. i,j, Hi-C pile-up analysis of meiotic SE loci (i) and meiotic SE-interacting loci (j) in WT PSs and A-myb-mutant PSs. The normalized contact strength in the central pixel is shown on the top left. k, Model of the preset of 3D chromatin at meiotic SE loci through CTCF in mitotic spermatogonia. A-MYB strengthens these 3D contacts in meiotic PSs. l, Left, schematic for CRISPR deletion. Two guide RNAs (gRNAs) were designed to target an adjacent region from a CTCF-binding site at a meiotic SE locus. Right, expression of the Tex101 gene, as measured by quantitative reverse transcription PCR (qRT–PCR). Two independent mutant ESC lines (1 and 2) were examined. Error bars represent mean ± s.e.m.: P values were determined by unpaired one-tailed t-tests. Four independent experiments were performed.

Fig. 6 |. SCML2 and A-MYB establish unique 3D chromatin of the meiotic sex chromosomes.

a, Hi-C maps of the X chromosome showing normalized Hi-C interaction frequencies (100-kb bins) in WT THY1⁺ spermatogonia, KIT⁺ spermatogonia, PSs and RSs. b, Left, heat maps of the X chromosome showing normalized Hi-C interaction frequencies (100-kb bins) in Scml2-KO PSs. Right, red and blue Hi-C maps represent a log₂ ratio comparison of Hi-C interaction frequencies between WT and Scml2-KO PSs. c, Left, heat maps of the X chromosome showing normalized Hi-C interaction frequencies (100-kb bins) in A-myb-mutant PSs. Right, red and blue Hi-C maps represent a log2 ratio comparison of Hi-C interaction frequencies between WT and A-myb-mutant PS (right). d, Heat maps showing normalized Hi-C interchromosomal interactions (250-kb bins, chromosomes 1 and X) for WT THY1⁺ spermatogonia, WT PSs, Scml2-KO PSs and A-myb-mutant PSs. e, Heat maps showing normalized Hi-C interchromosomal interactions (250-kb bins, chromosomes 1 and 2) for WT THY1⁺ spermatogonia, WT PSs, Scml2-KO PSs and A-myb-mutant PSs. f, Model for the establishment of a unique 3D chromatin in the XY body and segregation of XY from autosomes in PS. g, Model of interchromosomal interactions in PSs. h, Schematic of the molecular pathway that establishes an XY-unique 3D chromatin in PSs.

Fig. 7 |. Models of 3D chromatin dynamics and gene regulation on autosomes and sex chromosomes during spermatogenesis.

a, Model showing the changes in chromosome interactions from mitotic spermatogonia to meiotic pachytene spermatocytes on autosomes. b, Model of the functional relationship between SCML2 and A-MYB. c, Model of 3D chromatin dynamic on the sex chromosomes during meiotic prophase I. At the onset of MSCI, at the early pachytene stage, DDR initiates MSCI, and subsequently SCML2 and A-MYB establish unique 3D chromatin of the sex chromosomes and segregate the sex chromosomes from the autosomes at the mid-pachytene stage.