High resolution maps of chromatin reorganization through mouse meiosis reveal novel features of the 3D meiotic structure

When germ cells transition from the mitotic cycle into meiotic prophase I (MPI), chromosomes condense into an array of chromatin loops that are required to promote homolog pairing and genetic recombination. To identify the changes in chromosomal conformation, we isolated nuclei on a trajectory from spermatogonia to the end of MPI. At each stage along this trajectory, we built genomic interaction maps with the highest temporal and spatial resolution to date. The changes in chromatin folding coincided with a concurrent decline in mitotic cohesion and a rise in meiotic cohesin complexes. We found that the stereotypical large-scale A and B compartmentalization was lost during meiotic prophase I alongside the loss of topological associating domains (TADs). Still, local subcompartments were detected and maintained throughout meiosis. The enhanced Micro-C resolution revealed that, despite the loss of TADs, higher frequency contact sites between two loci were detectable during meiotic prophase I coinciding with CTCF bound sites. The pattern of interactions around these CTCF sites with their neighboring loci showed that CTCF sites were often anchoring the meiotic loops. Additionally, the localization of CTCF to the meiotic axes indicated that these anchors were at the base of loops. Strikingly, even in the face of the dramatic reconfiguration of interphase chromatin into a condensed loop-array, the interactions between regulatory elements remained well preserved. This establishes a potential mechanism for how the meiotic chromatin maintains active transcription within a highly structured genome. In summary, the high temporal and spatial resolution of these data revealed previously unappreciated aspects of mammalian meiotic chromatin organization.


Figure S1. Isolation of stage specific nuclei through spermatogenesis.
(A) A schematic of the strategy for isolating nuclei from each stage.Intra-nuclear markers used for sorting were labeled with different colors.The relative expression level of each protein is represented as the thickness of their corresponding lines.Nuclei from unDiff.SGA to the Meiotic G1 stage were isolated with a combination of DAPI, PLZF, DMRT1 and STRA8.Nuclei from the Meiotic S stage to Leptonema were isolated with a combination of DAPI, SCP3, DMRT1 and STRA8.Nuclei from the leptonema stage to diplonema were isolated with a combination of DAPI, SCP3, H1T1 and SCP1.
(B-C) Strategies for isolating nuclei from unDiff.SGA to the Meiotic G1 stage.
(B) Nuclei from unDiff.SGA to Meiotic G1 stage were isolated from the 2C population.
(C) A snapshot from flow cytometry illustrating the gating for each population.Undifferentiated spermatogonia (unDiff SGA) were identified as 2C nuclei expressing the PLZF (promyelocytic leukemia zinc finger) protein (Buaas et al. 2004;Costoya et al. 2004).Other "pre-meiotic" populations were identified using a combination of DMRT1, which suppresses meiotic entry and is expressed throughout spermatogonial development (Matson et al. 2010), and STRA8, which can trigger meiotic entry but is expressed both in spermatogonia (Zhou et al. 2008;Endo et al. 2015) and in cells that enter MPI (Endo et al. 2015).We identified three populations of nuclei along the continuum of expression of these proteins to represent a putative commitment trajectory of mitotic germ cells as they entered meiosis.The differentiating spermatogonia (Diff.SGA) were collected as nuclei that expressed high levels of DMRT1 and no STRA8.The mitotic-to-meiotic transition nuclei were identified as a population with reduced DMRT1 and slightly increased STRA8.To further examine the transition between mitotic and meiotic stages, the transition population of nuclei was separated into Transition I and Transition II stages based on the DMRT1 signal.Finally, the nuclei with increased STRA8 and reduced expression of DMRT1 were collected from the 2C population as Meiotic G1.PLZF+ = unDiff.SGA; P2= Diff.SGA; P3= Transition I; P4= Transition II; P5= MeioticG1.The reproducibility was assessed by HiCRep (Yang et al. 2017;Lin, Sanders, and Noble 2021).The stratum-adjusted correlation coefficient (SCC) was calculated with a smooth parameter h=11, bin size=25,000 and maximal genomic distance to include in the calculation dBPMax = 500000.

Figure S3. Contact probability decay at single nucleotide resolution within 5 Kb
The read pairs at different distances were counted at a resolution of 1 bp.The contact probabilities were calculated by dividing the pair number at each distance by total pairs.The reads were shifted to the nucleosome dyad by adding 73 bp to the forward reads and subtracting 73 bp from the reverse reads.(A) Compartmentalization was identified by spiking in interactions from the meiotic S stage.Various percentages of read pairs from the meiotic S stage were integrated with read pairs from the pachytene stage.Subsequently, matrices with a resolution of 100 Kb were constructed and balanced.The matrices for chromosome 2, which include different percentages of spike-in, were then plotted.Additionally, the first principal component (PC1) was calculated from the balanced matrices.
(B) Compartmentalization was identified by spiking in interactions from the leptotene stage.Various percentages of read pairs from the leptotene stage were integrated with read pairs from the pachytene stage.Subsequently, matrices with a resolution of 100 Kb were constructed and balanced.The matrices for chromosome 2, which include different percentages of spike-in, were then plotted.Additionally, the first principal component (PC1) was calculated from the balanced matrices.

