Main

Genome sequences of several pathogens have revealed a partitioning of chromosomes, with housekeeping genes often being located in the central core and antigen genes being located in subtelomeric regions5,6. These assemblies suggest that the linear organization of the genome may be important for restricting high levels of recombination to regions that code for antigens and for ensuring that all but one antigen is repressed.

Recently, genome-wide Hi-C analyses have begun to uncover the 3D organization of chromosomes at high resolution4, which has highlighted the critical role of spatial organization and compartmentalization of DNA in the regulation of gene expression and recombination2,3. In addition, microscopy-based analyses of the unicellular eukaryotic parasites Plasmodium falciparum and T. brucei have indicated that nuclear organization may be important for the mutually exclusive expression of antigens7,8,9. However, to our knowledge, the proteins that are involved in shaping genome architecture and controlling antigen expression have not yet been identified in any organism.

This study aimed to identify the process that restricts antigen expression. Specifically, we sought to identify proteins that are important for maintaining genome architecture and to determine whether global and/or local changes in chromatin conformation affect antigen expression.

In T. brucei—which is the causative agent of human sleeping sickness—the key antigens are the variant surface glycoproteins (VSGs). Most VSG genes—of which there are about 2,500—are found in long subtelomeric arrays of megabase chromosomes6. In addition, about 65 VSG genes are located on mini-chromosomes (50–150 kb in length)10 and a smaller subset of VSG genes is located in distinct telomere-proximal polycistronic transcription units, called expression sites11. Expression sites are grouped into metacyclic-form and bloodstream-form expression sites (MESs and BESs, respectively) on the basis of the life-cycle stage during which they can be activated. VSG genes are transcribed only when they are located within an expression site and only one of about 15 BESs is transcribed at a time, which ensures that the expression of VSG genes is mutually exclusive11. Therefore, a genome sequence that contains both subtelomeric VSG gene arrays and telomeric expression sites, which is lacking in the available T. brucei genome (isolate TREU 927)6, is required to elucidate the molecular link between genome architecture and antigenic variation.

Using PacBio single-molecule real-time (SMRT) sequencing technology, we generated an approximately 100-fold genome-sequence coverage of the T. brucei 427 Lister isolate (the most commonly used laboratory isolate) and assembled the reads into megabase chromosomes, of which there are 11 (96 contigs, Fig. 1, Extended Data Table 1). To order and orient contigs without relying on scaffolds of related parasite isolates (which may have undergone genome rearrangements), we took advantage of two ubiquitous features of chromosome organization: a distance-dependent decay of DNA–DNA interaction frequency and substantially higher interaction frequencies between DNA loci located on the same chromosome, compared to those on different chromosomes4. The high degree of subtelomeric heterozygosity enabled us to assemble the complete T. brucei genome with phased diploid subtelomeric regions (Extended Data Figs. 1, 2, Supplementary Data). In addition, RNA sequencing (RNA-seq) revealed a notable partitioning of the genome into a transcribed homozygous core and non-transcribed heterozygous subtelomeric regions, which encode the vast repertoire of antigens (Fig. 1).

Fig. 1: Long-read and Hi-C-based de novo assembly of the T. brucei Lister 427 genome.
figure 1

Only one of the two homologous chromosomes (chr.) is depicted for the homozygous chromosomal core regions (22.71 Mb). Both chromosomes are shown for the heterozygous subtelomeric regions (19.54 Mb). Relative transcript levels (window size, 5,001 bp; step size, 101 bp) are shown as a black line above each chromosome. BESs and MESs were assigned to the respective subtelomeric region if an unambiguous assignment based on DNA interaction data was possible (see Supplementary Information). Centromeres were assigned based on KKT2 ChIP–seq data30.

Analysis of the frequency of intra-chromosomal DNA–DNA interaction suggested a strong compartmentalization of the T. brucei genome: centromeres and junctions between the core and subtelomeres function as the most prominent boundaries of DNA compartments. In addition, the frequency of DNA–DNA contact was substantially higher across subtelomeric regions compared to core regions, which indicates that subtelomeres are more compact than the core region (Fig. 2a, Extended Data Fig. 3). Therefore, the partitioning of the genome into transcribed housekeeping genes and non-transcribed antigen genes that is observed in the genome assembly and transcriptome data is mirrored by the 3D organization of the genome. In T. brucei, RNA polymerase II transcription can occur in the absence of canonical promoter motifs12,13. Thus, the high degree of compaction across subtelomeric regions probably prevents the spurious initiation of transcription and ensures mutually exclusive expression of a single VSG gene from a BES. In addition, BES–BES interactions were much more frequent than interactions among randomly chosen genomic loci, suggesting a clustering of BESs (Fig. 2b). Taken together, the Hi-C data suggest a distinct compartmentalization of the T. brucei nucleus.

Fig. 2: Hi-C and ChIP–seq reveal partitioning of the T. brucei genome into distinct domains.
figure 2

a, Hi-C heat maps of chromosomes 3 and 6 at 20-kb resolution. Horizontal blue, black and red lines mark heterozygous subtelomeric, homozygous core regions and BESs, respectively. Centromeres are marked by asterisks. b, Scatter plot showing inter-chromosomal interaction frequencies among centromeres (cen) (n = 206 bins; P = 0.0029), VSG genes in silent expression sites (VSGs) (n = 54 bins; P = 1.63 × 10−6) and rRNA genes (n = 40 bins; P = 0.0177) compared to a matching background sample, which was randomly selected from the interaction matrix (50-kb bin size). The background sample (grey) matches the genomic feature (red) in size and number. Selected bins with zero values were removed from both the query and background sample. P values are based on Welch’s t-test (two-sided). Black lines represent the mean. c, ChIP–seq data showing the enrichment (compared to input material) of the cohesin subunit SCC1 (n = 3 biologically independent experiments) across representative tRNA and rRNA genes (window size, 501 bp; step size, 101 bp). Black, red and blue boxes represent protein coding, tRNA and rRNA genes, respectively. Tick marks on the x axis represent 5-kb intervals. 3′B refers to one of the two alternative subtelomeric ends (A or B) at the 3′end of chromosome 11. d, ChIP–seq data showing cohesin (n = 3 biologically independent experiments), H3.V and H4.V (each in n = 2 biologically independent experiments) enrichment across three transcriptionally repressed BESs (window size, 2,001 bp; step size, 501 bp). Red flags mark BES promoters and black boxes indicate the locations of VSG genes.

Higher-order genome structures are established and maintained by architectural proteins such as CCCTC-binding factor (CTCF) and cohesin14. Histone variants are also enriched at many compartment boundaries15, but the role of these variants in shaping genome architecture remains unknown. Although CTCF appears to be absent in non-metazoans16, the major subunit of cohesin (SCC1) is present in T. brucei and the depletion of this subunit causes deregulation of VSG expression17. However, it has remained unclear whether this is a direct effect because SCC1 depletion strongly affects cell-cycle progression and growth rate, leading to rapid parasite death18.

