Cell-cycle dynamics of chromosomal organisation at single-cell resolution

Improved multiplexed single-cell Hi-C

We redesigned our original single-cell Hi-C approach¹² to employ a 4-cutter restriction enzyme, Tn5-transposase tagmentation, flow cytometry sorting of single cells into 96-well plates and partially automated library preparation (Fig. 1a). We processed a total of 2924 single F1 hybrid 129 × Castaneus mouse ESCs grown in 2i media using 1.5 million reads per cell on average (Extended Data Fig. 1a-c). Analysis of various quality control (QC) metrics (Extended Data Fig. 1d-i) confirmed the robustness of the protocol over different batches, and in particular identified karyotypic imperfections within the single-cell sample (Extended Data Fig. 1j). Overall, 1992 cells (68.1%) passed stringent quality control filters for further analysis. This provided median coverage of 393,506 restriction fragments, and 127,233 distinct >1kb contacting pairs per cell (Fig. 1b). We found the incidence of trans-chromosomal contacts in the single-cell maps was low (median 5.87%, Fig. 1c), indicative of high library quality¹³ and reflecting specific and restricted interfaces of chromosome territories (Extended Data Fig. 2). This data provided us with the opportunity to study variation in chromosomal conformation across a large single-cell pool, on the basis of a high-resolution reference ensemble map (Extended Data Fig. 3)

Integrated cell-cycle phasing

Mammalian mitotic chromosomes are rod-shaped, and bulk Hi-C on sorted M phase cells has revealed a characteristic scale of cis-chromosomal contact distances peaking between 2-12 Mb¹¹. We examined the cells with the highest fraction of contacts at this mitotic range, and found their overall distribution of contact distances was similar to that previously observed for mitotic cells (Extended Data Fig. 4a). These single-cell mitotic maps were enriched for striking patterns of trans-chromosomal alignments (Fig. 1d, Extended Data Fig. 4b), suggesting that the chromosomes interacted along their entire length (head-to-head and tail-to-tail) (1CDX4_242 in Extended Data Fig. 4b) as though organized by the mitotic spindle in telophase (see, also 1CDX_74 and 1CDX_453 in Extended Data Fig. 4b for partial or disordered chromosomal alignments). Comparison of the frequency of contacts in the mitotic range to the frequency of contacts at shorter distances (short-range contacts < 2 Mb; Fig. 1e) shows a remarkable circular pattern evocative of a gradual chromosome remodelling process as cells enter and exit mitosis. These observations suggest chromosome conformation may be used to phase single cells at various stages of the cell cycle around mitosis.

Early- and late-replicating domains are defined by their copy number dynamics during S phase, and show distinct enrichment for active and inactive epigenetic marks respectively¹⁴. Indeed, we found normalized TAD coverage across the single-cell dataset reflects a strong correlation between domains previously annotated as early- or late-replicating in mouse ESCs¹⁵ (Extended Data Fig. 4c-d). We computed a repli-score for each cell based on the copy-number ratio of early-replicating regions to total coverage. G1 cells would be expected to have low repli-scores, whereas early S cells would be associated with increasing scores, peaking in mid-S before declining again in late S phase and returning to a low level in G2/M. Combining repli-score with the circular pattern exhibited by the frequencies of short-range and mitotic contacts yields exactly that expected trajectory (Fig. 1e). Together, the mitotic signatures and time of replication analysis define two major anchors to support phasing of the entire single-cell Hi-C dataset along the cell cycle.

