ARTICLE
                  https://doi.org/10.1038/s41467-021-23270-1                 OPEN

                  The impact of chromosomal fusions on 3D genome
                  folding and recombination in the germ line
                  Covadonga Vara 1,2, Andreu Paytuví-Gallart 1,2,3, Yasmina Cuartero4,5, Lucía Álvarez-González1,2,
                  Laia Marín-Gual 1,2, Francisca Garcia6, Beatriu Florit-Sabater1,2, Laia Capilla2,7, Rosa Ana Sanchéz-Guillén1,2,11,
                  Zaida Sarrate 1, Riccardo Aiese Cigliano3, Walter Sanseverino 3, Jeremy B. Searle8, Jacint Ventura7,
                  Marc A. Marti-Renom 4,5,9,10, François Le Dily 4,5 & Aurora Ruiz-Herrera 1,2 ✉
1234567890():,;




                  The spatial folding of chromosomes inside the nucleus has regulatory effects on gene
                  expression, yet the impact of genome reshuffling on this organization remains unclear. Here,
                  we take advantage of chromosome conformation capture in combination with single-
                  nucleotide polymorphism (SNP) genotyping and analysis of crossover events to study how
                  the higher-order chromatin organization and recombination landscapes are affected by
                  chromosomal fusions in the mammalian germ line. We demonstrate that chromosomal
                  fusions alter the nuclear architecture during meiosis, including an increased rate of hetero-
                  logous interactions in primary spermatocytes, and alterations in both chromosome synapsis
                  and axis length. These disturbances in topology were associated with changes in genomic
                  landscapes of recombination, resulting in detectable genomic footprints. Overall, we show
                  that chromosomal fusions impact the dynamic genome topology of germ cells in two ways:
                  (i) altering chromosomal nuclear occupancy and synapsis, and (ii) reshaping landscapes of
                  recombination.




                  1 Departament de Biologia Cellular, 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 Sequentia Biotech,
                  Barcelona, Spain. 4 Centre for Genomic Regulation, The Barcelona Institute for Science and Technology, Barcelona, Spain. 5 CNAG-CRG, Centre for Genomic
                  Regulation, The Barcelona Institute of Science and Technology, Barcelona, Spain. 6 Servei de Cultius Cel.lulars, Universitat Autònoma de Barcelona,
                  Cerdanyola del Vallès, Spain. 7 Departament de Biologia Animal, Biologia Vegetal i Ecologia, Universitat Autònoma de Barcelona, Cerdanyola del
                  Vallès, Spain. 8 Department of Ecology and Evolutionary Biology, Corson Hall, Cornell University, Ithaca, NY, USA. 9 Pompeu Fabra University,
                  Barcelona, Spain. 10 ICREA, Pg. Lluís Companys 23, Barcelona, Spain. 11Present address: Instituto de Ecología AC (INECOL), Red de Biología Evolutiva, Xalapa,
                  Veracruz, Mexico. ✉email: aurora.ruizherrera@uab.cat

                  NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                                   1
ARTICLE                                                                  NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1




H
         igher-order chromatin structure demarcates the limits                postmeiotic cells (i.e., round spermatids), and discuss the impli-
         of gene-regulatory domains1–3. Thus, disturbances of                 cations for evolution and fertility.
         domain architecture due to genome reshuffling (i.e.,
inversions, fusions, or indels) represent a nongradual change in
gene regulation because shifting of domain boundaries exposes                 Results
genes to novel regulatory environments4. Models5,6 and growing                Variation in recombination rates in natural populations of
experimental evidence7,8 suggest that indels and inversions can               house mice. We analyzed the variation in recombination rates in
alter interactions between contiguous topological associated                  wild-caught BRbS mice and the potential impact of Rb fusions on
domains (TADs), which can lead to oncogene activation, mor-                   these patterns. We conducted an integrative approach that
phological alterations, and novel gene functions. However, the                combined the cytological mapping of CO events directly in male
impact of balanced chromosomal rearrangements, such as                        germ cells (reflecting recombination at the Mbp scale) (Figs. 1, 2
Robertsonian (Rb) fusions9, on genome architecture and its                    and Supplementary Fig. 1) together with estimates of linkage
heritability are less explored. This is of particular relevance since         disequilibrium based on SNP genotyping (recombination at the
Rb fusions represent the most common chromosomal rearran-                     kbp scale) (Supplementary Fig. 2).
gement in nature (from plants to mammals)10, and are linked to                   We first experimentally determined the number (frequency) of
recurrent miscarriages, infertility, and aneuploid offspring in               COs along chromosomal axes in 45 wild mice by the
humans11. In fact, it has long been suggested, although not yet               immunodetection of the recombination protein MLH1 (a marker
empirically demonstrated at the genome level, that the presence               of COs) on pachytene chromosomes (Fig. 1a–c and Supplemen-
of chromosomal fusions in the germ line can alter segregation                 tary Table 1) (see “Methods”). We also included three laboratory
patterns (the so-called interchromosomal effect12).                           mice (strain C57BL/6J, BL6) for comparison. The CO survey in
   Germ cells are a unique cell model to test the genome-wide                 wild-caught mice included 15 all-acrocentric individuals (two
impact of chromosomal fusions—they have sequential develop-                   standard populations, 2n = 40) and 30 mice with Rb fusions
mental stages that involve dramatic and tightly regulated chro-               (three populations, 2n = 39–28), allowing for microscopic
mosomal movements and chromatin remodeling. These include                     visualization of a total of 1468 spermatocytes (Fig. 1b–d and
changes in intra-/interchromosomal interaction ratios, distance-              Supplementary Table 1). Of Rb mice, the population Sant Sadurní
dependent interaction frequencies, genomic compartments,                      d’Anoia (SS) was characterized by a low number (from one to
TADs, occupancy of insulator proteins (CTCF and cohesins), and                three) of Rb fusions (2n = 39–37), whereas the Castelldefels
gene expression13–16. The delicate fine-tuning between chromatin               (CAS) and Viladecans (VIL) populations have from four to six Rb
remodeling, architectural proteins, and cell-specific gene expres-             fusions (2n = 33–28) (Fig. 1d and Supplementary Table 1).
sion is crucial during the first meiotic prophase (prophase I)                    The observed average number of COs ranged from 20.16
when homologous chromosomes align, pair, synapse, and                         (±1.18) MLH1 foci per cell to 21.32 (±1.25) in wild-caught
recombine16.                                                                  standard mice, and from 18.13 (±1.78) to 21.82 (±2.21) in Rb mice
   Recombination has a dual role in sexual reproduction: (i) it               (Fig. 1d). A population-level analysis of COs showed that mice
assembles new combinations of allelic variants, contributing to               from Rb populations with four or more fusions (CAS and VIL
the maintenance of genetic diversity, and (ii) establishes physical           populations) showed greater interindividual variability regarding
associations between homologous chromosomes that enable                       the number of MLH1 foci numbers per cell (Fig. 1d). In fact, mean
faithful chromosomal segregation during meiosis. Importantly,                 numbers of COs per cell were positively correlated with diploid
recombination can be modulated not only by factors that control               numbers (Spearman, P < 0.001), and therefore negatively with the
the formation of meiotic crossovers (COs) during early meiosis                number of Rb fusions (Spearman, P < 0.0001) (Supplementary
(e.g., chromosome axis length is determined by chromatin loop                 Fig. 1a). Despite fewer CO events, Rb mice showed significantly
length17–19) but also by large-scale structural reorganizations that          more interindividual CO covariation than standard mice
can dramatically alter the genomic landscape20–22. Yet, the                   (Mann–Whitney test, P < 0.0001) (Supplementary Fig. 1b). We
impact of large-scale genome reshuffling (e.g., chromosomal                    detected statistical differences in mean numbers of COs per cell
fusions) on the three-dimensional genome topology in germ cells               among populations, including the wild-caught standard popula-
and its implications for recombination remain unknown.                        tions MON (Caldes de Montbuí) and BOI (Castellfollit del Boix)
   Here, we take advantage of chromosome conformation capture                 and the different Rb populations (Fig. 1d). Overall, mice with Rb
followed by deep sequencing (Hi-C) in combination with cyto-                  fusions (SS, VIL, and CAS populations) showed significantly lower
logical analysis of CO events and SNP genotyping to study how                 mean numbers of COs per cell (20.09 ± 1.88, 20.12 ± 2.16, and
genome folding and recombination landscapes are affected by                   19.74 ± 1.87) when compared to wild-caught standard (MON and
chromosomal fusions in the mammalian germ line. We analyzed                   BOI populations) (20.95 ± 1.21 and 20.56 ± 1.37) and BL6 mice
wild Western European house mice (Mus musculus domesticus)                    (21.67 ± 1.52) (Wilcoxon test, P < 0.05) (Fig. 1d), hence confirm-
from the northeast of the Iberian Peninsula, belonging to the so-             ing previous surveys conducted on different European Rb systems
called Barcelona Rb system (BRbS) characterized by a recent                   of the house mice21,25.
evolutionary origin of chromosomal fusions21,23,24. The standard                 Mirroring descriptions from other mammals18,19, both the
karyotype of M. m. domesticus consists of 40 acrocentric chro-                mean numbers of COs and meiotic double-strand breaks (DSBs)
mosomes, in contrast to BRbS mice, which have a variety of                    (here exemplified as RAD51 foci) were correlated in both
diploid numbers (2n) ranging from 2n = 27 to 2n = 4023. This                  standard (BL6 and wild-caught mice) and Rb mice (Spearman,
system is characterized by the presence of different Rb fusions               P < 0.05, Supplementary Fig. 1c). The variation in RAD51 foci
distributed in nongeographically coincident clines, leading to a              numbers between laboratory mice and BRbS mice is expected
progressive reduction in diploid numbers toward the center of the             based on CO counts, as these two processes are interrelated18,19.
range23. This natural model permits interrogation of the impact               As such, differences in DSB numbers between standard and Rb
that Rb fusions have on chromatin remodeling and fine-scale                    mice, together with the observation of a significantly high CO/
recombination in the germ line. In particular, we studied how                 DSB ratio in Rb mice (Supplementary Fig. 1d), suggest that a
chromosomal fusions alter the nuclear architecture at different               higher proportion of DSBs are resolved as COs when compared
hierarchical levels in meiotic (i.e., primary spermatocytes) and              to standard mice (Supplementary Fig. 1e).


2                 NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                      ARTICLE




