The circulating cell-free DNA landscape in sepsis is dominated by impaired liver clearance

Abstract

Circulating cell-free DNA (cfDNA) is a promising molecular biomarker, but its role in severe infection is unclear. Here, we profile cfDNA from sepsis patients and controls, demonstrating a 41-fold increase during disease. Methylation-based deconvolution revealed similar cfDNA compositions in the two groups, suggesting that cfDNA accumulation during disease is due not to excess cell death but to impaired hepatic clearance. Fragmentation and end-motif patterns both support this hypothesis, suggesting prolonged exposure of cfDNA to circulating nucleases. In addition, we show that cfDNA retains nucleosome footprints informative of gene activity. By developing a novel method to quantify these footprints and integrate them with single-cell data, we report an increase in cfDNA from Kupffer cells and liver parenchyma in patients with liver dysfunction. Finally, we show that cfDNA contains pathogen-derived material, highlighting its diagnostic potential. This high-throughput, multimodal study provides a reference for understanding cfDNA's role in sepsis and critical illness.

Keywords: DNA methylation; cell-free DNA; critical illness; epigenetics; epigenomics; fragmentomics; infection; liquid biopsies; metagenomics; sepsis.

Graphical abstract

Figure 1

Levels of cfDNA in circulation reflect sepsis severity (A) cfDNA concentration in the plasma of sepsis patients and healthy controls (top) and in sepsis patients within different hospital care settings (bottom). p values were derived using a Wilcoxon rank-sum test. Boxplots represent median and interquartile ranges (IQRs) of cfDNA concentration. (B) cfDNA concentration in sepsis patients stratified by level of respiratory support. Categories are ordered by increasing invasiveness. p values were derived using a Wilcoxon rank-sum test. Boxplots represent median and IQRs of cfDNA concentration. (C) cfDNA concentration in sepsis patients stratified by requirement for treatment with vasopressors or inotropes. p values were derived using a Wilcoxon rank-sum test. Boxplots represent median and IQRs of cfDNA concentration. (D) Correlation between mean arterial pressure (x axis) and cfDNA concentration (y axis). p values and correlation estimates were computed with a Pearson correlation test.

Figure 2

An experimental approach for multimodal profiling of circulating cfDNA (A) Schematic of possible mechanisms contributing and clearing DNA from the circulation. Solid arrows represent putative causal relationships between disease processes. Dashed red and blue arrows represent possible mechanisms of cfDNA production and clearance during sepsis, respectively. (B) Experimental approach followed in this study. Blood samples were collected from sepsis patients and used for cfDNA isolation. cfDNA was subsequently used for library preparation and sequencing with the TAPS protocol, which enabled simultaneous characterization of its methylation status, its fragmentation patterns, and evidence of nucleosome positioning. This information was then related back to clinical measures of illness severity and used to infer the likely tissues of origin of cfDNA.

Figure 3

Variation in the cfDNA methylome reflects organ function and disease processes (A) mCpG conversion rates in positive and negative spike-in controls as estimated based on our full cohort. Each dot represents a sample, with colors indicating library preparation batch. (B) Two-dimensional-density plot showing the variability (standard deviation, x axis) and average proportion (%, y axis) of cfDNA methylation at 1-kb genomic windows as estimated based on our full cohort. Shades of color are proportional to the number of genomic regions located within a given range. (C) Histogram of cfDNA methylation variability (standard deviation) across all 1-kb windows in the genome (top) as estimated based on our full cohort. The proportional overlap between 1-kb windows and known genomic annotations is shown for the top 10% most variable and the bottom 10% least variable genomic windows (bottom). (D) Two-dimensional-density plots showing the variability and average cfDNA methylation at first intron (top) and promoter (bottom) regions as estimated based on our full cohort. Shades of color are proportional to the number of regions located within a given range. (E) The proportion of methylation variance explained by different variables was estimated at intronic regions using variance partitioning analysis. This analysis was conducted separately for the full cohort (left) and for sepsis patients only (right). Violin plots show the distribution of variance explained by each variable, with dots representing estimates for individual introns. (F) Principal-component analysis based on methylation at first introns. Dots represent samples, with colors indicating either disease status (top) or ALT concentration (bottom). Crosses indicate samples for which ALT measurements were not available. (G) Volcano plot showing GSEA-derived enrichment scores (x axis) and FDR-adjusted p values (y axis) for genomic regions ranked by strength of association between cfDNA methylation and disease status (left). Enrichment was estimated based on our full cohort. A subset of gene sets is highlighted. Enrichment score distributions for the topmost enriched gene sets are also shown (right). (H) Volcano plot showing GSEA-derived enrichment scores (x axis) and FDR-adjusted p values (y axis) for genomic regions ranked by strength of association between cfDNA methylation and ALT levels (left). Enrichment was estimated based on sepsis samples only. A subset of gene sets is highlighted. Enrichment score distributions for the topmost enriched gene sets are also shown (right).