We studied the distributions of contact distance scales in single cells using non-linear dimensionality reduction (NLDR) (Extended Data Fig. 4e) and clustering (Extended Data Fig. 5). Both approaches indicate the data can be globally represented by a cyclic, gradual process of coordinated change that links the highly distinct mitotic and replicating states characterized above. To capture this variation precisely we used a supervised approach that groups cells to phases and then orders cells within phases. We first identified cells approaching mitosis (pre-M regime, orange) and coming out of mitosis (post-M regime, red) by thresholding short-range and mitotic contact frequencies (based on Fig 1e, see Methods). Cells in both regimes can be ordered by the ascending and descending frequency of mitotic contacts for the pre-M and post-M regimes, respectively. Next, we observed that cells with less than 63% short-range contacts do not include the replication group and likely represent a G1 regime (blue). The mean distance of long-range contacts in these cells is correlated to the overall frequency of short-range contacts, allowing ordering of cells within this group (Extended Data Fig. 6a). Cells with over 63% short-range contacts show a continuum of repli-score and short-range contact enrichments (Extended Data Fig. 6b). These were modelled using two additional regimes, one ordered based on a gradual increase in repli-score and short-range contacts (the early-S regime, green), and one ordered based on a gradual decrease in repli-score and further increase in short-range contacts (the late-S/G2 regime, purple). Using this phasing strategy, we observed a smooth, cyclical transition of the distribution of genomic contact distances (Fig. 1f), which is insensitive to batch effects (Extended Data Fig. 6c-d) and robust to subsampling of the data or analysis of specific chromosomes (Extended Data Fig. 6e-f). In summary, we developed an in-silico cell-cycle phasing approach that organises single cells along a smooth trajectory of highly varied chromosome conformations (Fig 1g, Extended Data Fig. 5) leading from the mitotic state through replication and back into mitosis.

To confirm our Hi-C based cell-cycle phasing, we analysed 1169 cells sorted based on Geminin and DNA content using fluorescence-activated cell sorting (FACS) and found that the ordering of cells based on their contacts is remarkably consistent with their stage as determined by FACS (Fig. 2a-b). We further hypothesise that the refined ordering we impose within each coarse-grained stage is likely to represent cell-cycle progression at high resolution. This was initially demonstrated by correlating Hoechst/Geminin FACS indices to the inferred phasing rank of 342 single-cell Hi-C profiles (Fig. 2c, Extended Data Fig. 6g), deriving correlation even within classical sorting gates. To further support this idea we analysed the distribution of key statistics for single cells ordered according to their inferred cell-cycle stage. For reference, we projected the repli-score on the cell-cycle progression (Fig. 2d), validating that it was low during G1 and G2, showed consistent gradual increase during early S phase, and decreased during late S phase. The total number of contacts (Fig. 2e), which unlike the repli-score is a-priori independent of the inferred phasing, was shown to increase gradually during S-phase, peaking at G2 at almost a 2-fold higher level than G1, before dropping sharply after mitosis. The fraction of trans-chromosomal contacts exhibits an unforeseen cell-cycle dynamic, with a smooth curve from its lowest point of around 3% in G2 and mitosis, to ∼8% in early S phase (Fig. 2f). Next, we quantified the trans-alignment of chromosomes, which is a strong characteristic of mitotic cells. Remarkably, despite being a-priori unrelated to the inferred G1 ordering scheme, we found the trans-alignment score declined smoothly from its peak at mitosis to its nadir at the end of the G1 phase, immediately before the start of S phase (Fig. 2g; see Extended Data Fig. 6h for examples). Collectively, these results validate our phasing strategy and open the way to analyse cell-cycle dynamics of Hi-C architectural features at high resolution.

Cell-cycle chromosomal dynamics

We identified 3434 TAD borders in the pooled Hi-C map using analysis of contact insulation signatures⁶ (Extended Data Fig. 3). We created sub-pools of G1, early-S and late-S/G2 cells in silico to compare insulation patterns between the three major cell-cycle phases. We found that the locations of TAD borders are generally unchanged in these pools, but that the average intensity of the borders is dynamic (Extended Data Fig. 7a-c). We calculated the mean insulation at borders for each cell (Fig. 3a) and found that the intensity of insulation varies markedly over the cell cycle. Insulation is not detectable in mitotic maps (Extended Data Fig. 7d-e), and emerges rapidly upon exit from mitosis in early G1. Insulation increases to its maximum during G1, but declines when replication begins, plateauing at mid S-phase at its lowest interphase level and remaining fixed through G2 until mitosis. We note that since insulation is defined by the depletion of short-range contacts crossing TAD borders, it may be affected by changes in the frequency of contacts in the 25 kb-2 Mb distance range. However, this overall frequency only partially predicts single-cell insulation (Extended Data Fig. 7f).