Fig. 1 Genomic landscapes of recombination. a BRbS populations sampled. See Supplementary Table 1 for the population’s number assignment. The
diagram represents different types of chromosomal fusions. Chromosome type legend: Acr, all-acrocentric chromosomes of standard mice; Met Het, Rb
chromosomes in heterozygous state of Rb mice; Met Hom, Rb chromosomes in homozygous state of Rb mice. b Immunofluorescence of a spermatocyte at
pachytene stage from a Rb mouse (2n = 32): SYCP3 (green), centromeres (red), and DNA (DAPI, blue). Asterisks indicate Rb chromosomal fusions. Sex
chromosomes (XY) are indicated. Scale bar = 10 µm. Immunofluorescence replicates, n = 3. c Immunofluorescence of a spermatocyte at pachytene stage:
SYCP3 (green), MLH1 (red), and centromeres (blue). Asterisks indicate Rb chromosomal fusions. Sex chromosomes (XY) are indicated. Scale bar = 10 µm.
Immunofluorescence replicates, n = 3. d Boxplots depicting the number of MLH1 foci/cell per specimen represented (i) individually (colors correspond to
panel a) and (ii) per population. Boxplots are presented as mean values ± SD; center line, median; center diamond, mean. (i) Three laboratory mice (BL6)
are included for comparison. Boxplots indicating mean values and standard deviations are shown for each individual. BRbS populations: CAS Castelldefels,
BOI Castellfollit del Boix, MON Caldes de Montbuí, SS Sant Sadurní d’Anoia, VIL Viladecans. The diploid number of each mouse specimen is indicated on
top of each boxplot. Source data are provided as a Source Data file. (ii) Mean numbers of MLH1 foci/cell in laboratory mice (BL6) and wild-caught standard
(St) and Rb mice represented per population. P values (Wilcoxon’s rank-sum test followed by Dunn’s tests (two-sided) adjusted by Bonferroni, *P < 0.05)
represent differences between populations (BL6, laboratory mice (n = 55 cells); CAS (n = 296 cells), BOI (n = 174 cells), MON (n = 174 cells), SS (n = 206
cells), VIL (n = 605 cells)). Source data are provided as a Source Data file. e Percentage of chromosomal arms showing the different number (0, 1, 2, or 3)
of MLH1 foci per-chromosome (χ2 test, one-sided, *P < 0.05). The number of chromosome arms analyzed per arm type: n = 1046 for Acr BL6; n = 1140 for
Acr St (MON); n = 2233 for Acr Rb (VIL); n = 483 for Met Het (VIL); n = 1517 for Met Hom. Source data are provided as a Source Data file. f Distribution
of MLH1 foci along individual chromosomal arms with one (left panel) or two (right panel) MLH1 foci from panel e. The X axis represents the relative
positions on the chromosomal axes from the centromere (0%) to the distal telomere (100%) (χ2 test, one-sided, P < 0.05). The Y axis indicates the
frequency of MLH1 foci in each interval of chromosome arm length. Chromosome type legend: Acr St, all-acrocentric chromosomes of standard mice from
MON population (n = 1014 chromosomal arms with 1 CO; n = 120 chromosomal arms with 2 COs); Acr Rb, acrocentric chromosomes of Rb mice from VIL
population (n = 1232 chromosomal arms with 1 CO; n = 271 chromosomal arms with 2 COs); Met, arms of metacentric chromosomes of Rb mice from VIL
population (n = 1272 chromosomal arms with 1 CO; n = 112 chromosomal arms with 2 COs). Source data are provided as a Source Data file.




NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                                3
ARTICLE                                                                  NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1




Chromosomal fusions reshape genomic landscapes of recom-                      (iv) Rb chromosomes in heterozygous state of wild Rb mice (Met
bination. Since each pair of homologous chromosomes requires                  Het), and (v) Rb chromosomes in homozygous state of wild Rb
one CO to ensure proper chromosomal segregation (and hence                    mice (Met Hom) (Fig. 1e). In Rb mice from the VIL population,
fertility) during the first meiotic division (the so-called obligatory         our analysis indicated that irrespective of the chromosomal
CO26), we quantitatively assessed whether the number of COs per               complement (2n = 31 or 2n = 32), the CO distribution per-
chromosomal arm (0, 1, 2, or 3) was altered in Rb mice when                   chromosome arm had the same pattern (Supplementary Fig. 1f).
compared to wild-caught standard and laboratory mice. To that                 As a result, data from VIL mice were pooled together for sub-
aim and in order to reduce the intrinsic variation of the BRbS, we            sequent downstream analysis of CO distribution.
analyzed mice from one standard population (six mice from                        When considering the number of COs per-chromosome arm,
MON population, 2n = 40) and mice from one Rb population                      we detected significant differences between arm types (Acr BL6,
with a high number of chromosomal fusions (seven Rb mice from                 Acr St, Acr Rb, Met Het, and Met Hom) (χ2 test, P < 0.05)
VIL population, 2n = 31–32) (Fig. 1e, Supplementary Fig. 1f                   (Fig. 1e). The frequency of chromosome arms with zero CO was
and Supplementary Table 1). We also included three control                    low (≤3%) in standard mice (BL6 mice and wild-caught standard
BL6 mice. The chromosomes were assessed in five categories:                    mice), mirroring previous observations21. However, in Rb mice
(i) acrocentric chromosomes of laboratory mice (Acr BL6),                     with fusions in a homozygous state, we detected a statistically
(ii) acrocentric chromosomes of wild-caught standard mice (Acr                significant increased proportion of chromosomal arms with zero
St), (iii) acrocentric chromosomes of wild Rb mice (Acr Rb),                  CO (7.45%) (χ2 test, P < 0.05; Fig. 1e and Supplementary Fig. 1f).

4                 NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                     ARTICLE

Fig. 2 Effect of Rb fusions on recombination and synaptonemal complex length. a (i) Axis length (expressed as μm) analysis in standard (St) and
Robertsonian (Rb) mice from BRbS (see text and Supplementary Table 1 for further details) according to chromosome types (Dunn’s test, two-sided; P <
0.001; NS nonsignificant). Boxplots are presented as mean values ± SD; center line, median; center diamond, mean. Source data are provided as a Source
Data file. (ii) Analysis of CO density in the different chromosome types (Mann–Whitney test, two-sided; P < 0.0001). Boxplots are presented as mean
values ± SD; center line, median; center diamond, mean. Source data are provided as a Source Data file. Chromosome type legend for panels (i) and (ii):
Acr, all-acrocentric chromosomes of standard mice (n = 1140 chromosomal arms); Acr, Rb all-acrocentric chromosomes of Rb mice (n = 1711
chromosomal arms) Met Het, Rb chromosomes in heterozygous state of Rb mice (n = 396 chromosomal arms); Met Hom, Rb chromosomes in
homozygous state of Rb mice (n = 1294 chromosomal arms). b Immunofluorescence of primary spermatocytes of Rb mice at the pachynema stage,
labeling the synaptonemal complex with SYCP3 (green), the centromeres with CEN (blue), and MLH1 (red). Immunofluorescence replicates, n = 3. White
dashed circles: centromeric signals in heterozygous fusions. Red dashed: double-centromere signals in homozygous fusions. White arrowheads: Met Het
chromosomes. Red arrowhead: Met Hom with double-centromeric signals. Yellow arrowhead: Met Hom with a single centromeric signal. Scale bar =
10 µm. c Synapsis and recombination patterns found in mice with Rb fusions a heterozygous state (see Supplementary Table 1 for mouse codes). Left panel:
Percentage of heterozygous chromosomes according to synapsis pattern (synapsed, asynapsed, and open) for each of the Rb mice analyzed. Numbers on
the left refer to mouse ID (see Supplementary Table 1). Middle panel: Representation of the mean number of MLH1 per arm corresponding to the mice
analyzed. N = number of cells analyzed per individual. Right panel: Cumulative frequencies of MLH1 distributions in heterozygous metacentrics according
to synapsis pattern (synapsed, asynapsed, and open). Source data are provided as a Source Data file. d Double-centromeric signals and recombination. Left
panel: Distribution of MLH1 foci along individual chromosomal arms in homozygous metacentrics (Met Hom) with a single (1 CEN, blue) or double (2 CEN,
red) centromeric signal. The X axis represents the positions on the chromosomal axes from the centromeric end (black dot) to the distal telomere. The
Y axis indicates the frequency of MLH1 foci for each 10% interval of chromosomal length. Right panel: Cumulative frequencies of MLH1 distributions in
homozygous metacentrics with a single- or a double-centromeric signal. Source data are provided as a Source Data file. e Diagram depicting chromosomal
axis (SC) length, DNA loops, and crossovers (COs) observed in fused chromosomes at pachytene in Rb mice. Heterozygous metacentrics (Met Het)
chromosomes can occur in three states: synapsed, open, and asynapsed. COs are closer to the centromere in asynapsed chromosomes. Open
chromosomes have distal COs, whereas synapsed chromosomes can have more than two COs per arm and in interstitial-to-distal positions. Regardless of
the synapsis state, heterozygous metacentrics (Met Het) chromosomes have longer axes than homozygous metacentrics (Met Hom) chromosomes.
Conversely, DNA loop lengths in Met Het chromosomes are shorter than in Met Hom chromosomes.



Moreover, in Rb mice with homozygous fusions, we detected a                        As perturbations in obligatory CO frequencies and overall CO
significant decrease in the frequency of chromosomal arms with                   distribution can impair meiosis and affect reproductive fitness, we
two COs (Met Homo, 6.98%) (χ2 test, P < 0.05) when compared                     analyzed sperm viability. Rb mice had a significantly higher
with their nonfused counterparts (i.e., acrocentric chromosomes)                fraction of immobile sperm (66.1% ± 21.2 vs. 31.7% ± 8.8,
(Acr Rb, 18.23%) or Rb mice with heterozygous fusions (Met Het,                 Mann–Whitney test, P = 0.006) and increased sperm mortality
11.85%) (χ2 test, P < 0.05, Fig. 1e). These observations suggest that           (72.9% ± 18.5 vs. 32.1% ± 16.2, Mann–Whitney test, P = 0.005)
the reduced CO frequency observed in Rb mice is due to a                        when compared with standard mice.
reduction of chromosome arms with two COs, especially when                         Moreover, and in order to assess the fine-scale genomic impact
the Rb fusions are homozygous (Fig. 1e).                                        of both the overall reduction and chromosomal redistribution of
   We next examined how chromosomal CO distribution was                         COs, we analyzed the landscape of genomic diversity (i.e., number
affected by Rb fusions (Fig. 1f). We found significant differences               of alleles, observed/expected heterozygosity, and nucleotide diver-
in CO distribution between standard acrocentrics, Rb acro-                      sity, Supplementary Table 2) and divergence (expressed as FST
centrics, and Rb metacentrics (χ2 test, P < 0.05, Fig. 1f). CO                  values, Supplementary Tables 3 and 4) between populations
distribution along acrocentric chromosomes (i.e., nonfused) in                  together with estimates of recombination rates based on linkage
both standard (MON) and Rb mice (VIL) was distinct in distal                    disequilibrium (expressed as 4Ner/kbp) between standard and Rb
regions (90–100% of the chromosome length) (χ2 test, P< 0.05,                   mice (see “Methods”). Both principal component analysis (PCA)
Supplementary Fig. 1g). In the case of fused chromosomes, single                (Supplementary Fig. 2a) and estimations of population structure
CO events were slightly displaced (although not significant; χ2                  (Supplementary Fig. 2b) revealed that BRbS mice clustered
tests, P = 0.15) in homozygous fusions toward distal regions                    according to the geographical distribution of populations. More-
(80–90% of the chromosome arm length) when compared to                          over, and consistent with the cytological analysis, we observed an
heterozygous fusions (Supplementary Fig. 1g). On chromosome                     overall reduction of recombination rates (expressed as 4Ner/kbp)
arms with two COs, a bimodal distribution was detected in both                  in Rb mice when compared to standard mice (Mann–Whitney
standard and Rb mice, with differences between standard                         test, P < 0.001; Supplementary Fig. 2c, d). Consistent with the
acrocentrics, Rb acrocentrics, and Rb metacentrics (χ2 test, P <                variation in recombination rates observed in Rb mice, we also
0.05, Fig. 1f). In the case of acrocentrics (in both standard and Rb            detected differences in molecular diversity between populations
mice), we also detected differences in distal regions (90–100% of               (Supplementary Tables 2, 3, and 4). Although estimates of
the chromosome length) (χ2 test, P < 0.05, Supplementary Fig. 1g).              observed and expected heterozygosity were similar between
On fused chromosomes in a heterozygous state, the majority of                   populations, standard mice showed higher allelic richness than
COs were localized either in proximal (from 15 to 55% of the                    in Rb mice (Supplementary Table 2). Concomitant with the
chromosome arm length) or very distal (90–95% of the                            overall reduction of recombination rates, Rb mice showed higher
chromosome length) regions (Supplementary Fig. 1g). For                         FST estimates than standard mice, both when comparing
chromosomes in a homozygous state, proximal COs appeared                        divergence between standard and Rb populations and between
more displaced from the centromere (25–45% of the chromosome                    Rb populations (Tukey–Kramer test, P < 0.005) (Supplementary
arm length). Our results indicate that chromosomal fusions                      Tables 3 and 4).
reshape recombination landscapes by reducing both the overall
number of COs per-chromosome arm and the distribution of                        Modulation of chromosomal lengths and CO distribution.
recombination events along chromosomes.                                         Variation in both chromatin loop length and chromosomal axis


NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                              5
ARTICLE                                                                 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1


length can alter CO frequencies17–19; therefore, we sought to                reduction of COs per arm on homozygous Rb chromosomes
understand how these chromosomal features were affected                      when compared to heterozygous (Fig. 1e and Supplementary
by chromosomal fusions. We did not observe differences in the                Fig. 1f).
total axis length per cell in Rb mice when compared to standard                 Overall, we observed that Rb fusions had an effect on CO
mice (Supplementary Fig. 3a). However, when analyzing chro-                  chromosomal distribution and chromosomal axis length, depend-
mosomes in Rb mice according to chromosomal arm type, fused                  ing on whether they are in the heterozygous or homozygous state.
chromosome arm length was longer than the acrocentrics                       With consistent genome size, the longer chromosomal axes in
(Mann–Whitney test, P < 0.0001) (Fig. 2a and Supplementary                   heterozygous fusions would be attributed to shorter loops
Fig. 3a). Accordingly, fused chromosomes had significantly lower              (Fig. 2e). Despite the presence of different levels of asynapsis,
CO density than acrocentrics (Mann–Whitney test, P < 0.0001),                heterozygous fusions have longer axes allowing the formation of
irrespective of their state (homozygous or heterozygous) (Fig. 2a).          an obligatory CO per arm necessary for faithful chromosomal
Remarkably, heterozygous fused chromosomes were significantly                 segregation. The contrasting pattern found in homozygous
longer than when in a homozygous state (Wilcoxon test, P <                   fusions (short chromosomal axis, high frequency of arms with
0.0001) (Fig. 2a and Supplementary Fig. 3b). These observations              zero CO, and centromere interference) suggests that chromosome
suggest that Rb fusions affect both chromosomal axis length and              architecture could play a major role in reshaping CO distribution.
CO chromosome distribution. Likewise, CO formation is estab-
lished in a chromosome-specific manner, mirroring previous
observations in all-acrocentric mice17.                                      Chromosomal fusions reorganize spatial chromosome occu-
                                                                             pancy. We then analyzed whether Rb fusions impact the three-
                                                                             dimensional genome folding in germ cells. Using Fluorescence-
CO distribution in homozygous and heterozygous fusions.                      Activated Cell Sorting (FACS), we isolated highly enriched
Meiotic DSBs are repaired in the context of the chromosomal                  meiotic (primary spermatocytes at the pachytene/diplotene stage,
axes as homologous chromosomes pair and synapse27,28. Thus,                  P/D) and postmeiotic cells (round spermatids, RS) from wild-
disturbances of homologous pairing during prophase I are                     caught Rb mice (derived from the same population) with a high
expected to affect CO patterns. Concomitant with this view,                  number of fusions (Supplementary Table 1, Supplementary
perturbed pairing in Rb mice had an effect on CO distribution                Fig. 4a, b and “Methods” section). For each germ cell fraction, we
(Fig. 2b–e). We analyzed five Rb mice from the VIL population                 performed in situ Hi-C3,16 to generate genome-wide Hi-C maps
from which material was available, categorizing three different              for primary spermatocytes and round spermatids. The Rb somatic
pairing states observed in heterozygous fusions: (i) fully synapsed,         profile was derived from a wild-caught Rb male primary fibro-
(ii) open, and (iii) asynapsed chromosomes (Fig. 2b, c and Sup-              blast cell line belonging to the BRbS. Germ cell and somatic
plementary Fig. 3c). Disturbed synapsis influenced CO distribu-               data for all-acrocentric (standard) mice have been published
tion per-chromosome arm, while the average number of CO                      previously16. After filtering, an average of 254.9 million valid
events per arm was similar between individuals (Fig. 2c and                  interactions per cell type was detected (Supplementary Tables 5
Supplementary Fig. 3c). That is, while fully synapsed hetero-                and 6) with high correlation values between biological replicates
zygous fusions generally had one, homogeneously distributed, CO              (from 0.96 to 0.92) (Supplementary Fig. 4c and “Methods”
event per arm, asynapsed heterozygous fusions presented COs                  section).
either in intermediate positions (40–50% of the chromosome                      Genome-wide interaction maps confirmed the presence of six
length) or in distal regions (80–90% of the chromosome length)               Rb fusions (3.8, 4.14, 5.15, 6.10, 9.11, and 12.13) in all cell types
(Fig. 2c and Supplementary Fig. 3c). However, synaptic dis-                  analyzed from Rb mice (Fig. 3a and Supplementary Table 1). All
turbances in heterozygous chromosomes did neither significantly               chromosomes involved in Rb fusions showed higher interchro-
affect CO distribution nor reduce the number of COs per arm                  mosomal interaction values (interaction ratio ~3.0) in Rb mice,
when compared to acrocentrics (Fig. 1e), most probably due to                compared to standard mice, in all cell types (Fig. 3b and
variation in chromosomal axis length (Fig. 2a). Interestingly, we            Supplementary Table 1). However, we also detected different
observed that 28.57% of the CO events detected in asynapsed                  patterns of interchromosomal interaction ratios between non-
chromosomes were located at the border of synapsed and unsy-                 fused chromosomes. In fibroblasts, the few nonfused chromo-
napsed regions (40% of the chromosome length, Supplementary                  somes (e.g., chromosomes 1 and 2) had higher values of
Fig. 3c). These COs could prevent further asynapsis, as previously           heterologous interactions in Rb mice when compared to standard
suggested for hexavalent meiotic configurations in Rb mice29.                 mice (Fig. 3b, c and Supplementary Fig. 5). However, the spatial
   Both the reduced number of COs per arm and the                            genome architecture was highly reorganized in P/D, affecting
displacement away from the centromere observed in homozygous                 nearly all chromosomes (fused and nonfused) (Fig. 3d and
fusions compared to heterozygous (Fig. 1e and Supplementary                  Supplementary Fig. 6).
Fig. 1f) were unexpected based on the chromosomal speciation                    In P/D, all chromosomes showed high values of heterologous
theory30, thus, we sought to understand the mechanistic factors              interactions in Rb mice (interaction ratio ≥1.5) for all chromo-
behind this pattern. Remarkably, we detected that homozygous                 somes when compared to standard mice (Fig. 3b and Supple-
chromosomal fusions were associated with the presence of                     mentary Fig. 5), suggesting a genome-wide redistribution of the
double-centromeric signals at variable frequency among indivi-               spatial disposition of chromosomes inside nuclei. The presence of
duals (from 8.33 to 34.45%, Fig. 2d and Supplementary Fig. 3d).              heterologous associations was further demonstrated by the
We then investigated whether these double-centromeric signals                immunodetection of the centromeric constitutive heterochroma-
correlated with CO chromosomal distribution by analyzing mice                tin (exemplified as H3K9Me3 signals) on pachytene chromo-
from the VIL population with the same number of fusions (four                somes of Rb mice (Fig. 4a). In fact, heterologous interactions were
fusions in a homozygous state and one fusion in a heterozygous               dependent on the pairing state of Rb fusions (i.e., fully synapsed,
state). Notwithstanding cell-to-cell variability, the presence of a          open, or asynapsed) (Fig. 4b). Specifically, we observed more
double-centromere signal was associated with displacement of                 centromeric associations of acrocentric chromosomes when Rb
COs toward distal regions of the chromosome (Fig. 2d). This                  fusions showed open configurations (Fig. 4c). In addition, the sex
extended centromeric recombination suppression, coupled                      body (the X and the Y chromosome) was generally isolated from
with shorter chromosomal axes, likely explains the significant                autosomes (Fig. 4b). When heterozygous Rb fusions failed to fully

6                NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                               ARTICLE




synapse, there was a general disruption of interchromosome                      chromosomes (Fig. 4b). Interestingly, we detected full hetero-
associations, with an increased number of chromosomal associa-                  chromatinization of the sex body in >90% of cells with asynapsed
tions (Fig. 4c). Moreover, the sex body showed abnormal signals                 heterozygous Rb fusions (Fig. 4b), which could contribute to the
of heterochromatization being frequently associated with fused                  increased X/autosome heterologous interactions in P/D Rb mice


NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                      7
ARTICLE                                                                    NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1


Fig. 3 Effect of Rb fusions on the higher-order chromatin structure. a Genome-wide ICE-corrected heatmaps (500 kbp) for the cell types analyzed
(fibroblasts, pachynema/diplonema—P/D, and round spermatids—RS) in Rb mice. Chromosomes involved in Rb fusions emerge as regions with high
interchromosomal interaction in all cell types (arrowheads). Rb fusions are the following: 3.8, 4.14, 5.15, 6.10, 9.11, and 12.13. b Heatmaps depicting
genome-wide interchromosomal interaction ratio between standard (St) and Rb mice. Chromosomes in red indicate higher interactions in Rb than in
standard mice, whereas chromosomes in blue indicate higher interactions in standard than in Rb mice. As expected, all chromosomes involved in fusions
show high interaction ratios in all the three cell types analyzed. Source data are provided as a Source Data file. c Inter-/intrachromosome interaction ratio
for the cell types analyzed (fibroblasts, P/D, and RS) in St and Rb mice. In fibroblasts (upper panel), values are the same for chromosome 19 in both St and
Rb mice. Source data are provided as a Source Data file. d Examples of interaction patterns between chromosomes involved and not involved in fusions.
(i) Interaction heatmaps representing chromosomes 1 and 2 (not involved in fusions) and chromosomes 3 and 8 (fused in a heterozygous state) in both St
and Rb mice. In chromosomes 3 and 8, the fusion becomes evident in interaction maps from Rb mice, with high interaction in the pericentromeric region of
the chromosomes (0–3 Mbp from the centromere). The observed scaling is consistent across all chromosomes. Source data are provided as a Source Data
file. (ii) Interaction plots for fibroblasts, P/D, and RS for chromosomes not involved (1 and 2) and involved in fusions (3 and 8). The observed scaling is
consistent across all chromosomes. Source data are provided as a Source Data file. e Boxplot showing the number of genome-wide interactions between
the X chromosome and autosomes detected in St and Rb mice per each cell type (Mann–Whitney test, ****P < 0.0001, two-sided). Boxplots are presented
as mean values (center line) ± SD. Source data are provided as a Source Data file. f Boxplots depicting genome-wide interchromosomal interactions per
million at pericentromeric regions (from the centromere up to 3.5 Mbp) between St and Rb mice (Mann–Whitney test, ****P < 0.0001, two-sided).
Boxplots are presented as mean values (center line) ± SD. For each cell type, two groups of chromosomes were compared: chromosomes involved in Rb
fusions (3.8, 4.14, 5.15, 6.10, 9.11, and 12.13) and chromosomes not fused (1, 2, 7, 16, 17, 18, 19, and X). Source data are provided as a Source Data file.



