Maternal-to-zygotic transition (MZT) is essential for the formation of a new individual, but is still poorly understood despite recent progress in analysis of gene expression and DNA methylation in early embryogenesis1,2,3,4,5,6,7,8,9. Dynamic histone modifications may have important roles in MZT10,11,12,13, but direct measurements of chromatin states have been hindered by technical difficulties in profiling histone modifications from small quantities of cells. Recent improvements allow for 500 cell-equivalents of chromatin per reaction, but require 10,000 cells for initial steps14 or require a highly specialized microfluidics device that is not readily available15. We developed a micro-scale chromatin immunoprecipitation and sequencing (μChIP–seq) method, which we used to profile genome-wide histone H3 lysine methylation (H3K4me3) and acetylation (H3K27ac) in mouse immature and metaphase II oocytes and in 2-cell and 8-cell embryos. Notably, we show that ~22% of the oocyte genome is associated with broad H3K4me3 domains that are anti-correlated with DNA methylation. The H3K4me3 signal becomes confined to transcriptional-start-site regions in 2-cell embryos, concomitant with the onset of major zygotic genome activation. Active removal of broad H3K4me3 domains by the lysine demethylases KDM5A and KDM5B is required for normal zygotic genome activation and is essential for early embryo development. Our results provide insight into the onset of the developmental program in mouse embryos and demonstrate a role for broad H3K4me3 domains in MZT.
This is a preview of subscription content
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Get time limited or full article access on ReadCube.
All prices are NET prices.
Smith, Z. D. et al. A unique regulatory phase of DNA methylation in the early mammalian embryo. Nature 484, 339–344 (2012)
Wang, L. et al. Programming and inheritance of parental DNA methylomes in mammals. Cell 157, 979–991 (2014)
Shen, L. et al. Tet3 and DNA replication mediate demethylation of both the maternal and paternal genomes in mouse zygotes. Cell Stem Cell 15, 459–470 (2014)
Guo, F. et al. Active and passive demethylation of male and female pronuclear DNA in the mammalian zygote. Cell Stem Cell 15, 447–458 (2014)
Park, S.-J. et al. Inferring the choreography of parental genomes during fertilization from ultralarge-scale whole-transcriptome analysis. Genes Dev . 27, 2736–2748 (2013)
Smith, Z. D. et al. DNA methylation dynamics of the human preimplantation embryo. Nature 511, 611–615 (2014)
Deng, Q., Ramsköld, D., Reinius, B. & Sandberg, R. Single-cell RNA-seq reveals dynamic, random monoallelic gene expression in mammalian cells. Science 343, 193–196 (2014)
Xue, Z. et al. Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing. Nature 500, 593–597 (2013)
Guo, H. et al. The DNA methylation landscape of human early embryos. Nature 511, 606–610 (2014)
Aoshima, K., Inoue, E., Sawa, H. & Okada, Y. Paternal H3K4 methylation is required for minor zygotic gene activation and early mouse embryonic development. EMBO Rep . 16, 803–812 (2015)
Morgan, H. D., Santos, F., Green, K., Dean, W. & Reik, W. Epigenetic reprogramming in mammals. Hum. Mol. Genet. 14, R47–R58 (2005)
Santos, F., Peters, A. H., Otte, A. P., Reik, W. & Dean, W. Dynamic chromatin modifications characterise the first cell cycle in mouse embryos. Dev. Biol. 280, 225–236 (2005)
Adenot, P. G., Mercier, Y., Renard, J. P. & Thompson, E. M. Differential H4 acetylation of paternal and maternal chromatin precedes DNA replication and differential transcriptional activity in pronuclei of 1-cell mouse embryos. Development 124, 4615–4625 (1997)
Lara-Astiaso, D. et al. Immunogenetics. Chromatin state dynamics during blood formation. Science 345, 943–949 (2014)
Shen, J. et al. H3K4me3 epigenomic landscape derived from ChIP–seq of 1,000 mouse early embryonic cells. Cell Res . 25, 143–147 (2015)
Kues, W. A. et al. Genome-wide expression profiling reveals distinct clusters of transcriptional regulation during bovine preimplantation development in vivo. Proc. Natl Acad. Sci. USA 105, 19768–19773 (2008)
Ruthenburg, A. J., Allis, C. D. & Wysocka, J. Methylation of lysine 4 on histone H3: intricacy of writing and reading a single epigenetic mark. Mol. Cell 25, 15–30 (2007)
Bannister, A. J. & Kouzarides, T. Regulation of chromatin by histone modifications. Cell Res . 21, 381–395 (2011)
Smallwood, S. A. et al. Single-cell genome-wide bisulfite sequencing for assessing epigenetic heterogeneity. Nat. Methods 11, 817–820 (2014)
Ooi, S. K. T. et al. DNMT3L connects unmethylated lysine 4 of histone H3 to de novo methylation of DNA. Nature 448, 714–717 (2007)
Otani, J. et al. Structural basis for recognition of H3K4 methylation status by the DNA methyltransferase 3A ATRX-DNMT3-DNMT3L domain. EMBO Rep . 10, 1235–1241 (2009)
Shirane, K. et al. Mouse oocyte methylomes at base resolution reveal genome-wide accumulation of non-CpG methylation and role of DNA methyltransferases. PLoS Genet . 9, e1003439 (2013)
Erkek, S. et al. Molecular determinants of nucleosome retention at CpG-rich sequences in mouse spermatozoa. Nat. Struct. Mol. Biol. 20, 868–875 (2013)
Creyghton, M. P. et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl Acad. Sci. USA 107, 21931–21936 (2010)
Nord, A. S. et al. Rapid and pervasive changes in genome-wide enhancer usage during mammalian development. Cell 155, 1521–1531 (2013)
Shen, Y. et al. A map of the cis-regulatory sequences in the mouse genome. Nature 488, 116–120 (2012)
Visel, A. et al. ChIP–seq accurately predicts tissue-specific activity of enhancers. Nature 457, 854–858 (2009)
Andreu-Vieyra, C. V. et al. MLL2 is required in oocytes for bulk histone 3 lysine 4 trimethylation and transcriptional silencing. PLoS Biol . 8, e1000453 (2010)
Leung, D. et al. Integrative analysis of haplotype-resolved epigenomes across human tissues. Nature 518, 350–354 (2015)
Benayoun, B. A. et al. H3K4me3 breadth is linked to cell identity and transcriptional consistency. Cell 158, 673–688 (2014)
Chiang, T. & Lampson, M. A. Counting chromosomes in intact eggs. Methods Mol. Biol. 957, 249–253 (2013)
Dahl, J. A. & Collas, P. A rapid micro chromatin immunoprecipitation assay (microChIP). Nat. Protocols 3, 1032–1045 (2008)
Dahl, J. A. & Klungland, A. Micro chromatin immunoprecipitation (μChIP) from early mammalian embryos. Methods Mol. Biol. 1222, 227–245 (2015)
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009)
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014)
Krueger, F. & Andrews, S. R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572 (2011)
Trapnell, C. et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515 (2010)
Lerdrup, M., Johansen, J. V., Agrawal-Singh, S. & Hansen, K. An interactive environment for agile analysis and visualization of ChIP–sequencing data. Nat. Struct. Mol. Biol. 23, 349–357 (2016)
McLean, C. Y. et al. GREAT improves functional interpretation of cis-regulatory regions. Nat. Biotechnol. 28, 495–501 (2010)
Schones, D. E., Smith, A. D. & Zhang, M. Q. Statistical significance of cis-regulatory modules. BMC Bioinformatics 8, http://dx.doi.org/10.1186/1471-2105-8-19 (2007)
Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589 (2010)
We thank the members of the Ren laboratory for support and critical suggestions throughout the course of this work. We are thankful to L. Ellevog, G. Flor Lien, T. Catterall, S. Khan, I. Johnson, the Norwegian Transgenic Center, the Animal Facility OUS and Genoway for help with embryo collection, animal care and morpholino injections. We are grateful to A. Local for the gift of recombinant histone octameres and to R. Ahmed, M. Indahl and E. Skarpen for assistance with staining and imaging of embryos. We are also thankful to K. Shirane (Kyushu University) for sharing data. I. Jung would like to give special thanks to Y. Lee. This work was funded by the Ludwig Institute for Cancer Research, U54HG006997 (to B.R.), American Heart Association Postdoctoral Fellowship (to I. Jung), the Oslo University Hospital Early Career Award (to J.A.D.), the Norwegian Cancer Society (to A.K., J.A.D.), the Anders Jahre Foundation (to J.A.D.) and the Norwegian Research council (to A.K.).
Inven2 (on behalf of J.A.D. and A.K.) have registered a patent application entitled ‘ChIP–seq assays’.
Reviewer Information Nature thanks R. Schultz and the other anonymous reviewer(s) for their contribution to the peer review of this work.
Extended data figures and tables
a, Genome browser snapshot for H3K4me3 μChIP–seq with different numbers of human NCCIT cells. The number of cells used is shown at the beginning of each genome browser track. Input ChIP–seq and RefSeq gene annotation are shown in the two last rows. b, The percentage of overlapping H3K4me3 peaks between the top-10,000 H3K4me3 peaks detected from multiple μChIP–seq experiments with different numbers of human NCCIT cells. c, A genome browser snapshot for H3K27ac and H3K4me3 μChIP–seq in mouse E14 ES cells. ENCODE ChIP–seq is also shown. Cell number in μChIP–seq experiments ranged from 100,000 to 500 cells. d, The percentage of overlapping H3K4me3 peaks between top 10,000 H3K4me3 peaks when comparing multiple μChIP–seq experiments performed with different numbers of mouse ES cells. Applying μChIP–seq to 1,000 cells allowed detection of 90% of ENCODE peaks and 94% of peaks from large scale ChIP–seq from the same ES-cell culture batch. e, The percentage of overlapping H3K27ac peaks between multiple experiments using different numbers of mouse ES cells. We considered the top 10,000 H3K27ac peaks in each sample. The sample named ENCODE are H3K27ac peaks from the ENCODE dataset. f, A scatter plot comparing RPKM values obtained from 1,000-cell μChIP–seq and ENCODE ChIP–seq results using an H3K4me3 antibody in mouse ES cells. Each dot represents an H3K4me3 peak identified using the H3K4me3 μChIP–seq results (n = 26,397). Pearson correlation coefficient is also shown. g, Receiver operating characteristic (ROC) curve plots for H3K4me3 peaks identified from multiple μChIP–seq experiments. True H3K4me3 peaks were defined using ENCODE H3K4me3 ChIP–seq result and false H3K4me3 peaks were created by randomly selecting genomic regions. For the 1,000-cell μChIP–seq result we also calculated ROC scores for two additional subsets of true H3K4me3 peaks; the top 30,000 and 40,000 H3K4me3 ENCODE peaks. h, Box plots of GC content comparing 500-cell H3K4me3 ChIP–seq detected and undetected ENCODE H3K4me3 peaks (Kolmogorov–Smirnov test P = 0.81). i, Box plots of ENCODE H3K4me3 signal comparing 500-cell H3K4me3 ChIP–seq detected and undetected H3K4me3 peaks (Kolmogorov–Smirnov test P < 2.2−16). Whiskers correspond to the highest and lowest points within the 1.5× interquartile range (h, i).
a, Genome browser snapshots of H3K4me3 and H3K27ac μChIP–seq results in MII oocytes, 2- and 8-cell embryos and ES cells. For each histone-modification mark, we generated two biological replicates. b, Genome browser snapshots of H3K4me3 μChIP–seq results from P12 and P15 oocytes. We generated two biological replicates. c, d, Bar plots show the Pearson correlation coefficient between two biological replicates for H3K4me3 (c) and H3K27ac (d). Pearson correlation coefficients were calculated using RPKM values for 1-kb-binned regions of the whole genome.
a, Box plots for H3K4me3 RPKM values at promoter regions after grouping promoters into 10 equally sized groups according to gene expression levels. b, Box plots for gene expression levels for each group are shown. Whiskers correspond to the highest and lowest points within the 1.5× interquartile range (a, b). c, Fraction of the mouse genome covered by H3K4me3 peaks at different developmental stages and in various adult tissues26. H3K4me3 peaks were defined using default MACS2 peak call settings except for oocyte. Oocyte H3K4me3 peaks were defined using MACS2 broad peak call settings. d, Scatter plot of H3K4me3 signals in 1kb bins between oocyte and 2-cell-stage embryos. Genomic regions marked with a strong H3K4m3 signal in the oocyte show reduced signal in 2-cell-stage embryos. e, Fraction of the mouse genome covered by H3K27ac peaks using default MACS2 peak call settings for different developmental stages.
a, Panel from Fig. 1d with consistently adjusted brightness and contrast to allow visual observation of H3K4me3 signal following rapid loss at the late 2-cell stage. The signal of the strongly staining stages appears saturated due to this. b, c, Quantification of H3K4me3 and DNA staining by integrated H3K4me3-specific and DAPI-specific fluorescence intensity, respectively, over pronuclei and nuclei was determined using an epifluorescence microscope and the ImageJ software. The nuclear DNA content of ‘early’ and ‘late’ 2-cell embryos reveals that DNA synthesis is almost completed at the ‘early’ time point. Therefore suggesting an active mechanism of H3K4me3 removal rather than passive dilution over DNA replication (n = 18, 18, 20, 41, 18, 74 for respective time points; error bars show s.e.m.).
a, Genome browser snapshot of H3K4me3 μChIP–seq and RNA-seq signals. The H3K4me3 μChIP–seq signal in oocytes shows enrichment over broader regions compared to sperm, 2- and 8-cell embryos and mouse ES cells (mESCs). The broad H3K4me3 domain call is shown as a green bar below the oocyte H3K4me3 track. The TSS of Foxa1 (ZGA gene) is located in a broad H3K4me3 domain. The Foxa1 gene becomes expressed in 2- and 8-cell embryos. b, Decreasing fraction of H3K4me3 domains in oocytes when increasing the distance threshold to merge adjacent H3K4me3 peaks. At 5 kb (asterisk), this fraction becomes stable, which is used as the threshold to call broad H3K4me3 domains in oocytes. c, Histogram showing fraction of H3K4me3 domains of various breadth for TSS-containing (n = 15,608) and non-TSS-containing (n = 47,934) domains. TSS-containing H3K4me3 domains tend to be broader (median, 12 kb) than non-TSS-containing domains (median, 3.7 kb) (P value < 2.2−16, Kolmogorov–Smirnov test). d, An anti-correlated relationship between genomic regions covered by H3K4me3 (x axis) and median-input-normalized H3K4me3 ChIP–seq RPKM values at 5,000 top-ranked promoters (y axis). The samples more widely covered by H3K4me3 tend to show relatively lower H3K4me3 RPKM values. e, Line plots are shown for median input normalized H3K4me3 RPKM values (y axis) with corresponding number of top-ranked promoters (x axis). When the genomic coverage of H3K4me3 is similar between samples (liver, cerebellum, and heart) the median of input normalized H3K4me3 RPKM values are highly consistent between samples at a given number of top-ranked promoters. f, Robustness of H3K4me3 RPKM adjustment factor. y axis indicates H3K4me3 RPKM adjustment factor for each cell type. Mouse ES cell data are used as a reference (that is, adjustment factor = 1). x axis indicates number of top-ranked promoters that was used to calculate H3K4me3 RPKM adjustment factor. The adjustment factors are very robust regardless of number of top-ranked promoters considered. g–i, Genome browser tracks for several loci show H3K4me3 μChIP–seq for oocytes, 2-cell and 8-cell embryos and mouse ES cells. y axis indicates auto scale for H3K4me3 raw reads (g), fixed scale for H3K4me3 RPKM (h), and fixed scale for H3K4me3 signal after applying the adjustment factor (i).
Extended Data Figure 6 Establishment of broad H3K4me3 domains in oocytes and removal in the early embryo.
a, b, Heat maps of H3K4me3 raw read counts normalized by total read numbers in indicated stages for non-TSS (a) and TSS-containing H3K4me3 domains (b) ordered according to increasing size. Horizontal positions reflect spatial distribution relative to the domain centre, and densities correspond to H3K4me3 signal intensity. c, Genome browser snapshot of H3K4me3 μChIP–seq from P12, P15 and MII oocytes. d, Superimposed tracks of H3K4me3 signal in P12, P15 and MII oocytes at the boundaries of non-TSS-containing H3K4me3 domains. The y axes show H3K4me3 signal, colouring reflects the total count of domains with a given signal intensity, and horizontal positions reflect spatial distribution relative to the domain border.
a, Percentage of broad domains displaying hypomethylated (<0.25), hypermethylated (>0.75) and intermediate (>0.25 and <0.75) levels within TSS-containing domains, non-TSS-containing domains, and outside broad domains. b, Histograms of average DNA methylation levels within TSS-containing (n = 15,135) and non-TSS-containing domains (n = 35,567), and outside of domains (n = 61,674). Domains containing at least 5 covered CpGs were included. c, Bar plots displaying the fraction of ZGA genes5 and genes not expressed that reside in broad H3K4me3 domains. (***P < 10−3, Fisher’s exact test). d, e, Heat maps of H3K4me3 and DNA methylation signal for indicated oocyte stages for non-TSS-containing H3K4me3 domains ordered according to increasing sizes. Horizontal positions reflect spatial distribution relative to the domain centre, and densities correspond to H3K4me3 signal intensity (d) or CpG methylation level (e). NGO, non-growing oocyte; GVO, germinal vesicle oocytes. f, Bar plots displaying the number of genes expressed from the maternal allele for genes that reside in broad H3K4me3 domains. The maternally expressed genes are significantly enriched within broad H3K4me3 domains compared to what would be expected by chance (***P < 10−3, Fisher’s exact test). g, k-means clustered heat maps of H3K4me3 and DNA methylation signal for P12, P15 and MII oocytes at H3K4me3 domain boundaries. Boundaries of domains larger than 5 kbp were clustered according to H3K4me3 and DNA methylation signals at the 2-kbp regions flanking the boundaries, and densities correspond to H3K4me3 signal or CpG methylation level. C1–C10, Cluster 1–10. h, i, Box plots for distance from domain boundary to the nearest CpG island (h) and bar plots for fraction of TSS-containing broad domains (i) in each cluster. Whiskers correspond to the highest and lowest points within the 1.5× interquartile range. j–l, Line plots are shown for average DNA methylation profiles at both 5′ and 3′ end of ZGA gene TSS-containing broad domains (n = 2,313) (j) and non-TSS-containing broad domains (n = 12,691) (k) and other TSS-containing domains (remaining TSS-containing domains after excluding ZGA gene TSS-containing domains, n = 6,390) (l). Domains larger than 10 kb were included. y axis indicates DNA methylation levels and x axis indicates distance relative to the domain boundaries. Black, green, blue, and orange colour indicates sperm, oocyte, P15, and P12, respectively (P < 2.2−16, Kolmogorov–Smirnov test of sperm DNA methylation within broad domains between ZGA gene- and non-TSS/other TSS-containing domains).
a, Box plots of the shortest distance between 2-/8-cell stage restricted cREs (n = 33,609) and between randomly selected cREs (n = 33,609). 2-/8-cell common cREs are significantly closer to each other than other cREs (P < 2.2−16, Kolmogorov–Smirnov test). b, H3K27ac levels at putative cREs in the vicinity of non-ZGA genes (n = 13,334). c, Enriched motifs found in different groups of stage restricted cREs. d, Box plots of gene expression values (FPKM) for each developmental stage for genes in the vicinity of stage restricted cREs (n = 1,612 for oocyte, n = 4,166 for 2-cell, n = 2,005 for 8-cell, n = 6,236 for 2/8-cell common, and n = 1,357 for mouse ES cells). Whiskers correspond to the highest and lowest points within the 1.5× interquartile range (a, b, d). e, Accumulated fraction curve of cREs according to distance between cREs and nearest broad H3K4me3 domains. Putative cREs in the vicinity of ZGA genes are positioned closer to broad H3K4me3 domains (P < 2.2−16, Kolmogorov–Smirnov test).
Extended Data Figure 9 Histone H3 lysine 4 methyltransferases and demethylases in relation to establishment and removal of broad H3K4me3 domains.
a, Gene expression levels for multiple H3K4 methyltransferases and demethylases in the oocyte, 2- and 8-cell-stage embryos and mouse ES cells. b, Genome browser snapshot of the KMT2B dependent Cdkn1a gene. The Cdkn1a gene is not expressed in oocytes, but highly expressed in 2- and 8-cell-stage embryos. c, Average GC content in oocyte H3K4me3 covered regions maintained or demethylated in the 8-cell-stage embryo. Regions where H3K4me3 is specifically maintained show higher GC content (P < 2.2−16, Kolmogorov–Smirnov test) d, Simple western immunoblots show Kdm5a and Kdm5b morpholino knock-down of KDM5A and KDM5B proteins in 2-cell embryos. β-actin is shown as loading control. For gel source data, see Supplementary Fig. 1. e, Quantification of H3K4me3 immunofluorescence signal in late 2-cell stage embryos after injection with morpholinos targeting Kdm5a and Kdm5b (n = 37), or 5-base mismatch-control morpholinos (n = 23). (H3K4me3 signal intensity, arbitrary units. P = 2.2−4, Kolmogorov–Smirnov test.). Error bars show the s.d.
a, Scatter plots show FPKM values of GENCODE annotated genes (each dot, n = 32,287) between two biological replicates for four different conditions. Pearson correlation coefficient values are also shown. b, Box plots show the effect of α-amanitin treatment on transcriptional inhibition (Kolmogorov–Smirnov test P < 2.2−16). c, Scatter plots show FPKM values of GENCODE annotated genes (each dot, n = 32,287) between Kdm5a and Kdm5b morpholino (Kdm5a+b MO)-injected embryos (x axis) and α-amanitin-treated embryos (y axis). We present log2-transformed average FPKM values from two biological replicates. Red lines indicate twofold up- or downregulation. d, Scatter plots show FPKM values of GENCODE annotated genes between Kdm5a+b MO injected embryos (x axis) and control MO injected embryos (y axis). We present log2-transformed average FPKM values from two biological replicates. Red lines indicate 1.5-fold up- or downregulation. Grey, all genes; blue, downregulated ZGA genes. e, Box plots are shown for FPKM values of maternally expressed genes in oocyte, Kdm5a and Kdm5b morpholino replicates, and control morpholino replicates. NS indicates non-significant P value from one-sided paired t-test. In 2-cell embryos maternally expressed genes are significantly down regulated in both Kdm5a and Kdm5b morpholino and control morpholino injected embryos (P values <2.2−16, Kolmogorov–Smirnov test). Whiskers correspond to the highest and lowest points within the 1.5× interquartile range (b, e). f, Rate of development to various embryonic stages for Kdm5a and Kdm5b morpholino injected and control morpholino injected embryos resulting from natural mating (control MO, n = 15; Kdm5a+b MO, n = 26). g, Number of embryos at the 8-cell and 4-cell or abnormal stage in Kdm5a+b MO (n = 46) and control MO (n = 42) (Fisher’s exact test P = 0.028, odds ratio = 3.5), for in vitro fertilized embryos for the experiment shown in Fig. 5c.
This file contains Supplementary Tables 1-9 and Supplementary Figure 1. (PDF 174 kb)
A 3D video composed of confocal laser scanning micrographs showing H3K4me3 levels in MII oocyte. H3K4me3 (green) and DNA (DAPI, blue). (MP4 8892 kb)
About this article
Cite this article
Dahl, J., Jung, I., Aanes, H. et al. Broad histone H3K4me3 domains in mouse oocytes modulate maternal-to-zygotic transition. Nature 537, 548–552 (2016). https://doi.org/10.1038/nature19360
Chromatin alterations during the epididymal maturation of mouse sperm refine the paternally inherited epigenome
Epigenetics & Chromatin (2022)
Histochemistry and Cell Biology (2022)
Transcription factor chromatin profiling genome-wide using uliCUT&RUN in single cells and individual blastocysts
Nature Protocols (2021)
Nature Genetics (2021)
Nature Metabolism (2021)