Chromatin immunoprecipitation with sequencing (ChIP–seq) revealed that in T. brucei SCC1 is enriched across tRNA and rRNA genes, termination sites of RNA polymerase II transcription and most of the 3′ ends of BESs (Fig. 2c, d, Extended Data Fig. 4). This pattern of cohesin enrichment is reminiscent of its distribution in humans and yeast, in which cohesin is found at insulator and boundary elements such as tRNA genes19,20. The observed distribution of SCC1 is also similar to that of histone variants H3.V and, to a lesser extent, H4.V in T. brucei (Fig. 2d, Extended Data Fig. 4; also see ref. 21). This raised the possibility that these two histone variants function together with SCC1 in shaping genome organization and the regulation of antigen expression.

To investigate a possible link between these histone variants, genome architecture and antigen expression, we determined the expression of VSG genes and genome architecture in ΔH3.V, ΔH4.V and ΔH3.VΔH4.V cells. No cell cycle defect was observed in these cell lines (Extended Data Fig. 5).

Laboratory-adapted isolates, such as the one used here, switch their expression of VSG isoforms at very low frequency (about 10−6 per population doubling), and homogenously express VSG-2 (Fig. 3a; also see ref. 22). Thus, an increase in heterogeneity of VSG gene expression can be caused by a loss of mutually exclusive expression of VSG genes in individual cells—that is, heterogeneity in antigen expression at the single-cell level—or an increased switching frequency in expression of VSG genes in different parasites (heterogeneity at the population level).

Fig. 3: Deletion of histone variants H3.V and H4.V leads to a switch in expression of VSG isoforms.
figure 3

a, scRNA-seq analysis of wild-type (n = 40) and ΔH3.VΔH4.V (first time point, n = 44) cells. Each row represents data from one cell. The number of sequencing reads was normalized to account for differences in library size, gene length and uniqueness of VSG gene sequence. Only uniquely mapping reads were considered. The total number of VSG transcripts per cell is set to 100% (for details, see Methods, Extended Data Fig. 7). The colour code indicates the contribution of individual VSG transcripts to the pool of VSG transcripts in a single cell. The dominant VSG isoform is depicted with an orange border. For selected cells, the read coverage is shown across VSG-2, VSG-8 and VSG-11 (with 500 bp of surrounding sequence). CDS, coding sequence. b, Outline of the VSG switching mechanisms described for T. brucei. Green and red flags mark the active and repressed promoters, respectively. Green lines and grey bars indicate regions of expected transcription for the two different scenarios. c, Sequencing coverage across BES1 (left, top), BES15 (left, bottom) and a hybrid BES consisting of the 5′-BES1 and 3′-BES15 (right). Coverage is based on SMRT sequencing reads >10 kb from ΔH3.VΔH4.V gDNA that map to BES1 and BES15. The cross represents the site of recombination. Boxes represent expression-site-associated genes and ψ denotes a pseudogene. d, scRNA-seq-based analysis of ΔH3.VΔH4.V cells that exclusively express VSG-2 (n = 42) or VSG-11 (n = 82). The average transcript levels (counts per billion, cpb) based on uniquely mapping reads across BES1 and BES15 are shown. Grey bars represent degree of uniqueness. e, 4C-like inter-chromosomal interaction profiles (based on Hi-C-data, 20-kb bin size) showing the average interaction frequencies of BES1 (top) and BES15 (middle) with chromosomes 3, 6 and 8 in wild-type cells and the fold change (log2) in interactions of BES15 (bottom) with chromosomes 3, 6 and 8 after deletion of H3.V and H4.V.

To distinguish between these possibilities and to identify the VSG genes that are expressed, we performed single-cell RNA-seq (scRNA-seq) of individual T. brucei cells. scRNA-seq data from a total of 40 wild-type and 378 ΔH3.VΔH4.V cells revealed that—whereas all wild-type cells expressed VSG-2—in 74% of the ΔH3.VΔH4.V cells, VSG-2 transcript levels contributed less than 20% of the total VSG mRNA; this indicates a switch in expression of VSG genes (Fig. 3a, Extended Data Figs. 6, 7). Activation of new VSG genes was not random, with VSG-11 being the dominant newly activated VSG gene in 230 out of 378 cells. In addition, several cells contained transcripts from multiple VSG genes, which points to a partial loss of mutually exclusive expression. To determine the stability of VSG-2 expression, we analysed ΔH3.VΔH4.V cells at two time points that were about 50 population doublings apart. Although the overall pattern remained the same (Fig. 3a, Extended Data Fig. 6), the percentage of cells that expressed only VSG-2 mRNA, or multiple VSG mRNAs, had declined by the second time point. This suggests that the process of VSG-2 deactivation had progressed further, and that the simultaneous expression of multiple VSG genes may have been a transient intermediate state. Analyses based on immunofluorescence and flow cytometry confirmed that the loss of VSG-2 mRNA resulted in a loss of VSG-2 expression (Extended Data Fig. 8). No major effect on the expression of VSG genes was observed upon deletion of H3.V or H4.V alone (Extended Data Fig. 8).

In T. brucei, the switching of expression of VSG genes occurs by two distinct mechanisms11: either by switching transcription from one BES to another (transcriptional switch) or by a recombination-based event that leads to the replacement of the previously active VSG gene with a new VSG gene from a different genomic location (recombinational switch, Fig. 3b).

To gain insight into the mechanism by which histone variants affect antigen expression, we sequenced ΔH3.VΔH4.V genomic DNA using SMRT sequencing technology. The SMRT data indicated that, in most cells, recombination had occurred between an expression-site-associated gene 8 (ESAG8) gene pair that was present in both BES1 and BES15. The data also revealed that the new chimeric BES contained three copies of ESAG8, one from BES1 and two from BES15 (Fig. 3c). scRNA-seq and Hi-C data support a recombination event (Fig. 3d, e). Hi-C data revealed that, upon deletion of H3.V and H4.V, the interaction frequency between VSG-11 and the 5′ end of chromosome 6—where VSG-2 is located in wild-type cells—increased, indicating that VSG-11 had relocated to chromosome 6.

Studies in different organisms have shown that the frequency of recombination is affected by spatial proximity and DNA accessibility23,24. Thus, to determine whether histone variants contribute to genome architecture and/or local DNA accessibility, we performed Hi-C and assays for transposase-accessible chromatin using sequencing (ATAC-seq). Hi-C data from ΔH3.V cells revealed marked changes in inter-chromosomal interactions (Fig. 4a, top) and a significant increase in interactions among repressed BESs (Fig. 4b), pointing to a loss of constraints that may have ‘anchored’ the BESs to specific nuclear sites. In support of these Hi-C data, fluorescence in situ hybridization (FISH) data revealed a strong clustering of telomeric repeats upon deletion of H3.V (Fig. 4c, d). By contrast, deletion of H4.V affected genome architecture only modestly (Fig. 4a, bottom). Unlike the Hi-C data, our ATAC-seq data indicated that promoter-proximal DNA accessibility increased upon H3.V or H4.V deletion (Fig. 4e). However, only ΔH3.VΔH4.V cells exhibited high DNA accessibility across the entire length of transcriptionally repressed BESs (Fig. 4e bottom, Extended Data Fig. 9).