(Fig. 3e). In contrast, heterologous interactions were reduced in               dependent on genomic distance, was switched at genomic dis-
RS for all chromosomes not involved in fusions (ratios < 1, Fig. 3b             tances of ~6 Mbp in fibroblasts. This was concomitant with slight
and Supplementary Fig. 5), a pattern confirmed by the analysis of                changes in A/B compartments (Fig. 5b), with Rb mice showing a
inter-/intrachromosomal interactions (Fig. 3c). In fibroblasts and               smaller mean compartment size (0.85 Mbp) than standard mice
P/D, Rb mice showed higher inter-/intrachromosomal interaction                  (1 Mbp).
ratios than all-acrocentric mice (Mann–Whitney test, P < 2.2e-                     Mirroring previous studies14,16, our genome-wide analysis
16), a pattern that was clearly reversed in RS with higher                      showed that most compartments were mostly lost in primary
interaction ratios in standard mice than in Rb mice for all                     spermatocytes (P/D) in Rb mice (Fig. 5b), when homologous
chromosomes (Mann–Whitney test, P < 2.2e-16) (Fig. 3c).                         chromosomes condense, align, pair, synapse, and recombine.
   Chromosome-specific interaction maps also revealed altered                    Consistent with this absence of compartments during prophase I,
three-dimensional chromosome folding in Rb mice, especially in                  eigenvector values were close to 0 (Fig. 5b) and inter-/
prophase I (Fig. 3d). For nonfused chromosomes (e.g., chromo-                   intrachromosomal interaction ratios reached a minimum for all
somes 1 and 2, which represented a pattern consistent across                    chromosomes (Fig. 3c). Importantly, the analysis of differential
nonfused chromosomes), interaction patterns were disrupted                      Hi-C matrices and intrachromosomal interaction ratios revealed
during P/D in Rb mice (Fig. 3d and Supplementary Fig. 6), losing                that, on average, probabilities of interactions were higher in Rb
the interchromosome interaction pattern previously described for                mice at genomic distances larger than 10 Mbp. Local adjustments
standard mice14,16. These results suggest the presence of genome-               in chromatin packing density due to differences in loop
wide conformational changes triggered by the presence of Rb                     positioning and size, together with variations in chromosomal
fusions, which was supported by the cytological analysis of                     axis length (as revealed by our cytological observations), likely
chromosome associations (Fig. 4c). For fused chromosomes (e.g.,                 give rise to the differences between standard and Rb mice Hi-
chromosomes 3 and 8, which have a pattern consistent across                     C maps.
fused chromosomes), we not only detected changes in mid-                           Robertsonian mice showed a distinctive interaction pattern in
chromosome interactions, but also high interaction values in                    postmeiotic cells, with a high incidence of interactions at shorter
pericentromeric regions (Fig. 3d–f and Supplementary Fig. 6).                   genomic distances (≤5 Mbp) (Fig. 5a). The effect of highly
   Overall, our data suggest that the presence of Rb fusions                    condensed chromatin in RS of Rb mice would increase short-
reorganize chromosomal nuclear occupancy genome-wide                            range contacts (relative to long-range contacts), as we observe.
(Fig. 4d, e), increasing the rate of heterologous interactions in               Although A/B compartments reappeared in RS, they were present
primary spermatocytes. Moreover, our results suggest that pairing               as a blurry plaid pattern of larger mean size (0.93 Mbp in Rb
disturbances of heterozygous fusions during prophase I have a                   mice) than in fibroblasts. Interestingly, the proportion of genomic
direct impact in chromosome spatial distribution (Fig. 4e). As                  bins with the same compartment status (A or B) was higher
such, when heterozygous fusions are synapsed, more permissive-                  between standard and Rb mice in RS (95% bins conserved, r2 =
ness for centromeric associations between acrocentric chromo-                   0.94) than in fibroblasts (88.24% bins conserved, r2 = 0.80)
somes is observed and the sex body is generally isolated at the                 (Fig. 5b). These results, together with the low inter-/intrachro-
periphery of nuclei. Conversely, when heterozygous Rb fusions                   mosomal interaction values detected in Rb mice (Fig. 3c), suggest
fail to fully synapse, there is a general disruption of interchromo-            that Rb fusions induce physical constraints in RS most likely due
some associations and the sex body shows abnormal signals of                    to the spatial organization of chromosomes around the
heterochromatinization and frequently associates with fused                     chromocenter31,32, thus favoring local interactions.
chromosomes (Fig. 4e).

                                                                                Chromosomal fusions induce TAD reorganization. Chromo-
Higher-order chromatin remodeling. We further investigated                      some fusions also affected chromatin remodeling at the sub-
whether the presence of Rb fusions had an effect on chromatin                   megabase scale as revealed by the number of TADs detected and
remodeling at a finer scale (e.g., compartments). In somatic cells,              the robustness of their boundaries at a 50-kbp resolution
the genome-wide analysis of differential Hi-C matrices showed                   (Fig. 5c–f and Supplementary Fig. 7). Similar to previous obser-
higher interactions at shorter genomic distances in standard mice               vations in standard somatic cells16, TADs were well defined in Rb
than in Rb mice (Fig. 5a). This interaction pattern, which is                   fibroblasts (Fig. 5c–f), with a total of 2391 detected. With an

8                   NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                   ARTICLE




Fig. 4 Interchromosomal associations. a Examples of immunofluorescence on primary spermatocytes at pachytene stage, labeling the synaptonemal
complex with SYCP3 (green), the centromeres (CEN, blue), and H3K9Me3 (red). Immunofluorescence replicates, n = 3. SYCP3 staining allowed for the
detection of the different heterozygous Rb fusion according to synapsis pattern: synapsed, open, and asynapsed (white arrows). H3K9Me3 shows
associations between different chromosomes (yellow dashed lines) and differential distribution in the XY sex body (orange dashed lines): in both arms
(as shown in the synapsed example), in one arm end (as shown in the open example), or fully around the sex body (as shown in the asynapsed example).
The sex body is indicated as XY. Scale bar = 10 µm. b Analysis of H3K9Me3 associations according to synapsis states (synapsed, open, and asynapsed)
of heterozygous metacentrics (left panel), sex chromosome/autosome associations (central panel), and the sex body on its own (right panel). Acr
acrocentric chromosomes, Met metacentric chromosomes, assoc. association. N = number of cells analyzed. Source data are provided as a Source Data
file. c Number of associated chromosomes (metacentrics (n = 60)) or acrocentrics (n = 60)) detected per cell depending on the synapsis state of
heterozygous metacentrics (Kruskal–Wallis, **P ≤ 0.005, one-sided). Boxplots are presented as mean values (center line) ± SD. syn synapsed, op open,
asyn asynapsed, ns nonsignificant. Source data are provided as a Source Data file. d Schematic representation of chromosome organization in P/D
according to the presence of Rb fusions. In standard mice, all chromosomes are acrocentric and are attached to the nuclear lamina. When Rb fusions are
present, chromosome organization is disrupted, affecting chromosome disposition inside the nucleus, either in a homozygous or a heterozygous state.
e Schematic representation of the centromeric associations detected with the H3K9Me3 signal in addition to the XY disposition according to synapsis
pattern of heterozygous Rb fusions (i.e., synapsed, asynapsed, or open). het heterozygous metacentric.



average length of 1.14 Mbps, TADs in Rb fibroblasts had a higher                 the TAD meta-border plot (Fig. 5e) that shows the specific border
variance of insulation scores than in standard fibroblasts (Fig. 5c).            interactions of loop domain TADs3,33.
Nevertheless, the majority (70%) of TADs were stable (Fig. 5d),                    In primary spermatocytes, there was a substantial reduction in
with general preservation of border conformation, as revealed by                the variance of TAD insulation score, and reduction of total


NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                            9
ARTICLE                                                                NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1




TADs detected, when compared to somatic cells (n = 288 in Rb                partially recovered in RS, being higher in Rb mice than standard.
and n = 294 in standard). The TAD insulation score variance was             In RS, Rb mice had fewer but larger (n = 3805, 0.71-Mbp average
higher in Rb mice (Fig. 5c). Moreover, meta-border plots show               length) TADs than standards (n = 4649, 0.59-Mbp average
asymmetric TAD borders (Fig. 5g), more characteristic of stripe             length). Only 40% of TADs were stable in RS (Fig. 5d). The
TAD domains33. The variance of TAD insulation scores was                    presence of weaker TAD borders in both standard and Rb RS

10              NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                        ARTICLE

Fig. 5 Variance in fine-scale compartmentalization. a Differential Hi-C matrices (log2 of fold change using Rb mice as a reference when compared with
standards) for chromosome 3 in cell types analyzed (fibroblasts, pachynema/diplonema—P/D, and round spermatids—RS), at a 500-kbp resolution. The
observed scaling is consistent across all chromosomes. Red indicates higher number of interactions in standard (St) mice when compared to Rb mice,
whereas blue represents higher number of interactions in Rb mice. b Chromosome (Chr) 3 region-specific ICE-corrected heatmaps at 50 kbp (from 30 to
55 Mbp), depicting compartment signal (1st eigenvector) for all cell types. Source data are provided as a Source Data file. c Variance of TAD insulator score
between St and Rb mice in all cell types (Mann–Whitney test, ****P < = 0.0001, ns P > 0.05, two-sided). Boxplots are presented as mean values (center
line) ± SD. ns nonsignificant. d Frequency of TAD reorganizations between standard and Rb for fibroblasts, P/D, and RS. Source data are provided as a
Source Data file. e Example of TAD border alignments along chromosome 1 of fibroblasts (from 75 to 86 Mbp). Examples of merged, split, and stable TADs
are indicated. TAD border scores are also shown, informing of the TAD boundary strength. Source data are provided as a Source Data file. f Schematic
representation of TAD reorganization. Merged TADs are the result of fusing two different TADs. Split TADs are those in which one TAD is divided into two
TADs. TADs are considered stable when there is an overlap above 75%. When TADs are found in a different organizations, they are considered
rearranged. g Metaplots for all TAD boundaries detected in Rb mice: fibroblasts (n = 2378), P/D (n = 288), and RS (n = 3798). Data on standard mice
were extracted from Vara and colleagues16.



(mean TAD insulator score = 5.46), when compared to fibro-                       occupancy occur genome-wide, affecting not only chromosomes
blasts (mean TAD insulator score = 8.8), could explain the                      involved in Rb fusions but also unfused chromosomes. The pat-
difference in TADs observed between these cell types, with TADs                 tern observed in primary spermatocytes is especially relevant as
more prone to reorganize in RS (Fig. 5d). In fact, rearranged,                  fundamental cellular processes take place, such as synapsis and
split, and merged TADs showed significantly lower TAD border                     pairing of homologous chromosomes and the formation and
scores than stable TADs (Mann–Whitney test, P = 2e-14). Meta-                   repair of DSBs. It is also in primary spermatocytes where DNA
border plots showed a blurred pattern, consistent with loop                     loops are organized along the chromosomal axis, with loop size
domain TAD borders and weak TAD insulator scores (Fig. 5g).                     and axis length inversely correlated17,18, which can covary within
                                                                                gametes19. We detected that the presence of Rb fusions can dis-
                                                                                rupt the standard nuclear architecture of peripherally attached
