Article Three-Dimensional Genomic Structure and Cohesin Occupancy Correlate with Transcriptional Activity during Spermatogenesis Graphical Abstract Authors Covadonga Vara, Andreu Paytuvı́-Gallart, Yasmina Cuartero, ..., Paul D. Waters, Marc A. Marti-Renom, Aurora Ruiz-Herrera Correspondence martirenom@cnag.crg.eu (M.A.M.-R.), aurora.ruizherrera@uab.cat (A.R.-H.) In Brief The formation of mammalian germ cells involves dramatic chromosomal movements and chromatin remodeling, whose regulatory pathways are far from understood. Vara et al. show how the dynamics of insulator proteins’ occupancy and transcriptional activity are coupled during the 3D genome re- organization that takes place during mouse spermatogenesis. Highlights d Meiotic and post-meiotic cells differ in their hierarchal genome organization d Meiotic cohesin loading correlates with gene expression and local insulation d Sperm cells contain compartments with active transcription but no clear TADs d Genes that escape X chromosome silencing associate with promoter cohesin occupancy Vara et al., 2019, Cell Reports 28, 352–367 July 9, 2019 ª 2019 The Author(s). https://doi.org/10.1016/j.celrep.2019.06.037 Cell Reports Article Three-Dimensional Genomic Structure and Cohesin Occupancy Correlate with Transcriptional Activity during Spermatogenesis Covadonga Vara,1,2,12 Andreu Paytuvı́-Gallart,1,2,3,12 Yasmina Cuartero,4,5,12 François Le Dily,4,5 Francisca Garcia,6 Judit Salvà-Castro,1,2 Laura Gómez-H,7 Eva Julià,4,8 Catia Moutinho,5 Riccardo Aiese Cigliano,3 Walter Sanseverino,3 Oscar Fornas,4,9 Alberto M. Pendás,7 Holger Heyn,5,9 Paul D. Waters,10 Marc A. Marti-Renom,4,5,9,11,* and Aurora Ruiz-Herrera1,2,13,* 1Departament de Biologia Cel.lular, Fisiologia i Immunologia, Universitat Autònoma de Barcelona (UAB), Cerdanyola del Vallès 08193, Spain 2Genome Integrity and Instability Group, Institut de Biotecnologia i Biomedicina (IBB), Universitat Autònoma de Barcelona (UAB), Cerdanyola del Vallès 08193, Spain 3Sequentia Biotech, Carrer Comte D’Urgell 240, Barcelona 08036, Spain 4Centre for Genomic Regulation (CRG), The Barcelona Institute for Science and Technology (BIST), Carrer del Doctor Aiguader 88, Barcelona 08003, Spain 5CNAG-CRG, Centre for Genomic Regulation, The Barcelona Institute of Science and Technology (BIST), Baldiri Reixac 4, Barcelona 08028, Spain 6Unitat de Cultius Cel.lulars, Universitat Autònoma de Barcelona (UAB), Cerdanyola del Vallès 08193, Spain 7Molecular Mechanisms Program, Centro de Investigación del Cáncer and Instituto de Biologı́a Molecular y Celular del Cáncer (CSIC-Universidad de Salamanca), Salamanca 37007, Spain 8Institut Hospital del Mar d’Investigacions Mèdiques (IMIM), Carrer del Doctor Aiguader 88, PRBB Building, Barcelona 08003, Spain 9Pompeu Fabra University (UPF), Doctor Aiguader 88, Barcelona 08003, Spain 10School of Biotechnology and Biomolecular Sciences, Faculty of Science, UNSW Sydney, NSW 2052, Australia 11ICREA, Pg. Lluı́s Companys 23, Barcelona 08010, Spain 12These authors contributed equally 13Lead Contact *Correspondence: martirenom@cnag.crg.eu (M.A.M.-R.), aurora.ruizherrera@uab.cat (A.R.-H.) https://doi.org/10.1016/j.celrep.2019.06.037 SUMMARY with three-dimensional (3D) remodeling during sper- matogenesis, ultimately reprogramming the genome Mammalian gametogenesis involves dramatic and for the next generation. tightly regulated chromatin remodeling, whose regu- latory pathways remain largely unexplored. Here, we generate a comprehensive high-resolution structural INTRODUCTION and functional atlas of mouse spermatogenesis by combining in situ chromosome conformation capture Mammalian genomes are packaged into a tailored chromatin structure, the regulation of which depends on several superim- sequencing (Hi-C), RNA sequencing (RNA-seq), and posed layers of organization, including epigenetic modifications chromatin immunoprecipitation sequencing (ChIP- (of both the DNA and nucleosomes) and the higher-order organi- seq) of CCCTC-binding factor (CTCF) and meiotic co- zation of chromatin compartments inside the nucleus. This orga- hesins, coupled with confocal and super-resolution nization is achieved by chromatins folding into loops, topologi- microscopy. Spermatogonia presents well-defined cally associating domains (TADs), and compartments (A and compartment patterns and topological domains. B), which can ultimately influence transcriptional activity (Dixon However, chromosome occupancy and compart- et al., 2012; Lieberman-Aiden et al., 2009; Rao et al., 2014). mentalization are highly re-arranged during pro- How these different levels of chromatin organization interact dur- phase I, with cohesins bound to active promoters in ing the cell cycle has just begun to be elucidated (Dekker et al., DNA loops out of the chromosomal axes. Compart- 2013). In somatic cells, the highly compartmentalized folding of ment patterns re-emerge in round spermatids, where the genome in interphase is lost during mitosis, when chromo- somes are linearly organized in consecutive chromatin loops cohesin occupancy correlates with transcriptional (Gibcus et al., 2018; Naumova et al., 2013). Recent studies in activity of key developmental genes. The compact mice have suggested remarkable chromatin architecture reprog- sperm genome contains compartments with actively ramming during the formation of germ cells (Alavattam et al., transcribed genes but no fine-scale topological do- 2019; Patel et al., 2019; Wang et al., 2019) and early develop- mains, concomitant with the presence of protamines. ment (Du et al., 2017; Flyamer et al., 2017). However, how the Overall, we demonstrate how genome-wide cohesin higher-order chromatin organization is configured during all occupancy and transcriptional activity is associated stages of spermatogenesis, and how insulator proteins and 352 Cell Reports 28, 352–367, July 9, 2019 ª 2019 The Author(s). This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). cohesins determine this organization to regulate transcription sequencing (RNA-seq) and chromatin immunoprecipitation activity, remains largely unexplored. sequencing (ChIP-seq) of CCCTC-binding factor (CTCF) and Germ cells represent a unique cell model, where unipotent meiotic cohesins (Figure 1A). These data have permitted the diploid cells (gonia) undergo extensive cellular differentiation comprehensive study of the close interplay between chromatin (gametogenesis) to form highly differentiated cells (oocytes higher-order organization dynamics and function during mouse and sperm) that ultimately form a totipotent embryo. In the spermatogenesis. case of mammalian males, this complex process is divided into three stages: (1) proliferation and differentiation of sper- RESULTS matogonia (Spg); (2) meiosis, a reductional division that pro- duces haploid cells through two consecutive cell divisions Dynamic Overall Chromatin Structure Reorganization (meiosis I and meiosis II); and (3) spermiogenesis, where round during Spermatogenesis spermatids (RSs) are transformed into densely compacted sper- To unveil changes in chromosome conformation during sper- matozoa. These sequential developmental stages involve dra- matogenesis, we developed a reproducible fluorescence-acti- matic and tightly regulated chromosomal re-organization and vated cell sorting (FACS) protocol to obtain, based on DNA chromatin remodeling. It is during the first meiotic prophase (pro- content and chromatin complexity, highly enriched (90.4% phase I) that homologous chromosomes align, pair, synapse, average enrichment) cell fractions for Spg, primary spermato- and recombine. All these processes are interconnected through cytes at the leptonema-zygonema (L/Z) and pachynema-diplo- four sequential stages: leptonema, zygonema, pachynema, and nema (P/D) stages, RSs, and sperm (Figures 1B, 1C, and S1) diplonema (Handel and Schimenti, 2010). (STAR Methods). For each germ cell fraction, as well as for a At leptonema, chromosomes cluster by their telomeres to the mouse primary fibroblast cell line as a somatic profile, we per- nuclear envelope in the bouquet (Reig-Viader et al., 2016). This formed in situ Hi-C (Rao et al., 2014) (Figures 1D1F and S1). structure promotes the pairing of homologous chromosomes After filtering the raw Hi-C interactions, an average of 237.86 by the formation of proteinaceous structures along chromo- million valid interactions were obtained per cell type (Tables somes formed by cohesins (i.e., REC8 and RAD21L; Gutiérrez- S1 and S2). The comparison between biological replicates re- Caballero et al., 2011; Llano et al., 2012) and proteins of the sulted in highly reproducible Hi-C maps (Figure S2). synaptonemal complex (SC). Meiotic recombination is triggered Genome organization changed during spermatogenesis (Fig- by the formation of double strand breaks (DSBs), caused by the ures 1F and 2A2F), as reflected by the analysis of distance- endonuclease protein SPO11 (Keeney et al., 1997). DSBs dependent interaction frequencies (Figure 2C) and inter- and are then initiated at zygonema, leading to synapsis between intra-chromosomal interaction ratios (Figure 2D). Fibroblasts homologous chromosomes. Subsequently, at pachynema, SCs showed high inter- and intra-chromosomal interaction ratios are completely established, creating bivalent structures with (>0.6) inversely correlated with chromosomal size (p < 0.001), resolved recombination producing crossover events. At diplo- suggesting distinct chromosomal compartmentalization within nema, homologous chromosomes start to segregate by the nuclei (Figure 2D). In contrast, inter- and intra-chromosomal disassembly of SCs to produce spermatocytes II, which undergo interaction ratios decreased 2-fold to about 0.3 for all chromo- a second meiotic division resulting in RSs. Finally, spermatids somes in Spg, suggesting that a commitment to enter meiosis become sperm ready for fecundation through spermiogenesis, is accompanied by a drastic remodeling of chromosomal terri- a differentiation stage that includes changes in cell morphology tories within the nucleus. Importantly, this decrease in inter- and DNA packaging via the replacement of histones by prot- and intra-chromosomal interactions was concomitant with the amines (testis-specific histone variants). dynamic changes of the so-called A-B compartments (Lieber- Despite recent analysis of genome conformational changes in man-Aiden et al., 2009). Likewise, fibroblasts and Spg, both in male germ cells (Alavattam et al., 2019; Patel et al., 2019; Wang interphase, shared similar contact probability patterns at short et al., 2019), how the higher hierarchal level of genome organiza- distances (from 0.5 to 7 Mbp) and had the lowest intermediate tion is related to gene expression and insulator proteins during interactions (between 1 and 10 Mb) (Figure 2C). However, at spermatogenesis remains unknown. It has been generally larger genomic distances (>10 Mb), fibroblasts showed a slight accepted that there are two waves of active transcription: one change in the slope at 10 Mbp with fewer interactions, whereas before entering meiosis and a second in the transition from Spg maintained the same trend up to 100 Mbp. RSs to sperm (Sassone-Corsi, 2002). However, recent transcrip- As meiosis progressed, compartments were mostly lost in tome analysis in germ cells suggests that these two transcrip- primary spermatocytes (L/Z and P/D), coinciding with pro- tional waves might take place earlier (da Cruz et al., 2016), which phase I, when homologous chromosomes condensate, align, evidences a finely tuned regulation of chromatin remodeling and pair, synapse, and recombine (Figures 2B, 3A, and 3B). Consis- active transcription. Here, we implement a reproducible flow tent with this absence of compartments during prophase I, cytometry protocol to isolate enriched male mouse germ cell eigenvector values were close to 0 (Figures 2B, 3A, and 3B), populations representing all stages of spermatogenesis: pre- and inter- and intra-chromosomal interaction ratios reached a meiotic (Spg), meiotic (leptonema, zygonema, pachynema, and minimum for all chromosomes (Figure 2D). An exception was diplonema), and post-meiotic cells (RSs and sperm). On these the sex chromosomes in P/D, which had detectable variations sorted germ cell populations, we performed genome-wide chro- in the inter- and intra-chromosomal interactions ratio. This un- mosome conformation capture analysis (in situ Hi-C [chromo- usual pattern most likely reflects meiotic sex chromosome inac- some conformation capture sequencing]), coupled with RNA tivation (MSCI), the process by which the X and Y chromosomes Cell Reports 28, 352–367, July 9, 2019 353 Figure 1. Higher-Order Chromatin Structure during Spermatogenesis (A) Experimental workflow: Hi-C, ChIP-seq, and RNA-seq on fluorescence-activated cell sorting (FACS)-enriched mouse germ cells. (B) Flow cytometry Hoechst Blue (UV355-460/50) and Hoechst Red (UV355-670/30) plot showing spermatogonia (Spg), leptonema-zygonema (L/Z), pachynema- diplonema (P/D), secondary spermatocytes (2c), round spermatids (RSs), and sperm. Recovered germ cell fractions presented the following average enrichment: 91% for Spg, 88.7% for L/Z, 89.2% for P/D, 92.9% for RS, and 90% for sperm. (C) Histogram of the differential DNA content showing cell events for each FACS-isolated germ cell population (Hoechst Red, UV355-670/30). (D) Overview of the spermatogenesis process (adapted from Reig-Viader et al., 2016). Numbers between parentheses indicate the diploid (2n) haploid (n) number for each cell type and the number of chromatids per chromosome (4c, 2c, or c). (E) Representative immunofluorescence images showing DAPI-stained DNA (gray/blue) and specific meiotic proteins for the different cell populations included in the present study. Fibroblasts and Spg have DAPI-stained DNA in gray. For L/Z and P/D, DAPI is shown in blue, SYCP3 in green, and gH2AX in red. The image represents a mosaic of two individual captured cells. In RSs, DAPI is blue, and H3K9me3 (marker for the constitutive heterochromatin at centromeres- chro- mocenters) is red. Scale bars, 10 um. (F) Iterative correction and eigenvector decomposition (ICE)-corrected Hi-C matrices for chromosome 15 at a 50-kbp resolution for the cell types analyzed. Deep blue lines indicate non-mapped bins. See also Figure S1. are transcriptionally inactivated in primary spermatocytes, primary spermatocytes, being more prominent in L/Z (Figures forming the sex body at the periphery of the nucleus (Turner, 2E and 2F). Both L/Z and P/D had an enrichment of counts 2007). Another structural feature characteristic of early meiosis at the most proximal sub-telomeric bins, but cells at L/Z had that was revealed in our analyses was the bouquet structure. We significantly higher interactions (Wilcoxon test, p < 0.05) than detected high inter-chromosomal contact between telomeres in any other cell type (Figures 2E and 2F). The analysis of the 354 Cell Reports 28, 352–367, July 9, 2019 (legend on next page) Cell Reports 28, 352–367, July 9, 2019 355 distance-dependent interaction frequencies revealed a distinct Two Rounds of TAD Reorganization in Primary chromosome organization for primary spermatocytes, with two Spermatocytes and Sperm abrupt changes in slope: the first at 2.54.5 Mbp and the sec- To further investigate the dynamics of the higher-order chro- ond at 40 Mbp (Figure 2C). matin structure at the sub-megabase scale, we identified After meiosis, haploid cells (RSs and sperm) had a distinctive TADs and examined the robustness of their boundaries using higher-order chromatin structure. Although A-B compartments TADbit (Serra et al., 2017) at 50-Kbp resolution. Similar to the re-appeared after being lost in previous stages, they presented A-B compartment patterns, TADs were well defined in both as a blurry plaid pattern of larger mean size (0.86 Mbp in RSs fibroblasts and Spg (Figure 3C). In primary spermatocytes, and 0.93 Mbp in sperm) compared to Spg (Figures 2F, 3B, S3A, however, there was a substantial reduction in the variance of and S3B; Table S3). Interestingly, the proportion of genomic TAD insulation score as well as an increase in TAD size, espe- bins with the same compartment status (A or B) was higher be- cially in L/Z (Figures 3C3E). TAD insulation scores were tween haploid cells (r2 = 0.80) than between haploid cells and partially recovered in RSs but, in contrast to previous observa- Spg (RSs versus Spg, r2 = 0.58; sperm vs Spg, r2 = 0.48). This tions (Jung et al., 2017; Wang et al., 2019), were drastically pattern of differential compartmentalization in RSs was also vali- reduced in sperm cells (close to 0) (Figures 3D, 3E, and S5). dated by three-dimensional fluorescence in situ hybridization A total of 2002 TADs, with an average length of 1.3 Mbp, (3D-FISH), where physical distances between pairs of loci were identified in fibroblasts, which was more than in Spg increased compared to Spg and fibroblasts (Figure S4; Table (834 TADs, mean size of 3.26 Mbp) (Figure 3C; Table S3). S4). Moreover, the inter- and intra-chromosomal interaction ratio Although fewer TADs were identified in primary spermato- values were higher in RSs for all chromosomes, with the excep- cytes (305 TADs in L/Z and 294 TADs in P/D), their bound- tion of sex chromosomes (Figure 2D), which are known to form aries had high strength scores (74.25% of TADs in L/Z, post-meiotic sex chromatin (PMSC) attached to the chromocen- and 79.59% in P/D, had scores >9; Table S3). This pattern ters (Namekawa et al., 2006). Thus, the formation of PMSC corre- contrasted with RSs, which had a large number of small lates with chromatin remodeling that results in low inter- and TADs (n = 4,649) with low border strength scores (Figure 3C; intra-chromosomal interactions. Remarkably, in sperm, inter- Table S3). Meta-border plots confirmed this dynamic (Fig- chromosomal interactions were greater than intra-chromosomal ure 3F). These results demonstrate that TADs also underwent interactions, compared to previous stages, and inversely corre- a large reorganization in genome structure during spermato- lated with chromosomal size (p < 0.001) (Figure 2D). Since ratio genesis, concomitant to larger changes in chromosome terri- values were higher than in fibroblasts, the higher-order chromatin tories and compartments. structure is likely densely packed in sperm but remains in chro- In summary, our data suggest at least four distinct patterns of mosome territories. The sex chromosomes did not follow the chromatin interactions during spermatogenesis progression: (1) autosomal pattern; the interactions ratio decreased by 2-fold interphase-like organization (e.g., Spg); (2) a condensed pattern (Figure 2D). In both haploid cell types, we detected higher inter- in prophase I (i.e., L/Z and P/D) where A-B compartments and chromosomal contact between centromeres, suggestive of the TADs are largely lost; (3) RSs with a blurry compartment plaid presence of the chromocenters (Figures 2E and 2F). In fact, it is pattern; and (4) sperm where TADs are reduced, but A-B com- known that in mouse sperm, centromeres are located at the partments are observed. center of the nucleus, while telomeres attach to the nuclear enve- lope (Haaf and Ward, 1995). This particular compartmentalization Functional Compartment Switching during in both cell types was also reflected by a decrease in interaction Spermatogenesis frequencies as genomic distance increased (Figure 2C). To assess whether A-B compartmentalization changes during These observations indicate that the genome suffers a major spermatogenesis correlated with differential expression of resi- structural re-organization during spermatogenesis with dy- dent genes, we investigated changes of compartment type be- namic and dramatic changes in chromosome occupancy and tween the stages where compartments could be clearly compartments. observed. The proportion of the genome organized in the Figure 2. Chromosomal Organization in Interphase, Pre-meiotic, Meiotic, and Post-meiotic Cells (A) Genome-wide ICE-corrected heatmaps at 500 kbp for the cell types analyzed. (B) Chromosome 18 region-specific ICE-corrected heatmaps at 50 kbp (from 20 Mbp to 60 Mbp), depicting compartment signal (first eigenvector) for all cell types. (C) Contact probability P(s) as a function of genomic distance in all cell types for autosomes (left panel) and the X chromosome (right panel). Discontinuous straight lines correspond to the fractal (green) and equilibrium (red) models (Lieberman-Aiden et al., 2009). Gray-shadowed area expands the genomic region from 0.5 to 7 Mbp. (D) Inter- and intra-chromosomal interaction ratios for each chromosome and cell type. Correlation values (**p < 0.001) between the inter- and intra-chromosomal ratio and chromosomal size (autosomes only) are shown for each cell type. (E) Heatmaps showing normalized inter-chromosomal interactions between chromosomes 1 and 2 in all cell types. Red circles represent high-contact regions. (F) Left panel: boxplots depicting inter-chromosomal interactions per million at sub-centromeric regions (from the centromere up to 3.5 Mbp) for all cell types (Wilcoxon test, **** < 0.0001; n.s., not statistically significant when compared to fibroblasts). Right panel: schematic representation of chromosomes and centromeres and telomeres in L/Z, P/D, RSs, and sperm. Dots represent centromeres (pink) and telomeres (green). In RSs and sperm, all centromeres associate in the center of the cell forming the chromocenter. Fib, fibroblast; Spg, spermatogonia; L/Z, leptonema-zygonema; P/D, pachynema-diplonema; RS, round spermatids; cen, centromeres of acrocentric chromosomes; tel, telomeres of non-centromeric ends. See also Figures S2–S5. 356 Cell Reports 28, 352–367, July 9, 2019 Figure 3. Chromosome-specific A-B Compartment Profiles and TAD Signals (A) Compartment signal (first eigenvector) across chromosome 18. (B) Density plots of eigenvector values considering autosomes. (C) TAD border alignments along chromosome 18. Dark gray arches represent TADs with higher intra-TAD interactions than expected. TAD border robustness (from 1 to 10) is represented by a color gradient. (D) Representation of TAD insulator score in mouse chromosome 18. (E) Variance of the TAD insulation scores for autosomes (left panel) and the X chromosome (right panel). (F) Meta-plots for all TAD boundaries detected in fibroblasts (n = 2002), Spg (n = 834), L/Z (n = 305), P/D (n = 294), RSs (n = 4649), and sperm (n = 1042). Fib, fibroblast; Spg, spermatogonia; L/Z, leptonema-zygonema; P/D, pachynema-diplonema; SpII, spermatocytes II; RS, round spermatids. See also Figures S3–S5. Cell Reports 28, 352–367, July 9, 2019 357 Figure 4. Dynamics of Gene Expression during Gametogenesis (A) Two extensive waves of transcription take place during spermatogenesis. (B) Differentially expressed genes (DEGs) for each pairwise comparison, including known nuclear protein-coding RNA (RefSeq), long noncoding RNA, antisense RNA, processed transcripts, and mitochondrial nuclear protein-coding RNA. (C) Boxplots representing genome-wide expression (as fragments per kilobase of transcript per million mapped reads [FPKM] values) according to A-B compartment assignment (N/A, not assigned compartments). (D) Gene Ontology (GO) analysis of expressed genes in cell-specific A compartments. Bubble size represents the number of genes assigned to each GO. Only GO terms with two or more genes are represented. Spg, spermatogonia; P/D, pachynema-diplonema; RS, round spermatids. A compartment was 45.7% in fibroblast and reduced to 39.4% in gene expression analysis between cell types generated lists of Spg, before a rise in RSs (46.9%) and sperm (48.6%). Since A differentially expressed genes (DEGs), which were classified as compartments correlate with open chromatin and active genes (1) protein-coding RNA; (2) long non-coding RNA (lncRNA); (3) (Lieberman-Aiden et al., 2009), the switching of compartment antisense RNA (asRNA); (4) processed and unprocessed pseu- types can provide insight into genome function during spermato- dogenes; or (5) unannotated and unconfirmed transcripts genesis. We therefore performed low-input RNA-seq (four bio- (Figure 4; Table S6). The majority (76.55%) of DEGs were pro- logical replicates) on a selected group of the enriched germ tein-coding, with their abundance reducing as spermatogenesis cell populations: Spg, primary spermatocytes (P/D), RSs, and progressed. This was coupled with an increased expression of sperm. After filtering, an average of 35.37 million paired-end lncRNAs and non-coding asRNAs and pseudogenes (Figure 4; reads was obtained per cell type (Table S5). Table S6). Consistent with the global shutdown in gene expres- The total number of expressed genes decreased as spermato- sion that occurs during meiosis (Sassone-Corsi, 2002), the net genesis progressed, with 19,145 expressed in Spg, 15,480 in balance between pairwise comparisons (Spg versus P/D, RSs P/D, 14,706 in RSs, and 13,646 in sperm. Pairwise differential versus sperm, and Spg versus sperm) was negative, which 358 Cell Reports 28, 352–367, July 9, 2019 suggests higher transcriptional activity in RSs than in primary between cohesins and DSBs (SPO11-oligos hotspots; Lange spermatocytes (Figures 4A and 4B). et al., 2016) and H3K4me3 (Brick et al., 2012). Interestingly, Consistent with a correlation between chromatin remodeling both SPO11 hotspots and H3K4me3 marks correlated with co- and active transcription, genes in the A compartments had hesin occupancy (STAR Methods; p < 0.001), concomitant with significantly higher expression than those in the B compartments open chromatin states. Since DSBs are known to occur (Figures 4C and 4D). As expected, genes related to spermato- genome-wide before being recruited at the chromosomal axes genesis (e.g., morphogenesis and cell differentiation) were en- to be repaired (Lange et al., 2016), the observed correlation be- riched in the A compartments (Figure 4D). In RSs, during the sec- tween SPO11 hotspots and cohesion peaks also suggests the ond wave of transcription (Figure 4A), the most representative presence of cohesins in DNA loops out of the axes. Gene Ontology (GO) term was ‘‘system process,’’ which Remarkably, we detected a correlation of cohesin genomic included 27 genes involved in sensory perception, including distribution with gene expression and local insulation in both olfactory receptors. Likewise, genes with important roles in fertil- primary spermatocytes and RSs (Figures 5 and 6). In P/D, ization were also expressed in RSs, including the acrosome most cohesin peaks (80.7% of RAD21L and 83.3% of REC8) reaction (e.g., Plcz1 and Smcp). Therefore, the transformation were located in promoter regions (less than 2 kbp from the tran- of RSs into spermatozoa was accompanied by the transcription scriptional start site [TSS]) of genes with significantly (p < 0.01) of genes related to spermiogenesis and sperm function located higher expression than genes without promoter-associated in newly created A compartments. cohesin peaks (Figures 5C5E). In RSs, cohesin peaks associ- ated with TSSs were reduced to 77.7% for REC8 and 45.5% Differential CTCF and Cohesin Loading Correlates with for RAD21L, suggesting an unequal re-distribution of meiotic co- Gene Expression and Chromatin Remodeling hesins later in spermatogenesis (Figures 5C and 5D). Regardless Western blots confirmed the presence of CTCF and meiotic of RAD21L re-distribution, the overlapping of genes with cohesins (REC8 and RAD21L) in whole-testis protein extracts. RAD21L in their promoters between P/D and RSs was statisti- Immunofluorescence (IF) then revealed a previously unreported cally significant (p = 0e+00, Fischer’s exact test). Despite this pattern for CTCF, with signal along all chromosome axes in pri- reduction of peaks in RSs, genes with cohesin peaks in their mary spermatocytes that were more intense on autosomals than TSSs still had significantly higher expression than genes without on the X (Figure S6). There was a weak CTCF signal in RSs and a peaks (Figures 5E and 5F). In contrast, only 25.9% of the CTCF clear, cloud-like signal for cohesins (Figure S6). ChIP-seq anal- peaks were within promoter regions (i.e., 2 kbp upstream of ysis confirmed these differences in IF signal intensity and density TSSs) in primary spermatocytes, with the majority located in in- (Figure S6). In fact, we detected 19,347 CTCF ChIP-seq peaks in tergenic regions (Figure 5C). In all cases, CTCF and cohesin primary spermatocytes (P/D), with the vast majority (97.1%) peaks were preferentially located (nearly 80% of peaks) in A being lost in RSs (Figures 5A and 5B). compartments (permutation test, STAR Methods; p < 0.001), We then performed ChIP-seq in P/D and RSs on the meiotic consistent with a correlation between chromatin remodeling cohesins REC8 and RAD21L. In primary spermatocytes (P/D), and active transcription. Although the TAD insulation score we detected 11,618 REC8 and 9639 RAD21L peaks distributed was reduced in primary spermatocytes (Figures 3D and 3E), across the genome, with substantial overlap (55.8% for some TADs still remained, with boundaries enriched for CTCF RAD21L and 46.3% for REC8) (Figures 5A and 5B). Of the total and cohesins (STAR Methods; p < 0.01) (i.e., genomic regions peaks, 3999 were common to all CTCF (20.6%), RAD21L with lower TAD insulation scores; Figure 6B) and associated (41.5%), and REC8 (34.4%). In RSs, we found 11,559 REC8 with gene expression (Figure 6C). Meta-border plots of CTCF peaks and 12,507 RAD21L peaks, of which 6383 (51% of and cohesin peaks confirmed this trend (Figure 6D). The same RAD21L and 55.2% of REC8) overlapped (Figures 5A and 5B). pattern was observed in RSs (Figures 6B6D), suggesting that The vast majority (90.5%) of REC8 peaks observed in primary despite global chromatin remodeling, local insulation is main- spermatocytes were maintained in RSs, which contrasted tained by insulator proteins and can affect gene expression. RAD21L, where 40% of the peaks were specific to RSs. Expressed genes with CTCF, RAD21L, and REC8 in their pro- Despite the overlap of RAD21L and REC8 peaks detected moters in P/D had GO term enrichments related to protein regu- genome-wide in P/D, the distribution of REC8 and RAD21L im- lation, modification, and polymerization (e.g., Usp42), as well as munolabeling along the axes was not continuous or significantly DNA repair and cellular response to DNA damage stimulus (e.g., correlated, as measured by super-resolution microscopy (r2 = Herc2), suggesting a regulatory role in spermatogenesis pro- 0.15; Figure S6). The close proximity of meiotic cohesins gression. Interestingly, genes with CTCF at their promoters genome-wide (on average, peaks are scattered every 264.3 were involved in the transcriptional machinery (e.g., Nsa2), Kbp for RAD21L and 219.2 Kbp for REC8) was much lower whereas genes with both cohesins but not CTCF were essential than the estimated DNA loop size at the axes for mouse pachy- for posterior neural development of the embryo (e.g., Cdh2). In nema (1.52 Mbp; Patel et al., 2019). This pattern suggests that RSs, CTCF was almost absent, but remarkably, when it co-local- the peaks correspond not only to cohesins loaded at the chro- ized with promoters, it was associated with genes implicated in mosomal axes (expected every 1.52 Mbp), but also to DNA key pathways of embryo development (e.g., Nanog). Genes with loops out of the axes. Given the close interplay among the orga- both cohesins in their promoters were related to the regulation of nization of chromosomal axes, DNA loops, and the formation of cell growth and the Wnt signaling pathway (e.g., Amer3) and gen- DSBs during early prophase I (Kleckner et al., 2003; Ruiz-Herrera eral nervous system development (e.g., Ccd88a; Figure 5F). et al., 2017; Wang et al., 2015), we examined the correlation Genes with RS-specific RAD21L peaks had functions involved Cell Reports 28, 352–367, July 9, 2019 359 Figure 5. CTCF and Cohesin Profiles in Primary Spermatocytes and RSs (A) Venn diagrams for CTCF and cohesins in P/D and RSs considering peak overlaps per cell type and peak overlaps per protein. (B) Representative examples of CTCF and cohesins’ genomic distribution along chromosome 1 in P/D and RSs. For each cell type, A-B compartments, gene expression (represented as log FPKM), CTCF peaks, and cohesin peaks are displayed. (C) Genome-wide distribution of CTCF and cohesin-occupied sites in relation to TSSs and other genomic features in P/D and RSs. (D) Insulator peak frequencies relative to TSSs of genes in P/D and RSs. (E) Boxplots representing expression (FPKM values) of genes with CTCF and cohesin peaks located at the TSS. Asterisks represent statistically significant differential gene expression when compared with all genes in the mouse genome (p < 0.01). (F) Examples of CTCF and cohesion-occupied sites in P/D and RSs for the expressed gene Cdc88a. See also Figure S6. in oligondendrocyte differentiation and cardiac ventricle devel- Sex Chromosome Silencing Is Coupled with Higher- opment (e.g., Notch1), whereas RS-specific REC8 peaks were Order Chromatin Restructuring found to be involved in nervous system development (e.g., The eutherian mammal X chromosome is composed of evolu- Kif3c). tionary strata that were isolated from recombination with the 360 Cell Reports 28, 352–367, July 9, 2019 Figure 6. Local Insulation, Cohesin Occupancy, and Gene Expression (A) Representative examples of CTCF and cohesins’ genomic distributions across a specific region of chromosome 19 (from 7 to 9.3 Mbp) in P/D and RSs. For each cell type, Hi-C interaction maps (50-Kbp bins), gene expression (represented as log FPKM), CTCF peaks, cohesin peaks, and genes from NCBI Ref Seq annotation are displayed. Green and orange highlights indicate differentially expressed genes (DEGs) in each cell type: Gm17227 and Uqcc3 for P/D and Atl3; Hrasls5, 1700092M07Rik, Eef1g, and Asrgl1 for RSs. (B) Distribution of CTCF and cohesin peaks at TAD borders. The y axes represent the TAD insulation Z-score relative to random genomic regions (based on 10,000 permutation tests with randomization, p < 0.01). (C) Distribution of gene expression for CTCF and cohesin peaks located at TAD borders. The y axes represent the TAD insulation Z-score relative to random genomic regions (based on 10,000 permutation tests with randomization, p < 0.01). (D) Meta-plots for all peaks detected in P/D and RSs. P/D, pachynema-diplonema; RS, round spermatids. See also Figure S6. Y chromosome at different evolutionary times (Lahn and Page, 7A and 7B). In fact, distance-dependent interaction frequencies 1999). In mice, the X chromosome evolutionary strata are highly revealed that genomic interactions at medium genomic dis- rearranged, and there are 22 ampliconic regions known to tances (110 Mbp) were higher in the X chromosome in post- escape MSCI (Mueller et al., 2008) (Figure S7; Table S7). The meiotic cells than in autosomes of the same cell type (Figure 2C). X chromosome lost its A-B compartment pattern once meiosis For higher-order chromatin structures at the sub-megabase was initiated, from L/Z onward (Figures 7A and 7B), consistent scale, all chromosomes except X had equivalent TAD dynamics with MSCI chromatin remodeling. This translated into an overall during spermatogenesis. In the X of meiotic and post-meiotic reduction in gene expression in P/D compared to Spg, although cells, there was a trend for fewer and larger TADs compared to ampliconic regions were still expressed in P/D, RSs, and sperm those observed in Spg (Figure 7D; Table S8). The genomic posi- (Figure 7C). Importantly, this organization of the X chromosome tions of TADs were not equally distributed along the X; in RSs, was maintained in post-meiotic cells (RSs and sperm) (Figures there were more short TADs detected in the distal regions of Cell Reports 28, 352–367, July 9, 2019 361 (legend on next page) 362 Cell Reports 28, 352–367, July 9, 2019 the X (Figure 7D), consistent with the increase of autosomal vari- ization of meiotic chromosomes in both early and late stages ance of TAD insulation scores in RSs compared to L/Z and P/D of spermatogenesis, which was reflected at different levels: (1) (Figure 3E). However, the remainder of the X in RSs maintained inter- and intra-chromosomal interaction ratios; (2) distance- the fewer and larger TADs observed in L/Z and P/D. In sperm, dependent interaction frequencies; (3) genomic compart- the high density of short TADs was only observed in proximal ments; (4) topological domains; and (5) occupancy of insulator X chromosome long arm (Xq). proteins. We provide evidence of a delicate fine-tuning among chromatin remodeling, architectural proteins, and cell-specific Differential CTCF and Cohesin Loading in Autosomes gene expression. and the X Chromosome Our analyses complement and extend recently published The proportion of cohesin peaks was lower for the X chromo- works (Alavattam et al., 2019; Patel et al., 2019; Wang et al., some than for any autosome in both P/D and RSs (Figure 7E). 2019) with a comprehensive view of the sequential develop- The X represents 4.9% of the mouse genome, but only ac- mental stages during mouse spermatogenesis. Pre-meiotic, counts for 1.82% of the total RAD21L peaks in P/D and 1.54% meiotic, and post-meiotic cells all presented differences in inter- of peaks in RSs. Likewise, REC8 peaks on the X chromosome and intra-chromosomal interaction ratios that were distinct from were underrepresented, with 2.12% of the total peaks in P/D somatic cells (i.e., fibroblasts). Although Spg maintained the and 2.03% of peaks in RSs. The underrepresentation of REC8 equivalent proportion of A-B compartments as fibroblasts, they ChIP-seq peaks was confirmed by super-resolution microscopy showed a drastic remodeling of chromosomal compartmentali- in pachynema, with decreased labeling of REC8 (44.7% ± zation. As the first meiotic prophase begins, additional chromatin 10.7%) on the X chromosome axis relative to those of autosomes remodeling appears. The inter- and intra-chromosomal interac- (Figure S6). In contrast, RAD21L labeling actually increased tion ratio reaches a minimum in primary spermatocytes, with (25.6% ± 10.6%) on the X (Figure S6). the A-B compartments almost lost in both L/Z and P/D, in Overall, 75% of cohesin peaks co-localized with promoter re- contrast to previous observations (Patel et al., 2019). Thus, chro- gions (less than 2 kbp upstream of TSSs) in both P/D and RSs in mosomal occupancy and compartmentalization inside nuclei the X chromosome (Figure 7F). These genes had significantly were re-arranged (i.e., higher-chromatin structure relaxation) higher expression than genes without promoter-associated during prophase I, permitting DNA-scaffold assembly and the cohesin peaks, despite MSCI and PMSC (Figures 7G and 7H), formation and repair of DSBs, with no distinction between suggesting that genes that escape MSCI correlate with promoter autosomes. cohesin occupancy. As with the cohesins, CTCF was depleted The distance-dependent interaction we detect suggests dif- on the X chromosome relative to the autosomes (1.26% of total ferences in previously reported mitotic and meiotic chromosome peaks were on the X in P/D). However, CTCF was almost folding (Gibcus et al., 2018; Naumova et al., 2013; Wang et al., completely depleted from the X chromosome in RSs (just two 2019). Prophase I cells display two changes in contact probabil- peaks were detected). This was not observed for individual auto- ity: the first between 2.5 and 4.5 Mbp and the second at 40 Mbp. somes, on which the proportion of CTCF peaks were generally This organization can result from two features of chromosome stable between P/D and RSs. Collectively, these observations assembly during prophase I. First, the chromatin is anchored point to a novel role for meiotic cohesins in genome organization as long DNA loops in a protein scaffold composed of specific and function during meiotic prophase I and spermiogenesis. meiotic cohesins (e.g., REC8 and RAD21L) (Gutiérrez-Caballero et al., 2011; Llano et al., 2012) and proteins of the SC (e.g., DISCUSSION SYCP3) (Henderson and Keeney, 2005), preventing interactions below 40 Mbp. Second, the association of recombination hot- Here, we provide a high-resolution structural and functional atlas spots and cohesins in primary spermatocytes suggests that of mouse spermatogenesis. Our data reveal the compartmental- cohesion-mediated transcription in genomic regions out of the Figure 7. Higher-Order Chromatin Structure and Gene Expression in the X Chromosome (A) ICE-corrected Hi-C matrices for the X chromosome in mouse germ cells, at a 50-Kbp resolution. Plaid blue regions correspond to non-mapped bins. (B) Representation of compartment signal (first eigenvector) along the mouse’s X chromosome. (C) Upper panel: Overview of spermatogenesis with two pairs of autosomes (black and gray lines) and the sex chromosomes (X and Y as red lines). Meiotic sex chromosome inactivation (MSCI) characterizes prophase I (shown as a pink cloud; also known as the sex body). In RSs, all centromeres associate in the center of the cell forming the chromocenter, with the X adjacent, forming the post-meiotic sex chromatin (PMSC). Lower panel: Boxplots of expression (FPKM values) of X genes binned according to A-B compartment or ampliconic. N/A are not assigned to a compartment or amplicon. (D) X chromosome TAD alignments. Arches represent TADs with higher (darker gray) and lower (lighter gray) than expected intra-TAD interactions. (E) CTCF and cohesin distribution on the X in P/D (green) and RSs (black). Ampliconic regions are green, and evolutionary strata are displayed, with blue rep- resenting older strata and red representing newer strata (see Figure S7 for further details). (F) X chromosome distribution of CTCF and cohesins’ occupancy relative to TSSs and other genomic features in P/D and RS. (G) Expression (FPKM values) of genes with CTCF and cohesin peaks located at promoter regions (2 kbp from TSS). Asterisks represent statistically significant differential gene expression when compared with all genes in the mouse genome (p < 0.01). Boxes represent first and third quartiles, whereas black bars in boxes represent the median values. (H) Expression changes (versus Spg) of representative X genes that reduce expression in P/D and maintain low levels during spermatogenesis (e.g., Cybb), genes that reduce expression in P/D and increase in RSs and sperm (i.e., Cdk16), genes that increase expression in RSs (e.g., Actr1), and genes that increase expression in P/D and then reduce (i.e., Cldn34c4). Fib, fibroblast; Spg, spermatogonia; L/Z, leptonema-zygonema; P/D, pachynema-diplonema; RS, round spermatids. See also Figure S7. Cell Reports 28, 352–367, July 9, 2019 363 axes can provide an environment conducive to both gene loops emerge (Haering and Jessberger, 2012). This structure expression and the formation of DSBs. This is reflected by the in- could block REC8 and RAD21L access to chromatin, preventing teractions observed at shorter distances (2.5–4.5 Mbp). detection at the axis by ChIP-seq. Our data suggest that cohe- Compartments re-emerged in post-meiotic cells with a sins associate with active promoters most probably located in remarkable high-level organizational difference between RSs DNA loops out of the axes, hinting at a functional role that and sperm. Chromosomal territories were re-established in adds to the well-known structural role of these cohesins in the both cell types but with chromatin less densely packed in RSs, formation of meiotic chromosomal axes (Llano et al., 2012). In as revealed by the interaction ratio and 3D-FISH. In RSs, intra- fact, mice with REC8 and RAD21L deficiencies are infertile; chromosomal interactions were greater than inter-chromosomal meiosis is arrested in early prophase I (Bannister et al., 2004; interactions, the opposite to what was observed in sperm. In Herrán et al., 2011), precluding analysis of potential regulatory both cell types, A compartments (which correlate with gene roles for cohesins during spermiogenesis. expression) were larger than in pre-meiotic cells and were over- Our results also point to the role of meiotic cohesins in regu- represented when compared to fibroblasts and Spg. Moreover, lating gene expression after meiosis, which was characterized the transformation of RSs into spermatozoa was accompanied by the occupancy of RAD21L at the TSSs of expressed genes by the transcription of genes related to spermiogenesis and in RSs. We suggest that CTCF and cohesins could have a sperm function located in newly created A compartments. This synergistic role in establishing transcriptional hubs for early agrees with a second burst of gene expression in RSs (da Cruz embryonic development in meiotic prophase I (where the three et al., 2016) and supports the more recent idea that sperm are proteins co-localize), while also fine-tuning subsequent sper- not inactive cells (Jodar et al., 2016; Jung et al., 2017). In fact, matogenesis progression. Moreover, cohesins in RSs appear it was shown that a large number of sperm promoters are in an to correlate with the expression of genes implicated in early active epigenetic state, suggesting that this genetic information embryonic development (e.g., Nanog), providing immediate can influence embryo development upon fertilization (Jung et al., access to the major molecular pathways involved in organogen- 2017). esis upon fertilization. These results suggest that some cellular There was a remarkable reduction of TADs in sperm. One functional activity is predetermined in early spermiogenesis, strength of our approach that differs from others (Jung et al., before RSs are differentiated into sperm, which raises the 2017; Wang et al., 2019) was the cell-sorting strategy, which intriguing possibility that cohesins play a role in such functional permitted the isolation of sperm from the rest of the germ cell predetermination. populations. Coupled with the Hi-C simulations, this permitted Here, we extend initial observations (Alavattam et al., 2019; us to attribute higher-order chromatin structures specifically Patel et al., 2019; Wang et al., 2019) of X chromosome remodel- to sperm. Thus, the reduction of the TAD insulation score in ing in prophase I to pre-meiotic and post-meiotic cells. Our sperm can be linked to histone replacement by protamines. functional and structural analyses show that silencing of the Protamines replace most histones in sperm, with the help of X chromosomes (MSCI) is accompanied by distinct changes transition proteins and H2A histone variants, folding DNA into in the higher-order chromatin structure at different levels: (1) toroidal subunits at the kbp level (Barral et al., 2017). Because more intra- than inter-chromosomal interactions, a pattern of the structural constraints of protamines, it is tempting to already present in Spg; (2) strong compartmentalization; and speculate that the compact chromatin organization in sperm (3) a reduction in both the TAD number and signal once prophase is associated with the presence of A compartments at the I initiates. Remarkably, architectural proteins were remodeled, Mbp scale, which likely correlates with active histone marks in which included a reduction of cohesin peaks (compared to auto- sperm (Jung et al., 2017), but not with the formation of fine-scale somes) and an almost complete depletion of CTCF. These data TAD structures (i.e., kbp). demonstrate that a reduction of the cohesin load correlates with The genomic atlas for CTCF and meiotic-specific cohesins an absence of chromosome X compartmentalization (A/B). presented here provides an unprecedented view of the connec- Despite MSCI, gene expression still takes place, which corre- tion among chromatin remodeling, architectural proteins, and lates with meiotic cohesin occupancy but not with CTCF. gene expression. Strikingly, the vast majority of cohesin peaks Altogether, these observations support a role for cohesin in detected in primary spermatocytes were localized within pro- regulating gene transcription. moter regions of genes, with significantly higher expression In summary, we have implemented a robust workflow to than non-cohesin-associated genes. The conflicting patterns provide an integrated structural and functional framework of of RAD21L on the X chromosome, with reduced ChIP peak the 3D organization of the mouse genome in germ cells. We numbers but increased super-resolution microscopy labeling detected a fine-tuned balance among chromatin remodeling, intensity, are counterintuitive. The increased labeling intensity architectural proteins, and gene expression during spermato- at the axis suggests that RAD21L has an important X chromo- genesis. Overall, our results provide insights into how the some scaffolding function, but the reduced number of peaks structural organization of the genome influences cellular differ- (compared to the autosomes) indicates that the detected peaks entiation, especially in the context of the dramatic chromatin correspond to cohesins that do not interact directly with DNA at changes that take place during the formation and differentiation the axis. In fact, REC8 and RAD21L bind to both head domains of the mammalian male germline. Future functional perturbation of a structural maintenance of chromosome (SMC) heterodimer, analyses will help us understand the mechanism by which 3D forming a ring-like protein structure topologically encircling sis- genome folding changes shape transcriptional activity during ter chromatids to the SC, a core axis from which chromatin spermatogenesis. 364 Cell Reports 28, 352–367, July 9, 2019 STAR+METHODS (CLC–2017–01). C.V. is supported by a FPI predoctoral fellowship from the Ministry of Economy and Competitiveness (BES-2015-072924). A.P.-G. is supported by a ‘Doctorats Industrials’ predoctoral fellowship (AGAUR). H.H. Detailed methods are provided in the online version of this paper is a Miguel Servet (CP14/00229) researcher funded by the Spanish Institute and include the following: of Health Carlos III (ISCIII), the Agencia Estatal de Investigación (AEI), and FEDER (SAF2017-89109-P). C.M. is an Asociación Española Contra el Cáncer d KEY RESOURCES TABLE (AECC) postdoctoral fellow. M.A.M.-R. acknowledges support by the Euro- d LEAD CONTACT AND MATERIALS AVAILABILITY pean Research Council under the 7th Framework Program FP7/2007-2013 d EXPERIMENTAL MODEL AND SUBJECT DETAILS (ERC grant agreement 609989) and the European Union’s Horizon 2020 B Animals Research and Innovation Programme (grant agreement 676556). A.M.P. and B Cell lines A.R.-H. also acknowledge support from MeioNet (BFU2015-71786-REDT). d METHOD DETAILS B Fluorescence Activated Cell Sorting (FACS) of mouse AUTHOR CONTRIBUTIONS male germ cells A.R.-H. conceived and devised the study. C.V., A.P.-G., Y.C., F.L.D., H.H., B Fibroblast cell culture M.A.M.-R., and A.R.-H. designed experiments. C.V., Y.C., F.L.D., F.G., B In nuclei Hi-C J.S.-C., L.G.-H., E.J., C.M., O.F., and A.R.-H. performed experiments. C.V., B Hi-C data processing, binning and normalization A.P.-G., Y.C., F.L.D., L.G.-H., A.M.P., P.D.W., M.A.M.-R., and A.R.-H. B Correlation coefficient analysis analyzed the data. R.A.C., W.S., A.M.P., H.H., M.A.M.-R., and A.R.-H. contrib- B Averaged contact probability P(s) uted to reagents and data collection. C.V., A.P.-G., P.D.W., and A.R.-H. wrote B Inter-/intra-chromosome interaction ratio the first draft of the manuscript with input from all authors. All authors read and approved the final version of the manuscript. B Inter-subtelomeric interaction quantification B Sperm simulations DECLARATION OF INTERESTS B A-B compartments and TAD calling B Compartment switching The authors declare no competing interests. B 3D-FISH and confocal microscopy B Western blot analysis of CTCF and cohesins in germ Received: March 8, 2019 cells Revised: April 18, 2019 B Spermatocyte spreads and immunofluorescence anal- Accepted: June 7, 2019 ysis of CTCF and cohesins Published: July 9, 2019 B Stimulated emission depletion (STED) microscopy SUPPORTING CITATIONS B ChIP-sequencing B ChIP-seq peak calling and annotation The following reference appears in the Supplemental Information: Cortez et al. B RNA-sequencing (2014). B RNA profiling B Analysis of correspondence between compartments REFERENCES and gene expression B Local insulation, cohesin occupancy and gene expres- Alavattam, K.G., Maezawa, S., Sakashita, A., Khoury, H., Barski, A., Kaplan, sion N., and Namekawa, S.H. (2019). Attenuated chromatin compartmentalization B X chromosome evolutionary strata and ampliconic re- in meiosis and its maturation in sperm development. Nat. Struct. Mol. Biol. 26, 175–184. gions Bannister, L.A., Reinholdt, L.G., Munroe, R.J., and Schimenti, J.C. (2004). Po- d QUANTIFICATION AND STATISTICAL ANALYSIS sitional cloning and characterization of mouse mei8, a disrupted allelle of the d DATA AND CODE AVAILABILITY meiotic cohesin Rec8. Genesis 40, 184–194. Barral, S., Morozumi, Y., Tanaka, H., Montellier, E., Govin, J., de Dieuleveult, SUPPLEMENTAL INFORMATION M., Charbonnier, G., Couté, Y., Puthier, D., Buchou, T., et al. (2017). Histone Variant H2A.L.2 Guides Transition Protein-Dependent Protamine Assembly Supplemental Information can be found online at https://doi.org/10.1016/j. in Male Germ Cells. Mol. Cell 66, 89–101.e8. celrep.2019.06.037. Brick, K., Smagulova, F., Khil, P., Camerini-Otero, R.D., and Petukhova, G.V. (2012). Genetic recombination is directed away from functional genomic ACKNOWLEDGMENTS elements in mice. Nature 485, 642–645. Bushnell, B. (2014). BBMap: A Fast, Accurate, Splice-Aware Aligner. (Depart- This work was supported by the Ministry of Economy and Competitiveness ment of Energy), Joint Genome Institute. (BFU2017-89408-R to A.M.P.; BFU2013-47736-P and BFU2017-85926-P to Carbon, S., Ireland, A., Mungall, C.J., Shu, S., Marshall, B., and Lewis, S.; the M.A.M.-R.; and CGL2014-54317-P and CGL2017-83802-P to A.R.-H.) and AmiGO Hub; the Web Presence Working Group (2008). AmiGO: online access the Agència de Gestió d’Ajuts Universitaris i de Recerca (AGAUR) (DI2015 to to ontology and annotation data. Bioinformatics 25, 288–289. A.R.-H. and R.A.C., as well as SGR468 to M.A.M-R.). Work at CRG, BIST, and UPF was in part funded by the Spanish Ministry of Economy and Compet- Cortez, D., Marin, R., Toledo-Flores, D., Froidevaux, L., Liechti, A., Waters, itiveness, ‘Centro de Excelencia Severo Ochoa 2013-2017’ (SEV-2012-0208), €tzner, F., and Kaessmann, H. (2014). Origins and functional evolution P.D., Gru and ‘Centro de Excelencia Marı́a de Maeztu 2016-2019.’ CIC-IBMCC is sup- of Y chromosomes across mammals. Nature 508, 488–493. ported by the Programa de Apoyo a Planes Estratégicos de Investigación da Cruz, I., Rodrı́guez-Casuriaga, R., Santiñaque, F.F., Farı́as, J., Curti, G., Ca- de Estructuras de Investigación de Excelencia, co-funded by Junta de poano, C.A., Folle, G.A., Benavente, R., Sotelo-Silveira, J.R., and Geisinger, A. Castilla y León (CSI239P18) and the European Regional Development Fund (2016). Transcriptome analysis of highly purified mouse spermatogenic cell Cell Reports 28, 352–367, July 9, 2019 365 populations: gene expression signatures switch from meiotic-to postmeiotic- Leng, N., Dawson, J.A., Thomson, J.A., Ruotti, V., Rissman, A.I., Smits, related processes at pachytene stage. BMC Genomics 17, 294–313. B.M.G., Haag, J.D., Gould, M.N., Stewart, R.M., and Kendziorski, C. (2013). Dekker, J., Marti-Renom, M.A., and Mirny, L.A. (2013). Exploring the three- EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq exper- dimensional organization of genomes: interpreting chromatin interaction iments. Bioinformatics 29, 1035–1043. data. Nat. Rev. Genet. 14, 390–403. Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Dixon, J.R., Selvaraj, S., Yue, F., Kim, A., Li, Y., Shen, Y., Hu, M., Liu, J.S., and Abecasis, G., and Durbin, R.; 1000 Genome Project Data Processing Sub- Ren, B. (2012). Topological domains in mammalian genomes identified by group (2009). The sequence alignment/map format and SAMtools. Bioinfor- analysis of chromatin interactions. Nature 485, 376–380. matics 25, 2078–2079. Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, Liao, Y., Smyth, G.K., and Shi, W. (2014). featureCounts: an efficient general P., Chaisson, M., and Gingeras, T.R. (2013). STAR: ultrafast universal RNA-seq purpose program for assigning sequence reads to genomic features. Bioinfor- aligner. Bioinformatics 29, 15–21. matics 30, 923–930. Du, Z., Zheng, H., Huang, B., Ma, R., Wu, J., Zhang, X., He, J., Xiang, Y., Wang, Lieberman-Aiden, E., van Berkum, N.L., Williams, L., Imakaev, M., Ragoczy, Q., Li, Y., et al. (2017). Allelic reprogramming of 3D chromatin architecture dur- T., Telling, A., Amit, I., Lajoie, B.R., Sabo, P.J., Dorschner, M.O., et al. ing early mammalian development. Nature 547, 232–235. (2009). Comprehensive mapping of long-range interactions reveals folding Feng, J., Liu, T., Qin, B., Zhang, Y., and Liu, X.S. (2012). Identifying ChIP-seq principles of the human genome. Science 326, 289–293. enrichment using MACS. Nat. Protoc. 7, 1728–1740. Llano, E., Herrán, Y., Garcı́a-Tuñón, I., Gutiérrez-Caballero, C., de Álava, E., Flyamer, I.M., Gassler, J., Imakaev, M., Branda~o, H.B., Ulianov, S.V., Abden- Barbero, J.L., Schimenti, J., de Rooij, D.G., Sánchez-Martı́n, M., and Pendás, nur, N., Razin, S.V., Mirny, L.A., and Tachibana-Konwalski, K. (2017). Single- A.M. (2012). Meiotic cohesin complexes are essential for the formation of the nucleus Hi-C reveals unique chromatin reorganization at oocyte-to-zygote axial element in mice. J. Cell Biol. 197, 877–885. transition. Nature 544, 110–114. Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold Gel, B., and Serra, E. (2017). karyoploteR: an R/Bioconductor package to plot change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, customizable genomes displaying arbitrary data. Bioinformatics 33, 3088– 550–571. 3090. Marco-Sola, S., Sammeth, M., Guigó, R., and Ribeca, P. (2012). The GEM Gel, B., Dı́ez-Villanueva, A., Serra, E., Buschbeck, M., Peinado, M.A., and Ma- mapper: fast, accurate and versatile alignment by filtration. Nat. Methods 9, linverni, R. (2016). regioneR: an R/Bioconductor package for the association 1185–1188. analysis of genomic regions based on permutation tests. Bioinformatics 32, Mi, H., Huang, X., Muruganujan, A., Tang, H., Mills, C., Kang, D., and Thomas, 289–291. P.D. (2017). PANTHER version 11: expanded annotation data from Gene Gibcus, J.H., Samejima, K., Goloborodko, A., Samejima, I., Naumova, N., Nue- Ontology and Reactome pathways, and data analysis tool enhancements. Nu- bler, J., Kanemaki, M.T., Xie, L., Paulson, J.R., Earnshaw, W.C., et al. (2018). A cleic Acids Res. 45, D183–D189. pathway for mitotic chromosome formation. Science 359, 652. Mueller, J.L., Mahadevaiah, S.K., Park, P.J., Warburton, P.E., Page, D.C., and Gutiérrez-Caballero, C., Herrán, Y., Sánchez-Martı́n, M., Suja, J.Á., Barbero, Turner, J.M.A. (2008). The mouse X chromosome is enriched for multicopy J.L., Llano, E., and Pendás, A.M. (2011). Identification and molecular charac- testis genes showing postmeiotic expression. Nat. Genet. 40, 794–799. terization of the mammalian a-kleisin RAD21L. Cell Cycle 10, 1477–1487. Namekawa, S.H., Park, P.J., Zhang, L.-F., Shima, J.E., McCarrey, J.R., Gris- Haaf, T., and Ward, D.C. (1995). Higher order nuclear structure in mammalian wold, M.D., and Lee, J.T. (2006). Postmeiotic sex chromatin in the male germ- sperm revealed by in situ hybridization and extended chromatin fibers. Exp. line of mice. Curr. Biol. 16, 660–667. Cell Res. 219, 604–611. Naumova, N., Imakaev, M., Fudenberg, G., Zhan, Y., Lajoie, B.R., Mirny, L.A., Haering, C.H., and Jessberger, R. (2012). Cohesin in determining chromosome and Dekker, J. (2013). Organization of the mitotic chromosome. Science 342, architecture. Exp. Cell Res. 318, 1386–1393. 948–953. Handel, M.A., and Schimenti, J.C. (2010). Genetics of mammalian meiosis: Patel, L., Kang, R., Rosenberg, S.C., Qiu, Y., Raviram, R., Chee, S., Hu, R., regulation, dynamics and impact on fertility. Nat. Rev. Genet. 11, 124–136. Ren, B., Cole, F., and Corbett, K.D. (2019). Dynamic reorganization of the Henderson, K.A., and Keeney, S. (2005). Synaptonemal complex formation: genome shapes the recombination landscape in meiotic prophase. Nat. where does it start? BioEssays 27, 995–998. Struct. Mol. Biol. 26, 164–174. Herrán, Y., Gutiérrez-Caballero, C., Sánchez-Martı́n, M., Hernández, T., Viera, Picelli, S., Faridani, O.R., Björklund, A.K., Winberg, G., Sagasser, S., and A., Barbero, J.L., de Álava, E., de Rooij, D.G., Suja, J.Á., Llano, E., and Pendás, Sandberg, R. (2014). Full-length RNA-seq from single cells using Smart- A.M. (2011). The cohesin subunit RAD21L functions in meiotic synapsis and seq2. Nat. Protoc. 9, 171–181. exhibits sexual dimorphism in fertility. EMBO J. 30, 3091–3105. Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a flexible suite of utilities for Jodar, M., Sendler, E., and Krawetz, S.A. (2016). The protein and transcript comparing genomic features. Bioinformatics 26, 841–842. profiles of human semen. Cell Tissue Res. 363, 85–96. Ramı́rez, F., Ryan, D.P., Gru€ning, B., Bhardwaj, V., Kilpert, F., Richter, A.S., Jung, Y.H., Sauria, M.E.G., Lyu, X., Cheema, M.S., Ausio, J., Taylor, J., and Heyne, S., Du €ndar, F., and Manke, T. (2016). deepTools2: a next generation Corces, V.G. (2017). Chromatin states in mouse sperm correlate with embry- web server for deep-sequencing data analysis. Nucleic Acids Res. 44 (W1), onic and adult regulatory landscapes. Cell Rep. 18, 1366–1382. W160–W165. Keeney, S., Giroux, C.N., and Kleckner, N. (1997). Meiosis-specific DNA dou- €ning, B.A., Villaveces, J., Ramı́rez, F., Bhardwaj, V., Arrigoni, L., Lam, K.C., Gru ble-strand breaks are catalyzed by Spo11, a member of a widely conserved Habermann, B., Akhtar, A., and Manke, T. (2018). High-resolution TADs reveal protein family. Cell 88, 375–384. DNA sequences underlying genome organization in flies. Nature Communica- Kleckner, N., Storlazzi, A., and Zickler, D. (2003). Coordinate variation in tions 9, 189. meiotic pachytene SC length and total crossover/chiasma frequency under Rao, S.S.P., Huntley, M.H., Durand, N.C., Stamenova, E.K., Bochkov, I.D., conditions of constant DNA length. Trends Genet. 19, 623–628. Robinson, J.T., Sanborn, A.L., Machol, I., Omer, A.D., Lander, E.S., and Aiden, Lahn, B.T., and Page, D.C. (1999). Four evolutionary strata on the human E.L. (2014). A 3D map of the human genome at kilobase resolution reveals X chromosome. Science 286, 964–967. principles of chromatin looping. Cell 159, 1665–1680. Lange, J., Yamada, S., Tischfield, S.E., Pan, J., Kim, S., Zhu, X., Socci, N.D., Rau, A., Gallopin, M., Celeux, G., and Jaffrézic, F. (2013). Data-based filtering Jasin, M., and Keeney, S. (2016). The landscape of mouse meiotic double- for replicated high-throughput transcriptome sequencing experiments. Bioin- strand break formation, processing, and repair. Cell 167, 695–708.e16. formatics 29, 2146–2152. 366 Cell Reports 28, 352–367, July 9, 2019 Reig-Viader, R., Garcia-Caldés, M., and Ruiz-Herrera, A. (2016). Telomere ho- Tian, T., Liu, Y., Yan, H., You, Q., Yi, X., Du, Z., Xu, W., and Su, Z. (2017). meostasis in mammalian germ cells: a review. Chromosoma 125, 337–351. agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Robinson, M.D., McCarthy, D.J., and Smyth, G.K. (2010). edgeR: a Bio- Nucleic Acids Res. 45, W122–W129. conductor package for differential expression analysis of digital gene expres- Turner, J.M.A. (2007). Meiotic sex chromosome inactivation. Development sion data. Bioinformatics 26, 139–140. 134, 1823–1831. Ruiz-Herrera, A., Vozdova, M., Fernández, J., Sebestova, H., Capilla, L., Froh- lich, J., Vara, C., Hernández-Marsal, A., Sipek, J., Robinson, T.J., and Rubes, Wang, S., Zickler, D., Kleckner, N., and Zhang, L. (2015). Meiotic crossover J. (2017). Recombination correlates with synaptonemal complex length and patterns: obligatory crossover, interference and homeostasis in a single pro- chromatin loop size in bovids-insights into mammalian meiotic chromosomal cess. Cell Cycle 14, 305–314. organization. Chromosoma 126, 615–631. Sánchez-Guillén, R.A., Capilla, L., Reig-Viader, R., Martı́nez-Plana, M., Pardo- Wang, Y., Wang, H., Zhang, Y., Du, Z., Si, W., Fan, S., Qin, D., Wang, M., Duan, Camacho, C., Andrés-Nieto, M., Ventura, J., and Ruiz-Herrera, A. (2015). On Y., Li, L., et al. (2019). Reprogramming of meiotic chromatin architecture dur- the origin of Robertsonian fusions in nature: evidence of telomere shortening ing spermatogenesis. Mol. Cell 73, 547–561.e6. in wild house mice. J. Evol. Biol. 28, 241–249. Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis (Springer- Sassone-Corsi, P. (2002). Unique chromatin remodeling and transcriptional Verlag). regulation in spermatogenesis. Science 296, 2176–2178. Serra, F., Baù, D., Goodstadt, M., Castillo, D., Filion, G.J., and Marti-Renom, Yang, T., Zhang, F., Yardımcı, G.G., Song, F., Hardison, R.C., Noble, W.S., M.A. (2017). Automatic analysis and 3D-modelling of Hi-C data using TADbit Yue, F., and Li, Q. (2017). HiCRep: assessing the reproducibility of Hi-C data reveals structural features of the fly chromatin colors. PLoS Comput. Biol. using a stratum-adjusted correlation coefficient. Genome Res. 27, 1939–1949. 13, e1005665. Yu, G., Wang, L.G., and He, Q.Y. (2015). ChIPseeker: an R/Bioconductor pack- Tarazona, S., Garcı́a-Alcalde, F., Dopazo, J., Ferrer, A., and Conesa, A. (2011). age for ChIP peak annotation, comparison and visualization. Bioinformatics Differential expression in RNA-seq: a matter of depth. Genome Res. 21, 2213– 31, 2382–2383. 2223. Thorvaldsdóttir, H., Robinson, J.T., and Mesirov, J.P. (2013). Integrative Geno- Zhu, L.J., Gazin, C., Lawson, N.D., Pagès, H., Lin, S.M., Lapointe, D.S., and mics Viewer (IGV): high-performance genomics data visualization and explora- Green, M.R. (2010). ChIPpeakAnno: a Bioconductor package to annotate tion. Brief. Bioinform. 14, 178–192. ChIP-seq and ChIP-chip data. BMC Bioinformatics 11, 237–247. Cell Reports 28, 352–367, July 9, 2019 367 STAR+METHODS KEY RESOURCES TABLE Reagent or Resource Source Identifier Antibodies Anti-mouse SYCP3 Abcam Cat#ab97672; RRID:AB_10678841 Anti-Rabbit DMC1 Santa Cruz Biotechnologies Cat#sc-22768; RRID:AB_2277191 Anti-mouse Cy5 Jackson Immunoresearch Cat#115-175-146; RRID:AB_2338713 Anti-rabbit FITC Jackson Immunoresearch Cat#111-095-003; RRID:AB_2337972 Anti-CTCF Millipore Cat#07-729; RRID:AB_441965 Anti-REC8 Courtesy of A.M. Pendás N/A Anti-RAD21L Courtesy of A.M. Pendás N/A Anti-HRP-PO Bio-Rad Cat#1706515; RRID:AB_11125142 Chemicals, Peptides and Recombinant Proteins Collagenase Type II Life Technologies #17101015 DNase I Sigma-Aldrich #DN25-10MG Trypsin from bovine pancreas Sigma-Aldrich #T9935-100MG Hoechst 33342 Life Technologies #H3570-10ml Complete Protease Inhibitor Roche #1187358001 MboI New England Biolabs #R0147M Biotin-14-dATP Life Technologies #19524-016 DNA Polymerase I, large (Klenow) Fragment New England Biolabs #M0210M T4 DNA Ligase New England Biolabs #M0202M Dynabeads MyOne Streptavidin T1 Life Technologies #65001 T4 Polunucleotides Kinase New England Biolabs #M0201L T4 DNA Polymerase New England Biolabs #M0212M Klenow Fragment 30 - > 50 exo- New England Biolabs #M0212M AMPure XP Beads Beckman-Coulter #A63880 2x Laemmli Buffer Bio-Rad Laboratories #1610737 Trans-Blot Turbo Transfer Packs Bio-Rad Laboratories #1704158 Ponceau S solution Sigma-Aldrich #P7170-1L ClarityTM Western ECL Substrate, 500 ml Bio-Rad Laboratories #1705061 OneDay ChIP kit Diagenode C01010080 Unblocked protein A Diagenode C03020002 NEBNext End Repair Module New England Biolabs E6050S NEBNext End Repair Reaction Buffer New England Biolabs B6052S NEBNEXT dA-tailing Reaction Buffer New England Biolabs B6059S Deposited Data Hi-C data This paper GEO:GSE132054 ChIP-seq data This paper GEO:GSE132054 RNA-seq data This paper GEO:GSE132054 SPO11-oligos hotspots Lange et al., 2016 GEO:GSE84689 H3K4me3 data Brick et al., 2012 GEO:GSE35498 Hi-C data from sperm Jung et al., 2017 SRR3225862 and SRR3225863 Experimental Models: Organisms/Strains C57BL/6J Charles River Laboratories N/A (Continued on next page) e1 Cell Reports 28, 352–367.e1–e9, July 9, 2019 Continued Reagent or Resource Source Identifier Software and algorithms BBDuk (version 10/2015) Bushnell, 2014 https://sourceforge.net/projects/bbmap/ TADbit (version 0.2.0.23) Serra et al. (2017 https://github.com/3DGenomes/TADbit GEM (version 1.7.1) Marco-Sola et al., 2012 https://sourceforge.net/projects/gemlibrary/ HiCExplorer (version 1.8.1) Ramı́rez et al., 2018 https://github.com/deeptools/HiCExplorer HiCRep (version 1.4) Yang et al., 2017 https://github.com/MonkeyLB/hicrep BEDtools (version 2.17) Quinlan and Hall, 2010 https://github.com/arq5x/bedtools2 SAM tools Li et al., 2009 http://samtools.sourceforge.net/ Deep tools Ramı́rez et al., 2016 https://github.com/deeptools/deepTools MACS2 Feng et al., 2012 https://github.com/taoliu/MACS IGV tools Thorvaldsdóttir et al., 2013 https://igv.org/ ChIPseeker Yu et al., 2015 https://github.com/GuangchuangYu/ChIPseeker KaryoploteR Gel and Serra, 2017 http://bioconductor.org/packages/release/bioc/ html/karyoploteR.html RegioneR Gel et al., 2016 https://bioconductor.org/packages/release/bioc/ html/regioneR.html Panther Mi et al., 2017 http://www.pantherdb.org/ AmiGo Carbon et al., 2008 http://amigo.geneontology.org AIR This paper https://transcriptomics.cloud LEAD CONTACT AND MATERIALS AVAILABILITY Further information and requests for resources and reagents should be directed to, and will be fulfilled by, the Lead Contact: Aurora Ruiz-Herrera (aurora.ruizherrera@uab.cat). EXPERIMENTAL MODEL AND SUBJECT DETAILS Animals C57BL/6J (B6) male mice at 8-17 weeks of age were purchased from the Charles River Laboratories. Animal maintenance and exper- imental procedures were carried out according to the Ethics Committee on Animal and Human Experimentation (CEEAH) guidelines from Universitat Autònoma de Barcelona (UAB). Cell lines A primary fibroblast cell line derived form a male mouse (C57BL/6J strain) previously established in our lab (Sánchez-Guillén et al., 2015), and cultured in DMEM medium supplemented with 10% fetal bovine serum and 1% PenStrep at 37 C and 5% CO2. METHOD DETAILS Fluorescence Activated Cell Sorting (FACS) of mouse male germ cells Male mice were dissected by an inguinal incision to remove whole testis and cauda epididymis. Sperm cells were recovered from the epididymis by pressing the deferent tube gently. Both decapsulated testis and epididymis were mechanically disaggregated together with the ejaculate, incubated in GBSS with collagenase type II (Life Technologies) (0.5 mg/ml) and DNase I (Sigma Aldrich) (1 ng/ml) at 33 C for 15min with agitation. Fetal bovine serum (ThermoFisher Scientific) was added to the mix at a proportion of 5% in order to inactivate trypsin. The cell suspension was filtered through a 70 mm diameter cell strainer and subsequently centrifuged for 3 min at 1,800 xg. Cells were then dyed with 5 mg/ml of Hoechst 33342 (Life Technologies) for 30 min at 33 C with agitation, then held at 4 C until sorting. Germ cells were sorted using a BD InfluxTM (BD Biosciences) coupled with an ultraviolet laser (355 nm). Subsequently, three sorting strategies were followed in order to obtain a total of six highly enriched germ cell populations: (i) spermatogonia (2n, 2c), (ii) spermatocytes at leptonema/zygonema (L/Z) stage (2n, 4c), (iii) spermatocytes at pachynema/diplonema (P/D) stage (2n, 4c), (iv) secondary spermatocytes (1n, 2c), (v) round spermatids and (vi) sperm. Strategy #1: Four main germ cell populations (spermatogonia, P/D, secondary spermatocytes and round spermatids) were first isolated by plotting Hoechst Blue (UV355-460/50) versus Hoechst red (UV355-670/30) emissions to discriminate cells by both their DNA content and their complexity (Figure S1A). The spermatogonia fraction (2n, 2c) represented a heterogeneous population of undifferentiated and differentiating spermatogonia. Cell Reports 28, 352–367.e1–e9, July 9, 2019 e2 Strategy #2: In order to discriminate primary spermatocytes in L/Z and P/D stages, an ‘‘in solution’’ immunofluorescence using DMC1 (1:300, Santa Cruz Biotechnologies) and SYCP3 (1:1000, Abcam) antibodies was performed before sorting. Briefly, DMC1 (1:300, Santa Cruz Biotechnologies) and SYCP3 (1:1000, Abcam) antibodies were added to the cell suspension that resulted from testis disaggregation. Cells showing DMC1(+)/SYCP3(+) staining corresponded to primary spermatocytes in L/Z stage, whereas cells at P/D stage appeared with DMC1(-)/SYCP3(+) staining (Figure S1B). Strategy #3: A third approach was followed to obtain highly enriched populations of sperm and round spermatids as both cell types have the same DNA content and can only be discriminated by morphology. We followed a re-gating approach since cell morphology can be discriminated in the flow sorter by how cells project a side shadow (side scatter, SSC) and a forward shadow (forward scatter, FSC), a plot for each initial gate was drawn (SSC versus FSC plot). With this approach, and applying a restrictive re-gating approach, sperm and round spermatids were isolated with high enrichment (Figure S1C). This allowed us to exclude elongating spermatids from these samples. Likewise, isolation of sperm from cauda epididymis was conducted by re-gating of the c fraction, resulting in isolated sperm Figure S1D). Irrespective of the sorting strategy, all cell populations were collected after sorting in 1x PBS, centrifuged for 5 min at 1,800 xg. The supernatant was discarded, and cell pellets were flash-frozen at 80 C until use. Sorts were typically 3-6 hours to collect between 0.2 3 106 and 3.2 3 106 cells, depending on the germ cell population. Cell enrichment of each flow-sorted population was evaluated by immunofluorescence using specific meiotic proteins and DAPI morphology (Figure S1). Spermatogonia were distinguished using an anti-CD90 (Thy-1+) (1:10) antibody and DAPI morphology. For primary spermatocytes, prophase-I stages (leptonema, zygonema, pachynema and diplonema) were identified based on SYCP3 (1:400) and ɣH2AX patterns (1:300). For secondary spermatocytes, immunofluorescence using anti-H3K9me3 (1:500), anti-actin (1:400) and anti-tubulin (1:1000) was performed. Cell enrichment of round spermatids and sperm fractions was determined based on nucleus morphology and DAPI pattern. Cells were fixed on slides and then mounted with DAPI diluted in Vectashield (Vector Lab- oratories). In all cases, slides were analyzed using a fluorescence microscopy (Axiophot, Zeiss) coupled with a ProgRes  CS10plus, Jenoptik camera. Representative images were captured with ACO XY (A. Coloma, Open Microscopy). Between 50 and 100 cells were counted for each flow-sorted population. Only sorted populations with enrichment above 80% were considered for subsequent experiments. Fibroblast cell culture The mouse primary fibroblast cell line was cultured to 100% confluence. A cell fraction was kept for quality controls (see below) and the rest were fixated accordingly to the crosslink step of the Hi-C method (see in nuclei Hi-C section) for adherent cells. Quality controls consisted of generating a modal karyotype and flow cytometry cell cycle analysis. Chromosome spreads were obtained using standard protocols. Briefly, cells were arrested in metaphase by adding 80 mL of Colcemid (10 mg/ml) (Sigma) to 10ml of medium for 2h and then trypsinised. Cells were spun down at 600 xg for 7 minutes and resuspended in 5ml of hypotonic solution (0.075M KCl) for 20 minutes at 37 C. Chromosomes were then fixed by addition of fixative solution (3:1 methanol/acetic acid) and metaphase spreads were obtained by dropping 15 ml of cell suspension onto a cleaned dry slide. Slides were kept at 20 C until use. Metaphases were stained homogenously with Giemsa solution for analysis of the modal karyotype. An optical microscope (model Zeiss Axioskop) equipped with a charged coupled device camera (ProgResR CS10Plus, Jenoptik) was used for the microscope analysis. Good-quality metaphases were captured with the program Progress Capture 2.7.7 and analyzed for each specimen, obtaining the modal karyotype. For the cell cycle analysis, fibroblasts were cultured until cells reached 100% confluence. Next, cells were fixated with 1% form- aldehyde for 10’ at RT and incubated with glycine 0.125M for 50 at RT and for 15’ at 4 C to stop the crosslinking reaction. Subse- quently, cells were scrapped off the flask, pelleted and resuspended in 1x PBS to finally be stained with Hoechst 33342. Cell cycle was analyzed using a BD InfluxTM (BD Biosciences) coupled with an ultraviolet laser (355 nm) to reveal that over 74% of the cells were in G1-phase (Figure S1). In nuclei Hi-C The generation of Hi-C libraries was optimized from the protocol developed by Rao et al. (2014). Different approaches were used depending on the cell type analyzed (somatic and germ cells, see below). Two replicates for each FACS-sorted population were ob- tained (with the exception of primary spermatocytes at L/Z stage where one single experiment was performed) from a total of 7.7 3 106 spermatogonia, 1.7 3 106 primary spermatocytes at L/Z stage, 20.9 3 106 primary spermatocytes at P/D stage, 15.7 3 106 sec- ondary spermatocytes, 107.7 3 106 round spermatids and 113.8 3 106 sperm. Somatic cells: Fibroblasts were grown until reaching confluence in supplemented DMEM medium. Cells were then washed with PBS and fixed with 1% formaldehyde for 10’ at RT. Cells were then incubated with glycine 0.125M for 50 at RT and for 15’ at 4 C to stop the crosslinking reaction. Following, 2 mL Trypsin (Fischer Scientific) were added, and cells were incubated for 8’ at RT and then washed twice with PBS. Cells were scraped and collected in a tube and centrifuged for 5 minutes at a maximum speed of 1,800 xg. Germ cells: Germ cells at a concentration of 1 million per 500 ml were incubated in formaldehyde (1%) for 10 min with agitation prior to FACS. Glycine (Sigma Aldrich) was added at a final concentration of 0.125 M and incubated with agitation at room temperature for e3 Cell Reports 28, 352–367.e1–e9, July 9, 2019 5 min and then at 4 C for 15 min. Cells were then centrifuged for 10 min at 290 xg at 4 C and resuspended in 3 mL 1x PBS in case of direct Hoechst staining, or in the according volume of block solution if immunofluorescence was the following step. Each crosslinked cell aliquot was resuspended with lysis buffer and incubated on ice for 30’ and then centrifuged for 50 at 1,800 xg. Pellets were washed with 1x NEB2 buffer (twice) and resuspended with NEB2 buffer with 10% SDS at RT and incubated for 10’ at 65 C with agitation (300 rpm). NEB2 buffer with 10% Triton X-100 solution was added and cells were incubated for 30’ at 37 C. Following, cells were centrifuged for 50 at 1,800 xg (4 C) and washed with 1x NEB2 buffer twice. An aliquot for ND (Non-Digested) control was taken from the sample to be processed and incubated at 37 C together with the digested sample. 400 U of MboI were added to the rest of the samples and chromatin was digested O/N at 37 C with agitation. An aliquot from the digested samples was taken for digestion controls. The full samples were kept at 37 C while the digestion control was performed. Proteinase K (10 mg/ml) (ThermoFisher) was added and the aliquoted samples were incubated for 45’-60’ at 65 C followed by Phenol:Chloroform purification. The quality of the sample was checked by running on an 0.8% agarose gel. After digestion, a small aliquot from the samples was kept as a Non-Ligated control and another aliquot was directly ligated without reparation (ligation control). Samples were centrifuged for 50 at 1,800 xg. Samples were then washed with 1x NEB2 buffer, twice. After the second wash, samples were directly resuspended with the reparation mix (1x NEB2 buffer, 0.05 mM dCTP, 0.05 mM dTTP, 0.05 mM, 50 mM biotin-dATP (ThermoFisher), 50U Kleenow (NEB, M02010M), H2O). Samples were incubated for 45’ at 37 C and for 10’ at 65 C and centrifuged for 50 at 1,800 xg and then resuspended with ligation buffer [1x NEB T4 ligase buffer, 10% Triton x-100, 0.1 mg/ml BSA, 5 ml (10000 U cohesive) ligase (2000 U/ ml) (NEB, M0201M), 963 ml H2O]. Samples were incubated at 16 C for at least 4h or O/N with mixing, then centrifuged for 50 at 1,800 xg and resuspended in 1x NEB2 buffer. Samples were incubated with RNaseA (10 mg/ml) (ThermoFisher) for 15 minutes at 37 C. The mix was incubated with Proteinase K (10 mg/ml) at 65 C O/N to reverse the cross-link. Samples were cooled to RT and purified with Phenol/Chloroform/Isoamyl alcohol. DNA content was measured on a Qubit. Samples were sonicated: 20 s time ON, 60 s time OFF, 8 cycles. Samples were then loaded in an electrophoresis gel of 1.2% agarose to check fragment size. Samples were incubated for 30’ with rotation at RT with Dynabeads MyOne Streptavidin T1 beads (Life Technologies) and 2x Bind- ing Buffer (10 mM TrisHCl, 1mM EDTA, 2M NaCl). Beads were washed twice with Binding Buffer and resuspended in the end repair mix [1x NEB T4 DNA ligase buffer with 10 mM ATP, 25 mM dNTP mix, 10U/ml NEB T4 PNK (NEB M0201), 3U/ml T4 DNA polymerase I (NEB M0203), 5U/ml NEB DNA polymerase I (Klenow) (NEB 0210)]. Samples were incubated for 30’ at RT. Beads were washed with 1x Binding Buffer, twice. Beads were subsequently resuspended in the dATP attachment master mix (1x NEBuffer 2, 0.5mM dATP, 5U/ml NEB Klenow exo minus (NEB, M0212), H2O). Samples were incubated at 37 C for 30’. The beads were washed with 1x Binding Buffer, twice and resuspended with 1x NEB Quick ligation buffer (NEB, B6058). NEB T4 DNA ligase and Illumina indexed adapters were added and thoroughly mixed. Samples were incubated for 15’ at RT, beads were captured, and the supernatant was discarded. Beads were then washed with 1x Binding Buffer, twice. The beads were resus- pended in Tris buffer (Elution buffer). PCR (5ml beads + 5.5 ml H2O, 1.25 ml primer P5 25 mM, 1.25 ml primer P7 25 mM, 12,5 cocktail master mix NEB) was performed (98 C 30 s, 98 C 10 s, 60 C 30 s x 8 cycles, 72 C 30 s, 72 C 50 ). After PCR, beads were captured on a magnet and the PCR was transferred to a new tube. The library was quantified using a Qubit fluorometer. 5ml from the PCR products were run in a 1.2% agarose electrophoresis gel to confirm range of sizes. A 1:1 amount of AmpureBeads (Beckman Coulter) was added to the samples and incubated for 10’ at RT. The beads were washed twice with 70% ethanol without mixing. The beads were then eluded and incubated at RT, for 50 . AmpureBeads were separated on a magnet, and the solution was transferred to a fresh tube. 5ml of each sample was loaded to an electrophoresis gel of 1.2% agarose for a final size check. DNA quantity was measured with the Qubit fluorometer. Libraries were submitted for Illumina sequencing (paired-end 75bp each side on HiSeq 2500, v4). Hi-C data processing, binning and normalization Quality check and trimming were conducted using BBDuk (version 10/2015) (Bushnell, 2014). Setting a minimum read length of 35 bp and a minimum Phred quality score of 20, adapters and low-quality reads were removed while preserving their longest high-quality regions. Then, reads were processed with TADbit (version 0.2.0.23) (Serra et al., 2017), which makes use of the GEM (version 1.7.1) mapper (Marco-Sola et al., 2012) to iteratively map them against the mouse genome (version mm10). Reads were mapped from 15 bp toward using a step size of 5 bp. The filters used in order to remove possible artifacts were the following: ‘‘self-circle,’’ ‘‘dangling- end,’’ ‘‘error,’’ ‘‘extra dangling-end,’’ ‘‘too short,’’ ‘‘too large,’’ ‘‘duplicated,’’ and ‘‘random breaks.’’ The maximum molecule length parameter was set at 2 times the 99.9 percentile of the insert size distribution, returned by the ‘‘insert_size’’ from TADbit. The maximum distance of a read to a cleavage site was set to the 99.9 percentile of the insert size distribution. An in-house script was used for binning and data normalization. This script imported the ‘‘HiC_data’’ module of TADbit, read the map files generated after the artifacts filtering step, binned the reads into a square matrix of 50 Kbp, and stored the matrix into a file in NPZ format (raw matrix). Afterward, HiCExplorer (version 1.8.1) (Ramı́rez et al., 2018) was used to normalize with the ICE (Iterative Correction and Eigenvector decomposition). Correlation coefficient analysis Pairwise comparisons between biological replicates were performed using HiCRep (version 1.4) (Yang et al., 2017), under a smooth- ing parameter of 5 and a considered distance over 10 Mbp. Since HiCRep only handles intra-chromosome raw matrices, each pair- wise comparison yielded 20 correlation scores (19 autosome chromosomes + sex chromosome X). The Y chromosome was Cell Reports 28, 352–367.e1–e9, July 9, 2019 e4 excluded from the analysis due to the lower number of interactions detected in our analysis (less than 1% of the overall detected interactions) and the highly repetitive nature of this chromosome. The correlation between 2 replicates was defined as the mean of the 20 correlation scores. Averaged contact probability P(s) ICE-normalized matrices were scaled with a factor of 1/sum(matrix). The resulting matrices were then input to ‘‘hicPlotDistVsCounts’’ from the HiCExplorer package in order to obtain the contact probability P(s). Inter-/intra-chromosome interaction ratio ICE-normalized data stored in matrices were exported with HiCExplorer to the GInteractions format, which consists of 7 columns: chromosome, start and end from bin 1, chromosome, start and end from bin 2, and the amount of interaction. The GInteractions tables were imported in R for further quantification of intra-chromosome and intra-chromosome interactions and plotting. Inter-subtelomeric interaction quantification ICE-normalized matrices were scaled with a factor of 1,000,000/sum(matrix) and exported with HiCExplorer to GInteractions format. The GInteractions tables were imported in R for this inter-telomere interaction quantification. Since the telomeric and centromeric regions (annotated from the beginning of each chromosome to 2.9 Mbp according to the UCSC Table Browser) were masked due to the low-count filtering step prior to ICE normalization, we only considered inter-chromosome interactions between loci located within genomic positions 3 to 3.5 Mbp in each chromosome. Differences in the subtelomeric interaction frequencies between cell types were assessed with the Wilcoxon test. Sperm simulations In order to validate our enriched sperm population, we simulated six Hi-C sperm datasets of 100 3 106 reads with different proportions, from 0 to 100% by steps of 20%, of fibroblast reads. Both sperm and fibroblasts reads were derived from our generated libraries. Previously published data on sperm (SRR3225862 and SRR3225863 accessions from Jung et al. (2017) were also down- loaded. These datasets underwent a quality check, Hi-C data processing, binning and normalization steps as described above. The resulting raw Hi-C matrices were used for correlation coefficient analysis while the ICE-normalized matrices were used to calculate the averaged contact probability P(s) (Figure S5). A-B compartments and TAD calling Raw matrices were used for the definition of A-B compartments. Columns with a low number of counts were filtered out using TADbit, setting the parameter min_count to 10. Since TADbit fits the column count distribution into a polynomial distribution, columns with a number of counts smaller than the first antimode of the distribution, which cannot be smaller than the min_count parameter, are filtered out. Then, the genome-wide matrices were normalized by the expected interactions at a given distance and by visibility by means of one-iteration of the ICE method. The correlation analysis was also performed with TADbit. In-house scripts computed A-B compart- ments from the first eigenvector, using 0 as threshold to differentiate both compartments and the gene density to label them. TADs were identified using an in-house script that imported the ‘‘Chromosome’’ module of TADbit and added the raw and the ICE- normalized matrices of each chromosome separately. Filtered bins, due to low counts, were included in order to mask them when calling TADs. TAD insulation scores were obtained by first normalizing the different matrices for read depth in order for the scores to be comparable. Each matrix was then scaled to have 100M reads. Afterward, TAD insulation scores were obtained from the output given by the ‘‘hicFindTADs’’ program from HiCExplorer. Compartment switching BED files with a resolution of 50 Kbp were available from the compartments definition step. Each genomic bin of 50 Kbp had its cor- responding compartment attributed. Pairwise comparisons between cell types -genome-wide and per-chromosome- were per- formed; the ratio of compartment switching was calculated as the number of genomic bins with a compartment change (A > B or B > A) divided by the total number of bins. From these files, a matrix file was created with 50 kbp-binned genomic coordinates as rows and cell types as columns, filled by the corresponding compartment labeling in each bin and cell type. Cell-specific A compart- ments were defined as those bins being compartment A in a cell type and compartment B in the remaining cell types. 3D-FISH and confocal microscopy A total of 10 commercially available Bacterial Artificial Chromosomes (BACs) (Source BioScience) from mouse chromosomes 12 and 14 were selected according to their genomic location in A or B compartments (5 BACs for each compartment) (Table S4). The selected probes were separated by a range of genomic distances ranging from 66 kbp and 1.45 Mbp. Probes were selected based on their mappability and repetitive content in the mm10 genome assembly using the UCSC browser (https://genome.ucsc.edu/). Probes were labeled with either dUTP-Dig (Sigma Aldrich) or dUTP-Cy3 (Enzo LifeSciences) by Nick Translation (Abbot Molecular). The 3D fluorescence in situ hybridization (3D-FISH) protocol was performed on mouse fibroblasts as following. Briefly, fibroblasts were cultured on slides overnight, whereas the germ cell suspension obtained after testis disaggregation was placed on slides e5 Cell Reports 28, 352–367.e1–e9, July 9, 2019 and incubated 1h at 37 C for cells to adhere on the surface. Slides were then fixed with 4% paraformaldehyde, washed in 1X PBS and then incubated in 60% Glycerol for 30-60 min, then finally snap-frozen in liquid nitrogen and kept in 50% formamide/2x SSC. As required, slides were thawed, treated with pepsin (0.01N HCl with 0.005% pepsin), washed with saline solutions (2x SSC, 50 mM MgCl2, 1x PBS) and post-fixed with 1% paraformaldehyde. After 1h incubation in 50% formamide/2x SSC, the hybridization solution was placed on the slides. Both slides and probes were simultaneously denatured at 75 C for 2 min and then incubated for 48h at 37 C. Slides were washed in 2xSCC at 37 C and in 0.1x SSC at 60 C and dUTP-Dig was detected followed incubation with anti- Dig FITC (1:150) for 45 min at 37 C. Finally, slides were washed with 4x SSC with 0.2% Tween-20 and mounted with DAPI diluted in Vectashield. Germ cells from four males were subject to the protocol simultaneously. Physical distances between hybridization signals were evaluated using a confocal microscopy (Leica SP5) using a 63x objective lens. Stacks of 0.34 mm-wide slices were captured for each cell nucleus, with a mean number of 20 stacks per nuclei. Images were reconstructed and analyzed using IMARIS (IMARIS Image Analysis Software), establishing nucleus volume (mm3), physical distances between signals (mm) and relative physical distance between signals and nucleus surface (mm). For the data analysis, pairwise distances were computed for particular nuclei. In total, we imaged 10 different loci (5 located in A compartments and 5 loci located in B compartment) across chromosomes 12 and 14. For each pair of probes, we measured the correlation between physical distances (in mm) and genomic distance (kbp) in pre-meiotic (fibroblasts and spermatogonia) and post-meiotic cells (round spermatids) (Figure S4; Table S4). A total of 527 measure- ments between seven pairs of probes were performed (Figure S6). We observed a positive correlation between physical and genomic distance (Spearman p value < 0.05) (Figure S4). But more importantly, we detected a differential pattern for both types of compart- ments (A and B) in round spermatids when compared to spermatogonia and fibroblasts (Figure S4). Lineal regressions were almost identical in spermatogonia and fibroblasts. In round spermatids, however, physical distances between pairs of probes were greater, confirming different chromatin organization in late spermatogenesis. Western blot analysis of CTCF and cohesins in germ cells Proteins from adult mice testis and fibroblast were extracted with RIPA buffer (NaCl 150 mM, Triton X-100 1%, Sodium deoxycholate 0.5%, SDS 0.1% and Tris-HCl 50 mM pH = 8.0). After protein quantification with the PierceTM BCA Protein Assay Kit (ThermoFisher Scientific), 6 replicates of both fibroblast and testis extract (three replicates of 30 mg and three of 40 mg protein) were denaturalized with 2x Laemmli and loaded into an 8% polyacrylamide gel until proteins reached the end of the gel [determined by the Precision Plus ProteinTM Dual Color Standards (Bio-Rad)]. Proteins were then transferred to a nitrocellulose membrane using the Trans-Blot TurboTM Transfer System from Bio-Rad (10 minutes of transfer time). Next, membranes were stained with the Ponceau S solution (Sigma-Aldrich) to ensure that proteins had correctly transferred to the membrane. Then, the membranes were washed and incubated with blocking solution (TBS 1x with 0.1% Tween-20, 5% fat-free milk and 1x PBS) for at least an hour. Each set of testis and fibroblast protein extracts were incubated overnight at 4 C with anti-rabbit CTCF (1:2500), anti-rabbit RAD21L (1:2000) and anti-rabbit REC8 (1:2000) respectively, and all were simultaneously incubated with anti-rabbit ßTubulin (1:5000) as a control. Anti- bodies were detected the following day with 1-hour incubation with anti-rabbit HRP-PO (1:15000). Finally, membranes were detected with ClarityTM Western ECL Substrate and the results captured with the Molecular Imager VersadocTM (Bio-Rad). Spermatocyte spreads and immunofluorescence analysis of CTCF and cohesins Spermatocyte spreads were obtained from frozen mouse testis. Testes were mechanically disaggregated until obtaining a cell suspension in 1x PBS. The cell suspension was then distributed into different slides and incubated with 1% Lypsol for 16 minutes followed by a 20-minute incubation with 4% paraformaldehyde. Then slides were left to dry and washed with twice PhotoFlo 1% (Kodak) and then blocked with PBS-Tween-20 (0.005%). Slides were incubated overnight at 4 C with the following antibodies: anti-rabbit CTCF (1:50), anti-rabbit RAD21L (1:20) and anti-rabbit REC8 (1:20); one per slide and all of them combined with anti- mouse SYCP3 (1:400). Primary antibodies were detected with anti-rabbit Cy3 (1:200) combined with anti-mouse FITC (1:200). Slides were finally mounted with DAPI and analyzed with a fluorescence microscopy (Axiophot, Zeiss) coupled with a ProgRes  CS10plus, Jenoptik camera. Representative images were captured with ACO XY (A. Coloma, Open Microscopy). Stimulated emission depletion (STED) microscopy Fab fragments were used for blocking and double labeling of REC8 and RAD21L antibodies (both raised in rabbit). STED microscopy (SP8, Leica) was used to generate the super-resolution images of REC8 and RAD21L foci along the chromosome axes. Secondary antibodies for STED imaging were conjugated to Alexa 555 and 488 (Invitrogen). Slides were mounted in Prolong Antifade Gold without DAPI. Fluorescence signals (red to green ratio) were measured along the 19 autosomal and XY axial elements of pachytene cells using the LAS X software from Leyca. Signal intensities were standardized and the overlay profiles of RAD21L and REC8 were plotted. Regression analysis was performed to determine the correlation between their profiles. The values of the coefficients of determination R2 are shown in the scatterplots. ChIP-sequencing For chromatin immunoprecipitation, antibodies for CTCF (10 ml per sample) and cohesins RAD21L and REC8 (30 ml per sample) were used. Two biological replicates of ChIP-sequencing (ChIP-seq) were performed using FACS-sorted primary spermatocytes at P/D Cell Reports 28, 352–367.e1–e9, July 9, 2019 e6 stage (14 million of cells) and round spermatids (20 million of cells) from adult mice. In short, cells were incubated on ice in lysis buffer I (5 mM PIPES, 85 mM KCl, 0.5% NP-40, Protease Inhibitors) and lysis buffer II (1% SDS, 10 mM EDTA, 50 mM Tris-HCl, Protease Inhibitors). Samples were then sonicated with the Biorruptor pico (30 s ON, 30 s OFF, 10 cycles) to obtain fragments around 200 bp. Sonicated lysates were centrifuged, and the supernatant was diluted in cold IP buffer (Diagenode). A small aliquot of each sample was kept as input and the remaining was divided in three aliquots, where antibodies for CTCF, RAD21L and REC8 were added (one in each aliquot) and then incubated overnight at 4 C. The antibody-chromatin pull-down was performed with Unblocked Protein A beads (Diagenode) and beads were then eluded with elution buffer (1% SDS, 0.1M NaHCO3). After elution samples were centrifuged, and the supernatant was incubated with 200 mM NaCl at 65 C overnight to reverse the crosslink. Finally, proteins were digested and the DNA purified with phenol:chloroform:isoamyl alcohol. Libraries were prepared for each sample, repairing fragment ends by incubating the samples for 30 minutes with the NEBNext end repair mix (New England Biolabs). Then samples were purified with AMPure beads. Samples were then incubated with the A-tailing NEBNext mix (New England Biolabs) and purified again with AMPure beads. Then, adaptors were ligated to the sample, which was subsequently purified with AMPure beads (0.8X). Libraries were PCR-enriched using NEBNext indexed primers, 12 PCR cycles for RAD21L and CTCF samples and 14 cycles for REC8 samples. Libraries were finally purified with AMpure beads (0.7X). ChIP-seq peak calling and annotation Quality check and trimming were performed using BBDuk. Setting a minimum read length of 35 bp and a minimum Phred quality score of 20, adapters and low-quality reads were removed while preserving their longest high-quality regions. Single end reads were obtained for CTCF ChIP and paired end reads for the cohesin ChIP. Reads were trimmed and mapped to the mouse reference genome (mm10) using a Galaxy server. After mapping, non-unique mapped reads with a mapping quality less than 30 were filtered out with SAMtools (Li et al., 2009). Then, to assess read coverage distribution across the genome, bigWig files for each sample were generated with DeepTools (Ramı́rez et al., 2016). Following the server recommendations, filtered files were merged with Picard (Broad Institute) and split again with SAMtools before peak calling. Peak calling was performed using MACS2 (Feng et al., 2012). To find the optimal parameters for running the peak call function, a cross-correlation analysis between reads mapping to plus and minus strands was performed, and the d parameter was estimated for each sample. Peaks were called using as extension size the d estimated for each sample respectively. Coverage, reads and peaks were visualized with IGVtools (Thorvaldsdóttir et al., 2013). BED peak files were imported to Rstudio and annotated to the reference mouse genome (mm10) using the R/bioconductor package ChIPseeker (Yu et al., 2015) together with the knownGene table from the UCSC resources (TxDb.Mmusculu- s.UCSC.mm10.knownGene). The number of peaks overlapping between CTCF and cohesins and between cell types was obtained by using the bioconductor package ChIPpeakAnno (Zhu et al., 2010). ChIP-seq data was represented using the packages ggplot2 (Wickham, 2016) and karyploteR (Gel and Serra, 2017) in combination with ChIPseeker and ChIPpeakAnno. RNA-sequencing Full-length single-cell RNA sequencing libraries were prepared using the Smart-seq2 protocol (Picelli et al., 2014) with minor modifications. Pools containing between 20,000 and 40,000 cells were obtained by FACS-sorting for four cell types (spermatogonia, primary spermatocytes at P/D, round spermatids and sperm) from adult mice. Four independent biological replicated were included in the analysis. Briefly, cells were sorted into 1.5ml eppendorfs containing lysis buffer. Reverse transcription was performed using SuperScrpit II (Invitrogen) in the presence of oligo-dT30VN, template-switching oligonucleotides and betaine. The cDNA was amplified using the KAPA Hifi Hotstart ReadyMix (Kappa Biosystems), ISPCR primer and 16 cycles of amplification. Following purification with Agencourt Ampure XP beads (Beckmann Coulter), product size distribution and quantity were assessed on a Bioanalyzer using a High Sensitvity DNA Kit (Agilent Technologies). Amplified cDNA was fragmented using Nextera XT (Illumina) and amplified with indexed Nextera PCR primers. Products were purified twice with Agencourt Ampure XP beads and quantified again using a Bioanalyzer High Sensitivity DNA Kit. Sequencing of Nextera libraries was carried out on a HSeq2500 (Illumina) to obtain > 30 million pair ends reads per sample. RNA profiling The Artificial Intelligence RNA-seq Software as a Service (SaaS) platform (https://transcriptomics.cloud) was used to analyze RNA- seq data. AIR accepts raw next generation sequencing Illumina FastQ data as input. RNA-seq data was uploaded to the site and validated in order to automatically pair forward and reverse files (in case of paired-end samples) as well as to check its format and integrity. Quality was assessed using FastQC. A new analysis was defined from the ‘‘new analysis’’ screen, where the samples to be included were selected for analysis, along the reference genome (> 120,000 available genomes from NCBI, Ensembl and JGI). The analysis included quality trimming, Differential Gene Expression (DGE) followed by a Gene Ontology Enrichment Analysis (GOEA). Once the analysis was launched, bad quality reads were removed using BBDuk by setting a minimum length of 35 bp and a minimum Phred-quality score of 25. Afterward, high quality reads were mapped against the reference genome with STAR (Dobin et al., 2013) using the end-to-end alignment mode and gene expression quantification was performed with featureCounts (Liao et al., 2014). e7 Cell Reports 28, 352–367.e1–e9, July 9, 2019 The statistical analysis started by filtering lowly expressed genes using HTSFilter (Rau et al., 2013). Four statistical methods were then used for the identification of differentially expressed genes: DESeq2 (Love et al., 2014), edgeR (Robinson et al., 2010), EBSeq (Leng et al., 2013) and NOISeq (Tarazona et al., 2011). Data normalization was performed with the Trimmed Mean of M-values (TMM) method. Finally, GOEA was performed with in-house scripts based on hypergeometric tests (Tian et al., 2017). Multiple testing corrections controlling false positives from high-throughput experiments were also performed with the Benjamini-Hochberg method. The statistics section included: the Principal Component Analysis (PCA) clusterization of the samples, general plots for the interpre- tation of the experiment and several tabs in which the Differential Expressed Genes (DEG) and GOEA data is shown in different tables. Specifically, we selected the mouse genome (GRCm38) from the Ensembl release 89. NOISeq was used to explore the DEG due to the variability of biological replicates. Raw expression and FPKM (Fragments Per Kilobase Million) values were downloaded from AIR. The Principal Component Analysis (PCA) was performed using the noise correction function from the NOISeq package. The raw expression values were converted to Counts Per Million (CPM) using edgeR and, afterward, expression values for each cell type were averaged for further analyses. Analysis of correspondence between compartments and gene expression Cell-specific A-B compartments were intersected with BEDTools (version 2.26) against a BED file with the TSS of genes derived from the GRCm38 gene annotation from Ensembl (release 89). Genes in each compartment were grepped (Bash command) with the table of FPKM values downloaded from AIR (see above), producing the expression profiles represented as boxplots for each cell type and compartment. Statistical significance among pairwise comparisons was tested using the Wilcoxon test. For GOEA in cell-specific A-compartments we extracted gene Ontology (GO) terms from expressed genes (CPM > 1) located in A compartments. GOEA was performed following the approach used by AIR on the DEG. It consists of a series of hypergeometric tests carried out for each GO term thus identifying significant enriched GO terms relative to the expected genome background. P values were corrected with the Benjamini-Hochberg procedure to reduce false positives. We considered as significantly enriched the terms with corrected P values % 0.01. Results were displayed using bubble plots. Moreover, expressed genes with peaks on their promoters were classified according to cell type and protein (CTCF, RAD21L and REC8). Genes were then analyzed with PANTHER (Mi et al., 2017) using the online tool AmiGO (Carbon et al., 2008). Local insulation, cohesin occupancy and gene expression The GRCm38 gene annotation from Ensembl (release 89) was downloaded in GTF format. It was parsed with an in-house script to extract in BED format the promoter regions 2 kbp upstream the Transcriptional Start Site (TSS) for each gene. Afterward, ChIP-seq peaks were intersected with BEDTools to obtain gene lists with ChIP-seq peaks in their promoters. These gene lists were grepped (Bash command) with the downloaded table with the expression values from the AIR platform (see below). Protein peaks were associated with changes in the corresponding TAD insulation score or gene expression using the function meanInRegions of the R package regioneR (version 1.10) (Gel et al., 2016). This package generates an expected distribution from 10.000 permutations from the values (e.g., TAD insulation score or FPKM values) observed in random genomic locations, thus calculating the Z-score of the TAD insulation score or the gene expression observed on ChIP peaks. Z-scores from flanking regions (+/ 250 Kbp) of peaks were calculated with the function localZScore. Meta-border plots were created using ICE-normalized matrices normalsed by the number interactions expected at a given distance with ‘‘hicFindEnrichedContacts’’ program from the HiCExplorer package (parameters: ‘‘–method obs/exp,’’ ‘‘–perchr’’). Subsequently, sub-matrices of interaction counts 5 Mbp up and downstream of specific regions (e.g., TAD boundaries or protein peaks) were generated and averaged among them. The interaction counts underwent a log10 transformation and were plot with ‘‘hicPlotTADs’’ program from HiCExplorer. X chromosome evolutionary strata and ampliconic regions Mouse X strata were extrapolated from the human X chromosome using synteny information (https://www.ensembl.org/ Mus_musculus/Location/Synteny?r=X). Likewise, the boundaries of each amplicon in the current assembly were defined by flanking paralogs (Table S7). QUANTIFICATION AND STATISTICAL ANALYSIS All sequencing data was checked and trimmed with BBDuk (Bushnell, 2014). Hi-C data was processed with TADbit (version 0.2.0.23) to obtain raw interaction matrices, compartments and TADs. HiCExplorer was used to: (i) normalize raw interaction matrices, (ii) create interaction heatmaps, (iii) predict averaged contact probabilities P(s) and (iv) obtain TAD insulation scores. Correlation values between biological replicates were performed using HiCRep (Yang et al., 2017). Normalized interaction matrices were exported with HiCExplorer to the GInteractions format and imported in R to calculate the inter-intra-chromosome interaction ratio. Normalized interaction matrices were scaled with a factor of 1,000,000/sum(matrix) with a custom Python script, exported with HiCExplorer to GInteractions format and imported in R to create boxplots. Differences between cell types were assessed with the Wilcoxon test (p value < 0.05). Cell Reports 28, 352–367.e1–e9, July 9, 2019 e8 ChIP-seq data was processed using the Galaxy server. BigWig files for each sample were generated with DeepTools (Ramı́rez et al., 2016) and peaks were called using MACS2 (Feng et al., 2012). BED peak files were imported to R and annotated using the R/Bioconductor package ChIPseeker (Yu et al., 2015) and UCSC resources (TxDb.Mmusculus.UCSC.mm10.knownGene). ChIPpeakAnno (Zhu et al., 2010) was used for peak intersection. Changes in the insulation score on peaks were statistically assessed using the function meanInRegions of the R package regioneR (version 1.10) (Gel et al., 2016). RNA-seq data was processed with the Artificial Intelligence RNA-seq (AIR) Software as a Service (SaaS) platform (https:// transcriptomics.cloud). In this case, raw expression values were converted to Counts Per Million (CPM) using edgeR, and averaged for each cell type. FPKM values were averaged for each cell type and BEDTools (version 2.26) was used to intersect genomic coordinates to relate gene expression with compartments, insulation scores and ChIP-seq marks. Differences in gene expression in the intersected regions between cell types were assessed by means of Wilcoxon tests (p value < 0.05). DATA AND CODE AVAILABILITY The accession number for the Hi-C, RNA-seq and ChIP-seq data reported in this paper is GEO: GSE132054. e9 Cell Reports 28, 352–367.e1–e9, July 9, 2019