Fig. 4: Histone variants H3.V and H4.V influence global and local chromatin structures.
figure 4

a, Hi-C heat map showing the fold change (log2) in DNA–DNA interaction frequency between wild-type and ΔH3.V (top two panels) and ΔH4.V cells (bottom two panels). b, Fold change in DNA–DNA interaction frequency among VSG genes located in BESs compared to background DNA–DNA interactions in wild-type, in ΔH4.V, ΔH3.V and ΔH3.VΔH4.V cells. The ratios for each cell line were calculated for 100 randomly selected background regions. Mean ± s.d. Significance was determined using Welch’s t-test (two-sided). c, FISH with probes against telomeric repeats (Alexa Fluor488, green), n = 2 biologically independent experiments. Scale bar, 5 μm. d, Quantification of telomere signal to determine the fraction of cells containing large telomeric clusters (white arrows in c) was performed using Imaris 8 and is based on the analysis of 1,128 cells. Means ± s.d. of two replicates are shown (wild type: n = 116, 221; ΔH3.V: n = 140, 102; ΔH4.V: n = 146, 107 and ΔH3.VΔH4.V: n = 190, 106). e, ATAC-seq data (n = 2 biologically independent experiments) across BES15 (repressed in wild-type cells). gDNA read coverage (bottom) is shown to illustrate mappability of reads (window size, 501 bp; step size, 101 bp). Red flag and black box indicate the position of the promoter and the VSG gene, respectively. Tick marks on x axes represent 20-kb intervals. f, Model illustrating the influence of H3.V and H4.V on genome architecture and local DNA accessibility. H3.V and H4.V single knockouts alone mediate only partial opening of BESs (half open arrow) and H3.V knockout leads to a spatial rearrangement of BESs inside the nucleus, whereas deletion of both histone variants is required to obtain the fully opened BESs (open arrow) and spatial proximity of BESs that facilitate recombination (red cross) and lead to the expression of a new VSG isoform.

In summary, the Hi-C and ATAC-seq data indicate that although deletion of H3.V was responsible for the majority of genome architectural changes and increased BES clustering, this alone was not sufficient to induce a switch in expression of VSG genes. Only the concurrent deletion of H3.V and H4.V, which also strongly increased DNA accessibility across transcriptionally repressed BESs, enhanced the rate of recombination-based switching of VSG genes.

The depletion of histone H3 was previously shown to upregulate BES proximal-promoter activity—presumably via a general increase in DNA accessibility—but did not cause deregulation of VSG genes25. We hypothesize that the marked increase in switching frequency of VSG gene expression results from the combination of decreased spatial distance between BESs and increased local DNA accessibility (Fig. 4f).

The activation of new VSG genes did not occur at random; this non-random activation has previously been observed for infections of different hosts26,27. In a small number of cells, we detected transcripts from different VSG isoforms. This loss of mutually exclusive expression of VSG genes may be caused by increased DNA accessibility upon the deletion of histone variants, which may result in promiscuous RNA polymerase II transcription. Our observations that even in ΔH3.VΔH4.V cells not all expression sites are transcribed and that specific ‘pairs’ of VSG genes tend to be co-expressed, suggest that there are additional constraints imposed by genome organization or VSG protein structure28. At the genome level, co-activated VSG genes may have to be localized in close proximity to ensure sufficient levels of an activating factor8,29; alternatively, differences in VSG protein structure may make it impossible for the parasite to tolerate certain mosaic surface coats.

In this study, we have demonstrated how evolutionarily conserved features of genome architecture can be exploited for the de novo scaffolding of phased diploid genomes. The use of Hi-C, scRNA-seq and ATAC-seq—to our knowledge, all used here for the first time in T. brucei—opened opportunities for genome assembly and the characterization of the mechanism that underlies VSG switching in ΔH3.VΔH4.V cells. Our data reveal that histone variants can function as architectural proteins, and that changes in global genome architecture and local chromatin configuration can induce extensive switches in antigen expression.

Methods

No statistical methods were used to predetermine sample size. The experiments were not randomized and investigators were not blinded to allocation during experiments and outcome assessment.

Cell culture

All T. brucei strains used in this study are derivatives of the Lister 427 bloodstream-form isolate. Cells were cultured in HMI-11 medium (HMI-931 without serum plus) at 37 °C up to a density of 106 cells/ml. If required, drugs were used at standard concentrations.

Cell lines

ΔH3.V and ΔH4.V cells used in this study have previously been published21,32. After generation of a transgenic cell line, the correct tagging of a gene or the deletion of gene was verified by PCR. Cell lines were not tested for mycoplasma contamination.

Ty1-H3.V/ΔH3.V cells. To delete the first H3.V allele (Tb927.10.15350), the regions upstream and downstream of the H3.V CDS were PCR-amplified using the following primer pairs: H3.V_01_F, H3.V_02_R and H3.V_03_F, H3.V_04_R (see Supplementary Table 4 for a full list of oligonucleotides) and cloned into pyrFEKO-Puro33 using InFusion HD Cloning Plus reagents (Clontech) at PvuII/HindIII and BamHI/SbfI restriction sites. The resulting plasmid was linearized with PvuII and SbfI and stably transfected into the H3.V locus of T. brucei wild-type cells to generate H3.VH3.V cells. To add an N-terminal 2× Ty1 tag to the second H3.V allele, the sequence of 326-bp upstream of the H3.V CDS (H3.V_05_F, H3.V_06_R) was cloned into the ApaI/NotI linearized vector pPOTv3-2×Ty1 using InFusion HD Cloning Plus reagents (Clontech). Downstream of the blasticidin resistance marker and the Ty1-tag, a 417-bp DNA sequence homologous to the H3.V CDS 5′-end (H3.V_07_F, H3.V_08_R) was amplified (leaving out the ATG start codon) and likewise inserted using SacI and NheI restriction sites. The tag sequence was subsequently replaced by a codon-optimized version: oligonucleotides containing two codon-optimized Ty1 coding sequences (H3.V_09 and H3.V_10) were annealed, digested with HindIII and SacI and ligated into the HindII/SacI-linearized plasmid. Finally, the plasmid was linearized with ApaI and NheI restriction enzymes and stably transfected into H3.VH3.V cells to generate Ty1-H3.V/ΔH3.V.

ΔH3.VΔH4.V double-knockout cells. To delete H4.V, the upstream (H4.V_11_F, H4.V_12_R) and downstream (H4.V_13_F, H4.V_14_R) regions flanking the H4.V CDS (Tb927.2.2670) were amplified from bloodstream-form wild-type gDNA and purified using NucleoSpin Gel and PCR Clean-up kit (Macherey-Nagel). The PCR product of the upstream region was inserted into the plasmid pyrFEKO-Neo33 using InFusion HD Cloning Plus reagents (Clontech) between HindIII and AgeI restriction sites. The downstream region was integrated by ligation using BamHI and SbfI restriction sites. The neomycin resistance cassette was exchanged with a blasticidin or a phleomycin resistance marker, respectively, using the BglII and XbaI sites that flank the resistance marker. To this end, the blasticidin and phleomycin cassettes were excised from pyrFEKO-BSD or pyrFEKO-Phleo, purified and ligated into the target plasmid. The plasmids were linearized with NheI and SbfI and stably transfected into the previously published ΔH3.V cell line32.