Next, we studied the dynamics of compartmentalisation of the genome into active (A) and non-active (B) groups of TADs by analysing the cis-chromosomal compartment linkage. For each cell we estimated the depletion of A-B contacts given the observed A-A and B-B counts (Methods, Fig. 3b). Similar to insulation, compartmentalisation is weakest during mitosis, as expected from the condensed uniform structure at that stage¹¹. However, in contrast to TAD insulation, compartmentalisation is weak during G1 phase and increases upon entering S phase, with a further increase in late S phase, reaching a maximum in G2, just before a rapid and complete loss approaching mitosis. Compartmentalisation calculated from trans-chromosomal contacts shows exactly the same trend (Extended Data Fig. 7g). The distinct cell-cycle dynamics of TAD insulation and compartment structure suggest that independent mechanisms contribute to these phenomena.

To study the cell-cycle dynamics of insulation at individual TAD borders and compare it to the global trend above, we grouped single cells into 20 time-slots according to their inferred cell-cycle phase, and pooled data from each time-slot to compute an insulation time series for each border. We then clustered the borders based on their dynamics in these time-slots, and examined the mean insulation in each cluster per cell, along the inferred phasing. Borders demonstrate similar qualitative dynamics in all clusters. However, each cluster varies in timing and intensity of deterioration in insulation during S-phase (Fig. 3c). Cluster C4 includes borders showing the strongest G1 insulation and C2 borders show the earliest loss of insulation in S. Interestingly, borders in both clusters tend to be located between early-replicating TADs (Fig. 3d). Borders in C3 show a gradual increase in insulation well into S-phase, followed by abrupt loss of insulation in mid-S. This cluster is enriched for borders separating early-replicating from mid/late-replicating TADs (Fig. 3c). C1, which has the weakest G1 insulation, shows gradual and delayed loss of insulation, and is enriched for borders between late-replicating domains. These data suggest that loss of insulation in S-phase is temporally coupled to the replication process. Specifically, replication through a border or in both adjacent TADs is linked to loss of insulation (C1, C2, C4), whereas replication of only one adjacent TAD is linked to a temporary increase in insulation (C3). In all cases, insulation remains low post-replication, not re-emerging until the next cell cycle.

We identified pairs of looping loci by analysis of pooled matrices and shuffled pooled matrices (Extended Data Fig. 7h). We focused on 2036 top-scoring loops that were associated with convergent CTCF motifs and evidence for CTCF/Cohesin binding (Methods). We then grouped loops according to the replication time of their anchors (early-early vs. late-late) and computed for each group and each single cell the total number of looping contacts in 20×20 kb bins around the loop anchors, normalized by the number of contacts in the surrounding 60×60 kb bin (Fig. 3e). We found that loop enrichment increases after mitosis, is generally stable through G1, but then decreases slightly in early S phase for early-replicating loops and in late S phase for later-replicating loops. Spatial enrichment analysis shows conserved localized loop enrichment (Fig. 3f). We do not detect loops that are specific to a certain cell-cycle phase (Extended Data Fig. 7i), suggesting that even when decreased in average intensity during replication, loops are significantly observed in the data. However, the assay's resolution is not sufficient to allow detection of individual loop turnover at single-cell resolution, leaving unresolved regimes for rapid loop dynamics. Interestingly, loops are still apparent in cells approaching mitosis despite the fast chromosomal condensation – unlike cells after mitosis, which do not have any contact enrichment around loops (Fig. 3f). The notable lack of loop contact enrichment in mitotic maps does not necessarily imply loops are disassembled at this stage, since mitotic condensation may increase the local contact frequency around loop anchors and decrease Hi-C sensitivity.