Olfactory receptor genes and interchromosomal interactions.                     acrocentrics to center-localized metacentrics35, ultimately indu-
The increased rate of heterologous interactions detected in Rb                  cing ectopic heterologous interactions. Such alteration in nuclear
primary spermatocytes resulted in the emergence of new inter-                   occupancy genome-wide, together with the pairing disturbances
chromosomal interactions (n = 249, representing 0.5% of the                     observed with fusions in a heterozygous state, detected in a
mouse genome) that involved all chromosomes (e.g., chromo-                      chromosome-specific manner, would trigger stochastic physical
somes 2, 3, and 7, Fig. 6a). To explore the potential role of Rb-               interactions and expose chromosomal domains to novel reg-
specific interchromosomal interactions, we analyzed the gene                     ulatory environments. Remarkably, this disruptive pattern also
ontology of the 2000 genes contained in these regions. Impor-                   affected chromosome axis length and, as a result, higher-order
tantly, we detected enrichment for sensory perception genes                     chromatin remodeling. As chromatin is organized into DNA
(Fig. 6b), specifically, olfactory receptors (OR) (n = 118) and                  loops that emerge out of meiotic chromosomal axes36, variations
vomeronasal receptors (VR) (n = 92). These OR/VR genes with                     in axis length alter the number and size of these loops19, ulti-
heterologous interaction were located in 11 out of the 94 OR                    mately affecting medium- and long-range interactions. Given that
clusters previously described in the mouse genome34 (Fig. 6c).                  the assembly of chromatin loops into chromosomal axes, and the
These included mouse chromosomes 4, 6, 7, 9, 10, 14, 16, 18, and                formation and repair of meiotic DSBs (and subsequent COs) are
19. Interestingly, we detected that individual OR clusters inter-               interconnected18,19,37, the remodeling of chromosomal interac-
acted with a wide range of regions from different chromosomes,                  tions translated into a reorganized recombination landscape in Rb
with interactions being either chromosome-specific (e.g., the                    mice. Thus, Rb fusions pose mechanical constraints on the spatial
cluster on chromosome 19 interacts solely with chromosome 3)                    genome architecture that affect not only the hierarchical three-
or multichromosomal (e.g., interactions between two clusters on                 dimensional organization of the genome, but also the chromo-
chromosome 10 with multiple regions of chromosomes 1, 2, 3, 4,                  somal distribution of CO events. The molecular mechanisms
7, 8, and 9) (Fig. 6d). Furthermore, we detected significantly more              behind this phenomenon remain unknown.
repeats in OR gene clusters when compared with the Rb-specific                      Importantly, the redistribution of COs across chromosomal
interchromosomal interactions and with the rest of the mouse                    arms was consistent with the observation of low recombination
genome (Supplementary Fig. 8a) (Wilcoxon test, P ≤ 0.0001),                     rates at a fine scale (kbp), and high genetic divergence genome-
particularly LINE/L1, and LTR/ERV family retrotransposons                       wide. This provides evidence that disturbances in CO distribution
(Supplementary Fig. 8b).                                                        due to chromosome fusions result in detectable genomic foot-
                                                                                prints. According to the ‘suppressed recombination’ model20, a
Discussion                                                                      reduction in recombination is expected within reorganized
Here, we provide evidence that chromosomal fusions affect three-                regions in the heterokaryotype (e.g., heterozygous Rb fusions).
dimensional genome topology and meiotic recombination, high-                    Although this was consistent with our observations (the pro-
lighting the implications of large-scale genome reorganizations on              portion of chromosome arms with zero CO was slightly higher in
genome function and fertility. Our results show that chromoso-                  heterozygous Rb fusions than nonfused chromosomes), we also
mal fusions pose important mechanistic constraints in the nuclear               detected a strong reduction of recombination in homozygous Rb-
architecture of germ cells, affecting heterologous interactions,                fused chromosomes, mirroring early cytological observations38.
chromosomal synapsis, and meiotic recombination. This was                       Decreased recombination in homozygous Rb fusions is counter-
reflected at different hierarchical levels: (i) chromosomal nuclear              intuitive but could be explained by the presence of double-
occupancy, (ii) inter- and intrachromosomal interactions, (iii)                 centromeric signals that represented misaligned centromeres39.
chromosomal axis length, and (iv) topological domains.                          Since centromeres can reduce recombination rates (so-called
   Our high-throughput analyses in combination with cytological                 centromere interference40), the presence of double centromeres
observations show that disturbances in nuclear chromosomal                      could magnify this effect, interfering with both the formation of

NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                                 11
ARTICLE                                                                   NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1




Fig. 6 Interchromosomal interaction and olfactory receptors. a Interaction profiles of pairs of chromosomes (3 and 19, 7 and 9, and 2 and 14) in primary
spermatocytes (pachynema/diplonema, P/D) of Rb mice. In all three examples, interchromosomal interactions are depicted by asterisks. St standard mice,
Rb Rb mice, chr chromosome. Source data are provided as a Source Data file. b Gene Ontology Enrichment Analysis (GOEA) of genes in newly detected
interchromosomal interactions in pachynema/diplonema of Rb mice. Only significant gene ontology (GO) terms with an adjusted P value below 0.01 are
shown. The X axis represents the number of genes associated to each GO term. Source data are provided as a Source Data file. c Mouse ideogram showing
the localization of olfactory clusters (green) described in literature34, and the interchromosomal interactions detected in this work (red). d Circus plot
representing the interchromosomal interactions in the mouse genome related with sensory perception, which are mostly genes from the olfactory (Olfr)
and vomeronasal (Vmn) receptor family. The number of genes found in each region and the gene family are shown in parenthesis. Source data are provided
as a Source Data file.

COs during prophase and subsequently affecting chromosome                      have an effect on mouse fertility, although mild enough to allow
segregation in metaphase I. Previous studies in BRbS reported                  their persistence within populations.
disturbances in heterochromatization at the centromeres of fused                  Remarkably, we also detected that chromosomal axes of Rb
chromosomes21, and a higher frequency of apoptotic spermato-                   fusions in the heterozygous state were longer when compared to
genetic cells in homozygous mice than in heterozygous mice41. In               nonfused chromosomes. Considering that the presence of Rb
fact, the presence of Rb fusions can induce subfertility, as revealed          fusions can alter the nuclear disposition of chromosomes during
by our observations of the significant reduction in sperm motility              prophase I, this can lead to a delay in synapsis in heterozygous Rb
and vitality between Rb and standard mice. This, together with                 fusions, resulting in the detection of elongated chromosomal axes
previous reports on the variation of the sperm head morphology                 in a chromosome-specific manner. Given that chromosomal axes
and spermatogenic activity in Rb mice42, suggests that Rb fusions              elongate, the substrate for the formation of COs increases in


12                 NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                                    ARTICLE

heterozygous Rb fusions, as both features are correlated18,19.                     recombination analysis. Animals were housed and treated in strict accordance with
Conversely, the cytological data indicated the presence of shorter                 ethical guidelines approved by the Universitat Autònoma de Barcelona (Spain).
                                                                                      The BRbS system mice included in the present study consisted of three
chromosomal axes in homozygous Rb fusions (and therefore,                          populations without Rb fusions (2n = 40; MON, BOI, and Olost) and three
longer DNA loops). This was consistent with the Hi-C data,                         populations with Rb fusions (2n = 39–28; SS, CAS, and VIL). Mice from Rb
suggesting that Rb fusions induce variance in higher-order                         populations are characterized by having between one and six Rb fusions involving
chromatin organization. This adds to initial reports on global                     12 different chromosomes (Rb(3.8), Rb(4.14), Rb(5.15), Rb(6.10), Rb(9.11), and Rb
modulation of chromosomal CO frequencies by both chromatin                         (12.13)), either in heterozygous or homozygous states (Supplementary
                                                                                   Table 1)23,24. The BRbS is characterized by Rb fusions present as chromosomal