Ty1-SCC1/ΔSCC1 cells. To delete the first SCC1 allele (Tb927.7.6900), the flanking regions upstream (Scc1_15_F, Scc1_16_R) and downstream (Scc1_17_F, Scc1_18_R) of the SCC1 CDS were amplified, digested with PvuII/HindIII and BamHI/SbfI, respectively, and ligated into pyrFEKO-Hyg33 at PvuII/HindIII and BamHI/SbfI restriction sites. Wild-type cells were transfected with the linearized plasmid (PvuII/SbfI) to obtain SCC1/ΔSCC1 cells. For the N-terminal Ty1-tagging of the second SCC1 allele, the 3′ end of the SCC1 5′ UTR was amplified (Scc1_19_F, Scc1_20_R), digested with ApaI and NotI and ligated into the ApaI/NotI-linearized vector pPOTv3-2×Ty1. Next, the 5′ end of the SCC1 CDS was amplified (leaving out the ATG start codon) (Scc1_21_F, Scc1_22_R), digested with SacI and NheI and ligated into the likewise-digested plasmid. The Ty1-tag was exchanged by a codon-optimized version as described for N-terminal tagging of H3.V (see above). The ApaI/NheI linearized plasmid was stably transfected into SCC1SCC1 T. brucei cells to generate Ty1-SCC1/ΔSCC1.

The cell line in which both endogenous H4.V alleles are knocked out and ectopic overexpression of a Ty1-tagged version of H4.V can be induced has previously been published21.

In situ Hi-C

Because ΔH3.V32 and ΔH4.V21 cells had been generated in a ‘single marker’ background34, we generated the ΔH3.VΔH4.V cells in a single marker background and compared the Hi-C profiles of the transgenic cell lines to those generated from single marker cells. Thus, all ‘wild-type’ Hi-C data shown in Figs. 24 and Extended Data Figs. 3, 4 are generated from single marker cells. Hi-C data from ‘true’ wild-type cells (Lister 427, MiTat 1.2) were also generated, but used only for the genome assembly.