Given the relative stability of looping contacts that demarcate TAD structures and the dynamic cell-cycle rise and fall in TAD insulation, we next asked if chromatin structure within TADs changes during the cell cycle, or if TADs constitute stable chromatin blocks that are re-organised dynamically. We computed average intra-TAD contact distances and used these as a measure of TAD condensation per cell, reasoning that condensed TADs will be characterised by higher mean contact distance than de-condensed TADs^3,16. We compared condensation for groups of TADs while controlling for size and time of replication (Fig. 3g). In G1, early-replicating TADs are de-condensed compared to mid- or late-replicating TADs, consistent with open chromatin in active regions. Early-replicating TADs are further de-condensed in early S phase, before showing a remarkable increase in condensation at mid-S, likely after their replication is complete. Following that, and during late S and G2, early-replicating TADs increase their mean condensation to levels comparable to those of late-replicating TADs in G1. Condensation of mid-early replicating TADs follows a similar trend, with some delay and damped amplitude. Mid-late and late-replicating TADs, on the other hand, show constitutively high condensation in G1 and throughout S, with further marked increases observed for mid-late TADs in G2. We note that the potential involvement of trans sister chromatid contacts in driving the observed increase in mean contact distance cannot be determined unambiguously. However, the delay between replication of early domains at the start of S phase and their increased compaction in mid/late S suggests sister-chromatid contacts cannot fully account for the dynamics observed.

We next generated 1247 QC positive single-cell Hi-C maps from haploid mouse ESCs, to allow systematic analysis of conformation in single-copy chromosomes. We confirmed haploid ESC maps show comparable QC metrics to diploid cells, with less than half of the contacts per cell and some specific patterns of chromosome loss that we filtered in-silico (Extended Data Fig. 8a-c). We derived a pooled ensemble haploid Hi-C map and computed TAD structures, insulation and A/B compartments (compared to diploid in Extended Data Fig. 3f). We then used the same approach for cell-cycle phasing that was developed for diploid cells (Fig. 4a). Remarkably, we observe a similar cyclical trend, and analogous dynamics of insulation, loop enrichment and TAD condensation (Extended Data Fig. 8d-g). The universality of these processes shows that haploid ESCs are comparable to normal diploid cells, and a good model for studying chromosome conformation and dynamics, while providing the opportunity for unambiguous 3D modelling of single-copy chromosomes for G1 cells.

A spectrum of epigenetic compartments

We defined A-scores for each TAD as the fraction of trans-chromosomal contacts made with the A-compartment, and found a spectrum that generalizes the original discrete partitioning of the genome into two compartments (Fig. 4b). A-scores correlate with time of replication; high A-score TADs are early-replicating and low A-score TADs are late-replicating (Fig. 4c). These scores are conserved between diploid and haploid cells and growth conditions (Extended Data Fig. 9a), suggesting that a continuum of compartment affinities is a robust property of chromosomal domains. Such a continuum may be created if individual TADs switch compartments in individual cells, thereby creating a gradient of A-compartment affinities. Alternatively, it is possible that the A-score gradient is conserved at the single-cell level, representing a structural continuum, rather than a two-compartment structure that is heterogeneous across a cell population. To explore these questions at single-cell resolution we needed sufficient TAD coverage in each cell, so we computed a single-cell A-compartment association score for each TAD, using the mean A-score of its long-range (> 2 Mb) cis-contacting TADs. Cis A-association scores and trans A-scores are highly consistent on average (Fig. 4d), but the cis A-association score is defined by more contacts, while being constrained by the chromosomal larger-scale domain architecture. The standard deviation of Aassociation scores for individual TADs across the entire haploid single-cell cohort was less than 0.06 (Fig. 4e), suggesting limited variability of A-score between cells. Further analysis (Fig. 4f) showed that A-score variation decreases in early and late S-phase, consistent with the overall increase in A/B compartment strength above. Both sorting TADs by their mean A-association scores (Fig. 4g) and grouping these distributions by unsupervised clustering (Extended Data Fig. 9b) showed that TADs vary quantitatively in their A-association, but no significant group of TADs shows a bimodal compartment association distribution. The same analysis indicated that compartment association polarises in late S/G2, with many TADs shifting from broadly spread association in G1 or early S, to preferential contacts with TADs in a narrower range of A-scores following replication.