loop size and chromosomal axis length17–19, further suggesting                     polymorphisms, thus not fixed within populations.
that recombination landscapes can be altered within cells by Rb                       Moreover, a primary fibroblast cell line derived from a male mouse from the
fusions affecting both the loop/axes length ratio and the spatial                  BRbS system (#954, 2n = 30) previously established in our lab46 was used as a
chromosome occupancy.                                                              somatic control in the Hi-C experiments. Cells were cultured in DMEM medium
                                                                                   supplemented with 10% fetal bovine serum and 1% PenStrep at 37 °C and 5% CO2.
   Likewise, Rb fusions not only pose restrictions for chromoso-
mal interactions in primary spermatocytes but also in postmeiotic
cells. Chromosomes are highly condensed in round spermatids                        Spermatocyte spreads and immunofluorescence. Direct analyses of recombi-
                                                                                   nation are normally based on the detection of either COs or their final products,
with all centromeres associated in the center of the cell forming                  the chiasmata, visible cytogenetically in meiocytes in later stages of the first meiotic
the chromocenter16. As heterologous contacts were reduced in Rb                    division (i.e., pachytene chromosomes)18,37. Here, we analyzed the physical loca-
mice, we hypothesize that chromosomal fusions restrict interac-                    tion of COs along the axes of pachytene chromosomes using the immuno-
tions between nonfused chromosomes in round spermatids. This                       fluorescence staining technique to detect MLH1, a protein that localizes type I
                                                                                   (interfering) COs along with one of the proteins involved in synaptonemal com-
restricted and condensed pattern favors intrachromosomal con-                      plex (SC) formation (the synaptonemal complex protein 3, SYCP3). The chro-
tacts, resulting in the reorganization of TADs in postmeiotic cells.               mosomal distribution of MLH1 foci can be considered a proxy for the positioning
   Our observations have important evolutionary and develop-                       of COs18,47,48. The position of the centromeres along the chromosomal axes was
mental implications. The dynamic three-dimensional genome                          visualized by staining centromeric proteins using the sclerodactyly and tel-
topology of germ cells can be affected by chromosomal fusions in                   angiectasia (CREST) serum (Fig. 1a and Supplementary Table 1). In addition, we
                                                                                   also detected RAD51 (marker of programmed DSBs) at the early stages of prophase
two ways: (i) altering chromosomal nuclear occupancy, and (i)                      I and H3K9Me3 (a marker of centromeric constitutive heterochromatin). Sper-
reshaping landscapes of recombination. The redistribution of                       matocyte spreading and immunofluorescence were performed48. Testes were
chromosomal nuclear occupancy in spermatocytes that result                         mechanically disaggregated until obtaining a cell suspension in 1× PBS. The cell
from Rb fusions brings new genomic regions into close proximity,                   suspension was then distributed into different slides and incubated with 1% lypsol
                                                                                   for 16 min followed by a 20-min incubation with 4% paraformaldehyde. Then
predisposing to the occurrence of additional rearrangements, as                    slides were left to dry and washed twice with PhotoFlo 1% (Kodak) and then
previously suggested43,44. Moreover, rearranged nuclear occu-                      blocked with PBS–Tween-20 (0.05%). Slides were incubated overnight at 4 °C with
pancy would expose chromosomal domains to novel regulatory                         the following primary antibodies: anti-mouse MLH1 (BD Pharmingen Cat#
environments, potentially affecting gene expression and/or reg-                    551092, 1:50), anti-human CREST (courtesy of M. Fritzler, 1:100); anti-rabbit
                                                                                   SYCP3 (Abcam Cat# ab15093, 1:400), anti-rabbit H3K9me3 (Abcam Cat# ab8898,
ulation, as initially proposed by the integrative breakage model of                1:300), and anti-mouse RAD51 (Millipore Cat# PC-130, 1:100). Primary antibodies
genome architecture5,6. As such, new chromosomal interactions                      were detected with anti-rabbit Cy3 (Jackson ImmunoResearch Laboratories Cat#
resulting from chromosomal fusions may rewire or attenuate                         111-165-003, 1:200) combined with anti-mouse FITC (Jackson ImmunoResearch
gene networks, providing new grounds for evolutionary novelty                      Laboratories Cat# 115-095-003, 1:200) and anti-mouse Cy5 (Jackson ImmunoR-
                                                                                   esearch Laboratories Cat# 115-175-146, 1:200). Spermatocyte preparations were
in the long run. This was the case for olfactory receptor family                   visualized and captured using a Zeiss Axioskop epifluorescence microscope
clusters detected in meiotic-specific interchromosomal interac-                     equipped with the appropriate filters and a charged coupled device camera
tions in Rb mice. As this gene family is expressed in the                          (ProgRes® CS10plus, Jenoptik).
male germ line16, with a function in spermatogenesis and
fertilization42, altered regulation of their expression could have an              Sperm analysis. Analyses were conducted in a subset of 15 of male mice from the
adaptive role. In fact, interchromosomal interactions between OR                   BRbS. This included six standard mice from BOI population and nine mice with Rb
genes has been previously described in sensory neurons45, sug-                     fusions from VIL and CAS populations. Briefly, the right epididymis was obtained
gesting that associations between multiple chromosomes can                         for each specimen, and the caudal portion was drained to acquire epididymal
                                                                                   spermatozoa. The process was monitored with an Olympus SZ30 stereoscope
selectively regulate transcription of individual OR genes. Further                 microscope. Spermatozoa were processed following the protocols for the exam-
functional studies will be needed to support this hypothesis in the                ination and processing of human semen samples described by the World Health
germ line.                                                                         Organization49. Around 1000 spermatozoa were analyzed per individual using an
   Here, we demonstrate that chromosomal fusions affect the                        Olympus CH30 microscope.
                                                                                       For each sample, a slide with a drop of fresh spermatozoa suspension was
three-dimensional genome-wide topology in germ cells, ulti-                        analyzed. Depending on the characteristics of the spermatozoa movement, they
mately reshaping recombination landscapes. The modulation of                       were classified into the following categories: (i) fast progressive and linear motility,
recombination implies a close interplay between different factors                  (ii) slow progressive motility by swinging or doing circular movements, (iii)
that are involved in chromatin remodeling, centromere inter-                       nonprogressive motility but the movement of the head/tail, and (iv) no motility. In
                                                                                   order to determine sperm vitality, a drop of eosin Y (0.5%, Sigma-Aldrich) and a
ference, and chromosomal axis length. Documenting how such                         drop of sperm suspension were mixed on a slide. After 2 min, the preparation was
changes in genome organization affect gene expression and reg-                     assessed under the microscope. This method allows the identification of live (not
ulation is an important dimension for further understanding the                    stained) or dead spermatozoa (stained by the inability to expel eosin).
effect of genome reshuffling on evolution and fertility. Our results
can provide the impetus for the exploration of the functional and                  CO data analysis. For the analysis of COs and centromere position, only
structural basis of genomes in a broad context, reinforcing the                    pachytene spermatocytes were considered. Only axis-associated MLH1 foci were
link between the three-dimensional genome architecture, genome                     counted considering the number of MLH1 foci per arm, chromosome, and per cell.
integrity, and fertility.                                                          SC length (expressed in μm) was calculated as the mean length of all autosomal SCs
                                                                                   per cell in individual mice. Thus, we measured the SC length of each chromosome
                                                                                   and the physical distance between COs. For the analysis of RAD51, only leptotene
                                                                                   and zygotene spermatocytes were considered.
Methods                                                                                The Micromeasure 3.3 software50 was used for the analysis of chromosome-
Animals and cell lines. We sampled a total of 63 wild-caught house mice (M. m.     specific recombination maps based on the distances between adjacent MLH1
domesticus) from six populations, covering the extent of the BRbS (Fig. 1a and     foci18,48. For each chromosome, the position of each MLH1 focus was recorded as a
Supplementary Table 1). All animals included in the study were previously          relative position (percentage of the SC total length), identified by the
karyotyped21,23,46 confirming the presence of Rb fusions (Supplementary Table 1).   immunofluorescence signal of the centromeres in each cell. Discrimination of
Three males from the laboratory strain BL6 were also included in the               heterozygous chromosome states was accounted for, in parallel with the presence of

NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                                               13
ARTICLE                                                                             NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1


double-centromere signals in homozygous metacentrics. Detailed analyses of CO            two groups: standard and Rb mice. Mann–Whitney tests were used to analyze
chromosomal distribution, asynapsis, and double centromeres, were conducted in           differences between groups.
animals belonging to one standard population (MON) and one Rb
population (VIL).
    CO frequency plots were constructed for each chromosome type (acrocentric,           FACS of mouse male germ cells. Testis cell disaggregation and FACS were
Rb in heterozygous state, and Rb in homozygous state). As CO data were not               conducted16. Briefly, germ cells at a concentration of 1 million per 500 µl were
normally distributed, analysis of variation in the number and position of MLH1           incubated in formaldehyde (1%) for 10 min prior to FACS. Glycine (0.125 M) was
foci along chromosomes among different groups was assessed using nonparametric           added and incubated with agitation at room temperature for 5 min and then at 4 °C
tests, including Bartlett’s test (P = 0.1372), thus demonstrating equal variance.        for 15 min. Cells were then centrifuged for 10 min at 290 × g at 4 °C and resus-
Accordingly, the Mann–Whitney test was used in comparisons between groups,               pended in 3 ml of 1× PBS with Hoechst staining.
Kruskal–Wallis to test intragroup differences (either by population or type of               Germ cells were sorted using a BD InfluxTM (BD Biosciences) coupled with the
fusion), and Wilcoxon’s rank-sum test followed by Dunn’s tests adjusted by               BD FACSTM software (version 1.0) and an ultraviolet laser (355 nm). Subsequently,
Bonferroni for multiple comparisons. Data were expressed as the mean ± standard          two main germ cell populations (P/D and RS) were isolated by plotting Hoechst
deviation (SD). In addition, Pearson’s χ2 test was used to compare arm proportions       Blue (UV355–460/50) vs. Hoechst red (UV355–670/30) emissions to discriminate
of the different number of MLH1 foci per arm type and the differences in the             cells by both their DNA content and their complexity. Cell populations were
relative distribution of MLH1 foci along chromosomal arms and between                    collected after sorting in 1× PBS and centrifuged for 5 min at 1800 × g. The
chromosomal arms according to the presence of Rb fusions (acrocentric nonfused,          supernatant was discarded, and cell pellets were flash-frozen at −80 °C until use.
heterozygous metacentric, and homozygous metacentric). Moreover, Spearman                Sorting experiments lasted between 3 and 6 h to collect between 0.2 × 106 and 3.2 ×
correlations were calculated to test relationships among MLH1 and RAD51 foci,            106 cells, depending on the germ cell population.
MLH1 foci and diploid number, and MLH1 foci and the number and state of                      Cell enrichment of each flow-sorted population was evaluated by
chromosomal fusion. Spearman correlation analysis was also performed for MLH1            immunofluorescence using specific meiotic proteins and DAPI morphology. For
foci and chromosome arm length. In all statistical analyses, a P ≤ 0.05 was              primary spermatocytes, prophase-I stages (leptonema, zygonema, pachynema, and
considered statistically significant.                                                     diplonema) were identified based on SYCP3 (1:400) and ɣH2AX patterns (1:300).
                                                                                         Cell enrichment of round spermatids was determined based on nucleus
                                                                                         morphology and DAPI pattern16. Cells were fixed on slides and then mounted with
SNP genotyping. Analyses of genomic divergence and recombination rate at the             DAPI diluted in Vectashield (Vector Laboratories). Slides were analyzed using
Mbp scale were conducted using genotyping data from a subset of 34 mice from             fluorescence microscopy (Axiophot, Zeiss) coupled with a ProgRes® CS10plus,
two standard populations (2n = 40; BOI and Olost) and two Rb populations (2n =           Jenoptik camera. Representative images were captured with ACO XY (A. Coloma,
28–39; SS and CAS populations) retrieved from24 (Supplementary Table 1). Data            Open Microscopy). Between 50 and 100, cells were counted for each flow-sorted
consisted of the Mouse Universal Genotyping Array (MegaMUGA), which con-                 population. Only sorted populations with an enrichment above 80% were
sisted of 77,808 evenly distributed SNP markers built on the Illumina Infinium            considered for subsequent experiments.
platform51. SNPs were filtered to remove markers with missing values >5%
threshold using PLINK version 1.952. This resulted in a final data set of 63,344
informative SNPs distributed across all chromosomes, with the exception of
chromosomes 8 and Y. The final data set was considered for subsequent analyses of         In nuclei Hi-C. The generation of Hi-C libraries was conducted following Vara and
genome-wide screening of divergence and estimates of recombination rate.                 collaborators16. Rb mice included in the Hi-C analysis were selected based on their
    The ADMIXTURE software53 was used to estimate individual ancestry and                karyotype characteristics (high number of Rb fusions) and availability of testis
admixture proportions assuming K populations based on a maximum likelihood               material. All mice included in the Hi-C experiments belong to the same population
method. Analyses were run only for SNPs with a greater than 95% genotype call.           (VIL) and were included in the recombination analysis, showing similar patterns of
The numbers of clusters (K) were evaluated applying Evanno’s ΔK54 with three             CO distribution (Supplementary Fig. 1f and Supplementary Table 1). Two repli-
different K values (K = 3, 4, 5) showing the lowest likelihood values. ADMIXTURE         cates for cell type were obtained from a total of 3.4 × 106 primary spermatocytes at
analyses were plotted using the R package Pophelper v2.3.055. In addition, multiple      the P/D stage, 12.8 × 106 round spermatids previously isolated by FACs. In addi-
dimensional scaling analysis was performed using PLINK, by first generating a             tion, a total of 10 × 106 Rb fibroblasts (two biological replicates) was also included.
genome file (–genome flag) from the vcf file containing the SNPs, and then the mds          Libraries were submitted for Illumina sequencing (paired-end 75 bp each side on
file that was plotted in the R environment.                                               HiSeq 2500, v4).


Genome-wide screening of genetic divergence and diversity. We estimated the              Hi-C data processing, binning, and normalization. The quality check and
number of alleles (Na), allelic richness (Ar), observed heterozygosity (Ho), expected    trimming step of raw data was carried out using BBDuk (version 10/2015)64.
heterozygosity (He), inbreeding coefficient (FIS), and nucleotide diversity (pi). Na      Setting a minimum read length of 35 bp and a minimum Phred quality score of 20,
and Ar were estimated using the hierfstat package v0.04-2256 implemented in R.           adapters and low-quality reads were removed while preserving their longest high-
Allelic richness was refracted for a minimum of 22 alleles (or 11 diploid samples),      quality regions. After the quality check, the reads were processed with TADbit
which was the lowest observed sample size between the three groups. Ho and He            (version 0.2.0.23)65, which makes use of the GEM (version 1.7.1) mapper66 to
were calculated using PLINK v1.90b6.1252, and FIS and pi with VCFtools 162               iteratively map them against the mouse genome (version mm10). Reads were
v.0.1.16. One-thousand bootstraps were performed for pairwise FST57 estimations          mapped from 15 bp toward using a step size of 5 bp. The filters used to remove
with the StAMPP v1.6.1 package in R58.                                                   possible artifacts were the following: self-circle, dangling-end, error, extra dangling-
    Pairwise FST comparisons were conducted between populations genome-wide              end, too short, too large, duplicated, and random breaks. The maximum molecule-
and considering chromosomes involved in fusions (3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14,   length parameter was set at two times the 99.9 percentile of the insert-size dis-
and 15) and not involved in fusions (1, 2, 7, 16, 17, 18, 19, and X). Estimated FST      tribution, returned by the insert_size from TADbit. The maximum distance of a
values were adjusted with the Bonferroni correction to minimize type I errors.           read to a cleavage site was set to the 99.9 percentile of the insert-size distribution.
Tukey–Kramer tests (JMP package version 5.1.2, SAS Institute Inc.59) were used to            An in-house script was used for binning and data normalization. This script
analyze differences between groups.                                                      imported the HiC_data module of TADbit, read the map files generated after the
                                                                                         artifact 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
Estimates of recombination rates. The program LDhelmet60 was applied for the             3.3)67 was used to normalize with the ICE (Iterative Correction and Eigenvector
estimation of recombination rates at a fine scale (kbp). As LDhelmet has a 25-            decomposition). The normalized matrices of standard and Rb were then compared
diploid-sample limit (50 haplotypes), we sampled a random subset of 25 indivi-           by the log2 ratio method using hicCompareMatrix from HiCExplorer to obtain the
duals from the 34 individuals included in the SNP analysis, using the vcftools —         differential matrices.
max-indv option and by chromosome using —recode and —chr61. LDhelmet                         Pairwise comparisons between biological replicates derived from the Hi-C
estimates the recombination rates from phased chromosomes or haplotypes from a           experiments were performed using HiCRep (version 1.4), under a smoothing
population; thus, we phased our data using the software SHAPEIT62 by using the           parameter of 5 and a considered distance over 10 Mbp16,68. The Y chromosome
option –rho 0.001. Once phased, each chromosome vcf file was split into two               was excluded from the analysis due to the lower number of interactions detected in
groups, according to standard or Rb samples. Each file was then transformed to            our analysis (<1% of the overall detected interactions) and the highly repetitive
LDhelmet input snps and pos files with the —ldhelmet flag. The likelihood tables           DNA that characterizes this chromosome. The correlation between 2 replicates was
were generated using LDpop63, and then transformed to LDhelmet format fol-               defined as the mean of the 20 correlation scores.
lowing the software’s manual indications. We performed the analysis per-chro-                Normalized matrices of standard and Rb mice were transformed into observed/
mosome, based on the SNP data. Estimations of the population-scaled                      expected matrices using the tool hicTransform from HiCExplorer (version 3.3)67.
recombination rate ρ = 4Ner were obtained using the parameters recommended by            This tool computes the expected matrix as the sum of contacts per genomic
software’s developers, where Ne is the effective population size and r the genetic       distance divided by the maximal possible contacts. Then, the observed matrix was
map distance across the region analyzed. Using this approach, we established             divided by this calculated expected matrix and the observed/expected matrix is
recombination rates in windows of 50 SNPs across the mouse genome considering            obtained.

14                    NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                                                       ARTICLE

   The observed/expected matrices were compared by the log2 ratio method using        median, while the lower and upper bounds correspond to the first and third
hicCompareMatrix from HiCExplorer to obtain the differential matrices. Finally,       quartiles; the upper whisker extends the largest value within 1.5 times the inter-
these matrices were plotted by chromosome using a 500-kb resolution.                  quartile range above the 75th percentile and the lower whisker extends to the
                                                                                      smallest value within 1.5 times the interquartile range below the 25th percentile.
Interchromosome/intrachromosome interaction ratio. ICE-normalized data
stored in matrices were exported with HiCExplorer to the GInteractions format,        Reporting summary. Further information on research design is available in the Nature
which consists of seven columns: chromosome, start and end from bin 1, chro-          Research Reporting Summary linked to this article.
mosome, start and end from bin 2, and the amount of interaction. The GInter-
actions tables were imported in R for further quantification of interchromosome
and intrachromosome interactions and plotting.                                        Data availability
                                                                                      The data that support this study are available from the corresponding author upon
                                                                                      reasonable request. Raw and processed Hi-C data from Rb mice generated in the course
Interchromosomal interaction analysis. Using R (version 3.6.1), the GInteraction      of this study are available in the NCBI GEO repository under accession number
tables were subset by chromosome, so the analysis of each chromosome inter-           GSE145978. The Hi-C data set from standard mice was retrieved from Vara and
chromosomal interactions could be done individually. Then, the mean of inter-         collaborators16, which is available in the NCBI GEO repository, under the accession
actions of a given chromosome with others was calculated16. Finally, data were        number GSE132054. Source data are provided with this paper.
plotted as a heatmap where red stands for more interaction in Rb mice and blue in
standard mice. The GIinteraction tables subset by chromosome were then plotted
as interaction profiles considering the interaction for each genomic position.         Received: 4 August 2020; Accepted: 16 April 2021;
    The analysis of the intrachromosomal interaction ratio in RS (haploid cells)
allowed for a predicted classification of fusions into homozygous or heterozygous
by quantifying the interactions of fused chromosomes when compared to
chromosomes not involved in Rb fusions. In this manner, we established that           References
fusions 3.8, 6.10, and 5.15 were present in a heterozygous state (ratio between 1.2   1.    Dixon, J. R. et al. Topological domains in mammalian genomes identified by
and 1.6) while fusions 9.11, 12.13, and 4.14 were present in a homozygous state             analysis of chromatin interactions. Nature 485, 376–380 (2012).
(ratio around 2).                                                                     2.    Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions
                                                                                            reveals folding principles of the human genome. Science 326, 289–293 (2009).
Chromosomal-specific interchromosomal interaction analysis. Inter-                     3.    Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution
chromosomal regions were statistically defined as bins with a standard deviation             reveals principles of chromatin looping. Cell 159, 1665–1680 (2014).
(Z-score) higher than 2.58. Bin regions were also intersected with BEDTools           4.    Kraft, K. et al. Serial genomic inversions induce tissue-specific architectural
intersect (version 2.25) against promoter regions (-2 kbp to TSS) from the mouse            stripes, gene misexpression and congenital malformations. Nat. Cell Biol. 21,
GENCODE annotation vM14 to obtain the overlapping genes. These genes were                   305–310 (2019).
then included in a Gene Ontology Enrichment Analysis (GOEA). We used CIR-             5.    Deakin, J. E. et al. Chromosomics: Bridging the gap between genomes and
COS (version 0.69-8)69 to plot the genomic positions of the interchromosomal                chromosomes. Genes 10, 627 (2019).
interactions across the mouse genome.                                                 6.    Farré, M., Robinson, T. J. & Ruiz-Herrera, A. An integrative breakage model of
    For the analysis of repetitive sequences within interchromosomal regions,               genome architecture, reshuffling and evolution. BioEssays 37, 479–488 (2015).
repeatMasker70 annotation on the Mus musculus genome version mm10 was                 7.    Franke, M. & Gómez-Skarmeta, J. L. An evolutionary perspective of regulatory
downloaded and parsed to become a BED file. BEDtools (version 2.26) was used to              landscape dynamics in development and disease. Curr. Opin. Cell Biol. 55,
intersect equal-sized bins considered as interchromosomal interaction regions               24–29 (2018).
against the RepeatMasker annotation with option -wo, thus counting the number         8.    Bompadre, O. & Andrey, G. Chromatin topology in development and disease.
of overlapping bases. All equal-sized bins of the genome were also intersected, and         Curr. Opin. Genet. Dev. 55, 32–38 (2019).
the number of overlapping bases was also counted and considered as the                9.    Robertson, W. R. B. Chromosome studies. I. Taxonomic relationships shown
background-repeat profile of the mm10 genome. A Mann–Whitney test was                        in the chromosomes of Tettigidae and Acrididae: V-shaped chromosomes and
applied to test for significant differences.                                                 their significance in Acrididae, Locustidae, and Gryllidae: chromosomes and
                                                                                            variation. J. Morphol. 27, 179–331 (1916).
A/B compartments and TAD calling. Analyses were conducted at the genome-              10.   Robinson, T. J. & King, M. Species evolution: the role of chromosome change.
wide level. For A/B compartment calling, columns with a low number of counts                Syst. Biol. 44, 578 (1995).
were filtered out using TADbit, setting the parameter min_count to 10. Since           11.   Scriven, P. N. Robertsonian translocations–reproductive risks and indications
TADbit fits the column count distribution into a polynomial distribution, columns            for preimplantation genetic diagnosis. Hum. Reprod. 16, 2267–2273 (2001).
with a number of counts smaller than the first antimode of the distribution, which     12.   Lucchesi, J. C. & Suzuki, D. T. The interchromosomal control of
cannot be smaller than the min_count parameter, are filtered out. Then, the                  recombination. Annu. Rev. Genet. 2, 53–86 (1968).
genome-wide matrices were normalized by the expected interactions at a given          13.   Alavattam, K. G. et al. Attenuated chromatin compartmentalization in meiosis
distance and by visibility by means of one iteration of the ICE method. The                 and its maturation in sperm development. Nat. Struct. Mol. Biol. 26, 175–184
correlation analysis was also performed with TADbit. In-house scripts computed              (2019).
A/B compartments from the first eigenvector, using 0 as the threshold to differ-       14.   Patel, L. et al. Dynamic reorganization of the genome shapes the
entiate both compartments and the gene density to label them.                               recombination landscape in meiotic prophase. Nat. Struct. Mol. Biol. 26,
    TADs were identified using an in-house script that imported the Chromosome               164–174 (2019).
module of TADbit and added the raw and the ICE-normalized matrices of each            15.   Wang, Y. et al. Reprogramming of meiotic chromatin architecture during
chromosome separately. Filtered bins, due to low counts, were included to mask              spermatogenesis. Mol. Cell 73, 547–561.e6 (2019).
them when calling TADs. TAD insulation scores were obtained by first                   16.   Vara, C. et al. Three-dimensional genomic structure and cohesin occupancy
normalizing the different matrices for read depth in order for the scores to be             correlate with transcriptional activity during spermatogenesis. Cell Rep. 28,
comparable. Each matrix was then scaled to have 100 M reads. Afterward, TAD                 352–367.e9 (2019).
insulation scores were obtained from the output given by the hicFindTADs              17.   Wang, S., Zickler, D., Kleckner, N. & Zhang, L. Meiotic crossover patterns:
program from HiCExplorer.                                                                   obligatory crossover, interference and homeostasis in a single process. Cell
                                                                                            Cycle 14, 305–314 (2015).
Compartment switching. BED files with a resolution of 50 kbp were available            18.   Ruiz-Herrera, A. et al. Recombination correlates with synaptonemal complex
from the compartments definition step. Each genomic bin of 50 kbp had its cor-               length and chromatin loop size in bovids—insights into mammalian meiotic
responding compartment attributed. Pairwise comparisons between cell types—                 chromosomal organization. Chromosoma 126, 615–631 (2017).
genome-wide and per-chromosome—were performed; the ratio of compartment               19.   Wang, S. et al. Per-nucleus crossover covariation and implications for
switching was calculated as the number of genomic bins with a compartment                   evolution. Cell 177, 326–338.e16 (2019).
change (A > B or B > A) divided by the total number of bins. From these files, a       20.   Farré, M., Micheletti, D. & Ruiz-Herrera, A. Recombination rates and genomic
matrix file was created with 50-kbp-binned genomic coordinates as rows and cell              shuffling in human and chimpanzee—a new twist in the chromosomal
types as columns, filled by the corresponding compartment labeling in each bin               speciation theory. Mol. Biol. Evol. 30, 853–864 (2013).
and cell type. Cell-specific A compartments were defined as those bins being            21.   Capilla, L. et al. Genetic recombination variation in wild Robertsonian mice:
compartment A in a cell type and compartment B in the remaining cell types.                 on the role of chromosomal fusions and Prdm9 allelic background. Proc. Biol.
                                                                                            Sci. 281, 1–18 (2014).
Quantification and statistical analyses. The statistical analyses were performed       22.   Ullastres, A., Farré, M., Capilla, L. & Ruiz-Herrera, A. Unraveling the effect of
using R. Statistical parameters and tests are reported in the Figures and Figure            genomic structural changes in the rhesus macaque - implications for the
Legends. Boxplots are represented in the manner that the middle line is the                 adaptive role of inversions. BMC Genomics 15, 530 (2014).

NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications                                                  15
ARTICLE                                                                         NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1


23. Medarde, N., López-Fuster, M. J., Muñoz-Muñoz, F. & Ventura, J. Spatio-          54. Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of
    temporal variation in the structure of a chromosomal polymorphism zone in            individuals using the software STRUCTURE: a simulation study. Mol. Ecol.
    the house mouse. Heredity 109, 78–89 (2012).                                         14, 2611–2620 (2005).
24. Vara, C. et al. PRDM9 diversity at fine geographical scale reveals contrasting    55. Francis, R. M. pophelper: an R package and web app to analyse and visualize
    evolutionary patterns and functional constraints in natural populations of           population structure. Mol. Ecol. Resour. 17, 27–32 (2017).
    house mice. Mol. Biol. Evol. 36, 1686–1700 (2019).                               56. Goudet, J. hierfstat, a package for r to compute and test hierarchical F-
25. Merico, V., Pigozzi, M., Esposito, A., Merani, M. & Garagna, S. Meiotic              statistics. Mol. Ecol. Notes 5, 184–186 (2005).
    recombination and spermatogenic impairment in Mus musculus domesticus            57. Weir, B. S. & Cockerham, C. C. Estimating F-statistics for the analysis of
    carrying multiple simple Robertsonian translocations. Cytogenet. Genome Res.         population structure. Evolution 38, 1358–1370 (1984).
    103, 321–329 (2003).                                                             58. Pembleton, L. W., Cogan, N. O. I. & Forster, J. W. StAMPP: an R package for
26. Petes, T. D. Meiotic recombination hot spots and cold spots. Nat. Rev. Genet.        calculation of genetic differentiation and structure of mixed-ploidy level
    2, 360–369 (2001).                                                                   populations. Mol. Ecol. Resour. 13, 946–952 (2013).
27. Baudat, F., Imai, Y. & de Massy, B. Meiotic recombination in mammals:            59. JMP®, Version 14.1.0. SAS Institute Inc., Cary, NC, 1989–2019.
    localization and regulation. Nat. Rev. Genet. 14, 794–806 (2013).                60. Chan, A. H., Jenkins, P. A. & Song, Y. S. Genome-wide fine-scale recombination
28. Lam, I. & Keeney, S. Mechanism and regulation of meiotic recombination               rate variation in Drosophila melanogaster. PLoS Genet. 8, e1003090 (2012).
    initiation. Cold Spring Harb. Perspect. Biol. 7, a016634 (2015).                 61. Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27,
29. Ribagorda, M. et al. Meiotic behavior of a complex hexavalent in heterozygous        2156–2158 (2011).
    mice for Robertsonian translocations: insights for synapsis dynamics.            62. Delaneau, O., Marchini, J. & Zagury, J.-F. A linear complexity phasing method
    Chromosoma 128, 149–163 (2019).                                                      for thousands of genomes. Nat. Methods 9, 179–181 (2012).
30. Rieseberg, L. H. Chromosomal rearrangements and speciation. Trends Ecol.         63. Alexander, T. A. & Machiela, M. J. LDpop: an interactive online tool to calculate
    Evol. 16, 351–358 (2001).                                                            and visualize geographic LD patterns. BMC Bioinforma. 21, 1–4 (2020).
31. Brinkley, B. R. et al. Arrangements of kinetochores in mouse cells during        64. Bushnell, B. BBmap: a fast, accurate, splice-aware aligner. https://www.osti.
    meiosis and spermiogenesis. Chromosoma 94, 309–317 (1986).                           gov/biblio/1241166-bbmap-fast-accurate-splice-aware-aligner (2014).
32. Hoyer-Fender, S., Singh, P. B. & Motzkus, D. The murine heterochromatin          65. Serra, F. et al. Automatic analysis and 3D-modelling of Hi-C data using
    protein M31 is associated with the chromocenter in round spermatids and Is a         TADbit reveals structural features of the fly chromatin colors. PLoS Comput.
    component of mature spermatozoa. Exp. Cell Res. 254, 72–79 (2000).                   Biol. 13, e1005665 (2017).
33. Barrington, C. et al. Enhancer accessibility and CTCF occupancy underlie         66. Marco-Sola, S., Sammeth, M., Guigó, R. & Ribeca, P. The GEM mapper: fast,
    asymmetric TAD architecture and cell type specific genome topology. Nat.              accurate and versatile alignment by filtration. Nat. Methods 9, 1185–1188 (2012).
    Commun. 10, 1–14 (2019).                                                         67. Wolff, J. et al. Galaxy HiCExplorer: a web server for reproducible Hi-C data
34. Degl’Innocenti, A., Meloni, G., Mazzolai, B. & Ciofani, G. A purely                  analysis, quality control and visualization. Nucleic Acids Res. 46, W11–W16 (2018).
    bioinformatic pipeline for the prediction of mammalian odorant receptor gene     68. Yang, T. et al. HiCRep: assessing the reproducibility of Hi-C data using a
    enhancers. BMC Bioinforma. 20, 474 (2019).                                           stratum-adjusted correlation coefficient. Genome Res. 27, 1939–1949 (2017).
35. Berríos, S. et al. Robertsonian chromosomes and the nuclear architecture of      69. Connors, J. et al. Circos: an information aesthetic for comparative genomics.
    mouse meiotic prophase spermatocytes. Biol. Res. 47, 16 (2014).                      Genome Res. 19, 1639–1645 (2009).
36. Zickler, D. & Kleckner, N. Meiotic Chromosomes: integrating structure and        70. Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0. http://
    function. Annu. Rev. Genet. 33, 603–754 (1999).                                      repeatmasker.org (2015).
37. Capilla, L., Garcia Caldés, M. & Ruiz-Herrera, A. Mammalian meiotic
    recombination: a toolbox for genome evolution. Cytogenet. Genome Res. 150,
    1–16 (2016).                                                                     Acknowledgements
38. Bidau, C. J., Giménez, M. D., Palmer, C. L. & Searle, J. B. The effects of       This work was supported by the Ministry of Economy and Competitiveness (CGL2014-
    Robertsonian fusions on chiasma frequency and distribution in the house          54317-P and CGL2017-83802-P to A.R.-H., BFU2013-47736-P and BFU2017-85926-P to
    mouse (Mus musculus domesticus) from a hybrid zone in northern Scotland.         M.A.M.-R. and CGL2010-15243 to J.V.) and the Agència de Gestió d’Ajuts Universitaris i
    Heredity 87, 305–313 (2001).                                                     de Recerca, AGAUR (DI2015 and SGR1215 to A.R-H. as well as SGR468 to M.A.M.-R.).
39. Cappelletti, E. et al. CENP-A binding domains and recombination patterns in      C.V. and L.A.-G. were supported by FPI predoctoral fellowships from the Ministry of
    horse spermatocytes. Sci. Rep. 9, 15800 (2019).                                  Economy and Competitiveness (BES-2015-072924 and PRE-2018-083257). L.M.-G. was
40. Lynn, A., Ashley, T. & Hassold, T. Variation in human meiotic recombination.     supported by an FPU predoctoral fellowships from the Ministry of Science, Innovation
    Annu. Rev. Genomics Hum. Genet. 5, 317–349 (2004).                               and University (FPU18/03867). M.A.M.-R. acknowledges support by the European
41. Medarde, N. et al. Impact of the number of Robertsonian chromosomes on           Research Council under the 7th Framework Program FP7/2007-2013 (ERC grant
    germ cell death in wild male house mice. Chromosom. Res 23, 159–169 (2015).      agreement 609989) and the European Union’s Horizon 2020 research and innovation
42. Sans-Fuentes, M. A., García-Valero, J., Ventura, J. & López-Fuster, M. J.        program (grant agreement 676556). A.R.-H. also acknowledges support from MeioNet
    Spermatogenesis in house mouse in a Robertsonian polymorphism zone.              (BFU2015-71786-REDT).
    Reproduction 140, 569–581 (2010).
43. Berríos, S. et al. Model of chromosome associations in Mus domesticus            Author contributions
    spermatocytes. Biol. Res. 43, 275–286 (2010).                                    A.R.-H. conceived and devised the study. C.V., A.P.-G., F.L.D., Y.C., M.A.M.-R., and A.R.-H.
44. Branco, M. R. & Pombo, A. Intermingling of chromosome territories in             designed experiments. C.V., Y.C., F.L.D., F.G., B.F.-S., L.M.-G., L.C., Z.S., and A.R.-H.
    interphase suggests role in translocations and transcription-dependent           performed experiments. C.V., A.P.-G., L.A.-G., L.M.-G., R.A.S.-G., and A.R.-H. analyzed the
    associations. PLoS Biol. 4, e138 (2006).                                         data. R.A.-C., W.S., J.B.S., J.V., M.A.M.-R., and A.R.-H. contributed specimens, karyotypic
45. Monahan, K., Horta, A. & Lomvardas, S. LHX2- and LDB1-mediated trans             data, reagents, and resources. C.V. and A.R.-H. wrote the first draft of the paper with inputs
    interactions regulate olfactory receptor choice. Nature 565, 448–453 (2019).     from all authors. All authors read and approved the final version of the paper.
46. Sánchez-Guillén, R. A. et al. On the origin of Robertsonian fusions in nature:
    evidence of telomere shortening in wild house mice. J. Evol. Biol. 28, 241–249
    (2015).                                                                          Competing interests
47. Froenicke, L., Anderson, L. K., Wienberg, J. & Ashley, T. Male mouse             The authors declare no competing interests.
    recombination maps for each autosome identified by chromosome painting.
    Am. J. Hum. Genet. 71, 1353–1368 (2002).
48. Segura, J. et al. Evolution of recombination in eutherian mammals: insights
                                                                                     Additional information
                                                                                     Supplementary information The online version contains supplementary material
    into mechanisms that affect recombination rates and crossover interference.
                                                                                     available at https://doi.org/10.1038/s41467-021-23270-1.
    Proc. R. Soc. B Biol. Sci. 280, 20131945 (2013).
49. World Health Organization. WHO Laboratory Manual for the Examination
                                                                                     Correspondence and requests for materials should be addressed to A.R-H.
    and Processing of Human Semen (WHO Press, 2010).
50. Reeves, A. MicroMeasure: a new computer program for the collection and           Peer review information Nature Communications thanks Akihiko Okuda and the other,
    analysis of cytogenetic data. Genome 44, 439–443 (2001).                         anonymous, reviewers for their contribution to the peer review of this work.
51. Morgan, A. P. et al. The mouse universal genotyping array: from substrains to
    subspecies. Genes|Genomes|Genet. 6, 263–279 (2016).                              Reprints and permission information is available at http://www.nature.com/reprints
52. Purcell, S. et al. PLINK: a tool set for whole-genome association and
    population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).
                                                                                     Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
53. Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of
                                                                                     published maps and institutional affiliations.
    ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).


16                   NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23270-1                                                         ARTICLE

                  Open Access This article is licensed under a Creative Commons
                  Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made. The images or other third party
material in this article are included in the article’s Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not included in the
article’s Creative Commons license and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this license, visit http://creativecommons.org/
licenses/by/4.0/.


© The Author(s) 2021




NATURE COMMUNICATIONS | (2021)12:2981 | https://doi.org/10.1038/s41467-021-23270-1 | www.nature.com/naturecommunications        17