In situ Hi-C was performed based on previously published protocols35,36 and adapted to T. brucei: 2 × 108 cells (wild type, single marker, ΔH3.V, ΔH4.V and ΔH3.VΔH4.V) were collected and resuspended in 40 ml of 1× trypanosome dilution buffer (1× TDB; 0.005 M KCl, 0.08 M NaCl, 0.001 M MgSO4 ×7H2O, 0.02 M Na2HPO4, 0.002 M Na2HPO4 ×2H2O, 0.02 M glucose). Cells were fixed in the presence of 1% formaldehyde for 20 min at room temperature by addition of 4 ml of formaldehyde solution (50 mM Hepes-KOH pH 7.5, 100 mM NaCl, 1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, 11% formaldehyde). The reaction was stopped by addition of 3 ml of 2 M glycine and incubation for 5 min at room temperature and 15 min on ice. Cells were washed twice in 1× TDB and the cell pellet was snap-frozen in liquid nitrogen. Cells were resuspended in 1 ml of permeabilization buffer (100 mM KCl, 10 mM Tris pH 8.0, 25 mM EDTA) supplemented with protease inhibitors (1.5 mM pepstatin A, 4.25 mM leupeptin, 1.06 mM PMSF, 1.06 mM TLCK) and digitonin (200 μM final concentration) and incubated for 5 min at room temperature. Cells were washed twice in 1× NEBuffer3.1 (NEB, B7003S) and resuspended in 342 μl of 1× NEBuffer3.1. After addition of 38 μl of 1% SDS, and an incubation at 65 °C for 10 min, the SDS was quenched by addition of 43 μl of 10% Triton-X 100 (Sigma) and the incubation was continued at room temperature for 15 min. Another 35 μl of water, 13 μl of 10× NEBuffer3.1 and 100 units of MboI (NEB, R0147M) were added and the chromatin was digested at 37 °C overnight while shaking. To inactivate MboI, the sample was incubated at 65 °C for 20 min. Restriction fragments were biotinylated by supplementing the reaction with 60 μl of fill-in mix (0.25 mM biotin-14-dATP (Life Technologies, 19524016), 0.25 mM dCTP, 0.25 mM dGTP, 0.25 mM dTTP (Fermentas), 40 U of DNA polymerase I, large (Klenow) fragment (NEB, M0210)) and incubation at 23 °C for 4 h. The end-repaired chromatin was transferred to 665 μl of ligation mix (1.8% Triton-X 100, 0.18 mg BSA, 1.8× T4 DNA Ligase Buffer (Invitrogen, 46300018) and 5 μl of T4 DNA ligase (invitrogen, 15224025) were added. The ligation was performed for 4 h at 16 °C with interval shake. Crosslinks were reversed by adding 50 μl of 10 mg/ml proteinase K (65 °C for 4 h) and another addition of 50 μl of 10mg/ml proteinase K, 80 μl of 5M NaCl and 70 μl of 10% SDS (65 °C, overnight).

The DNA was precipitated with ethanol and resuspended in 257 μl of TLE (10 mM Tris-HCl, 0.1 mM EDTA, pH 8.0). SDS was added to a final concentration of 0.1% and the sample was split among two tubes for sonication (Covaris S220; microtubes, 175 W peak incident power, 10% duty factor, 200 cycles per burst, 240 s treatment). The samples were recombined and the volume was adjusted to 300 μl with TLE. Fragments between 100 and 400 bp in size were selected using Agencourt AMPure XP beads (Beckman Coulter), according to the manufacturer’s instructions. The DNA fragments were eluted off the beads in 55 μl of TLE.

For end-repair and biotin removal from unligated ends, 70 μl of end-repair mix was added (1× Ligation buffer (NEB), 357 μM dNTPs, 25U T4 PNK (NEB, M0201), 7.5U T4 DNA polymerase I (NEB, M0203), 2.5U DNA polymerase I, large (Klenow) fragment (NEB, M0210)) and incubated for 30 min at 20 °C and 20 min at 75 °C. To inactivate the enzymes, EDTA was added to a final concentration of 10 mM. To isolate biotin-labelled ligation junctions, 50 μl of 10 mg/ml Dynabeads MyOne Streptavidin C1 (Life Technologies, 65001) were washed with 400 μl of 1× Tween washing buffer (TWB; 5 mM Tris-HCl pH 7.5, 0.5 mM EDTA, 1 M NaCl, 0.05% Tween 20), collected with a magnet, resuspended in 400 μl of 2× binding buffer (10 mM Tris-HCl pH 7.5, 1 mM EDTA, 2 M NaCl) and added to the sample suspended in 330 μl TLE. Biotinylated DNA was bound to the beads by incubating the sample for 15 min at room temperature with slow rotation. Subsequently, the DNA-bound beads were captured with a magnet, washed twice with 400 μl of 1× binding buffer, washed once in 100 μl of 1× TLE T4 ligase buffer and resuspeded in 41 μl of TLE. For polyadenylation, 5 μl of 10× NEBuffer 2.1, 1 μl of 10 mM dATP and 3 μl of 5 U/μl of Klenow fragment (3′→ 5′ exo-) (NEB, M0212) and incubated for 30 min at 37 °C followed by deactivation for 20 min at 65 °C. Beads were reclaimed with a magnet, washed once with 400 μl 1× Quick ligation buffer (NEB, M2200) and resuspended in 46.5 μl of 1× Quick ligation buffer (NEB, M2200). 2.5 μl of DNA Quick ligase (NEB, M2200) and 0.5 μl of 50 μM annealed TruSeq adapters were added and incubated for 1 h at room temperature. Beads were separated on a magnet, resuspended in 400 μl of 1× TWB (5 mM Tris-HCl, 0.5 M EDTA, 1 M NaCl, 0.05% Tween-20) and washed for 5 min at room temperature with rotation. Beads were washed on the magnet with 200 μl 1× binding buffer and 200 μl of 1× NEBuffer 2.1 and resuspended in 20 μl of 1× NEBuffer 2.1. The library was amplified in 8 separate reactions of 50 μl. Per reaction, 1.5 μl of 25 μM TruSeq PCR primer cocktail (TruSeq PCR primer cocktail_F, TruSeq PCR primer cocktail_R; see Supplementary Table 4), 25 μl of 2× Kapa HiFi HotStart Ready Mix (Kapa Biosystems, KR0370) and 21.5 μl of water were added to 2 μl of library bound to the beads. Amplification was performed as follows: 3 min at 95 °C, 5 cycles of 20 s at 98 °C, 30 s at 63 °C and 30 s at 72 °C, 1 cycle of 1 min at 72 °C, hold at 4 °C. The PCR reactions were pooled and the beads were removed from the supernatant using a magnet. The library was purified by addition of 1.5 volumes of Agencourt AMPure XP beads (Beckman Coulter), according to the manufacturer’s instructions. The sample was eluted off beads using 25 μl of 1× TLE buffer, transferred to a fresh tube and the concentration was determined using Qubit (Qubit dsDNA HS Assay Kit, Thermo Fisher) and qPCR (KAPA SYBR FAST qPCR Master Mix, Kapa Biosystems), according to the manufacturer’s instructions. Library size distributions were determined on a 5% polyacrylamide gel. Paired-end 76-bp sequencing was carried out using the Illumina NextSeq 500 system with high and mid output NextSeq 500/550 kits according to the manufacturer’s instructions.

Mapping of Hi-C reads and generation of interaction matrices

Reads were trimmed at their ligation junction using the truncator of the HiCUP pipeline37 (version 0.5.9 devel), mapped using bwa mem (https://arxiv.org/abs/1303.3997) of the bwa kit version 0.7.12-r1039 and those with a quality score q > 0 were retained. Reads were mapped to the T. brucei Lister 427 genome version 9 (Tb427v9). For each chromosome, this genome contains the core region of two homologous chromosomes only once, whereas both of the the respective heterozygous subtelomeric regions belonging to the two homologous chromosomes are included. During the mapping of Hi-C reads, contigs that displayed alternative variations of an assembled allele were removed (Tb427v9_without_allelic_variants) to keep these loci visible in the Hi-C matrices. The primary analysis of Hi-C reads was performed with HiC-Pro38 (version 2.10.0) to visually inspect reproducibility among replicates. Raw matrices were normalized for differences in ploidy (for example, read counts at diploid regions were multiplied by 0.5), balanced by iterative correction using HiC-Pro (default settings) and converted into a homer compatible format39 using a custom Python script (see ‘Code availability’). To enable comparisons between different Hi-C experiments, each value in the balanced interaction matrix was divided by the respective column sum.

Distance-dependent decay of interaction frequencies

To visualize the distance-dependent decay, interaction frequencies between Hi-C bin pairs were grouped on the basis of the linear distance between the pairs, and the distribution of the median interaction frequency across distances was plotted.

Co-localization of genomic loci

To determine whether a region of interest interacted more or less than expected by chance, the median and mean interaction frequencies were calculated for bins overlapping with the regions of interest. In addition, a ‘background’ interaction frequency was determined by randomly selecting regions of identical size from the same matrix. Significance was determined using Welch’s t-test.

To identify changes in DNA–DNA interaction frequencies after deletion of histone variants, the ratio of ‘feature median interaction frequencies’/‘background median interaction frequency’ was determined for 100 randomly selected background regions in wild-type, ΔH3.V, ΔH4.V and ΔH3.VΔH4.V cells. Significance was determined using Welch’s t-test. Co-localization analyses were performed using balanced interaction matrices (50-kb resolution). Bins with zero values were excluded from the analyses.

4C-like analysis

To visualize interactions between one genomic region (bait) and all other genomic sites, relevant bins were extracted from a 20-kb Hi-C matrix. An average interaction value for every genomic bin was calculated if the bait regions spanned more than one bin.

Hi-C matrix visualization

Matrices were plotted based on the colour palettes provided by seaborn (https://seaborn.pydata.org). To generate differential heat maps, a pseudo-count of 0.000001 was added to each interaction value of the numerator and denominator matrix before division.

SMRT sequencing

Genomic DNA was isolated and precipitated from 3 × 108 cells of the T. brucei 427 17.13 P10 isolate40 using the Blood & Cell Culture DNA Midi Kit (Qiagen), according to the manufacturer’s instructions. In addition, the DNA was purified in a phenol chloroform extraction using Manual Phase Lock Gel 2-ml (heavy/light) tubes (5Prime) and suspended in 100 μl of TE buffer. SMRT library preparation and sequencing was performed at the Icahn School of Medicine at Mount Sinai.

Genomic DNA library preparation and sequencing was performed primarily using the manufacturer’s instructions for the P6-C4 sequencing enzyme and chemistries. In short, ~5 μg gDNA was quantified and diluted to 150 μl using elution buffer (Qiagen) at 33 ng/μl and then sheared to ~20 kb by centrifugation at 4,500 rpm for 50 s using a G-tube spin column (Covaris). The sheared DNA was then re-purified using Agencourt AMPure XP beads (Beckman Coulter) at 0.45×. Next, ~1.6 to 3.2 μg of DNA from each batch was taken into DNA damage and end repair. In brief, the DNA fragments were repaired by adding 21.1 μl of DNA damage repair solution (1× DNA damage repair buffer (1× NAD+, 1 mM ATP and 0.1 mM dNTP) and 1× DNA damage repair mix) and incubation at 37 °C for 20 min. DNA ends were repaired by adding 1× end repair mix to the solution and incubation at 25 °C for 5 min, followed by an additional 0.45× Agencourt AMPure XP purification step. Next, 0.75 μM of blunt adaptor was added to the DNA and 1× template preparation buffer (0.05 mM ATP and 0.75 U/μl T4 DNA ligase) was added to a final volume of 47.5 μl. This solution was incubated at 25 °C overnight, followed by incubation at 65 °C for 10 min to inactivate the ligase. To remove un-ligated DNA fragments, exonuclease cocktail (1.81 U/μl Exo III 18 and 0.18 U/μl Exo VII) was added to the library followed by a 60 min incubation at 37 °C. Two additional 0.45× Agencourt Ampure XP purification steps were performed to remove <2,000-bp molecular weight DNA and organic contaminants.

The size of the library was validated using an Agilent DNA 12000 chip. Before P6-C4-based sequencing, Blue Pippin size selection was applied to remove molecules <7,000 bp. This step was conducted using Sage Science Blue Pippin 0.75% agarose cassettes to select libraries in the range of 7,000–50,000 bp. Primers were annealed to the size-selected SMRTbell at a ratio of 20:1 with the full-length libraries by denaturation (80 °C for 2 min) and slow cooling (0.1 °C/s to 25 °C). The polymerase-template complex was bound to the P6 enzyme using a ratio of 10:1 polymerase to SMRTbell at 0.5 nM for 4 h at 30 °C and then held at 4 °C until ready for magnetic-bead loading. The magnetic-bead loading step was conducted at 4 °C for 60 min. The magnetic-bead-loaded, polymerase-bound SMRTbell libraries were placed onto the RSII machine at a sequencing concentration of 50 pM and configured for a 240-min continuous sequencing run.

Assembly, post-assembly improvements and genome scaffolding

The 642,583 sequencing reads from wild-type gDNA (seven SMRT cells) were assembled into contigs following the RS_HGAP_Assembly.3 workflow from SMRT Analysis v2.3.041 with default parameters. In brief, a sequence seeding dataset with the longest sequencing reads was pulled apart. The remaining reads were mapped onto them to obtain error-corrected reads through a consensus procedure. The error-corrected reads were assembled by traditional overlap layout consensus with a Celera Assembler42. Contig sequences were polished with Quiver and were further joined and extended using PBJelly 2 (PBSuite_15.2.20)43. Remaining sequence errors were corrected using iCORN244 with previously published and new gDNA Illumina data available under GSM2586510 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2586510) and ERS1503958 (http://www.ebi.ac.uk/ena/data/view/ERS1503958), respectively.

For genome scaffolding, all Hi-C reads were combined, mapped onto the error-corrected HGAPv3 contigs and a balanced 10-kb bin heat map matrix of DNA–DNA interactions was generated as previously described. For the scaffolding, only a subset of the matrix that contained contigs larger than 50 kb was considered (Extended Data Fig. 1).

Then, contigs that did not exhibit a gradual distance-dependent decay in DNA–DNA interactions (which suggests that they may have been incorrectly assembled) were broken at the site of the putative mis-assembly. Available pipelines for the scaffolding of contigs based on DNA–DNA interaction frequencies are not designed for genomes that contain large regions of heterozygosity45,46; therefore, contigs were manually rearranged and/or inverted so that long-distance DNA interactions—visible as signal away from the diagonal—were reduced. The process of Hi-C read-mapping, heat map generation and contig repositioning or inversion was repeated until no contigs could be identified, the repositioning of which would have further minimized the signal away from the diagonal (Extended Data Fig. 1). To validate our scaffolding approach and to further improve the genome, we ran PBJelly 2, which enabled us to join together several of the contigs that we had placed next to each other and to reduce the number of contigs from 139 to 91. In addition, we compared the obtained Lister 427 scaffold to that of the previously assembled TREU 927 genome and found that the core regions exhibited a high degree of similarity and synteny (Extended Data Fig. 2), which validates our scaffolding approach. The subtelomeric regions are known to be different between the genomes. Finally, we observed that 27 out of the 33 core–subtelomere boundaries in the assembly are spanned by PacBio reads and/or contigs, which supports their linear proximity (see Supplementary Information).

For a comparison of different assembly strategies and an assessment of the genome quality, see Supplementary Information and Supplementary Table 3. The genome-built Tb427v9, which was used for all analyses performed in this study, is available in the European Nucleotide Archive with the ENA study accession number PRJEB18945 (http://www.ebi.ac.uk/ena/data/view/PRJEB18945).

Genome annotation

Annotations from the T. brucei TREU 927 genome were transferred using Companion (https://companion.sanger.ac.uk; accessed 9 January 2017)47. In addition, VSG genes were annotated based on similarity (>90% coverage and >95% identity to a VSG gene) of the available VSGnome dataset from the Lister 427 strain10 with NCBI-BLAST+ (version 2.2.31+). Overlapping VSG-gene entries were merged and named after the entry with the best bit score. For the identification of putative novel MESs, the assignment of BESs to chromosome ends and the identification of centromeres, see Supplementary Information.

Chromatin immunoprecipitation with sequencing

Assays were performed as previously described48 with minor alterations. For the immunoprecipitation of cell lines with Ty1-tagged proteins (Ty1-H3.V/ΔH3.V, Ty1-H4.V ΔH4.V/ΔH4.V and Ty1-SCC1/ΔSCC1 cells), 50 μl of Dynabeads protein G (ThermoFisher) was separated on a magnet, the supernatant was removed and the beads were resuspended in 200 μl of PBS-Tween (0.02%) containing 10 μg of BB2 antibody49, and then incubated with slow rotation at 4 °C overnight. Antibody-coupled protein G beads were separated on a magnetic rack and washed three times with PBS-Tween (0.02%). Five hundred microlitres of ChIP sample was added, and incubated at 4 °C overnight with slow rotation. For ChIP assessment of the distribution of histone variants, the DNA was fragmented using micrococcal nuclease; for ChIP assessment of the distribution of SCC1, the DNA was sheared using a Covaris S220 instrument before target protein binding. Sequencing reads were mapped using bwa mem and coverage plots were generated using COVERnant (v.0.3.2) (https://github.com/konrad/COVERnant)13.

ATAC-seq

To ensure reproducibility of the assays independent of cell-number variations, all assays were performed with 1 × 106 and 2 × 106 cells. To this end 3 × 107 cells were collected and washed in 30 ml of cold 1× TDB. The pellet was resuspended in 300 μl permeabilization buffer with protease inhibitors, 3 μl of 4 mM digitonin was added and incubated for 5 min at room temperature. The cells were pelleted, resuspended in 600 μl isotonic buffer with protease inhibitors and split in two samples, containing 1 × 107 and 2 × 107 cells, respectively. The transposition reaction was performed by adding 50 μl of transposition mix to the pellet (25 μl TD (2× reaction buffer from Nextera kit), 25 μl TDE1 (Tn5 transposase from Nextera kit), 22.5 μl nuclease-free water) and incubation for 30 min at 37 °C. For the gDNA control, 200 ng of gDNA was treated in the same manner. The DNA samples were purified using Qiagen MinElute PCR Purification Kit and eluted in 10 μl EB (10 mM Tris-HCl, pH 8). The transposed DNA fragments were amplified using the NEBNext High-Fidelity 2× PCR Master Mix (M0541) supplied with 2.5 μl of 25 μM barcoded primers and amplification for 13 cycles. The libraries were purified using AMPure XP beads (Beckman Coulter) according to the manufacturer’s instructions. The library fragment sizes between 150 and 1,000 bp were purified from a 6% polyacrylamide gel. Paired-end 76-bp sequencing was carried out using the Illumina NextSeq 500 system with a high-output NextSeq 500/550 kit, according to the manufacturer’s instructions. Sequencing reads were mapped using bwa mem and coverage plots were generated using COVERnant (v.0.3.2) (https://github.com/konrad/COVERnant)13.

scRNA-seq

T. brucei wild-type and ΔH3.VΔH4.V cells were sorted (0 cell, 1 cell or 50 cells) using a FACSAria III (BD Biosciences; precision: single-cell; nozzle: 100 μm). A forward-scatter area versus side-scatter area plot was used to gate and sort the cells. T. brucei cells were sorted in 48-well plates (Brand) filled with 2.6 μl of 1× lysis buffer (Takara) supplemented with 0.01 μl of RNase inhibitor (40 U/μl; Takara). Immediately after sorting, cells were placed on ice for 5 min and stored at −80 °C.

Lysates from 50 trypanosomes and single trypanosomes were supplemented with 0.2 μl of a 1:2 × 106 (scRNA-seq I) or a 1:20 × 106 (scRNA-seq II) dilution of ERCC Spike-In Control Mix 1 (Thermo Fisher, 4456740). Libraries were prepared using SMART-Seq v.4 Ultra Low Input RNA Kit (Takara) using a quarter of the reagent volumes recommended by the manufacturer. PCR amplification was performed using 26 cycles (scRNA-seq I) or 22 cycles (scRNA-seq II), according to the supplier’s recommendations. cDNA was purified using Agencourt AMPure XP beads (Beckman Coulter) and recovered in 15 μl of elution buffer (Takara). Libraries were quantified using the Qubit 3 Fluorometer with dsDNA Hs Assay kit (Life Technologies) and the quality of the libraries was assessed using a 2100 Bioanalyzer with High Sensitivity DNA kit (Agilent) (Extended Data Fig. 7). Similar to previously published approaches50, 0.5 ng of cDNA was subjected to a tagmentation-based protocol (Nextera XT, Illumina) using a quarter of the recommended volumes, 10 min for tagmentation at 55 °C and 1 min extension time during PCR amplification. Libraries were pooled (96 libraries for NextSeq) and sequencing was performed in paired-end mode for 2 × 75 cycles using Illumina’s NextSeq 500. Details of the sequencing results are listed in Supplementary Table 2.

Analysis of scRNA-seq

The reads were mapped to the combination of the Tb427v9 genome with the ERCC spike-in sequences, using bwa mem, version 0.7.16. The mapped data were processed using samtools51 (version 1.8) and MarkDuplicates tool (2.18.3-SNAPSHOT) from Picard (http://broadinstitute.github.io/picard/), and the read counts-to-features were done with bedtools52 (version 2.26.0). To assess the quality of the data, for each scRNA-seq experiment the read counts to the following groups were determined: reads mapping to ‘rRNA genes’, reads mapping to ‘protein-coding genes’ (CDS plus 89 bp for the 5′ UTR and 400 bp for the 3′ UTR)53, reads mapping to ‘other regions of the genome’ and ‘unmapped reads’. Reads that overlapped ‘rRNA genes’ and ‘protein-coding genes’ features at the same time were excluded from both groups and counted as reads mapping to ‘other regions of the genome’.

To assess the library complexity of the scRNA-seq datasets, the number of genes with ≥ 10 read counts was determined for each cell. The counts per feature group, as well as the number of genes captured in each scRNA-seq experiment, are available in Supplementary Table 2. Only those scRNA-seq datasets with more than 500 genes with ≥ 10 read counts were considered for the quantification of VSG gene expression.

Many VSG genes share a high degree of homology with each other. Therefore, to determine the expression levels for each of the 2,846 VSG genes annotated in Tb427v9, the number of uniquely mapping sequence reads obtained for each VSG gene was normalized to account for differences in uniqueness. The uniqueness of the VSG genes was determined by alignment of an in silico-generated dataset that matched the scRNA-seq datasets in read size and fragment-length distribution. For each cell, the transcript level of an individual VSG gene is shown as a percentage of the total VSG gene transcript level in that cell. Raw and normalized read counts are available in Supplementary Table 2, and a diagram explaining the VSG count normalization procedure is shown in Extended Data Fig. 7.

Total RNA-seq

Triplicates of T. brucei wild-type, ΔH3.V, ΔH4.V and ΔH3.VΔH4.V cells were grown to a density of ~106 cells/ml. Cell concentration was determined nine times for each replicate using a Coulter Cell Counter (Beckman Coulter) and 4.5 × 107 cells were collected from each culture. Cells were washed with 1× TDB, resuspended in 350 μl of buffer RA1 (Macherey-Nagel, NucleoSpin RNA) and 38 μl of 0.1 M DTT and 1 μl of a 1:10 dilution of ERCC Spike-In Control Mix (Thermo Fisher, 4456740) were added. Total RNA was purified from the lysate according to the NucleoSpin RNA kit protocol, and eluted in 30 μl of nuclease-free water. To deplete rRNA, 3.5 μg of total RNA was mixed with 2.6 μl of 5× hybridization buffer (500 mM Tris-HCl pH 7.4, 1 M NaCl), 0.459 pmol of 131 50-bp anti-rRNA oligonucleotides (covering 18S, M1, M2, M3, M4, M5, M6, 28S alpha and 28S beta rRNAs, a kind gift of C. Clayton; for the full list, see Supplementary Table 4) and 2.2 μl of water, denatured at 95 °C for 2 min and slowly (0.1 °C per s) cooled to 37 °C. One microlitre of prewarmed 10× RNaseH digestion buffer (200 mM Tris-HCl, pH 7.4, 500 mM KCl, 40 mM MgCl2, 10 mM DTT) and 10 U of RNaseH (ThermoFisher, AM2292) were added and the volume was adjusted to 16 μl with nuclease-free water. The mixture was incubated at 37 °C for 20 min. Residual oligonucleotides were subsequently digested by addition of 2 U of Turbo DNaseI and 5 μl of 10× Turbo DNase reaction buffer (AM2238, Thermo Fisher) in a total reaction volume of 50 μl and incubation at 37 °C for 20 min. The DNase was inactivated by addition of EDTA (15 mM final concentration) and heating at 75 °C for 10 min. The RNA was purified using RNAeasy Minelute columns (Qiagen), according to the manufacturer’s instructions, and eluted in 14 μl of RNase-free water. Double-stranded cDNA was synthesized from 100 ng of rRNA-depleted RNA using the NEBNext Ultra RNA Library Prep Kit for Illumina. The double-stranded cDNA was purified using Agencourt AMPure XP beads (Beckman Coulter) and eluted in 60 or 30 μl of 0.1× TE Buffer, respectively. Libraries were prepared and sequenced as previously described48.

Analysis of RNA-seq data

After adaptor clipping and quality trimming using cutadapt54 (version 1.10), the RNA-seq reads were mapped against the T. brucei genome using bwa mem of the bwa kit version 0.7.16. Only reads with a quality score q > 0 were retained. Feature quantification was performed with bedtools multicov subcommand. Differential gene expression analysis was then conducted with DESeq2 (v.1.20.0)55. Features with an adjusted P value (calculated based on Wald test and adjusted for multiple testing using the procedure of Benjamini and Hochberg56) below 0.1 were considered as differentially expressed. Supplementary Table 1 contains the raw counts for each gene in each individual RNA-seq replicate, as well as the fold change (in log2 scale) and P value adjusted for each sample versus wild type.

Fluorescence in situ hybridization

For each FISH assay, 1 × 107 bloodstream-form trypanosomes grown to a density of up to 1 × 106 cells/ml were collected and the cell pellet was washed once with 1× TDB and fixed for 15 min in 1× TDB containing 4% formaldehyde. Cells were washed with 1× TDB once and resuspended in 50 μl of 1× TDB. Gene frames were placed onto microscopy slides to cover an aminopropyltriethoxysilane-coated coverslip. Cells were pipetted onto the framed coverslip and allowed to settle for 5 min by gravity. The sample was washed twice for 2 min with 90 μl of 1× TDB, incubated with quenching solution (1× TDB containing 1 mg/ml NaBH4) for 10 min, washed twice with 1× TDB, permeabilized with 70 μl of 1× TDB containing 0.1% NP-40 for 5 min and washed twice with 1× TDB. Next, cells were treated with 1× PBS containing 1 mg/ml RNaseA for 30 min, washed twice with 1× TDB and incubated with 50 μl of a 1:1 dilution of hybridization buffer with 2× SSC for 30 min at room temperature. The labelled probe (see Supplementary Table 4) was diluted to a final concentration of 400 nM in 25 μl of hybridization buffer (50 (v/v) formamide, 10% (w/v) dextran sulfate, 2× SSPE, 250 μg/ml herring sperm DNA). The hybridization buffer was removed from the sample, 25 μl of hybridization solution containing the probe was added and the frame was closed with a plastic lid. The sample was incubated using a thermal block at 90 °C for 5 min and at 37 °C overnight. Next, the samples were washed for 30 min in 30 ml of 50% (v/v) deionized formamide and 2× SSC at 37 °C, for 10 min in 30 ml of 1× SSC at 50 °C, for 10 min in 30 ml of 2× SSC at 50 °C and for 10 min in 30 ml of 4× SSC at room temperature. Subsequently, cells were blocked in P1 buffer (100 mM maleic acid, 150 mM sodium chloride, pH 7.5, 4% BSA, 1% milk) for 1 h, incubated with primary antibody (sheep anti-digoxigenin antigen-binding fragment, Roche, diluted 1:2,000 in P1) for 45–60 min and washed with 0.5% Tween-20 and PBS 4 times for 4 min each. Cells were then incubated with the secondary antibody (Alexa Fluor 488 conjugated donkey anti-sheep IgG (H + L), Life Technologies, A11015, diluted 1:2,000 in P1) for 30–60 min. For further signal amplification, the samples were washed with PBS containing 0.5% Tween-20 4 times for 4 min and incubated with a rabbit anti-donkey IgG (H + L) DyLight 488 (Invitrogen) (diluted 1:2,000 in P1) for 30-60 min. The samples were washed twice for 10 min in 30 ml of PBS containing 0.5% Tween-20 and another 10 min in 30 ml of PBS, each time in a falcon tube on a shaker. Samples were mounted with 36 μl of Vectashield Mounting Medium with DAPI (Biozol) on a microscopy slide and were sealed with nail polish.

Immunofluorescence

Immunofluorescence was performed as previously described57, with minor alterations. Ten million cells per ml (wild type, ΔH3.V, ΔH4.V and ΔH3.VΔH4.V) were suspended in HMI-11 containing 2% formaldehyde, incubated for 5 min at room temperature and washed with 1× TDB. α-Tubulin was stained using the mouse monoclonal antibody Tat158 (1:200) and a secondary Alexa Fluor 594-conjugated chicken anti-mouse IgG (1:350, Invitrogen). VSG-2 was stained using CRD-depleted rabbit anti-VSG-240 (1:1,000) and a secondary Alexa Fluor 488-conjugated donkey anti-rabbit IgG (1:350, Invitrogen).

Fluorescence microscopy and image analysis

For imaging, a wide-field fluorescence Leica DMI6000 microscope with a mercury metal halide lamp and a HCX PL APO CS 100×/1.47 OIL objective was used. Images were captured with a Leica DFC 360 FX camera. Stacks with 32 slices and 6.3232 μm in height (0.0645 × 0.0644 × 0.1976 μm3 voxel size) were captured using identic exposure times for all conditions.

Quantification of ‘large’ telomere clusters was carried out using Imaris 8 software (Oxford Instruments). After segmenting individual nuclei in the DAPI channel, surfaces were rendered for the telomere FISH signal, while setting the quality filter > 1,000. All FISH signals that generated surfaces with a volume > 0.3 μm3 were classified as large telomere clusters, and scored.

FACS flow cytometry

VSG expression on the cell surface was quantified according to a previously published protocol40. In brief, 1 × 106 cells were centrifuged in a chilled microcentrifuge at 1,500g for 4 min at 4 °C. Cells were resuspended in 100 μl of ice cold HMI-11 and a VSG-specific antibody (anti-VSG-2 or anti-VSG-13)40 was added. After 60 min of incubation at 4 °C with gentle shaking, cells were washed three times in 500 μl of ice cold HMI-11, resuspended in 100 μl of cold HMI-11 and incubated with an Alexa Fluor 488-conjugated secondary antibody for 20 min. The cells were washed twice with 500 μl of 1× TDB and finally resuspended in 400 μl of 1× TDB before analysis with a FACSort flow cytometer (Becton Dickinson Biosciences).

To determine the cell-cycle profiles, 5 × 106 cells were collected (10 min, 1,300g, 4°C) and washed once with ice-cold 1× TDB. The cells were resuspended in 1 ml ice-cold PBS and 2 mM EDTA and fixed by adding 2.5 ml ice-cold methanol. After a 1-h incubation, cells were washed with 1 ml PBS and EDTA at room temperature and resuspend in 1 ml PBS and EDTA. One microlitre RNaseA (10 μg/μl) and 10 μl propidium iodide (1 μg/μl) were added. The stained cells were analysed with a FACSCalibur (Becton Dickinson) after a 30-min incubation at 37 °C.

Mapping the site of recombination

A library from ΔH3.VΔH4.V gDNA was prepared as described for wild-type gDNA and sequenced using the PacBio Sequel system by diffusion at 5 pM, with v2.1 chemistry and a 10-h movie. Reads > 10 Kb were extracted, split into ~2,500-bp chunks and mapped independently to the genome. Reads that contained chunks that mapped to BES1 and BES15 were kept and mapped again to BES1 and BES15. Based on the observed mapping pattern, a BES1–BES15 hybrid was constructed.

Code availability

Workflows and custom-made Unix Shell, Python and R scripts have been deposited at Zenodo (https://doi.org/10.5281/zenodo.823671). Documentation to reproduce the data analysis is provided.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.