To further explore this effect we discretised the A-association score continuum, deriving 20 groups of TADs that corresponded to putative pseudo-compartments. We tested how such pseudo-compartments contact one another trans- and cis-chromosomally, and observed significant cis-chromosomal interaction between TADs within the same and neighbouring pseudo-compartments (Fig. 4h and Extended Data Fig. 9e, top). In trans we observed a more binary contact structure in G1 cells, transitioning toward a gradient of pseudo-compartment preferences following replication (Fig. 4i). A-association analysis of diploid cells shows similar behaviour (Extended Data Fig. 9c,d,f). Analysis of RNA-seq, histone marks¹⁵ and Lamin association¹⁷ suggested pseudo-compartments are characterized by a gradient of epigenetic signatures (Extended Data Fig. 9g-h). Overall, the ordering of TADs according to A-association score is significantly correlated with a linear model weighting H3K4me1 and H3K4me3 scores (r=0.79, p≪0.001). Together, these data suggest pseudo-compartments capture large-scale trans- and cis-chromosomal conformation at the population and single-cell level, generalising the binary A/B compartment scheme.

3D whole genome modelling and analysis

Using G1 haploid cell data we developed a structure sampling approach for inferring whole genome 3D conformation. We applied stringent contact filtering to eliminate contacts that are more likely to be technical artifacts (Methods), removing a median of 0.89% of the contacts from each cell (Extended Data Fig. 10a). We applied a simulated annealing protocol¹² to sample over 30,000 whole genome 3D models for 190 single-cell datasets, using contacts as distance restraints at 0.1, 0.5 and 1Mb resolutions. The resulting 3D models satisfied nearly all observed cis- and trans-chromosomal contacts, with a median of <0.012% violated constraints per cell (Extended Data Fig. 10b-c). Subsampling 10,000 contacts per cell produced consistent models (Extended Data Fig. 10d-e), while cross-validation tests showed that structural modelling is sufficiently robust to re-discover trans-chromosomal distances between a pair of contacting chromosomes, even when all contacts between those chromosomes are discarded (Extended Data Fig. 10f).

To study changes in chromosomal structure as cells emerge from mitosis into G1, we analysed over 20,000 models from 126 cells with > 10,000 contacts per cell and < 0.5% constraint violations. Analysis of large-scale cis-chromosomal shapes shows that individual chromosomes rapidly de-condense from a rod-like mitotic conformation to a more spherical structure as G1 progresses (Fig. 5a-b, Extended Data Fig. 10g-h). We monitored the expansion of chromosomes by measuring the distance between adjacent 500 kb chromatin segments in A or B compartments in each cell (Fig. 5c). A and B compartments are similarly compact upon mitotic exit (P > 0.05, two-sided MW test), but as cells progress through G1, A compartments expands more rapidly and to a larger extent than B compartments (P < 0.005, onesided MW test), consistent with decondensed euchromatin and compact heterochromatin. Remarkably, the twenty high-resolution pseudo-compartments form a gradient that generalises this observation, with the strongest B-like pseudocompartments most compact and the strongest A-like pseudo-compartments the most expanded (Fig. 5d).

We then computed the mean nuclear radial positioning of A and B compartments as cells exit mitosis. We observed B compartments occupying more peripheral locations and A compartments tending toward the nuclear interior for cells 9-16 and onward (P < 0.001, one-sided MW test, Fig. 5e). Pseudo-compartments generalise this trend, showing progressively more differentiated radial positioning through G1 phase (Fig. 5f). The structural models also confirm and refine the preferential contacts among pseudo-compartments (as shown in Fig. 4h-i), and show that long-range cis-contacts are re-established very early in G1 phase, while trans compartmentalisation is not re-established until late G1 (Fig. 5g-h), consistent with 4C observations¹⁸. Finally, distributions of chromosomal conformations also suggest preferential contacts of Nucleolar Organiser Regions¹⁹ (NOR, P < 0.01 one-sided MW test, Extended Data Fig. 10i).

Methods