Figure S8 Aggregated interactions of conserved TADs
The TAD boundaries were identified from both unDiff.SGA and mouse ES cells (Yan et al. 2018) using Cooltools at a resolution of 10 Kb with a window size of 100 Kb.In order to pinpoint conserved boundaries, we initially extended 10 Kb on both sides of the boundaries identified from unDiff.SGA.Boundaries identified in unDiff.SGA that overlapped with boundaries identified in mES cells were classified as conserved boundaries.TAD domains in unDiff.SGA were defined as regions bordered by the two nearest boundaries.Only the TADs in unDiff.SGA defined by the conserved boundaries were used for this analysis.(B) The percentage of dots with both anchors, one anchor, and no anchor overlapping with CTCF binding sites.CTCF motifs were identified from CTCF chip-seq peaks by FIMO.The motifs were extended +/-5 kb when overlapping the anchors; (C) Distance between anchors of the identified dots; (D) Percentage of dots with different oriented anchors.Anchors that overlap with a single CTCF motif were selected for this analysis.(B) Aggregated promoter-enhancer interactions were plotted for both Micro-C and Hi-C data.The enhancers were derived from H3K27ac Chip-seq data (Lam et al. 2019 ;Maezawa et al. 2020).Common peaks from both studies were selected and only the sites that did not overlap with H3K4me3 peaks and hotspots were used for the calculation.Promoters and enhancers at distances ranging from 5 Kb to 5 Mb were paired and aggregated at a resolution of 1000 bp with a flanking region of 50 Kb.Peaks from these datasets that do not overlap with H3K4me3 peaks were considered as enhancers.Furthermore, both promoters and enhancers that coincided with CTCF binding sites within a range of ± 5 Kb were excluded from the calculations.Two promoters, or promoters and enhancers, at distances ranging from 5 Kb to 5 Mb were paired and aggregated at a resolution of 1000 bp with a flanking region of 50 Kb.

(
D) Representative images of isolated nuclei from unDiff.SGA to the Meiotic G1 stage.

Figure S4 .
Figure S4.Genome-wide chromatin reorganization derived from Micro-C

Figure S5 .
Figure S5.Contact probability plotted as a function of contact distances and their corresponding derivatives

Figure S6 .
Figure S6.PC1 profiles derived from biological replicates of both Micro-C and Hi-C data (A) PC1 along chromosome 2 for each replicate of Hi-C data.Principal component analysis was performed at a resolution of 100 Kb.

Figure S7 .
Figure S7.Identification of compartmentalization by Spike-in.

Figure S9 .
Figure S9.Changes in factors that regulate chromatin folding at the scale of TAD (A) Changes in RAD21.The populations ranging from Diff.SGA to meiotic G1 were gated as shown in Figure S1.The meiotic S population was gated based on a 2-4C DNA content, very low DMRT1 expression, and very strong STRA8 expression.The profiles of RAD21 immunostaining signals from the different gated stages were plotted for two biological replicates.

Figure S10 .
Figure S10.Western blot analysis of REC8 and RAD21.(A)Western blot analysis of mitotic kleisin RAD21 and meiotic kleisin REC8 in different "pre-meiotic" stages.TBP (TATA box binding protein) was used as the loading control.

Figure S11 .
Figure S11.Dots visualization in Hi-C/Micro-C matrices generated from ours and published data.

(
Figure S12.A summary of all dots identified from each stage(A) The number of dots identified from each cell type.The anchors that are less than 50 Kb were removed because the calling at this range was not accurate enough;

Figure S13 .
Figure S13.Aggregated interactions of pairwise CTCF from different orientations.

Figure S15 .
Figure S15.The schematic illustrates the enriched interaction pattern at various distances between CTCF binding sites and cohesin processivity.

Figure S16 .
Figure S16.A zoomed-in view of Figure 4F illustrates the depletion of short-range interaction enrichment.Aggregation of interactions surrounding CTCF binding sites with unified motif orientation.The analysis was performed at a resolution of 1 Kb, with a 150 Kb flanking region.Only the interactions at loci with forward-oriented CTCF was shown.

Figure S17 .
Figure S17.Snapshots of Micro-C matrices of at a resolution of 400 bp Micro-C matrices from meiotic S through MPI were plotted for the region spanning 31 Mb to 31.4 Mb on chromosome 5 at a resolution of 400 bp.The Micro-C matrix of the same region from mES(Hsieh et al. 2020) was plotted as a comparison.Genes in the same region were annotated at the bottom.Bars with arrows indicate the orientation of the genes.The difference of visibility mainly results from sequencing depth.

Figure S18 .
Figure S18.Aggregated interactions surrounding TSSs.(A)Aggregated interactions centered around TSSs were calculated for both Micro-C and Hi-C data.TSSs were derived from USCS annotated genes.The aggregated interactions were calculated at a resolution of 1000 bp with a flanking region of 50 Kb.

Figure S19 .
Figure S19.Aggregated interactions between regulatory elements(A) Aggregated promoter-promoter interactions were plotted for both Micro-C and Hi-C data.The promoters were derived from H3K4me3 Chip-seq data(Lam et al. 2019).The H3K4me3 modification sites that overlapped with hotspots were excluded from calculations.Two promoters at distances ranging from 5 Kb to 5 Mb were paired and aggregated at a resolution of 1000 bp with a flanking region of 50 Kb.

Figure S20 .
Figure S20.Aggregation of P-P/P-E interactions that do not overlap with CTCFSimilar to Figure5and FigureS19, promoters were identified based on H3K4me3 ChIP-seq data(Lam et al. 2019).Enhancers were defined as the common peaks found in two separate H3K27ac ChIP-seq datasets(Lam et al.2019; Maezawa et al. 2020).Peaks from these datasets that do not overlap with H3K4me3 peaks were considered as enhancers.Furthermore, both promoters and enhancers that coincided with CTCF binding sites within a range of ± 5 Kb were excluded from the calculations.Two promoters, or promoters and enhancers, at distances ranging from 5 Kb to 5 Mb were paired and aggregated at a resolution of 1000 bp with a flanking region of 50 Kb.