ARTICLE https://doi.org/10.1038/s41467-022-30296-6 OPEN 3D chromatin remodelling in the germ line modulates genome evolutionary plasticity Lucía Álvarez-González 1,2,7, Frances Burden3,7, Dadakhalandar Doddamani3,7, Roberto Malinverni4, Emma Leach 3, Cristina Marín-García1,2, Laia Marín-Gual 1,2, Albert Gubern1, Covadonga Vara 1,2, Andreu Paytuví-Gallart 1,2,5, Marcus Buschbeck 4,6, Peter J. I. Ellis 3 ✉, Marta Farré 3 ✉ & Aurora Ruiz-Herrera 1,2 ✉ 1234567890():,; Chromosome folding has profound impacts on gene regulation, whose evolutionary con- sequences are far from being understood. Here we explore the relationship between 3D chromatin remodelling in mouse germ cells and evolutionary changes in genome structure. Using a comprehensive integrative computational analysis, we (i) reconstruct seven ancestral rodent genomes analysing whole-genome sequences of 14 species representatives of the major phylogroups, (ii) detect lineage-specific chromosome rearrangements and (iii) identify the dynamics of the structural and epigenetic properties of evolutionary breakpoint regions (EBRs) throughout mouse spermatogenesis. Our results show that EBRs are devoid of pro- grammed meiotic DNA double-strand breaks (DSBs) and meiotic cohesins in primary sper- matocytes, but are associated in post-meiotic cells with sites of DNA damage and functional long-range interaction regions that recapitulate ancestral chromosomal configurations. Overall, we propose a model that integrates evolutionary genome reshuffling with DNA damage response mechanisms and the dynamic spatial genome organisation of germ cells. 1 Departament de Biologia Cel·lular, Fisiologia i Immunologia, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Spain. 2 Genome Integrity and Instability Group, Institut de Biotecnologia i Biomedicina, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Spain. 3 School of Bioscience, University of Kent, Canterbury, UK. 4 Cancer and Leukemia Epigenetics and Biology Program, Josep Carreras Leukaemia Research Institute (IJC), Campus ICO-GTP- UAB, Badalona, Spain. 5 Sequentia Biotech, Barcelona, Spain. 6 Program for Predictive and Personalized Medicine of Cancer, Germans Trias i Pujol Research Institute (PMPPC-IGTP), Badalona, Spain. 7These authors contributed equally: Lucía Álvarez-González, Frances Burden, Dadakhalandar Doddamani. ✉email: p.j.i.ellis@kent.ac.uk; m.Farre-Belmonte@kent.ac.uk; aurora.ruizherrera@uab.cat NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 U nveiling the genomic basis of speciation is a main research breaks) and DSBs generated by oxidative damage in mature priority in biology fuelled by the availability of an sperm25. Moreover, there is fine-tuning between chromatin unprecedented large number of genomic resources. remodelling, architectural proteins and cell-specific gene expres- Comparative genomics of both closely and distantly related sion during spermatogenesis12,16,26,27 that can affect the potential mammalian species have revealed that genomic regions impli- outcomes of genetic damage in the germ line. It is not known, cated in structural evolutionary changes are clustered in regions however, which of these sources contribute most to the formation more prone to break and reorganise1–4. In this context, evolu- of transmissible evolutionary chromosomal reorganisations. tionary breakpoint regions (EBRs) are considered genomic Here we investigate how 3D genome folding relates to the regions implicated in structural evolutionary changes that disrupt functional and the epigenetic features of evolutionary chromo- genomic syntenic regions1–4. In searching for the origin and somal reorganisations in the mouse germline. To do so, we (i) functional implications of this evolutionary rearrangements, reconstruct ancestral rodent genomes analysing whole-genome research has pointed to repetitive elements as possible drivers1,5,6, sequences of 14 rodent species representatives of the major while changes in gene expression caused by genome reshuffling phylogroups, (ii) detect lineage-specific chromosome rearrange- may provide a selective advantage through the development of ments and (iii) identify the dynamics of the structural and epi- new adaptive characters specific to mammalian lineages3,4,7–9. genetic properties of EBRs through mouse spermatogenesis by Given the diversity of factors associated with EBRs it is most applying integrative computational analyses. unlikely that sequence composition of genomes is solely We find that EBRs are located in chromatin environments that responsible for genomic instability during evolution, and that the become more accessible as meiosis progresses, especially in post- regulation of 3D genome folding is also a critical factor8,10–12. meiotic stages. Moreover, our results show that EBRs are devoid Mammalian genomes are packaged into a chromatin structure, of programmed meiotic DSBs and meiotic cohesins in primary the regulation of which depends on several superimposed layers spermatocytes but associated with functional long-range inter- of organisation, including chromosome territories within which action regions and sites of DNA damage in post-meiotic cells. chromatin is organised into compartments (open/closed), which Overall, we propose a model that integrates evolutionary genome in turn consist of topologically associated domains (TADs) and reshuffling with DNA damage response mechanisms and the DNA loops10,13,14. The characterisation of how chromatin con- dynamic spatial genome organisation of germ cells. As such, we formation and DNA-protein interactions have evolved during detect the presence of long-range interactions in spermatids mammalian diversification is providing a new interpretive recapitulating evolutionary syntenic associations present in the hypothesis on the mechanism(s) responsible for the origin of Muridae ancestor. Our results suggest 3D genome organisation of genome architecture and plasticity10,15,16. Distant loci within the post-meiotic cells (i.e., spermatids) to be a major contributor to genome interact in a regulatory manner during the cell the formation of transmissible evolutionary chromosomal cycle12,14,15, affecting their ultimate function thus providing reorganisations. grounds for exploring the dynamics of genome composition, the evolutionary relationships between species and, in the long run, speciation. This view has been unified by the ‘Integrative Break- Results age Model’, an interpretative evolutionary hypothesis which Patterns of chromosome rearrangements in rodents. To assess posits that the permissiveness of genomic regions to reorganise the evolutionary reshuffling of rodent genomes, we used can be determined by the higher-order chromatin structure10,11. DESCHRAMBLER28 with 14 Rodentia chromosome-level gen- As with any evolutionary change of state, chromosomal reor- ome assemblies and two outgroups (human and rabbit) (Sup- ganisations that originate in the germ line before meiosis (pro- plementary Table 1). Rodent species included in the study were liferating primordial germ cells, spermatogonia and oogonia), representatives of the major phylogroups, with diploid numbers during meiotic division (spermatocytes and oocytes) or in post- ranging from 22 to 72 chromosomes. We first defined the meiotic stages (i.e., round spermatids) can be transmitted to reconstructed ancestral chromosome fragments (RACFs) for subsequent generations. In such cases, chromosomal reorganisa- seven ancestors in the rodent lineage (Muridae, Eumuroidea, tions can reduce gene flow and potentially contribute to specia- Muroidea, Myodonta, Mouse Clade, Mouse Clade + Ctenohy- tion by the suppression of recombination in the reorganised stricia, and Rodentia) (Fig. 1 and Supplementary Fig 1). After regions between chromosomally different, but contiguous manual curation (see Methods), the final number of RACFs populations16–19. Low levels of recombination could lead to a varied from 26 in the Eumuroidea ancestor to 35 in the Mouse high divergence and fixation of new mutations in these regions20 clade and Myodonta ancestors (Table 1). Coverage of the mouse which, combined with the presence of genes related to species- genome ranged from 92.92% in the Rodentia ancestor up to specific evolutionary pressures may reinforce the adaptive value 97.71% in the Muridae ancestor (Table 1). of EBRs in the germ line8. Theoretical work has suggested that Given that our approach predicts RACFs that in some cases heritable rearrangements would occur in genomic regions that are may not represent full ancestral chromosomes, we estimated the accessible in germ cells and/or early totipotent developmental number of ancestral chromosomes as those RACFs ≥ 26 Mbp, the stages10, highlighting the existence of a constraining role of EBRs length of the smallest chromosome in the assemblies used29 (see in the germ line that needs further investigation. Methods). Our results predicted the same number of ancestral Different sources of potential genomic structural alterations can chromosomes in Eumuroidea (haploid complement, n = 24), one arise during spermatogenesis: (i) formation and repair through more in the Rodentia (n = 26) and Muroidea (n = 27) ancestors, homologous recombination (HR) and non-allelic homologous and one less in Muridae (n = 24) than previously reported using recombination (NAHR) of meiotically programmed double strand cytogenetics comparisons (Table 1). Despite these differences, our breaks (DSBs) catalysed by SPO11 in early stages of prophase I results recover the large majority of the syntenic associations in (i.e., primary spermatocytes in leptotene and pachytene stages)21, the Muridae ancestor previously reported using the mouse (ii) non-disjunction and meiotic drive in meiosis I and II22, (iii) genome (Mus musculus domesticus, MMU) as a reference. formation and repair through non-homologous DNA end joining Particularly, we were able to detect seven syntenic associations (NHEJ) or microhomology-mediated end joining (MMEJ) of (MMU7/19, MMU5/11, MMU16/11, MMU13/15, MMU2/13, DSBs generated in later stages of spermatogenesis (i.e., round MMU1/17, and MMU17/10), with the exception of MMU12/17. spermatids)23,24 and (iv) zygotic repair of SSBs (single strand The high resolution provided by our genomic comparative 2 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE A 19,11 Muridae Mouse Eumuroidea 7,10 16,30 Muroidea Rat 22,21 Chinese hamster Myodonta 22,32 Blind mole rat Mouse-clade 12,13 Lesser Egyptian jerboa Mouse-clade + Ctenohystricia American beaver Ords kangaroo rat 13,12 Long-tailed chinchilla Common degu Guinea pig Naked mole-rat Damara mole-rat Rodentia Ground squirrel Yellow-bellied marmot European rabbit Human 90 80 70 60 50 40 30 20 10 0 MYA B Rodentia Muroidea Mouse 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9 10 10 10 11 11 12 11 12 13 12 13 13 14 14 15 14 15 16 15 16 17 16 17 18 17 18 19 18 19 20 19 20 21 X 21 22 22 23 0 25 50 75 100 125 150 175 200 23 24 Chromosome length (Mb) 24 25 25 26 Un Un 1 5 9 13 17 21 25 X X 2 6 10 14 18 22 Un 0 25 50 75 100 125 150 175 0 25 50 75 100 125 150 175 200 3 7 11 15 19 23 X Chromosome length (Mb) Chromosome length (Mb) 4 8 12 16 20 24 Reference C Eumuroidea 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Un X 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Un X Muridae approach, which permitted the identification of small rearrange- syntenic associations, indicating a high degree of agreement ments allowed us to detect three syntenic associations in the between methodologies (Supplementary Fig 1). Muridae ancestor (MMU5/6, MMU11/17 and MMU5/12) not Using the mouse genome as a reference, our genomic approach described previously using cytogenetics (Supplementary Fig 1). allowed us to identify a total of 232 EBRs considering all lineages For Eumuroidea and Muroidea ancestors, our data recovers 7/8 (Table 1). EBR sizes varied between 2 bp to 3.7 Mbp with a NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 Fig. 1 Phylogenetic tree of rodent species compared and reconstructed ancestors. A Phylogenetic tree constructed from the pairwise divergence times between the house mouse and each of the rodent species allows reconstructing the karyotypes for seven ancestors from Rodentia to mouse. The chromosomal rearrangements between the ancestors are shown for each node: in blue the number of inversions, in red the number of inter-chromosomal rearrangements. Green coloured dots denote ancestral lineages; blue dots represent recent ancestors (Muridae, Eumoroidea, Muroidea and Myodonta, respectively). The red depicts mouse. B Reconstructed ancestral karyotypes for the ancestors Rodentia, Muroidea and mouse, coloured according to Rodentia RACFs. The smaller RACFs (less than 26 Mb) are shown as an unplaced chromosome (Un). C Pairwise comparison between the Eumuroidea and Muridae ancestors. Each ribbon represents the syntenic fragments between the ancestors, tilted ribbons indicate inversions. Syntenic fragments are coloured according to Rodentia RACFs. Abbreviations – MYA million years ago, Un unplaced fragments. Table 1 Summary statistics of the reconstructed ancestral chromosome fragments (RACFs) and evolutionary breakpoint regions (EBRs) in the seven rodent ancestors and mouse. Ancestor No. RACFs No. ancestral chr. FISHa Coverage of mouse genome (%) No. of EBRs Rate (EBR/my) Raw Curated Mouse – – – 100 54 2.58 Muridae 62 28 24 25 97.71 5 0.42 Eumuroidea 34 26 24 24 95.52 13 1.06 Muroidea 60 30 27 26 95.36 75 7.50 Myodonta 78 35 30 – 94.52 41 2.75 Mouse Clade 49 35 27 – 94.17 30 27.27 Mouse Clade + Ctenohystricia 59 33 25 – 93.66 10 5.00 Rodentia 53 31 26 25 92.92 4 0.40 aReference: 54. median size of 21.5 Kbp (Supplementary Fig 2). These EBRs are consistent with the highest rearrangement rate in the clade located in gene-dense regions (permutation test based on 10,000 (Supplementary Fig 4). permutations, z-score 7.46, p value 0.001) and co-localised with Our ancestral reconstructions allowed us to characterise the transposable elements, particularly ERVs (permutation test based genomic regions that were maintained syntenic for 73 My on 10,000 permutations, z-score 4.57, p value 0.001, Supplemen- (million years) of rodent evolution (multispecies Homologous tary Table 2). A total of 2524 genes were located within EBRs, and Syntenic Blocks, msHSBs). We identified a total of 26 msHSBs enriched in Gene Ontology (GO) terms related to immune covering 11.1% of the mouse genome, with a size ranging from response [134 genes, false discovery rate (FDR) 7.18E−36] and 513 Kbp to 45.68 Mbp with an average length of 11.61 Mbp. At chromatin silencing (11 genes, FDR 0.00925) (Supplementary least one msHSB was present in each mouse chromosome, except Fig 3). in MMU4, MMU5, MMU7, MMU9, and MMU14, indicating that Subsequently, we phylogenetically classified EBRs based on the these chromosomes are highly rearranged. A total of 5282 genes ancestral lineage in which they occurred, ranging from 4 to 75 were located within msHSBs, enriched in GO terms related to EBRs in the Rodentia and Muroidea ancestor, respectively mitochondrial transport (190 genes, FDR 0.006) and detection of (Table 1). Using the identified EBRs, we then calculated rates of stimulus (1372 genes, FDR 9.7e−24) (Supplementary Fig. 5). rearrangement as the number of EBRs per branch length between ancestors in million years (My), ranging from 0.40 to 27.27 EBRs/ My in the ancestor of all rodents and the mouse clade ancestor, Epigenetic landscape dynamics during mouse spermatogenesis. respectively (Table 1 and Fig. 1). The mouse clade ancestor To investigate the relationship between EBRs and chromatin showed a higher chromosomal rearrangement rate than the structure of mouse male germ cells, we first analysed the land- expected average rate of chromosomal rearrangements (3.18 scape of the higher-order chromatin organisation during sper- EBRs/My, χ2 test, p value <0.001). EBRs were further classified matogenesis (Fig. 2A). To that aim we re-analysed available based on whether they demarcate inversion (N = 128) or non- chromatin epigenetic data30 for spermatogonia (pre-meiotic), inversion (N = 104) EBRs, if they are associated with a fusion or primary spermatocytes (meiotic) and round spermatids (post- fission, and finally grouped into ancestral (N = 44), recent meiotic) using three histone marks: active associated histone (N = 134) or mouse-specific EBRs (N = 54), whether they modifications H3K4me3 and H3K27ac, and one repressed asso- occurred before or after the split of Myodonta, or are uniquely ciated histone modification (H3K27me3) (Supplementary present in the mouse genome (Table 1). Table 3). In order to avoid methodological bias, we used available We then determined the number and type of chromosomal ChIP-seq datasets for all three cell types (spermatogonia, primary rearrangements present at each ancestor leading from the spermatocytes and round spermatids) retrieved from the same Rodentia ancestor to mouse. A total of 240 chromosomal experimental study (see Methods). The percentage of the genome rearrangements comprised of inter-chromosomal (fissions and covered by each histone varied slightly between the cell types. fusions) (N = 129) and inversions (N = 111) were identified H3K4me3 ranged from 1.3 to 4.6% in spermatogonia and round across seven different ancestors (Fig. 1 and Table 2). Both inter spermatids, H3K27me3 from 1.7% to 3.8% while H3K27ac from and intra-chromosomal rearrangements occurred in all ancestral 1.1% to 1.3% of the mouse genome, respectively. lineages. Seven inversions and 10 fusion/fissions occurred Using ChromHMM, we created a chromatin state model at between Eumuroidea and Muridae ancestor; while 22 inversions 200 bp resolution defined by eight different chromatin states and 32 fusion/fissions characterised the evolutionary reshuffling (from E1 to E8) in spermatogonia, primary spermatocytes and between the Mouse clade ancestor and Myodonta (Table 2), round spermatids (Fig. 2B and Supplementary Table 4). States E1, 4 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE Table 2 Number and type of chromosome rearrangements between reconstructed rodent ancestors. Ancestor Inversions Fission/Fusion Conserved chr. Time (my) Muridae → Mouse 19 11 13 20.9 Eumuroidea → Muridae 7 10 14 11.8 Muroidea → Eumuroidea 16 30 6 12.3 Myodonta → Muroidea 22 21 13 10.0 Mouse clade → Myodonta 22 32 9 14.9 Mouse + Ctenohystricia → Mouse clade 12 13 14 1.1 Rodentia → Mouse + Ctenohystricia 13 12 15 2.0 E3 and E7 were dominant, covering 48.9%, 78.5% and 82.3% of followed by E5-E0-E0 (6% coverage). During spermatogenesis, the mouse genome, in round spermatids, primary spermatocytes 25.8% of the genome became active, whereas only 4.49% and and spermatogonia, respectively. As all three chromatin states E1, 1.94% transitioned to repressed or poised states (Fig. 2C). In E3 and E7 had low coverage of histone marks they were classified contrast, a total of 2.56% of the mouse genome became trivalent as background state (E0) (Fig. 2B). Active state E2 (enriched in in spermatids. As expected, both X and Y chromosomes are H3K27ac), increased from a coverage of 1.1% of the genome in enriched in ‘closed’ chromatin states (Supplementary Fig 6) as spermatogonia to 26.4% in round spermatids (Fig. 2B). State E4, they are subjected to meiotic sex chromosome inactivation however, was dominated by the repressive chromatin mark (MSCI) during prophase I and post-meiotic sex chromatin H3K27me3, spanning between 0.31% and 5.11% of the genome in (PMSC) in round spermatids12,31. spermatogonia and round spermatids, respectively. As for poised chromatin (states E5 and E8 with both H3K27me3 and H3K4me3 marks) it covered from 13.7% in spermatogonia to 2.92% in EBRs associate with open chromatin environments activated round spermatids; while state E6 (labelled as trivalent chromatin during meiosis. We next integrated the genome positions of showing all three histone marks) covered 0.07% to 2.94% in EBRs with chromatin states and structural datasets including Hi- spermatogonia and spermatids (Fig. 2B). C data, CTCF, meiotic cohesins (REC8 and RAD21L), CpG To assess the dynamics of chromatin state transitions through- islands, transcription start sites (TSS), ATAC-seq and RNA-seq out spermatogenesis we compared the transition of chromatin (see Methods) (Fig. 2D). In particular, we analysed 3D genome states from spermatogonia to spermatocytes and then to round folding dynamics (A/B compartments and TADs) using pub- spermatids for a given genomic region (Fig. 2C). This allowed us lished Hi-C maps generated for spermatogonia, primary sper- to investigate which regions of the mouse genome get activated or matocytes and round spermatids12 and compared with the repressed during spermatogenesis. To do so, genomic regions at dynamics of the epigenetic landscape dynamics (previous sec- 200 bp resolution defined by six different chromatin states (E0, E2, tion). Moreover, CTCF and meiotic cohesins binding sites were E4, E5, E6 and E8) were compared between the three cell types included for primary spermatocytes and round spermatids12. (spermatogonia, primary spermatocytes and round spermatids). Genome-wide associations of different datasets were evaluated This gave us a total of 192 combinations of cell type and using permutation tests (see Methods). chromatin state, with 34 combinations covering overall >98% of We first studied the epigenetic landscape of EBRs, by assessing the genome with at least 0.1% each. The trivalent state with all the co-location of the 35 three-cell type states with the genomic three histone marks (E6-E6-E6) was also included, representing a positions of EBRs using a multi-association permutation test (see coverage of 0.029% of the mouse genome and making a total of 35 Methods). Remarkably, EBRs are negatively associated (permuta- combinations studied (Fig. 2C, Supplementary Fig 6 and tion test based on 10,000 permutation, normalised z-score −0.05, Supplementary Table 5). p < 0.05) with the background state (E0-E0-E0), but highly When comparing the percentage of each mouse chromosome associated (permutation test based on 10,000 permutation, coverage for each of the 34 most frequent chromatin state normalised z-score > 0.01, p < 0.05) with active or poised combinations, we detected that autosomes have different coverage chromatin (Fig. 3A and Supplementary Fig 7). This association of the three-cell type states than the sex chromosomes, with state was stronger with states that transition to E6 and E8 in E0-E0-E0 having the largest coverage ranging from 45.43% on spermatids (normalised z-score > 0.05, p < 0.05), particularly with chromosome 11 to 62.24% on chromosome 3 (Supplementary those EBRs that occurred in the mouse lineage, suggesting that Fig 6). For the sex chromosomes, state E0-E0-E2 had the largest EBRs occur in chromatin environments prone to rapid change coverage (44.71% on chromosome X and 50.76% on chromosome during spermatogenesis. Y). State E5-E0-E0 was lower on the sex chromosomes (X:1.43% To further investigate labile chromatin landscapes, we analysed and Y:0.31%) compared to an average of 6.51% on the autosomes. the gene content of those three cell type states associated with State E5-E0-E2 decreased on chromosome Y with 0.4% coverage EBRs. A total of 10,925 unique protein-coding genes were present compared to an average of 2.3% on the autosomes. Instead, state in E6 and E8 regions in spermatids. GO enrichment analysis E0-E0-E8 increased on chromosome Y at 3.44% compared to an (≥1.5-fold enrichment and FDR < 0.5) identified GO terms related average of 1.4% on the autosomes. to protein localisation to cell junction and protein depho- We detected that most of the genome (54.45%) remained in the sphorylation (1.79 and 1.55 fold) as well as dendrite development same background chromatin state (E0-E0-E0) throughout and regulation of organelle assembly (1.57 and 1.51 fold). spermatogenesis (Supplementary Fig 6). This contrasted with At the gross structural level, EBRs were associated with regions the small proportion of the genome that is maintained active that altered their state during spermatogenesis (all associations (0.084% in E2-E2-E2), poised (0.20%, E8-E8-E8 and E5-E5- based on multiple permutation test based on 10,000 permuta- E5 < 0.001%), trivalent (0.029% E6-E6-E6) or repressed (0.12%, tions, normalised z-score > 0.01, p < 0.05) (Fig. 3A and Supple- E4-E4-E4) in all three cell types. Remarkably, 41.09% of the mentary Fig 7). In particular, EBRs are associated with the ‘closed’ genome changed chromatin state during spermatogenesis, with B compartment in pre-meiotic spermatogonia, but with the ‘open’ E0-E0-E2 being the most common transition (21.9% coverage), A compartment in meiotic spermatocytes and post-meiotic NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 B C Chromatin States Genome coverage (%) 1 1.38 8.07 48.97 Background 2 1.1 2.23 26.37 Active 1x10 3 1.18 78.51 10.63 Background Frequency 4 0.31 4.44 5.11 Repressed 5 11.47 0.05 0.01 Poised 6 0.07 0.88 2.94 Trivalent 7 82.3 3.67 3.06 Background 5x10 8 2.19 2.15 2.91 Poised Chromatin States H3K27me3 H3K27ac H3K4me3 Spermatogonia Spermatocytes I Spermatids State 0 State 2 State 4 State 5 State 6 State 8 0 Spermatogonia Spermatocytes I Spermatids D Spermatogonia Spermatocytes I Spermatids Chromosome 13 Chromosome 13 Chromosome 13 56 58 60 62 64 Mbp 56 58 60 62 64 Mbp 56 58 60 62 64 Mbp 50 50 50 25 25 25 2 2 2 A/B A/B A/B State 0 State 4 State 6 State 0 State 4 State 6 State 0 State 4 State 6 Chromatin States State 2 State 5 State 8 Chromatin States State 2 State 5 State 8 Chromatin States State 2 State 5 State 8 H3K4me3 H3K4me3 H3K4me3 30 70 30 0 0 0 H3K27ac H3K27ac H3K27ac 10 10 10 0 0 0 H3K27me3 H3K27me3 H3K27me3 10 10 10 0 0 0 RNAseq-log FPKM RNAseq-log FPKM RNAseq-log FPKM 30 30 30 0 0 0 TSS TSS TSS CpG CpG CpG EBR EBR CTCF CTCF 30 30 0 0 REC8 REC8 30 30 0 0 RAD21L RAD21L 30 30 0 0 DMC1 DSB 10 10 0 0 PRDM9 ATAC-seq 10 20 0 0 EBR EBR ATAC-seq 25 0 EBR EBR spermatids (Fig. 3B and Supplementary Fig 7). Consistent with associated with TAD boundaries in spermatogonia and sperma- this, EBRs are associated with ‘closed’ chromatin environments tocytes, but located within TADs in round spermatids (Fig. 3B (E0, E4, E5) in spermatogonia and ‘open’ chromatin environ- and Supplementary Fig 7). Although EBRs were associated with ments (E2, E6, E8) in both primary spermatocytes and round TSS, these regions do not present high levels of expression in spermatids. At a finer structural level, too, EBRs are associated spermatogonia but are expressed more highly in round with regions that undergo structural remodelling, being spermatids (Mann–Whitney test, p < 2.2e−16) (Supplementary 6 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE Fig. 2 Epigenetic landscape dynamics during mouse spermatogenesis. A Schematic representation of mouse spermatogenesis. Adapted from49,70. Diploid (2n) and haploid (n) numbers are indicated for each cell type as well as the number of chromatids per chromosome (4c, 2c, or c). B ChromHMM chromatin states based on marks H3K27me3, H3K4me3 and H3K27ac. Numbers in the table indicate the percentage of genome coverage for chromatin states in the three cell types analysed (spermatogonia, primary spermatocytes and round spermatids). A total of 6 major chromatin states were found, including background (states 1, 3 and 7; grey), active (state 2; red), repressed (state 4; blue), poised (states 5; purple and 8; pink) and trivalent (state 6; yellow). C Alluvial plots representing chromatin states transitions from spermatogonia to primary spermatocytes and round spermatids. Chromatin states 1, 3 and 7 from panel (B) were merged into state 0 (background). D Chromosome 13 region-specific heatmaps at 50 kbp resolution (from 55 Mbp to 65 Mbp), for all three cell types depicting compartment signal (A, B), chromatin states, H3K4me3, H3K27ac, H3K27me3, RNA-seq (represented as log FPKM), CTCF and cohesin peaks (REC8 and RAD21L) and ATAC-seq. The genomic locations of EBRs are displayed (salmon highlight) in each cell type. Abbreviations – EBRs evolutionary breakpoint regions, FPKM fragments per kilobase of transcript per million fragments mapped. Fig 8). Overall, our results suggest that EBRs localise preferen- subjected to unique mutational pressures. These include DSBs tially in genomic regions that become accessible as spermatogen- formed (potentially as the result of topoisomerase II activity) to esis progresses. Furthermore, evolutionary rearrangements should relieve DNA helix torsion. As round spermatids lack a template not disrupt TAD structures in spermatogonia or spermatocytes for homology directed repair, DSBs must be repaired by error- but may do so in round spermatids. prone methods such as NHEJ or MMEJ23,37. Using publicly available data (see Methods), we identified a total of 151,732 DSB hotspots in spermatids (from now on post- Chromosome rearrangements do not disturb meiotic chro- meiotic DSBs), covering 1.49% of the mouse genome. Regions mosomal architecture in prophase I. Next, we investigated containing post-meiotic DSBs ranged from 146 bp to 6662 bp, whether the meiotic chromosomal architecture had an influence with a mean and a median of 267 bp and 213 bp, respectively. The on the distribution of EBRs in primary spermatocytes. To do so, it location of these post-meiotic DSBs in the mouse genome is important to consider that meiotic chromosomes are organised correlated with chromosome size (Pearson correlation, r2 = 0.68, into DNA loops anchored to chromosome axes formed by the p = 0.00074), with longer chromosomes accumulating more DSBs synaptonemal complex (SC) during prophase I, a proteinaceous than shorter ones, except in the chrY (3.51% coverage) structure with a zipper-like morphology that mediates the (Supplementary Fig 9). Moreover, post-meiotic DSBs were also synapsis of homologous chromosomes32. The SC establishes the associated with repeat content in each chromosome (Pearson context in which synapsis and recombination between homologs correlation, r2 = 0.78, p = 0.00003). Specifically, post-meiotic take place, joining sister chromatids within the lateral compo- DSBs co-localised with transposable elements L1Md_T and nents of the SC axis by meiotic cohesins, such as REC833 and L1Md_A (GAT analysis, 1000 permutations with 12.6- and RAD21L34. Importantly, while meiotic DSBs can occur at DNA 11.1-fold, p = 0.001). loops, they are repaired in the context of the chromosomal axes. At the epigenetic and structural level, we detected that post- As chromosomal axis length is inversely correlated with the size meiotic DSBs tend to occur in genomic regions depleted of of chromatin loops emerging from the SC35, the number and epigenetic modifications (state E0-E0-E0). That is, DSBs tend to distribution of DSBs per chromosome is linked to the structural be excluded from active or poised epigenetic landscapes in round organisation of the genome during meiosis22,36. spermatids (Supplementary Fig. 7). Also, post-meiotic DSBs were In searching for the evolutionary plasticity of meiotic negatively associated with A-compartments (multiple permuta- chromosomal architecture we conducted permutation tests (based tion test based on 10,000 permutations, normalised z-score on 10,000 permutations) to evaluate the association between −0.06, p < 0.05) and TAD boundaries in round spermatids EBRs and the genomic position of DMC1 and PRDM9 sites (normalised z-score −0.07, p < 0.05). Remarkably, post-meiotic together with different structural meiotic features such as meiotic DSB hotspots are positively correlated with EBRs (normalised cohesins (RAD21L and REC8) in primary spermatocytes. z-score 0.06, p < 0.05, Fig. 3B and Supplementary Fig. 7), Importantly, we detected that EBRs were negatively associated suggesting that EBRs are located in the subset of DSBs that with DMC1 and PRDM9 sites (Fig. 3B). EBRs were likewise occur in open chromatin in round spermatids. Therefore, our negatively associated with meiotic cohesins (RAD21L and REC8) results indicate that transmissible genomic rearrangements in both primary spermatocytes and round spermatids (Fig. 3A). preferentially occur within accessible genomic regions that suffer Even though EBRs correlated with ‘open’ chromatin environ- DNA damage in post-meiotic cells (Fig. 3D). ments in primary spermatocytes (i.e., E0-E6-E6 and E8-E8-E6) (Fig. 3A), these regions were not associated with cohesins (Supplementary Table 6). As cohesins are necessary structural Functional features of post-meiotic cell-specific long-range parts of the DNA loops attached to the chromosomal axes and genomic interactions. A refined analysis of chromosomal- the formation and repair of DSBs are tightly regulated during specific Hi-C maps revealed the existence of genomic regions meiosis our results suggest the presence of purifying selection for participating in cell-specific long-range intra-chromosomal the occurrence of large scale chromosomal reorganisations in the interactions in round spermatids that were absent in spermato- germ line (Fig. 3C, D). These results are in line with previous gonia and primary spermatocytes. We dubbed these regions observations of a reduction of recombination rates in EBRs3,8,18. ‘intra-LRIs’ (intra-chromosomal long-ranged interactions) and detected the presence of 36 intra-LRIs specific to round sper- matids distributed across chromosomes 1, 3, 4, 5, 7, 9, 12, 13, 17 EBRs associate with DNA damage in post-meiotic cells. As and 19 (Fig. 4 and Supplementary Table 7) ranging in size from evolutionary rearrangements were found preferentially associated 0.5 Mbp to 2.05 Mbp (mean size = 0.8 Mbp). The average dis- with genomic regions that become accessible as spermatogenesis tance between intra-LRIs was of 24.2 Mbp, and the interactions progresses, we subsequently analysed the occurrence of EBRs in were non-transitive, i.e. the presence of interactions A ↔ B and the context of post-meiotic cells (i.e., round spermatids). Sper- A ↔ C does not imply a direct interaction B ↔ C, and a single matids are haploid cells with highly compacted genome (most of “hub” intra-LRI could interact with two or more “spoke” intra- the histones are replaced by protamines), as such they are LRIs in the same chromosome (Supplementary Table 8). NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 A Spermatogonia Spermatocytes I Spermatids 1 EBRs Ancestral EBRs 0 Mouse EBRs Recent EBRs -1 E0-E0-E0 E0-E0-E4 E0-E0-E6 E0-E0-E8 E0-E2-E0 E0-E2-E2 E0-E2-E8 E0-E4-E0 E0-E6-E6 E0-E8-E6 E2-E0-E2 E4-E4-E4 E5-E0-E2 E5-E0-E4 E5-E0-E8 E5-E4-E0 E0-E0-E2 E0-E4-E2 E0-E4-E6 E0-E8-E8 E2-E0-E8 E2-E8-E6 E5-E0-E0 E5-E4-E4 E0-E4-E4 E5-E4-E2 E8-E0-E8 E8-E6-E6 E8-E8-E6 E8-E8-E8 E5-E4-E6 E5-E6-E6 E5-E8-E6 E6-E6-E6 E8-E0-E2 Chromatin state transitions B C Spermatogonia EBRs 0.04 Ancestral EBRs Ancestral EBRs Spermatogonia 0.02 Mouse EBRs TADs 0 Mouse EBRs -0.02 Recent EBRs Recent EBRs -0.04 TAD Boundaries TAD Boundaries - 5 Mbp 0 + 5 Mbp A Compartments A Compartments 0.04 Compartment Switch Compartment Switch 0.02 CpG Islands CpG Islands EBRs 0 TSS TSS -0.02 -0.04 - 5 Mbp 0 + 5 Mbp 1 0 -1 EBRs Spermatocytes I Ancestral EBRs Ancestral EBRs 0.04 Mouse EBRs Mouse EBRs 0.02 Spermatocytes I Recent EBRs Recent EBRs TADs 0 TAD Boundaries TAD Boundaries A Compartments A Compartments -0.02 Comparment Switch Comparment Switch -0.04 CpG Islands CpG Islands TSS TSS - 5 Mbp 0 + 5 Mbp CTCF CTCF RAD21L RAD21L REC8 REC8 0.04 Type I PRDM9 Type I PRDM9 0.02 Type II PRDM9 Type II PRDM9 EBRs 0 DMC1 DMC1 ATAC-seq ATAC-seq -0.02 -0.04 - 5 Mbp 0 + 5 Mbp 1 0 -1 EBRs Spermatids Ancestral EBRs Ancestral EBRs 0.04 Mouse EBRs Mouse EBRs 0.02 Spermatids Recent EBRs Recent EBRs TADs 0 TAD Boundaries TAD Boundaries -0.02 A Compartments A Compartments Compartments Switch Compartments Switch -0.04 CpG Islands CpG Islands - 5 Mbp 0 + 5 Mbp TSS TSS CTCF CTCF RAD21L RAD21L 0.04 REC8 REC8 0.02 DSB DSB EBRs 0 ATAC-seq ATAC-seq -0.02 -0.04 1 0 -1 - 5 Mbp 0 + 5 Mbp D Spermatocytes I Spermatids To explore the functional role of genes located within intra- contained in intra-LRIs were non-annotated transcripts. When LRIs in round spermatids, we analysed the GO terms contained analysing GO terms, we detected that intra-LRIs were enriched (p in these regions. Out of the total 691 genes located in intra-LRIs, value < 0.05 and enrichment score (ES) > 1.3, see Methods) for 41% (n = 283) were protein-coding genes, whereas non-coding genes related with sensory perception (ES = 13.6), lipid biosynth- RNAs represented 34% of the total, followed by pseudogenes esis and gluconeogenesis (ES = 6.27), oxido-reduction processes (20%) and lncRNAs (1.6%). The remaining 3.4% of the genes (ES = 6.97), response to cytokines (ES = 1.94) and phagocytosis 8 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE Fig. 3 Multi-comparison analysis. A Heatmaps obtained by regioneR (multicomparison) displaying correlations between different EBRs (ancestral, recent and mouse specific) and chromatin state transitions between chromatin states (E) in spermatogonia, primary spermatocytes and round spermatids. B Heatmaps obtained by regioneR (multicomparison) displaying correlations between different EBRs (ancestral, recent and mouse specific) and, TAD boundaries, A compartments, compartment switch (from A to B and viceversa), CpG islands, transcription start sites (TSS) in the three cell types. CTCF, cohesins (RAD21L and REC8) and ATAC-seq were included for both primary spermatocytes and round spermatids. Primary spermatocytes also included PRDM9 sites (Type I and II) and DMC1 sites. Round spermatids included post-meiotic DSBs. C Interaction metaplots depicting insulator scores for TADs and EBRs in each cell type. D Working model depicting the disposition of the genome folding (DNA loops and compartments) in relation to cohesins, CTCF, meiotic DSBs and EBRs. In the case of primary spermatocytes DNA loops protrude out of the chromosomal axes with meiotic DSBs occurring inside TADs in A compartments; EBRs are associated with TAD boundaries. In the case of round spermatids, EBRs are associated with post-meiotic DSBs inside TADs in A compartments. Abbreviations – EBRs evolutionary breakpoint regions, TADs topological associated domains, DSBs double strand breaks, FPKM fragments per kilobase of transcript per million fragments mapped. (ES = 1.74) (Supplementary Data 1). Additionally, 33% of these syntenies MMU7/19, MMU5/11, MMU17/11, MMU13/15, genes were members of relevant superfamilies, such as serpines, MMU2/13, MMU5/6, and MMU5/12 (Fig. 5). Similarly, seven zinc finger proteins, olfactory receptors and vomeronasal syntenic associations (MMU7/19, MMU5/11, MMU13/15, receptors. Considering differentially expressed genes (DEG) in MMU2/13, MMU12/17, MMU5/8 and MMU5/6) from Eumur- spermatids when compared to spermatogonia, primary sperma- oidea and 15 syntenic associations (MMU7/19, MMU2/13, tocytes and sperm, GO terms in intra-LRIs were narrowed down MMU11/17, MMU13/15, MMU12/17, MMU1/4, MMU2/4, to signal transduction, ubiquitinization and chemotaxis, all MMU1/8, MMU5/6, MMU3/5, MMU3/8, MMU17/18, MMU8/ important biological processes related to spermatogenesis (Fig. 4). 15, MMU5/11 and MMU12/14) from Muroidea were found to be Once intra-LRIs were characterised, we determined whether connected by LRIs in mouse spermatids (Supplementary Fig 10 these intra-LRIs also participated in long-range interactions and Supplementary Data 4). The percentage of overlap between genome-wide in round spermatids, denoting these as inter- LRIs and syntenic associations was reduced to 53 and 60% when chromosomal LRIs (inter-LRIs). This was done by selecting considering both Eumuroidea (7 out of 13) and Muroidea (15 out significant higher interactions (z-score interaction >3, see Methods) of 25) ancestors, respectively suggesting that long-range genomic of intra-LRIs with genomic regions located in other chromosomes interactions in round spermatids are associated with recent (genome-wide approach). As such, we detected the presence of ancestral chromosomal states (i.e., mouse-rat ancestor). 119 spermatid specific inter-LRIs involving different mouse When analysing into more detail LRIs present in Muridae chromosomes (i.e., multiple interactions) (Fig. 4 and Supplemen- syntenic associations, we distinguished different types of tary Data 2). Out of the total 864 genes contained in inter-LRIs, rearrangements, including single fusions, single fissions and 71% were protein-coding genes, 12% pseudogenes, 9% ncRNAs complex reorganisations (Fig. 5). Single fusions from the Muridae and 7% long non-coding RNAs. As for GO terms, we detected ancestor included LRIs in mouse chromosomes 1, 8, and 17, enrichment (p-value < 0.05 and ES > 1.3, see Methods) for genes whereas single fissions involved LRIs in ancestral syntenies associated with negative regulation of peptidase activity (ES = 2.8), MMU7/19, MMU5/6 and MMU5/11. On the other hand, LRIs in G-protein coupled receptor signalling pathway (ES = 2.4), phago- chromosomes 2, 13 and 15 were involved in complex reorganisa- cytosis (ES = 2.04), response to cytokines (ES = 1.81), complement tions, combining both fusions and fissions. Importantly, EBRs activation (ES = 1.68) and arachidonic acid metabolic process resulting from this evolutionary reshuffling were strongly (ES = 1.5) (Supplementary Data 3). In contrast to intra-LRI, where associated with LRIs (permutation test based on 10,000 11% (n = 77) of the containing genes were DEG, the percentage of permutations, z-score = 136.05, p < 0.001). Remarkably, the DEGs within inter- LRIs increased up to 48.5% (n = 419), being majority of intra-LRIs (75%) were located at the proximity of genes mainly related with cell surface receptor signalling pathway EBRs (<2 Mbp), with very few spanning the breakpoint. (ES = 1.66) and translation (ES = 1.39) (Fig. 4). Out of 11 fusions-fissions events detected in the Muridae In terms of structural configuration, intra-LRIs were signifi- ancestor, 8 were connected by LRIs in mouse round spermatids cantly positively associated with genomic regions that became (11 intra-LRIs and 7 inter-LRIs). These evolutionary LRIs accessible during meiosis (permutation test based on 10,000 contained a total of 144 genes (30% of them DEG expressed in permutations, z-score 78.5, p < 0.001) and were located inside round spermatids), including 83 (57.6%) protein-coding genes, 35 TADs (z-score 13.11, p < 0.001). Moreover, intra-LRIs were (24.3%) pseudogenes and 9 (6.3%) ncRNAs. The remaining 6.3% significantly negative associated with CpG islands (z-score were non-annotated transcripts. As for GO terms, we detected an −11.9, p < 0.001) and meiotic cohesins such as: REC8 (z-score enrichment (p-value < 0.05 and ES > 1.3) for genes associated −12.6, p < 0.001), RAD21 (z-score −8.8, p < 0.001) and CTCF (z- with response to pheromone (ES = 8.06), epoxygenase P450 score −3.5, p < 0.001). On the other hand, inter-LRIs were not pathway (ES = 7.12) and sensory perception of smell (ES = 2.75), significantly associated with any structural feature. including vomeronasal and olfactory receptors. Additionally, exocrine-gland-secreting peptide family genes (5%)38 and Pramel Long-range genomic interactions in round spermatids recapi- genes39 were only present in evolutionary LRIs when compared tulate ancestral chromosomal configurations. To understand to LRIs not involved in ancestral chromosomal configurations. the evolutionary implications of LRIs we analysed whether LRIs found in round spermatids were related to the evolutionary his- Discussion tory of ancestral chromosomes in rodents. To do so, we analysed Here, we provide evidence for a key role of the 3D chromatin the presence of LRIs in the syntenic associations present in recent organisation in the formation of transmissible chromosomal rodent ancestors, such as Muridae (mouse-rat), Eumuroidea reorganisations in the male germline. Using a comprehensive (mouse-rat-Chinese hamster) and Muroidea (mouse-rat-Chinese computational approach, we first reconstructed seven ancestors hamster-mole rat). From the ten syntenic associations found in in the rodent lineage with a high degree of accuracy (more than the Muridae ancestor, seven of them (70%) were now connected 92% of the mouse genome represented). We then show that by LRIs in mouse spermatids. These included the ancestral EBRs are associated with chromatin environments that become NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 accessible as meiosis progresses, especially in post-meiotic Our results highlight the importance of interpreting the dynamic stages of spermatogenesis (i.e., round spermatids), and which spatial genome organisation of germ cells in the context of germline are susceptible to DNA damage. Importantly, we also reveal the DNA damage responses (DDR) and meiotic checkpoints (Fig. 6). presence of post-meiotic cell specific long-range interactions, During gametogenesis, there are two key waves of DNA damage which not only recapitulate evolutionary syntenic associations induction that may trigger genomic rearrangements: (i) pro- present in the Muridae ancestor leading to fissions in mouse, grammed meiotic DSBs catalysed by SPO11 in early stages of male but also ancestral chromosome fusions that created new mouse and female prophase I21, and (ii) male-specific induction of DSBs in chromosomes. elongating spermatids to relieve torsional stress during chromatin 10 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE Fig. 4 Long-range interactions in post-meiotic cells. A Iterative correction and eigenvector decomposition (ICE)-corrected Hi-C matrices for chromosome 19 at a 100-kbp resolution for the cell types analysed. Deep blue lines indicate non-mapped bins. Long-rage interactions are highlighted (yellow square) in round spermatids. B Circos plot of the mouse genome depicting the different types of EBRs identified (Mouse, Muridae, Eumuroidea, Muroidea, Myodonta, Rodentia, Mouse Clade and Mouse Clade + Ctenohystricia specific) and long-rage interactions (LRIs). LRIs containing differentially expressed genes (DEG) are shown in green. C Representative example of intra-chromosomal long-range interactions across a specific region of chromosome 19 (from 14 to 45 Mbp) in mouse round spermatids. Different genomic features are displayed: A/B compartments (in spermatogonia and round spermatids), chromatin states, H3K4me3 peaks, H3K27ac peaks, H3K27me3 peaks, RNA-seq (represented as log FPKM), CTCF peaks and cohesin peaks (REC8 and RAD21L). The genomic locations of EBRs (salmon highlight) and intra-chromosomal long-range interactions (square) are displayed. D Zoom-in of intra-chromosomal long-range interactions detected in chromosome 19, displaying A/B compartments, chromatin states, RNA-seq (represented as log FPKM), and genes from NCBI Ref Seq annotation. Differentially expressed genes (DEG) are shown (red arrow head). E Zoom-in of intra-chromosomal long-range interactions detected in chromosome 7, displaying A/B compartments, chromatin states, RNA-seq (represented as log FPKM), and genes from NCBI Ref Seq annotation. Differentially expressed genes (DEG) are shown (red arrow head). F Gene Ontology Enrichment Analysis (GOEA) of genes found in intra- chromosomal and inter-chromosomal long-range interactions (LRIs) in round spermatids. Only significant gene ontology (GO) terms with an DAVID enrichment cluster score (ES) > 1.3 are shown, ES is based on a FDR adjusted P value (Fisher’s exact test, two-sided). Abbreviations – EBRs evolutionary breakpoint regions, LRIs long-ranged interactions, GO gene ontology, FPKM fragments per kilobase of transcript per million fragments mapped. Fig. 5 Long-range genomic interactions and ancestral chromosomal configurations. A Pairwise circus plot comparing the syntenic regions between the mouse karyotype and the Muridae ancestor karyotype (MAK). Inter-LRIs are depicted in black and intra-LRIs in red. Chromosomes are colour code according to the Muridae ancestor. B Upper panel - Schematic representation of the different types of evolutionary rearrangements from the Muridae ancestral karyotype including: single fusions, single fissions and complex reorganisations. Lower panel - Schematic representation of the 3D architecture of LRIs found in round spermatids according to the different types of evolutionary rearrangements. Chromosomes are colour code according to the Muridae ancestor. Abbreviations – LRIs long-ranged interactions, MAK Muridae ancestor karyotype, MMU Mus musculus domesticus, chr chromosome. condensation40. Here we used the genomic distribution of DMC1 strongly associated with male specific post-meiotic DNA damage together with PRDM9 binding sites as a proxy for meiotic DSB locations than with non sex specific meiotic DSB locations. locations (SPO11-oligos hotspots41). As both male and female The permissiveness of some genomic regions to undergo meiotic DSBs are mainly driven by the same mechanism (despite chromosomal breakage can also be associated with changes in differential efficiencies42), the DMC1 and PRDM9 data analysed, chromatin accessibility3,10,43. This has been recently demon- though generated in males, encompass the majority of female strated in somatic cells, where TAD boundaries experience meiotic DSB hotspots. Importantly, we detected that chromosome strengthening upon induced DSBs to accommodate accessibility rearrangements were not associated to meiotic DSBs and did not of the plethora of proteins involved in the DDR44. Likewise, it has disturb meiotic chromosomal architecture in prophase I. Although been shown that highly interacting regions are normally tran- it was previously shown that both SPO11 hotspots and H3K4me3 scriptionally active45 and that promoter-enhancer interactions marks positively correlated with cohesin occupancy12, A can be facilitated by chromatin remodelling46. This view aligns compartments12 and open chromatin states (this study) in primary with the growing evidence pointing to an association between spermatocytes, EBRs (irrespective of the rodent ancestor con- EBRs and TAD boundaries in mammalian and bird somatic sidered) were significantly devoid of DMC1 and PRDM9 sites. cells3,10,47,48. However, this poses a paradox in that rearrange- Conversely, EBRs were very strongly associated with the location of ments that occur at TAD boundaries are selectively favoured DSB hotpots in post-meiotic spermatids, consistent with the fact because they are less likely to disrupt gene regulation, but the that such DSBs must be repaired via a less accurate mechanism DNA strand breaks that initiate rearrangements are more likely to such as NHEJ or MMEJ since no sister chromatid is present23,37. occur at open chromatin within TADs. Here, we resolve this We therefore conclude that transmissible rearrangements are more paradox by showing that EBRs are selectively associated with NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 3D organization Meiosis cycle Sources of genomic instability Fig. 6 Working model on how chromatin structure of germs cells modulates genome evolutionary plasticity. Upper panel: 3D organisation of meiotic cells. Spermatogonia present a somatic-like organisation with the genome folded into compartments (A and B) and subsequent TADs. Subsequently, there is an attenuation of compartments and TADs in primary spermatocytes (here exemplified as leptotene, zygotene, pachytene and diplotene stages). Post- meiotic cells (round spermatids and spermatozoa) recover a somatic-like configuration, although with particularities such as that TAD borders are not clearly defined, and long-range interactions appear12,49. Middle panel: Meiotic progression and checkpoints, Self-renewing spermatogonia commit to meiosis. The primary spermatocytes meiotic prophase I, which is subdivided into four stages: leptotene, zygotene, pachytene and diplotene. The first and second meiotic divisions result in round spermatids, which differentiate into sperm70. Meiosis includes three meiotic checkpoints: (i) response to unrepaired double-strand breaks (DSBs), (ii) transcriptional repression called meiotic silencing of unsynapsed chromatin (MSUC), and (iii) the spindle assembly checkpoint (SAC). Lower panel: Sources of genomic instability. Spontaneous/de novo mutations and SINEs/LINEs activity predominate in spermatogonia. Meiotic programmed DSBs and its repair through non-allelic homologous recombination (NAHR) are sources of structural genomic changes during prophase I, together with non-disjunction during the subsequent meiotic divisions. The resulting structural changes are under strong selective constrains since errors can be detected by either of the meiotic checkpoints in play. In the case of post-meiotic cells (i.e., sperm) genomic instability can result from oxidative damage and DSBs produced by topoisomerase II DNA disentangling. The repair of DSBs is normally taken place by non- homologous end joining (NHEJ), an error-prone process that can result in structural genomic changes. As no checkpoints are activated in this stage (hence, relaxed selective constrains), there is a high likelihood that chromosomal reorganisations are transmitted to the offspring. Abbreviations – LRIs long-ranged interactions, TADs topological associated domains. regions that form ‘closed’ chromatin in pre-meiotic (spermato- range genomic interactions in mouse round spermatids con- gonia) or meiotic (primary spermatocytes) cell types and only taining genes relevant for spermatogenesis, fertilisation and become ‘open’ chromatin in spermatids. In particular, these sperm chemotaxis. Critically, relevant inter-chromosomal long- regions lie within TADs in spermatids but are near TAD range interactions now present in mouse round spermatids cor- boundaries in other cell types (Fig. 6). respond to genomic regions that were present in single chro- Crucially, our analysis of contact frequency maps during mosomes from the recent Muridae ancestor; while current intra- spermatogenesis detected the presence of cell-specific long-range chromosomal long-range interactions in the present day mouse interactions that can recapitulate ancestral chromosomal states. round spermatids correspond to regions that were present in two We have previously shown that chromosomal reorganisations or more chromosomes in the Muridae ancestor. have an impact on 3D genome topology of germ cells in two We therefore propose that the landscape of chromosomal ways: (i) altering chromosomal nuclear occupancy, and (i) reorganisations (fusions and fissions) that took place during reshaping landscapes of recombination16,49. Importantly, the genome evolution in rodents is linked to the chromatin context redistribution of chromosomal nuclear occupancy in spermato- now present in mouse round spermatids. When reorganisations cytes that result from chromosomal fusions brings new genomic occur within 3D contact hubs, this may either separate interacting regions into close proximity, predisposing to the occurrence of regions (converting intra- to inter-chromosomal LRIs) or con- additional rearrangements and expose chromosomal domains to versely bring them together. Alternatively, when chromosomal novel regulatory environments, potentially affecting gene fissions occur, new inter-chromosomal LRIs may be created in expression and/or regulation16. Here we demonstrate that this order to maintain critical associations between regions that were can be the case, as reflected by the presence of cell-specific long- formerly contiguous. This is in fact a crucial aspect that links the 12 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE dynamic spatial genome organisation of germ cells with evolu- boundary were labelled as unplaced. Plots of the RACFs were produced with the R tionary chromosomal reorganisations. package syntenyPlotteR4. Our model also integrates the role of meiotic checkpoints in genome evolutionary plasticity. Once a rearrangement has Detection of Evolutionary Breakpoint Regions (EBRs). The evolutionary breakpoint regions (EBRs) in each lineage were identified by considering the mouse occurred by breaking and rejoining of the DNA, it must evade genome as a reference. EBRs were counted from the coordinates of manually elimination by cell cycle checkpoints, and/or viability selection on merged RACFs. The minimum size shared across all lineages was considered as a the resulting offspring. Germ cells have a complex surveillance breakpoint. The EBRs were phylogenetically classified depending on the ancestral network including three major meiotic checkpoints: (i) response lineage in which they occurred. EBRs were further separated by the type of rear- to unrepaired DSBs, (ii) transcriptional repression called meiotic rangement they delimited into inversion EBRs or inter-chromosome EBRs, if they demarcated inversions or were the result of fusion or fission events, respectively. silencing of unsynapsed chromatin (MSUC), and (iii) the spindle assembly checkpoint (SAC)31. Therefore, any chromosomal rear- Histone, ATAC and DSBs profiles. A total of 27 fastq files were downloaded from rangement occurring before or during meiosis has a high prob- NCBI (Supplementary Table 3). Read quality was checked using FastQC ability of leading to meiotic arrest due to the activation of any of (v0.11.9)55. Reads were trimmed using Trimmomatic (v0.39)56, with SE/PE setting these three checkpoints (Fig. 6). Moreover, for inter-chromosomal used depending on the read type. Adapter sequences were removed using ILLU- translocations, even in the event that such a rearrangement occurs MINACLIP, as well as reads with an average Phred score of <30 or <20 with the AVGQUAL:20 or 30 for DSB and other libraries, respectively. MINLEN:30 was in pre-meiotic cells and is not eliminated by the meiotic check- also used for DSBs libraries. Trimmed fastq files were aligned to the mouse genome points, around 50% of the resulting gametes will be aneuploid and mm10 using bwa-mem (v0.7.17-r1188)57. Samtools merge was used to combine unlikely to lead to viable offspring. Thus, the vast majority of BAM files of the same histone mark. rearrangements occurring before or during meiosis will be elimi- Histone mark data was analysed using MACS2 (v.2.2.7.1)58 with either the default settings of callpeak to produce narrow peaks or with–broad–broad-cut-off nated. This view aligns with our observation of EBRs being devoid 0.05 to produce broad peaks. Histone marks were defined as narrow or broad based of programmed meiotic DSBs and meiotic cohesins in primary on the ENCODE project59. For marks not defined on ENCODE, the broad settings spermatocytes. Conversely, genomic rearrangements occurring in were used. ATAC-seq data was analysed using MACS2 (v.2.2.7.1)58 with the default post-meiotic cells are not subject to meiotic checkpoints and are settings of callpeak to produce narrow peaks. As for DSBs, data was analysed using the same parameters described in the original study24. Briefly, MACS2 with the moreover guaranteed to be euploid since each cell contains exactly default settings and–bw600 -q0.01–broad–broad-cutoff 0.1 were used. MACS2 files one haploid genome. Therefore, chromosomal reorganisations of the different round spermatid stages (stage 1–9 and stage 15–16) were merged occurring after meiosis have a higher probability to be transmitted for further analysis to facilitate cross-comparison with datasets not stratified into than genetic insults occurring before or during meiosis (Fig. 6). different developmental stages in spermatids. Longest DSBs peaks were visualised And this is in fact what we observe in our computational in IGV (v2.12)60 and co-locate with large stretches of alpha satellite regions. Because alpha satellites are not fully assembled in the centromeres of all mouse approach, as EBRs are strongly associated with genomic regions of chromosomes, to avoid bias towards the assembled stretches of satellite regions, we DNA damage in round spermatids. excluded these DSBs peaks from further analysis. Summarising, our observations suggest that chromatin remo- ChromHMM (v1.22) was used for chromatin state analysis applying the default delling during spermatogenesis represents an emerging frame- binsize of 200 bp with the concatenated strategy61. The corresponding cell type specific input was used as a control to adjust the binarization threshold locally. work where evolutionary genomic variation can be generated and Once binarization was completed the model was learned with varying numbers of transmitted to the offspring. Understanding how genome states, with 8 states (from E1 to E8) chosen as the optimum model. instability affects gene expression and regulation in the germ line To identify how chromatin states change during spermatogenesis, genomic will allow to further determine the effect of genome reshuffling on locations of states in each cell type (spermatogonia, spermatocytes and spermatids) were compared using bedtools intersect62. Regions of the genome missing a evolution and reproduction. chromHMM state in any cell type were removed. The dominant states in any of the three cell types (E7, E3 and E1) were merged into a joint state named E0. Genomic locations were then labelled according to the states in each cell type, and the Methods transitions from one chromatin sate to another were plotted using ggalluvial with Ancestral karyotype reconstruction. The genome assemblies of 14 Rodentia ggplot2 in R. Consecutive 200 bp regions with the same three-cell type state species representatives of the major phylogroups, assembled at chromosome-level combination were merged, and 34 combinations with more than 0.1% genomic or with scaffold N50 > 3 Mbp (Supplementary Table 1), and two mammalian coverage were identified. This cut off was chosen because the total coverage of all outgroup species (human and rabbit) were used to generate pairwise alignments these regions represented >98% of the mm10 genome. with the mouse genome (mm10) using LASTZ with default parameters. LASTZ alignments were converted into chain and net files using Kent toolbox utilities Hi-C data. Hi-C data from mouse spermatogonia, primary spermatocytes and using the parameters -minScore = 1000, -linearGap = medium, C = 0, E = 30, round spermatids was obtained from GEO:GSE13205412 and processed with K = 3000, L = 3000, O = 400. The Y chromosome was omitted due to the difficulty TADbit (v0.2.0.23)63 and HiCExplorer (v3.6)12,64. Contact matrices were built at in assembling it to a sufficient degree of quality, enrichment of repeats and 50 kb resolution and normalised to 100 M reads. palindromes in the chromosome50. The coverage of nets of each species was cal- culated against mm10 to minimise the potential fragmentation introduced into the Long-range interactions. Hi-C interactions matrices at 50 Kbp resolution were reconstruction of the ancestral karyotypes. exported to Ginteractions format using HiCExplorer (v3.6). Intra-chromosomal Reconstructed ancestral chromosome fragments (RACFs) were generated by Hi-C interactions matrices were used as input to define long-range interacting DESCHRAMBLER algorithm28 using a syntenic fragment resolution of 300Kbp (LRI) regions for each mouse chromosome. This was done by selecting 50 Kbp bins and a minimum adjacency score of 0.0001. Pairwise divergence times between the with interactions higher than 100 that occurred at bins separated by at least house mouse genome and each of the studied species were found using 10 Mbp. Out of the resulted bins false positives were eliminated by removing low TimeTree51. The divergence times between species were then used to write a mappability regions. The mappability of each region was calculated using GenMap phylogenetic tree in Newick format and visualised using FigTree52. We (v1.3)65, removing regions with mappability below 0.5. Inter-chromosomal long- reconstructed seven different ancestors in the rodent lineage: Muridae, range interactions were identified based on a z-score > 3. Eumuroidea, Muroidea, Myodonta, the ancestor for the mouse-related lineage, the To evaluate whether intra-chromosomal LRI regions significantly interact more ancestor of all rodents except squirrels, and the Rodentia ancestor. with other regions in the genome, we created 1000 sets of regions with the same The number of RACFs produced from DESCHRAMBLER was higher than the length and distribution of the intra-chromosomal LRI regions, using bedtools number of chromosomes suggested by previous studies53,54. Consequently, the shuffle v2.29.262. Genomic regions contained in the random 1000 sets were adjacent RACFs in each of the reconstructed ancestors were manually merged extracted automatically using and in-house R script from the 50 kb inter- using both reference genome and other reconstructed ancestors, which were most chromosomal Hi-C interactions spermatids matrix. The interaction value of all closely related. This process was started on the Muridae ancestor using the house these regions was then transformed into Z-score to address the significant mouse genome as a point of reference, before working back in evolutionary time differences. using the closest related ancestors as a point of reference. To determine the final number of ancestral chromosomes, we followed the model suggesting that there is a conserved boundary of chromosome-size variation in all mammals29. Average Gene ontology enrichment. Using BioMart v2.46.1 in R v4.0.3, unique protein chromosome length was calculated for all genomes assembled at chromosome level coding gene IDs in mm10 for each classification were identified. These were input and minimum size boundary was estimated to be 26 Mbp. RACFs smaller than this into the PANTHER db66 for the gene ontology (GO) enrichment. Statistical NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 overrepresentation test was selected with either GO biological process complete or 2. Larkin, D. M. et al. Breakpoint regions and homologous synteny blocks in PANTHER (v17.0) pathways. Only GO terms with ≥1.5-fold enrichment and chromosomes have different evolutionary histories. Genome Res. 19, 770–777 FDR < 0.5 were considered statistically significant. Plots were created with ggplot2 (2009). v3.3.5 in R. In the case of LRIs, where non-coding transcripts were detected, 3. Capilla, L. et al. Mammalian comparative genomics reveals genetic and DAVID (v6.9) was used for the GO enrichment considering an enrichment cluster epigenetic features associated with genome reshuffling in rodentia. Genome score (ES) > 1.3 as default parameter67. Biol. Evol. 8, 3703–3717 (2016). 4. Farré, M. et al. Evolution of gene regulation in ruminants differs between Differentially expressed genes. RNA-seq data from mouse spermatogonia, pri- evolutionary breakpoint regions and homologous synteny blocks. Genome Res. mary spermatocytes and round spermatids was obtained from GEO:GSE132054 29, 576–589 (2019). and processed using AIR (https://transcriptomics.sequentiabiotech.com/)12. The 5. Longo, M. S., Carone, D. M., Green, E. D., O’Neill, M. J. & O’Neill, R. J. resulting expression file was then analysed following the manual of the R package Distinct retroelement classes define evolutionary breakpoints demarcating for empirical analysis of gene expression data: edgeR (v3.13)68. Then, the edgeR sites of evolutionary novelty. BMC Genom. 10, 1–14 (2009). algorithm glmQLFTest with default parameters (p-value < 0.05; log2Fold- 6. Farré, M., Bosch, M., López-Giráldez, F., Ponsà, M. & Ruiz-Herrera, A. Change = ± 1) was used to determine the differentially expressed genes, using Assessing the role of tandem repeats in shaping the genomic architecture of spermatids as a reference. great apes. PLoS One 6, e27239 (2011). 7. Groenen, M. A. M. et al. Analyses of pig genomes provide insight into porcine Multi-association analyses. Statistical association between different genomic demography and evolution. Nature 491, 393–398 (2012). features was evaluated using the RegioneR R package version. 1.2669. Based on the 8. Ullastres, A., Farré, M., Capilla, L. & Ruiz-Herrera, A. Unraveling the effect of regioneR package we have created a series of functions to allow the calculation of genomic structural changes in the rhesus macaque - implications for the associations between multiple regionsets. Due to the implementation of multiple adaptive role of inversions. BMC Genom. 15, 1–13 (2014). comparisons, the p-value calculation was adjusted using Benjamini–Hochberg 9. Farré, M. et al. Novel insights into chromosome evolution in birds, archosaurs, procedure. The value of the z-score for association with an adjusted p-value greater and reptiles. Genome Biol. Evol. 8, 2442–2451 (2016). than 0.005, was considered as 0. The z-score calculated was subsequently normalised 10. Farré, M., Robinson, T. J. & Ruiz-Herrera, A. An Integrative Breakage Model by dividing it by the square root of n where n is the number of regions present in the of genome architecture, reshuffling and evolution: The Integrative Breakage permuted regionset. All the permutations were performed using randomizeRegions Model of genome evolution, a novel multidisciplinary hypothesis for the study and NumOverlaps respectively as randomisation and evaluation function. The of genome plasticity. BioEssays 37, 479–488 (2015). genomic positions of gene, LINEs, LTRs, and ERVs for mm10 genome were 11. Deakin, J. E. et al. Chromosomics: bridging the gap between genomes and downloaded from the UCSC browser and used as input for our analysis. chromosomes. Genes 10, 627 (2019). 12. Vara, C. et al. Three-dimensional genomic structure and cohesin occupancy correlate with transcriptional activity during spermatogenesis. Cell Rep. 28, e9 Quantification and statistical analysis. The statistical analyses were performed using R. Statistical parameters and tests are reported in the Figures and Figure (2019). Legends when necessary. 13. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 326, 289–293 (2009). 14. Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution Reporting summary. Further information on research design is available in the Nature reveals principles of chromatin looping. Cell 159, 1665–1680 (2014). Research Reporting Summary linked to this article. 15. Vietri Rudan, M. et al. Comparative Hi-C reveals that CTCF underlies evolution of chromosomal domain architecture. Cell Rep. 10, 1297–1309 (2015). 16. Vara, C. et al. The impact of chromosomal fusions on 3D genome folding and Data availability recombination in the germ line. Nat. Commun. 12, 1–17 (2021). The data that support this study are available from the corresponding authors upon 17. Noor, M. A. F., Grams, K. L., Bertucci, L. A. & Reiland, J. Chromosomal reasonable request. The Hi-C, RNA-seq and cohesion ChIP-seq datasets from mice used inversions and the reproductive isolation of species. Proc. Natl Acad. Sci. USA in this study were retrieved from the NCBI GEO repository: GSE132054 (Spermatogonia 98, 12084–12088 (2001). Hi-C: GSM3840080; Spermatocytes I Hi-C: GSM3840082 and Spermatids Hi-C: 18. Farré, M., Micheletti, D. & Ruiz-Herrera, A. Recombination rates and genomic GSM3840083; Spermatogonia RNA-seq: GSM3840094; Spermatocytes I RNA-seq: shuffling in human and chimpanzee - A new twist in the chromosomal GSM3840095; Spermatids RNA-seq: GSM3840096; Spermatocytes I CTCF, RAD21L, speciation theory. Mol. Biol. Evol. 30, 853–864 (2013). REC8: GSM3840086, GSM3840087, GSM3840088 and Spermatids CTCF, RAD21L, 19. Capilla, L. et al. Genetic recombination variation in wild Robertsonian mice: REC8: GSM3840089, GSM3840090, GSM3840091). on the role of chromosomal fusions and Prdm9 allelic background. Proc. R. DSBs data from mice used in this study were retrieved from the ENA repository under Soc. B Biol. Sci. 281, 1786 (2014). the accession code PRJEB20038. ChIP-Seq data for H3K4me3, H3K27me3 and H3K27ac 20. Faria, R., Johannesson, K., Butlin, R. K. & Westram, A. M. Evolving modifications from mice was retrieved from the NCBI GEO repository: GSE49624 inversions. Trends Ecol. Evol. 34, 239–248 (2019). (Spermatogonia H3K4me3, H3K27me3 and H3K27ac: GSM1202705, GSM1202708, 21. Keeney, S., Giroux, C. N. & Kleckner, N. Meiosis-specific DNA double-strand GSM1202713; Spermatocytes I H3K4me3, H3K27me3 and H3K27ac: GSM1202706, breaks are catalyzed by Spo11, a member of a widely conserved protein family. GSM1202709, GSM1202714; Spermatids H3K4me3, H3K27me3 and H3K27ac: Cell 88, 375–384 (1997). GSM1202707, GSM1202710, GSM1202715; Spermatogonia input DNA: GSM1202723; 22. Wang, S. et al. Inefficient crossover maturation underlies elevated aneuploidy Spermatocytes I input DNA: GSM1202724 and Spermatids input DNA: GSM1202725). in human female meiosis. Cell 168, e17 (2017). ATAC-seq data were retrieved from the NCBI GEO repository: GSE102954 23. Ahmed, E. A., Scherthan, H. & de Rooij, D. G. DNA double strand break (Spermatocytes I replicate 1 and 2: GSM2751129, GSM2751130 and Spermatids replicate response and limited repair capacity in mouse elongated spermatids. Int. J. 1 and 2: GSM2751133, GSM2751134). Source data are provided with this paper. Mol. Sci. 16, 29923–29935 (2015). 24. Grégoire, M. C. et al. The DNA double-strand “breakome” of mouse spermatids. Cell. Mol. Life Sci. 75, 2859–2872 (2018). Code availability 25. Wright, C., Milne, S. & Leeson, H. Sperm DNA damage caused by oxidative All functions for the multi-association analyses can be found at https://github.com/ stress: modifiable clinical, lifestyle and nutritional factors in male infertility. RMalinverni/GenoMatriXeR/ (v1.0.0, https://doi.org/10.5281/zenodo.6379726). The rest of Reprod. Biomed. Online 28, 684–703 (2014). functions used can be found at https://github.com/Farre-lab/EBRs_HiC_Spermatogenesis_ 26. Alavattam, K. G. et al. Attenuated chromatin compartmentalization in meiosis paper (v1, https://doi.org/10.5281/zenodo.6376838). and its maturation in sperm development. Nat. Struct. Mol. Biol. 26, 175–184 (2019). 27. Patel, L. et al. Dynamic reorganization of the genome shapes the Received: 22 November 2021; Accepted: 21 April 2022; recombination landscape in meiotic prophase. Nat. Struct. Mol. Biol. 26, 164–174 (2019). 28. Kim, J. et al. Reconstruction and evolutionary history of eutherian chromosomes. Proc. Natl Acad. Sci. USA 114, E5379–E5388 (2017). 29. Li, X. et al. Chromosome size in diploid eukaryotic species centers on the average length with a conserved boundary. Mol. Biol. Evol. 28, 1901–1911 (2011). References 30. Hammoud, S. S. et al. Chromatin and transcription transitions of mammalian 1. Ruiz-Herrera, A., Castresana, J. & Robinson, T. J. Is mammalian chromosomal adult germline stem cells and spermatogenesis. Cell Stem Cell 15, 239–253 evolution driven by regions of genome fragility? Genome Biol. 7, R115 (2006). (2014). 14 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30296-6 ARTICLE 31. Subramanian, V. V. & Hochwagen, A. The meiotic checkpoint network: step- 64. Wolff, J. et al. Galaxy HiCExplorer 3: a web server for reproducible Hi-C, by-step through meiotic prophase. Cold Spring Harb. Perspect. Biol. 6, a016675 capture Hi-C and single-cell Hi-C data analysis, quality control and (2014). visualization. Nucleic Acids Res. 48, W177–W184 (2020). 32. Zickler, D. & Kleckner, N. Recombination, pairing, and synapsis of homologs 65. Pockrandt, C., Alzamel, M., Iliopoulos, C. S. & Reinert, K. GenMap: ultra- during meiosis. Cold Spring Harb. Perspect. Biol. 7, 1–28 (2015). fast computation of genome mappability. Bioinformatics 36, 3687–3692 33. Xu, H., Beasley, M. D., Warren, W. D., van der Horst, G. T. J. & McKay, M. J. (2020). Absence of mouse REC8 cohesin promotes synapsis of sister chromatids in 66. Mi, H. et al. PANTHER version 16: a revised family classification, tree-based meiosis. Dev. Cell 8, 949–961 (2005). classification tool, enhancer regions and extensive API. Nucleic Acids Res. 49, 34. Llano, E. et al. Meiotic cohesin complexes are essential for the formation of the D394–D403 (2021). axial element in mice. J. Cell Biol. 197, 877–885 (2012). 67. Huang, D. W. et al. The DAVID gene functional classification tool: a novel 35. Zickler, D. & Kleckner, N. Meiotic chromosomes: integrating structure and biological module-centric algorithm to functionally analyze large gene lists. function. Annu. Rev. Genet. 33, 603–754 (1999). Genome Biol. 8, R183 (2007). 36. Ruiz-Herrera, A. et al. Recombination correlates with synaptonemal complex 68. McCarthy, D. J., Chen, Y. & Smyth, G. K. Differential expression analysis of length and chromatin loop size in bovids—insights into mammalian meiotic multifactor RNA-Seq experiments with respect to biological variation. Nucleic chromosomal organization. Chromosoma 126, 615–631 (2017). Acids Res. 40, 4288–4297 (2012). 37. Cavé, T., Desmarais, R., Lacombe-Burgoyne, C. & Boissonneault, G. Genetic 69. Gel, B. et al. RegioneR: an R/Bioconductor package for the association analysis instability and chromatin remodeling in spermatids. Genes 10, 40 (2019). of genomic regions based on permutation tests. Bioinformatics 32, 289–291 38. Kimoto, H. et al. Sex- and strain-specific expression and vomeronasal activity (2016). of mouse ESP family peptides. Curr. Biol. 17, 1879–1884 (2007). 70. Reig-Viader, R., Garcia-Caldés, M. & Ruiz-Herrera, A. Telomere homeostasis 39. Church, D. M. et al. Lineage-specific biology revealed by a finished genome in mammalian germ cells: a review. Chromosoma 125, 337–351 (2016). assembly of the mouse. PLoS Biol. 7, e1000112 (2009). 40. Laberge, R. M. & Boissonneault, G. Chromatin remodeling in spermatids: a sensitive step for the genetic integrity of the male gamete. Arch. Androl. 51, Acknowledgements 125–133 (2005). This work was supported by the Ministry of Economy, Industry and Competitiveness 41. Lange, J. et al. The landscape of mouse meiotic double-strand break formation, (CGL2017-83802-P to A.R-H.), the Spanish Ministry of Science and Innovation processing, and repair. Cell 167, e16 (2016). (PID2020-112557GB-I00 to A.R-H. and RTI2018-094005-B-I00 to M.B.), the Agència de 42. Brick, K. et al. Extensive sex differences at the initiation of genetic Gestió d’Ajuts Universitaris i de Recerca, AGAUR (SGR1215 to A.R-H.) and the recombination. Nature 561, 338–342 (2018). Leverhulme Trust (RPG-2019-414 194 to P.J.I.E. and M.F.). C.V. and L.A.-G. were 43. Carbone, L. et al. Evolutionary breakpoints in the gibbon suggest association supported by FPI predoctoral fellowships from the Ministry of Economy and Compe- between cytosine methylation and karyotype evolution. PLoS Genet 5, titiveness (BES-2015-072924 to C.V. and PRE-2018-083257 to L.A.-G.). L.M.-G. was e1000538 (2009). supported by an FPU predoctoral fellowship from the Ministry of Science, Innovation 44. Sanders, J. T. et al. Radiation-induced DNA damage and repair effects on 3D and University (FPU18/03867 to L.M.-G.). Part of the analyses was conducted using the genome organization. Nat. Commun. 11, 6178 (2020). specialist High Performance Computing systems provided by Information Services at the 45. Sobhy, H., Kumar, R., Lewerentz, J., Lizana, L. & Stenberg, P. Highly University of Kent. interacting regions of the human genome are enriched with enhancers and bound by DNA repair proteins. Sci. Rep. 9, 1–10 (2019). Author contributions 46. Cho, W. K. et al. Mediator and RNA polymerase II clusters associate in P.J.I.E., M.F. and A.R.-H. designed experiments. L.A.-G., F.B., D. D., R.M., E.L., C.M.-G., transcription-dependent condensates. Science 361, 412–415 (2018). L.M.-G., C.V., A. P.-G., A.G., P.J.I.E., M.F. and A.R.-H. analysed the data. M.B., P.J.I.E., 47. Lazar, N. H. et al. Epigenetic maintenance of topological domains in the M.F. and A.R.-H. contributed reagents, and resources. L.A.-G., F.B., D. D., P.J.I. E., M.F. highly rearranged gibbon genome. Genome Res. 28, 983–997 (2018). and A.R.-H. wrote the first draft of the paper with inputs from all authors. All authors 48. Li, J. et al. A new duck genome reveals conserved and convergently evolved read and approved the final version of the paper. chromosome architectures of birds and mammals. Gigascience 10, 1–15 (2021). 49. Vara, C. & Ruiz-Herrera, A. Unpacking chromatin remodelling in germ cells: implications for development and evolution. Trends Genet. 38, 422–425 Competing interests (2022). The authors declare no competing interests. 50. Tomaszkiewicz, M. et al. A time- and cost-effective strategy to sequence mammalian Y Chromosomes: an application to the de novo assembly of gorilla Y. Genome Res. 26, 530–540 (2016). Additional information Supplementary information The online version contains supplementary material 51. Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. TimeTree: a resource for available at https://doi.org/10.1038/s41467-022-30296-6. timelines, timetrees, and divergence. Times 34, 1812–1819 (2017). 52. Rambaut, A. FigTree, a graphical viewer of phylogenetic trees (University of Correspondence and requests for materials should be addressed to Peter J. I. Ellis, Edinburgh, 2009). Marta Farré or Aurora Ruiz-Herrera. 53. Engelbrecht, A., Dobigny, G. & Robinson, T. J. Further insights into the ancestral murine karyotype: the contribution of the Otomys – Mus Peer review information Nature Communications thanks the anonymous, reviewer(s) comparison using chromosome painting. Cytogenet. Genome Res. 112, for their contribution to the peer review of this work. 126–130 (2006). 54. Romanenko, S. A., Perelman, P. L., Trifonov, V. A. & Graphodatsky, A. S. Reprints and permission information is available at http://www.nature.com/reprints Chromosomal evolution in Rodentia. Heredity 108, 4–16 (2012). 55. Andrews, S. FastQC: a quality control tool for high throughput sequence data. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Babraham Bioinformatics. http://www.bioinformatics.babraham.ac.uk/ published maps and institutional affiliations. projects/fastqc/. 56. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014). 57. Li, H. Aligning sequence reads, clone sequences and assembly contigs with Open Access This article is licensed under a Creative Commons BWA-MEM. Preprint at https://arxiv.org/abs/1303.3997 (2013). Attribution 4.0 International License, which permits use, sharing, 58. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, adaptation, distribution and reproduction in any medium or format, as long as you give R137 (2008). appropriate credit to the original author(s) and the source, provide a link to the Creative 59. The ENCODE Project Consortium. An integrated encyclopedia of DNA Commons license, and indicate if changes were made. The images or other third party elements in the human genome. Nature 489, 57–74 (2012). material in this article are included in the article’s Creative Commons license, unless 60. Robinson, J. Integrated genomics viewer. Nat. Biotechnol. 29, 24–26 (2011). indicated otherwise in a credit line to the material. If material is not included in the 61. Ernst, J. & Kellis, M. Chromatin-state discovery and genome annotation with article’s Creative Commons license and your intended use is not permitted by statutory ChromHMM. Nat. Protoc. 12, 2478–2492 (2017). regulation or exceeds the permitted use, you will need to obtain permission directly from 62. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for the copyright holder. To view a copy of this license, visit http://creativecommons.org/ comparing genomic features. Bioinformatics 26, 841–842 (2010). licenses/by/4.0/. 63. Serra, F. et al. Automatic analysis and 3D-modelling of Hi-C data using TADbit reveals structural features of the fly chromatin colors. PLoS Comput. Biol. 13, e1005665 (2017). © The Author(s) 2022 NATURE COMMUNICATIONS | (2022)13:2608 | https://doi.org/10.1038/s41467-022-30296-6 | www.nature.com/naturecommunications 15