• 1.Paweletz N. Walther Flemming: pioneer of mitosis research. Nat Rev Mol Cell Biol. 2001;2:72–75. doi: 10.1038/35048077. [DOI] [PubMed] [Google Scholar]
• 2.Lieberman-Aiden E, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009;326:289–93. doi: 10.1126/science.1181369. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 3.Sexton T, et al. Three-dimensional folding and functional organization principles of the Drosophila genome. Cell. 2012;148:458–72. doi: 10.1016/j.cell.2012.01.010. [DOI] [PubMed] [Google Scholar]
• 4.Nora EP, et al. Spatial partitioning of the regulatory landscape of the Xinactivation centre. Nature. 2012;485:381–5. doi: 10.1038/nature11049. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 5.Dixon JR, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485:376–80. doi: 10.1038/nature11082. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 6.Sofueva S, et al. Cohesin-mediated interactions organize chromosomal domain architecture. EMBO J. 2013;32:3119–29. doi: 10.1038/emboj.2013.237. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 7.Rao SSPSP, et al. A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping. Cell. 2014;159:1665–1680. doi: 10.1016/j.cell.2014.11.021. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 8.Zuin J, et al. Cohesin and CTCF differentially affect chromatin architecture and gene expression in human cells. Proc Natl Acad Sci U S A. 2014;111:996–1001. doi: 10.1073/pnas.1317788111. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 9.de Wit E, et al. CTCF Binding Polarity Determines Chromatin Looping. Mol Cell. 2015;60:676–684. doi: 10.1016/j.molcel.2015.09.023. [DOI] [PubMed] [Google Scholar]
• 10.Flavahan WA, et al. Insulator dysfunction and oncogene activation in IDH mutant gliomas. Nature. 2016;529:110–114. doi: 10.1038/nature16490. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 11.Naumova N, et al. Organization of the mitotic chromosome. Science. 2013;342:948–53. doi: 10.1126/science.1236083. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 12.Nagano T, et al. Single-cell Hi-C reveals cell-to-cell variability in chromosome structure. Nature. 2013;502:59–64. doi: 10.1038/nature12593. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 13.Nagano T, et al. Comparison of Hi-C results using in-solution versus in-nucleus ligation. Genome Biology. 2015;16:175. doi: 10.1186/s13059-015-0753-7. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 14.Pope BD, et al. Topologically associating domains are stable units of replication-timing regulation. Nature. 2014;515:402–405. doi: 10.1038/nature13986. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 15.Stamatoyannopoulos JA, et al. An encyclopedia of mouse DNA elements (Mouse ENCODE) Genome Biol. 2012;13:418. doi: 10.1186/gb-2012-13-8-418. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 16.Boettiger AN, et al. Super-resolution imaging reveals distinct chromatin folding for different epigenetic states - nature16496. pdf Nature. 2016;529:418–422. doi: 10.1038/nature16496. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 17.Peric-Hupkes D, et al. Molecular maps of the reorganization of genome-nuclear lamina interactions during differentiation. Mol Cell. 2010;38:603–13. doi: 10.1016/j.molcel.2010.03.016. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 18.Dileep V, et al. Topologically-associating domains and their long-range contacts are established during early G1 coincident with the establishment of the replication timing program Authors. Genome Res. 2015;25:1104–1113. doi: 10.1101/gr.183699.114. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 19.Britton-Davidian J, Cazaux B, Catalan J. Chromosomal dynamics of nucleolar organizer regions (NORs) in the house mouse: micro-evolutionary insights. Heredity. 2012;108:68–74. doi: 10.1038/hdy.2011.105. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 20.Ramani V, et al. Massively multiplex single-cell Hi-C. Nat Methods. 2017;14:263–266. doi: 10.1038/nmeth.4155. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 21.Stevens TJ, et al. 3D structures of individual mammalian genomes studied by single-cell Hi-C. Nature. 2017;544:59–64. doi: 10.1038/nature21429. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 22.Flyamer IM, et al. Single-nucleus Hi-C reveals unique chromatin reorganization at oocyte-to-zygote transition. Nature. 2017;544:110–114. doi: 10.1038/nature21711. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials