PSMC-FAC: A Statistical Framework for Correcting Loss of Heterozygosity in Low-Coverage Genomic Demographic Inference

1. Introduction

Reconstructing species’ demographic histories is a main goal in evolutionary biology and conservation genetics, providing insight into how populations responded to past climatic fluctuations, geological events, and anthropogenic pressures, and helping predict future responses. Classical genetic markers, such as mitochondrial DNA and microsatellites, yield valuable information about recent population dynamics, but their limited temporal resolution restrict demographic inference largely to the most recent thousands of years [1,2,3,4]. In contrast, whole-genome sequencing (WGS) has driven the emergence of population genomic approaches capable of inferring demographic change across hundreds of thousands of years. Among these, the Pairwise Sequentially Markovian Coalescent (PSMC) method represents a landmark advance, allowing the reconstruction of effective population size (Ne) trajectories from a genome across deep evolutionary timescales [5].

PSMC infers changes in Ne through time by modeling genome-wide patterns of heterozygosity, which reflect variation in the time to the most recent common ancestor (TMRCA) between homologous genomic regions. Because demographic contractions, expansions, and periods of stability leave distinct signatures in these patterns, PSMC provides a powerful framework for investigating long-term population history under a coalescent-based model. The method is particularly informative for timescales spanning approximately 10,000 to 1–3 million years ago and has been widely applied across diverse taxa, including humans and both wild and domesticated species [5,6,7,8,9].

The theoretical foundation of PSMC lies in coalescent theory, which provides a retrospective description of genealogical relationships by tracing sampled lineages backward in time until they merge at common ancestors [10,11]. In recombining genomes, genealogies vary along the sequence, forming a complex structure known as the ancestral recombination graph (ARG). Although the ARG provides a complete representation of ancestry, its high dimensionality makes it computationally intractable for most practical applications [12]. To overcome this limitation, approximations based on Hidden Markov Models have been developed, leading to the Sequentially Markov Coalescent (SMC), which assumes that genealogical changes along the genome follow a Markov process [13]. Building on this framework, PSMC models coalescent time at each genomic position as a hidden state and the observed homozygous or heterozygous genotype as the emission, enabling inference of Ne trajectories from a single genome that, under ideal assumptions, can approximate the history of an entire population [14,15].

Despite its broad utility, largely due to its non-parametric nature, PSMC is sensitive to data quality, demographic complexity, and analytical choices. PSMC further assumes neutrally evolving regions within a panmictic population. Therefore, background selection or the inclusion of overlooked constrained genomic elements can bias coalescent rate estimates and inferences of effective population size trajectories [16,17]. Population structure can be another source of bias. Past coalescent rate depends on migration and deme configuration in addition to changes in Ne [18,19], making interpretation of inference non-trivial. Beyond this change in expected curve due to population structure, PSMC may generate spurious, abrupt peaks under structured scenarios that do not reflect the underlying coalescent dynamics [20]. Similar instabilities can also arise from misspecification of the PSMC time-interval discretization (“time vector”) [21]. Although extensions such as MSMC and SMC++ improve resolution by incorporating multiple genomes [22,23], they retain the panmixia assumption, along with biases linked to it [20]; only more recent developments explicitly accommodate structured demographic scenarios [24].

Inferring population-wide demographic history in natural populations is especially difficult when population structure is ubiquitous [25,26]. Typically, multiple genomes are required to characterize underlying demographic processes and disentangle biological signals from methodological artefacts [20]. Consequently, empirical studies frequently adopt low- to medium-coverage sequencing strategies and prioritize sampling large numbers of individuals to better capture population-level variation. When considering genome data quality, sequencing errors, false-negative heterozygous calls, and insufficient coverage can bias heterozygosity estimates and distort inferred demographic histories [6,9]. A particularly important limitation of the PSMC method arises when using low-coverage genomes. Reduced coverage decreases detection of heterozygous sites, leading to systematic underestimation of individual heterozygosity and consequently biased inference of Ne and coalescent times. PSMC trajectories therefore become displaced downward and toward more recent times, while sharp demographic changes appear attenuated and temporally shifted [4]. One strategy to mitigate this bias during plotting of the PSMC curve is to estimate the false-negative rate (FNR), a statistical correction for lost heterozygosity through rescaling of inferred coordinates from the first estimation of genetic diversity. Traditionally, FNR has been estimated by downsampling a high-coverage genome to lower coverage, plotting both together, and visually adjusting FNR values of the low-coverage curve until it partially overlaps with the high-coverage reference curve [7,8,9,27]. Although this approach allows incorporation of more genomes and improves PSMC-based inferences, visual corrections are subjective, time-consuming, and difficult to reproduce. Because generating large numbers of high-coverage genomes remains cost-prohibitive for many laboratories, population genomics will increasingly rely on medium- and low-coverage datasets. Consequently, there is a growing need for fast, objective, and mathematically grounded approaches to estimate and evaluate FNR corrections and enable robust comparisons of demographic trajectories across samples.

Here, we present PSMC FNR-Automatized Correction (PSMC-FAC), a novel method that introduces a statistical framework for optimizing FNR correction of low-coverage PSMC trajectories by minimizing graphical distances between candidate FNR-corrected low-coverage curves and their corresponding high-coverage reference curves. We assessed this approach using downsampled high-coverage publicly available WGS data from multiple populations of humans (Homo sapiens), a wild species (European gray wolves, Canis lupus), and a domesticated species (cattle, Bos taurus), representing diverse demographic histories and varying levels of genome-wide heterozygosity. We show that PSMC-FAC provides reproducible and robust FNR estimates through an automated pipeline applicable to any WGS data in BAM format, independent of genome-wide heterozygosity and demographic history. By enabling accurate demographic inference from low-coverage genomes, PSMC-FAC substantially lowers sequencing cost requirements and broadens access to population genomic analyses, thereby facilitating large-scale and comparative studies of demographic history across a wide range of taxa.

2. Materials and Methods

2.1. Dataset Preparation

We downloaded previously published WGS data for multiple populations of interest from publicly available repositories, including humans, cattle, and wolves. Wolf genomes were obtained in FASTQ format [28] and aligned to the CanFam3.1 reference genome [29] using bwa-mem v0.7.17 [30]. Reads were filtered to remove spurious alignments and PCR duplicates and were sorted following the protocol of Sarabia et al. [27]. Human WGS data were downloaded in CRAM format [31] from the 1000 Genomes Project [32] and included nine individuals from three different geographical origins: three Han Chinese (CHB), three Yoruba from Ibadan, Nigeria (YRI), and three Tuscan from Italy (TSI). Files in CRAM format were decompressed and transformed to BAM format using samtools [33]. Cattle genomes were obtained in BAM format from the public Agricultural Research Service of the United States Department of Agriculture [34,35]. Genome-wide sequencing depth for all samples was calculated using the samtools depth function [36]. Only samples with mean coverage >15× were retained to enable downsampling. Metadata for all downloaded samples are provided in Table A1.

Downsampling was performed on all BAM files to multiple target depths (5×–15×) using the samtools view -s function. Genome-wide depth of coverage for each downsampled BAM file was verified using the samtools depth function (Table A1). Following the PSMC manual specifications [5], downsampled BAM files were converted to variant call format (VCF) using bcftools mpileup and call functions [36]. To run this step, the CanFam3.1 reference genome [29] was used for wolves, the GRCh38 (hg38) reference genome [37] for humans, and the bosTau9 reference genome [38] for cattle. Variant Call Format (VCF) files [39] were then converted to FASTQ format using the bcftools view function and the vcfutils.pl vcf2fq utility [5], filtering out variants with coverage lower than 5 (-d 5) and higher than twice the average depth of coverage (-D), retaining most heterozygous sites while filtering out gene duplicates. Finally, FASTQ files were converted to PSMC input format (PSMCFA) using fq2psmcfa [5] with a minimum phred-scaled base quality threshold of 20.

We obtained a PSMC inference from each individual sample and its downsamples and followed the literature to define time interval patterns. For wolves, we applied a customized time interval pattern of 64 atomic intervals arranged as six intervals of size one followed by 58 intervals of size one (16+581), as described in Freedman et al. [40]. For humans and cattle, we used the default PSMC time interval pattern (4+25*2+4+6), following Li and Durbin [5] for humans and Mei et al. [41] and Liu et al. [42] for cattle (see Appendix A for specific commands and the pipeline in Figure 1).

For each PSMC trajectory inferred from a downsampled genome, we applied False Negative Rate (FNR) correction across values ranging from 0 (no correction) to 0.99 in increments of 0.01, generating 100 corrected PSMC trajectories per downsample. Each corrected trajectory was then compared to the PSMC curve inferred from the original (non-downsampled) sequence. Because PSMC trajectories are defined over discretized time intervals that differ between inferences, direct pointwise comparison is not straightforward. To standardize comparisons, we defined a custom vector of 60 logarithmically spaced time points (in years) and projected all trajectories onto these shared temporal coordinates.

For cattle and wolves, the custom vector spanned 10 kya to 1.5 Mya. For humans, we repeated the analysis restricting the range to 50 kya to 1.5 Mya to avoid recent-time instabilities in PSMC. Early analyses attributed these peaks to Out-of-Africa demographic history [5], whereas more recent work has shown that they can arise from population structure [20] and from sensitivity to time-interval discretization [21]. Restricting the lower bound from 10 to 50kya was done to mitigate the influence of these artifacts.

2.2. A Mathematical Approach to Compute FNR for Low Coverage Samples

To quantify differences between PSMC trajectories, we employed two complementary distance-based metrics: (i) Hausdorff distance [43] and (ii) the discrete Fréchet distance [44]. Each PSMC trajectory is represented as an ordered set of points in the plane, where the x-coordinate corresponds to log-scaled time and the y-coordinate corresponds to inferred effective population size (Ne). Let

$$P=\left\{p_1,p_2,...,p_m\right\};Q=\left\{q_1,q_2,...,q_m\right\}$$

(1)

denote two PSMC curves, where P represents the reference trajectory obtained from a high-coverage genome (for example, 20×) and Q represents a trajectory derived from a downsampled genome after FNR correction. Before distance computation, trajectories are projected onto a shared temporal grid so that corresponding demographic histories can be compared consistently while preserving the piecewise-constant nature of PSMC inference.

2.2.1. Hausdorff Distance

The undirected Hausdorff distance measures the maximal local discrepancy between two curves. It is defined as:

$$H\left(P,Q\right)=max\left\{\begin{array}{c}max\\ p_i\in P\end{array}\begin{array}{c}min\\ q_j\in Q\end{array}‖p_i-q_j‖,\begin{array}{c}max\\ q_j\in Q\end{array}\begin{array}{c}min\\ p_i\in P\end{array}‖q_j-p_i‖\right\}$$

(2)

where ||.|| denotes the Euclidean norm [43].

For each point in one curve, the minimum distance to the other curve is computed. The Hausdorff distance corresponds to the largest of these minimal distances and therefore captures the worst-case disagreement between trajectories. This metric is sensitive to localized deviations, such as abrupt expansions or bottlenecks, and is useful for identifying regions where curves diverge strongly.

2.2.2. Discrete Fréchet Distance

While the Hausdorff distance ignores the sequential ordering of points, demographic trajectories represent ordered processes through time. To account for temporal progression, PSMC-FAC additionally uses the discrete Fréchet distance, which preserves ordering along both curves and measures similarity in overall shape. The discrete Fréchet distance between two polygonal curves is defined as:

$${\delta }_F\left(P,Q\right)=\begin{array}{c}min\\ \delta ,\tau \end{array}\begin{array}{c}max\\ \kappa \end{array}‖P_{\sigma \left(k\right)}-Q_{\tau \left(k\right)}‖$$

(3)

where σ and τ are non-decreasing index sequences that traverse the points of P and Q, respectively, from start to end, preserving ordering [44].

Intuitively, the Fréchet distance measures the minimum leash length required for two entities to traverse both curves without backtracking. Because temporal ordering is maintained, this metric captures global similarity in demographic trajectory shape rather than isolated local differences.

2.2.3. A Combination of Both Methods

By using both metrics, we can detect localized deviations and global differences in curve geometry, providing a more comprehensive and robust comparison than either metric alone (Appendix A). A custom distance-calculation script computed both the discrete Fréchet distance and the Hausdorff distance between the original path and each FNR-corrected path.

After evaluating the accuracy of both the Hausdorff and Fréchet distances for our samples, we defined the optimal FNR value for each downsampled genome as the FNR corresponding to the minimum discrete Fréchet distance relative to the original (high-coverage) trajectory. The Fréchet metric was selected as the primary optimization criterion because it preserves point ordering and thus reflects similarity in overall demographic trajectory shape and temporal progression. Hausdorff distances were retained to evaluate localized maximal deviations.

The set of optimal FNR values obtained across coverage levels was subsequently modeled using a polynomial regression of degree 2 to generate continuous correction curves. These fitted polynomial functions describe the relationship between sequencing depth and the FNR value that minimizes trajectory divergence, thereby enabling depth-dependent calibration of Ne estimates.

2.3. FNR and Heterozygosity as a Function of Coverage

We evaluated whether the optimal FNR corrections depend not only on genome-wide heterozygosity and sequencing coverage, but also on the demographic history of the population under study. To assess this, we systematically examined the relationship between sequencing depth and the most optimal FNR factor inferred for each individual. To model the relationship between sequencing depth and optimal FNR, we fitted a second-degree polynomial regression in which coverage was treated as the independent variable, and the most optimal FNR value as the dependent variable. This polynomial model was chosen to capture non-linear scaling of FNR with sequencing depth. The regression was fitted independently for each individual.

Because optimal FNR values were determined using two alternative distance metrics, we repeated the polynomial fitting procedure separately for FNR estimates obtained through minimization of the Hausdorff distance [43] and for those obtained through minimization of the Fréchet distance [44], allowing comparison between both metrics. Following Nadachowska-Brzyska et al. [9], we assumed that the optimal FNR is equal to zero for samples with mean coverage greater than 15×. Accordingly, for regression purposes, high-coverage genomes exceeding this threshold were treated as having FNR = 0.

2.4. Sum-of-Least-Squares Assessment of Goodness-of-Fit

To evaluate the effectiveness of FNR-based correction across coverage levels, we quantified the similarity between PSMC trajectories using an independent sum-of-least-squares criterion. For each downsampled individual and each coverage level (5×–15×), we first identified the most optimal FNR value by minimizing either the Hausdorff distance [43] or the Fréchet distance [44] between the FNR-corrected trajectory and the corresponding high-coverage (non-downsampled) reference trajectory. This procedure yielded, for each downsample and for each distance metric, a single corrected curve with an optimal FNR value.

Subsequently, and independently of the optimization step, we computed the Sum of Squared Errors (SSE) between the best FNR-corrected trajectory and the high-coverage reference trajectory. SSE was therefore used exclusively as an evaluation metric and was not involved in the selection of the optimal FNR value. To calculate SSE, at each time point of the custom vector, we obtained the difference in effective population size (Ne) between the low-coverage and high-coverage trajectories, squared, and summed across all time points. This procedure was repeated for each downsampled coverage level and for each individual separately, generating a set of SSE values corresponding to the best FNR correction at each depth. These SSE values were then used to assess how closely the optimally corrected trajectories approximated the high-coverage reference across coverage levels. Comparisons were based on direct inspection of SSE trends across depths and individuals.

3. Results

3.1. PSMC-FAC Enables Accurate FNR-Based Correction Across Species and Coverages

Across all three species (cattle, wolves, and humans) and populations analyzed, PSMC-FAC produced mathematically consistent corrections of the false negative rate (FNR), substantially improving concordance between downsampled and original (~20×) PSMC trajectories (Figure 2, Fig. A1). After projection onto a common logarithmic time grid with a custom time vector and comparison to the corresponding high-coverage inference, FNR-corrected trajectories generally converged toward the reference curve across the evaluated temporal window (Fig. A1).

Correction performance varied depending on both sequencing depth and demographic profile. In general, reconstructions from genomes of higher coverages (10–15×) showed close agreement with the original trajectory after FNR correction. In contrast, medium-to-low coverage datasets (5–6×) showed greater instability and reduced correction accuracy (Fig. A1). At 5× coverage, several corrected trajectories of gray wolves displayed slight increases in population size when the original 20× curve showed a slight decrease (Fig. A1 panels D2, E2, G2). This pattern, although apparently paradoxical, is consistent with expectations from the variant-calling pipeline, since a minimum depth filter of 5 (-d 5) was applied during PSMC preprocessing, causing sites with depth lower than 5× to be ignored. In addition, the reference genome used for wolves (CanFam3.1) is known to exhibit relatively low heterozygosity [29].

PSMC-FAC was also able to adjust downsampled curves toward high-coverage trajectories in demographic scenarios characterized by abrupt population size changes, such as strong bottlenecks or rapid expansions (e.g., Fig. A1 panels A1–A2, D1–D2, F1–F2, G1–G2). In such cases, corrected curves tended to smooth extreme transitions, suggesting that strong local curvature in inferred coalescent rates is especially sensitive to heterozygosity loss caused by downsampling.

3.2. Appropriate FNR-Based Correction Depends on Recent Demographic History

Correction performance differed among human populations and appeared to depend on recent demographic history. While Yoruba (YRI) trajectories were consistently well corrected across coverages (Fig. A1 panels M–O), PSMC-FAC showed reduced performance for Han (CHB; Figure S1 panels J–L) and Toscani (TSI; Fig. A1 panels P–R) genomes. This discrepancy coincides with a pronounced population size increase inferred in non-African samples between ~50–40 kya, where PSMC trajectories display a large, sharp elevation in effective population size during this interval. Because FNR optimization relies on minimizing Hausdorff and Fréchet distances between corrected and original trajectories [42,43], extreme local deviations can disproportionately influence the fitting procedure. If the last fragment of the demographic trajectory presents a sharp elevation, the algorithm will be biased to optimize correction parameters to better match this last maximum rather than the overall shape of the curve across the remaining temporal range.

To evaluate how PSMC-FAC adjusts downsampled curves to the original if the last peak is ignored, we repeated the FNR correction restricting the custom vector for the comparison to 50 kya–1.5 Mya for all genomes (Fig. A2). This modification produced no substantial change for wolves or cattle (Fig. A2 panels A–I), nor for Yoruba genomes (Fig. A2 panels M2–O2), where the discrepancy was not found previously. In contrast, correction performance for Han (Fig. A2 panels J2–L2) and Toscani (Fig. A2 panels P2–R2) genomes improved markedly. When the recent sharp elevation in the curve was excluded from the fitting interval, corrected trajectories showed substantially better global concordance with the original high-coverage curves.

3.3. FNR Corrections Are Robust Across Diverse Demographic Histories

Across species and populations, optimization of FNR values produced highly concordant results between the two distance metrics used. When corrections were computed over the 10 kya–1.5 Mya window, Hausdorff and Fréchet distances were strongly correlated across the full panel (Figure 3 and Figure A3, Table A2), indicating that both metrics identified nearly identical optimal FNR values. The only clear exceptions were the non-African human populations (Han and Toscani), consistent with the influence of the pronounced recent-time Ne peak described above. For all genomes, the relationship between sequencing coverage and inferred optimal FNR followed a highly regular polynomial trend. In every case, the fitted models yielded R² > 0.99, demonstrating a very close relationship between depth of coverage and the magnitude of FNR correction required. This pattern was consistent across cattle, wolves, and humans, despite their distinct demographic histories and differences in reference genomes and time interval parametrization. Importantly, Hausdorff and Fréchet distances performed nearly equivalently in identifying optimal FNR values. Differences in the selected FNR were minimal and resulted in only marginal variation in goodness-of-fit statistics for the coverage-FNR regression curves. Thus, correction accuracy does not depend strongly on the specific choice of distance metric, reinforcing the robustness of the PSMC-FAC optimization framework.

After restricting the optimization window to 50 kya–1.5 Mya and re-estimating FNR values, the previously observed discrepancy in non-African human populations largely disappeared (Fig. A4; Table A3). Han and Toscani genomes then showed near-perfect correlations between genomic coverage and optimal FNR, comparable to those observed in Yoruba, wolves, and cattle. In contrast, only negligible changes were observed for the other genomes, whose trajectories lack the extreme recent-time Ne increase. These results strongly suggest that deviations observed in specific populations were not due to intrinsic limitations of the correction framework but rather to localized, high-amplitude features in the recent portion of the trajectory that disproportionately influenced distance-based optimization. Once these features were excluded, the coverage-FNR relationship became uniformly stable and highly predictable across all species and populations analyzed, making it possible to use it as a mathematical tool to predict the FNR needed for other coverages.

Although FNR optimization behaved consistently within species and populations, the shape of the polynomial relationship between sequencing coverage and optimal FNR differed among taxa (Figs. A3-A5). These differences highlight the distinct demographic histories and genome-wide heterozygosity levels observed across cattle, wolves, and humans (Table A1). Thus, while the coverage-FNR relationship is highly predictable within a population (R² > 0.99 in all cases), it is not universally transferable across populations with divergent demographic backgrounds, so running a PSMC-FAC pipeline per species or per population is highly advisable.

4. Discussion

In conservation and population genomics studies, there is an inherent trade-off between sequencing depth and sampling breadth [45,46]. High-coverage whole-genome sequencing remains costly, and generating large numbers of genomes at ~20× depth is often financially prohibitive for many research groups. Conversely, sequencing a small number of individuals at high coverage risks failing to capture the demographic heterogeneity of structured populations [47]. Although low-coverage genomes are frequently considered unsuitable for demographic inference, previous work has shown that meaningful demographic signals can be retained even under reduced representation approaches such as RAD sequencing under specific conditions [48]. Nadachowska-Brzyska et al. [9] demonstrated that PSMC performance declines substantially below approximately 10× coverage and recommended higher depths for stable inference. Building on these observations, we show that low- and medium-coverage genomes can still yield informative demographic inference when the loss of heterozygosity is quantitatively corrected using a False Negative Rate (FNR) factor. By deriving population-specific calibration curves linking FNR and sequencing coverage, the PSMC-FAC framework reduces uncertainty and allows demographic analyses to rely on a combination of high- and low-coverage genomes, thereby lowering overall sequencing costs for sound population-based demographic inference.

4.1. FNR Correction in Low- and Mid-Depth Genomes: Reference Genome Effect

At low coverage (e.g., ~5×), error rates remain substantial and FNR correction cannot fully compensate for stochastic loss of information. In contrast, between approximately 8× and 15× coverage, corrected trajectories converge closely toward high-coverage references, substantially expanding the practical usability of medium-depth sequencing data. The reduced performance observed at very low coverage reflects limitations of the variant-calling pipeline rather than a failure of the correction framework itself. During the PSMC preprocessing, a minimum depth threshold of five reads is applied (parameter -d 5). When mean coverage is ~5×, a substantial proportion of true heterozygous sites inevitably falls below this threshold and is therefore excluded. These sites are subsequently treated as homozygous reference during consensus generation, leading to systematic underestimation of heterozygosity. Because PSMC infers demographic history from the spatial distribution of heterozygous tracts along the genome, such undercalling directly alters the inferred coalescent rate. The resulting distortion is therefore not random noise but a predictable shift in trajectory shape, affecting inferred effective population size and potentially generating spurious demographic features. Similar coverage-dependent biases have been documented previously, particularly by Nadachowska-Brzyska et al. [9], who showed that decreasing depth disproportionately impairs heterozygote detection and destabilizes demographic inference. Our results are consistent with these observations. Importantly, while PSMC-FAC corrects for heterozygosity loss attributable to false negatives, at sufficiently low coverage both stochastic and systematic variation introduce an irreducible source of error. This is reflected in the elevated sum-of-squared errors observed even at the optimal FNR value at the 5× downsamplings. Thus, while the framework substantially mitigates coverage-related bias, it cannot completely overcome the intrinsic limitations of very shallow sequencing data.

These effects may become particularly pronounced when the reference genome does not adequately represent the genome-wide diversity of the population under study. For example, the human reference assembly (hg19/GRCh37) was constructed from a limited number of individuals and does not capture global human genetic diversity [49,50]. Similarly, the domestic dog reference genome, although widely used in canid studies, does not fully represent the diversity present across wolf and dog populations [51,52]. Under such circumstances, mapping bias toward the reference allele can further exacerbate heterozygote undercalling, amplifying the apparent loss of heterozygosity at low coverage. Consequently, the accuracy of FNR correction depends not only on sequencing depth but also on the appropriateness of the chosen reference genome. This highlights the importance of considering reference bias when interpreting demographic reconstructions derived from low-coverage data.

4.2. Effects of Biases Introduced by PSMC Assumptions on Optimal FNR Calculation

The pronounced peak observed in non-African human populations between approximately 50 and 40 thousand years ago provides a revealing example of how features of a PSMC trajectory can influence FNR optimization. When a broad temporal window starting at 10 kya was used for distance minimization, this sharp recent-time elevation disproportionately influenced both Fréchet and Hausdorff distance metrics. As a consequence, optimization tended to favor FNR values that improved agreement around this localized feature rather than across the trajectory as a whole. Similar peaks were already noted in the original PSMC publication by Li and Durbin (2011), and their biological interpretation has remained controversial. More recent work has demonstrated that population structure alone can generate such sharp biasing peaks in SMC-based inference due to coalescent signatures (Nieto et al. 2025). Tournebize et al. (2025) further demonstrated that structure and admixture can produce nearly indistinguishable coalescent patterns. Admixture with Neanderthals between approximately 50 and 43 kya remains a plausible biological contributor to this pattern, as interspecific gene flow is known to generate transient increases in inferred effective population size under SMC-based frameworks (Cahill et al 2017; Sarabia et al. 2021). However, because these peaks are also sensitive to the temporal discretization parameters used in PSMC (Hilgers et al., 2025), this feature is interpreted conservatively as a methodological artefact lacking straightforward demographic interpretation. When the optimization window was restricted to exclude this recent interval (starting at 50kya), FNR–coverage relationships stabilized markedly across Toscani and Han populations, whereas little change was observed for other populations and species. This indicates that the instability did not originate from the correction framework itself but rather from biases present in the reference trajectory. In practice, adjustment of PSMC hyperparameters and careful selection of the temporal window should therefore precede FNR optimization in order to reduce artefact-driven bias in the reference curve.

This observation leads to a central conceptual clarification: PSMC-FAC optimizes trajectory similarity rather than recovering true demography. The method minimizes geometric distance between a low-coverage trajectory and a chosen high-coverage reference; consequently, any artefacts contained in the reference trajectory are inherently propagated through the correction process. If the reference curve contains artefacts arising from structure (Nieto et. al., 2025), background selection (Cousins et al 2024), time discretization (Hilgers et al 2025), inclusion of functional genomic elements in the problem sequence (Sellinger et al 2021), or model misspecification (Cousins et al 2025), the FNR correction will reproduce them. PSMC-FAC therefore corrects coverage-dependent heterozygosity loss but does not address the theoretical or methodological limitations intrinsic to the sequentially Markovian coalescent framework. Sellinger et al. (2021) demonstrated that SMC-based methods have intrinsic convergence properties sensitive to mutation–recombination ratio assumptions, potentially limiting applicability in some populations. The present framework does not resolve such limitations; instead, it isolates and corrects one specific and pervasive source of bias: false-negative heterozygous calls caused by limited sequencing depth.

4.3. Polynomial Relationship Between Coverage and Optimal FNR

The strong and consistent polynomial relationship observed between sequencing coverage and optimal FNR indicates that heterozygosity loss scales predictably with depth under stable analytical conditions. Within individual populations, this relationship is highly regular, allowing coverage-dependent correction curves to be estimated with high confidence. However, the shape and parameters of these polynomial functions differ across species and populations, reflecting variation in genome-wide heterozygosity levels, recombination landscapes, and demographic history.

Demographic features themselves also influence optimization behavior. Sharp events, such as rapid bottlenecks or expansions, generate regions of high local curvature in PSMC trajectories. These features can disproportionately influence distance-based comparisons and, consequently, shift the inferred optimal FNR values. Thus, FNR estimation is determined not only by sequencing coverage but also by the geometric properties of the underlying demographic trajectory. The combined use of Hausdorff and Fréchet distances [42,43] provides complementary perspectives on this problem. The Hausdorff distance is sensitive to localized maximum deviations, whereas the Fréchet distance captures overall trajectory similarity while preserving temporal ordering. Together, these metrics offer a robust framework for estimating coverage-dependent correction parameters while accounting for variation in demographic trajectory shape.

4.4. Empirical Applications and Future Implications

In practical terms, this framework enables demographic analyses in conservation and population genomics contexts where sequencing large numbers of individuals at high coverage is economically unfeasible. By calibrating coverage-dependent bias using a high-coverage reference genome, additional individuals sequenced at moderate depth can be incorporated into demographic analyses without relying on subjective or visually determined FNR adjustments. This substantially expands the range of sampling designs compatible with PSMC-based inference while maintaining methodological consistency across individuals.

These results should also be interpreted in light of the known temporal limits of PSMC inference. Patton et al. [55] demonstrated that SMC-based methods achieve their highest resolution at intermediate timescales, whereas accuracy declines for very recent demographic events. More recently, Peede et al. [56] emphasized that temporal resolution depends strongly on the distribution of coalescent and recombination events along the genome. To accommodate these constraints, PSMC-FAC allows the specification of custom temporal windows during FNR optimization, enabling the exclusion of intervals known to contain unstable or artefactual signals. Consequently, low-coverage genomes corrected using PSMC-FAC should not be interpreted at very recent or fine-scale temporal resolutions. Still, they can provide robust information for intermediate and deep-time demographic inference where PSMC performs most reliably.

PSMC-FAC provides a reproducible and mathematically grounded procedure for calibrating coverage-dependent bias. By replacing subjective visual adjustment of FNR values with distance-based optimization metrics, the framework improves reproducibility and allows explicit quantification of correction performance. Although currently implemented as a correction layer for PSMC trajectories, the underlying conceptual approach could potentially be extended to other SMC-based methods that rely on heterozygosity patterns, provided that suitable high-coverage reference data are available.

Overall, these results demonstrate that low-coverage genomes are not inherently unsuitable for demographic inference. When appropriately calibrated, they can approximate high-coverage demographic trajectories within predictable limits. By correcting the component of distortion attributable specifically to sequencing depth, PSMC-FAC lowers the practical barrier to demographic reconstruction and enables broader sampling strategies, making population-scale PSMC analyses more accessible to laboratories for which large-scale high-coverage sequencing remains economically unfeasible.

5. Conclusions

This study introduces PSMC-FAC, an automated and reproducible framework that corrects coverage-dependent biases in demographic reconstructions derived from low-coverage genomes. By replacing subjective visual adjustments with objective, distance-based optimization, the method enables consistent estimation of false-negative rate (FNR) corrections and substantially improves concordance between low- and high-coverage demographic trajectories across diverse species and demographic scenarios.

Our results show that low- and medium-coverage genomes can provide reliable demographic information when appropriately calibrated using population-specific correction curves. Although the framework does not address inherent theoretical limitations of Sequentially Markovian Coalescent models, it effectively corrects a major practical source of bias associated with sequencing depth. By reducing reliance on high-coverage datasets, PSMC-FAC lowers sequencing cost constraints, supports broader sampling designs, and expands the practical applicability of demographic inference in evolutionary biology, conservation genomics, and comparative population studies.