Figure 4

Analysis of tissues of origin of cfDNA during sepsis (A) Proportion of cfDNA estimated to arise from different tissues based on methylome deconvolution. Boxplots show median and interquartile ranges (IQRs) of estimated proportions in sepsis patients (red; n = 49 samples from 31 patients) and healthy controls (blue; n = 7 independent samples). (B) Correlations between time since hospital admission and proportion of MEP-derived cfDNA (left), as well as ALT levels and proportion of hepatocyte-derived cfDNA (right). Correlation coefficients and p values were estimated using Pearson correlation tests. (C) Correlation between proportions of cells in circulation (x axis) and their cellular contributions to cfDNA (y axis). Colors indicate the three main immune cell lineages. The identity line is shown for reference. (D) Cellular turnover rates estimated based on cfDNA GEQs and circulating cell proportions. Boxplots represent median and IQRs of cellular turnover proportions. The dashed red line indicates a turnover rate of 1, corresponding to cell death and cell generation being equal. (E) Correlation between markers of liver dysfunction (x axis) and hepatocyte-derived cfDNA GEQs (y axis). Blue lines and shaded regions indicate linear fits and confidence intervals. Correlation coefficients and p values were estimated using Pearson correlation tests. (F) Correlation between clinical markers of erythropoiesis (x axis) and MEP-derived cfDNA GEQs (y axis). Blue lines and shaded regions indicate linear fits and confidence intervals. Correlation coefficients and p values were estimated using Pearson correlation tests.

Figure 5

cfDNA fragmentation reveals impaired hepatic clearance during sepsis (A) The estimated percentage of cfDNA fragments (y axis) of different lengths (x axis) in linear (top) and logarithmic (bottom) scales. Each line represents a sample, with colors indicating control and sepsis patient samples (n = 51 samples from 31 patients and 7 control samples). (B) Percentage of cfDNA fragments classified as subnucleosomal, mononucleosomal, dinucleosomal, or polynucleosomal based on their observed length. Proportions are shown in both linear (top) and logarithmic (bottom) scales. Each dot represents a sample, with colors indicating disease status. Boxplots show median and IQR values for each sample group. p values were estimated using Wilcoxon rank-sum tests. (C) cfDNA fragmentation indices (y axis, in logarithmic scale) stratified by disease status (x axis). Each dot represents a sample, with colors indicating disease status. Boxplots show median and IQR values for each sample group. p values were estimated using Wilcoxon rank-sum tests. (D) cfDNA fragmentation indices (y axis, in logarithmic scale) stratified by time since admission (x axis). Each dot represents a sample, with colors indicating disease status. Boxplots show median and IQR values for each group. p values were estimated using a Kruskal-Wallis test. For sepsis samples, correlation coefficients and p values between fragmentation indices and time were computed using a Pearson correlation test and are also shown. (E) Comparison of average observed (y axis) and expected (x axis) cfDNA end-motif frequencies in our cohort. Each dot represents a 4-bp 5′ end motif, with a subset of motifs off the diagonal being highlighted. The identity line is shown for reference. (F) Principal-component analysis plot based on the frequencies of all 4-bp 5′ end motifs. Each dot represents a sample, with colors indicating disease status. (G) Volcano plot showing the correlation between end-motif frequency and bilirubin. Each dot represents a motif, with its correlation coefficient (x axis) and FDR-adjusted p value (y axis) from Pearson correlation tests shown. Significantly positively and negatively correlated motifs are highlighted in red and blue, respectively. (H) Heatmap of 5′ cfDNA end-motif frequencies across all samples in our cohort. Shades of color represent end-motif frequencies, with marginal color bars indicating the first nucleotide in the motif (horizontal axis) as well as disease status and liver-specific sequential organ failure assessment (SOFA) scores (vertical axis). Samples were grouped using hierarchical clustering of motif frequencies. A subset of motifs showing appreciable differences in frequency is highlighted. (I) Correlation plot between all pairwise combinations of liver function tests and cfDNA features. Each square indicates results from a pairwise correlation test, with shades of color indicating the estimated Pearson correlation coefficient. Variable names are shown, with colors indicating whether they were derived from liver function testing or from cfDNA sequencing. Squares were grouped using hierarchical clustering.

Figure 6

Cell-free DNA retains cell-type-specific nucleosome positioning signatures (A) Average WPS values and their associated 95% confidence intervals were estimated for all known gene regions using a 5-kb window centered at the TSS. Blue and red lines indicate estimates for sepsis patients and healthy controls, respectively (n = 51 samples from 31 patients and 7 control samples). Dotted vertical lines indicate the expected position of each nucleosome and are spaced 187 bp from one another. A regional plot focusing on a 2-kb region around the TSS is also shown (top). (B) Average WPS values at the TSS regions of active (right) and inactive genes (left), based on gene expression measurements from the GTEx study. Colors indicate disease status, and dotted vertical lines indicate the expected position of each nucleosome. (C) Average WPS values at the TSS regions of increasingly more active genes, based on gene expression measurements from the GAinS study. Colors indicate disease status. (D) Observed WPS values at the TSS regions of all known genes (solid lines) are shown alongside predictions derived from a dampened harmonic oscillator model (dotted lines). Colors indicate disease status. The equation used for model fitting is shown for reference. (E) Dampened harmonic oscillator models were fitted separately to each sample in our study. The coefficient of variation of estimates for each model parameter is shown as a bar plot. (F) Heatmap of correlation estimates computed between all model parameter estimates. Colors indicate Pearson correlation coefficients. Rows and columns were ordered by similarity using hierarchical clustering. (G) Dampened harmonic oscillator models were fitted to the TSS region of tissue-specific gene sets. Bar plots show the Pearson correlation coefficient (y axis) between each model parameter estimate and the proportion of hepatocyte-derived cfDNA as estimated based on DNA methylation. Each plot shows estimates for a different model parameter, with colors indicating the gene set used during model fitting. (H) Volcano plot of correlations between hepatocyte-derived cfDNA proportions estimated from methylome deconvolution and WPS exponential decay rates estimated for different liver cell types. The dotted line indicates the statistical significance threshold of Benjamini-Hochberg (BH) adjusted p < 0.05. Colors indicate the gene set used during model fitting.

Figure 7

Cell-free DNA contains information on infecting pathogens (A) Volcano plot showing microbial genera that are differentially abundant between sepsis patients and healthy controls. Microbial abundance in circulation (x axis) and statistical evidence (y axis) are shown, with each dot representing a microbial genus. Red dots indicate known contaminants reported in previous metagenomic studies. Marginal density plots indicate the distribution of known contaminant species. (B) Principal-component analysis was performed on a matrix of microbial abundances from sepsis patient samples. Each dot represents a sample, with colors highlighting microbiological results from blood cultures. Marginal plots show the contribution of different genera to each principal component. (C) Histogram showing the distribution of microbial abundances in cfDNA. The dotted line indicates the position of an outlier filter, which we used as a threshold to separate true signals from unspecific background. (D) Bar plot of abundances for the 11 genera passing our outlier filter threshold. Each bar represents a sample that tested positive for the genus in question, with abundance estimates shown on the x axis. (E) Stacked bar plot showing the proportions of microbial genera detected in cfDNA that are also detected in blood or bodily fluid cultures. These proportions were compared to those of the cfDNA background using a Fisher’s exact test.