I spent time on the plane home yesterday looking at the details of this p-value calculation, and here's what I'd say.
1) I don't think the authors erred in a literal sense in reporting this number. A p-value of p < 1×10⁻³⁰⁰ isn't narrow sense insane. With ~300,000 individuals, ~9,000 CAD events, a hazards ratio of 1.59 and a very small confidence interval [1.55–1.62], gets you a test statistic of 42 ish, which if you use a normal approximation produces an absurdly small tail probability.
2) However, that calculation is conditioned on a number of assumptions - proportional hazards, independence of errors, correct variance estimation, absence of residual confounding, absence of residual population structure, that the polygenic score is known without error, and others. Even minor departures from these assumptions, coupled with uncertainty in the construction of the PRS, would swamp the nominal tail probability. (This paper makes this point
arxiv.org/pdf/2412.20611).
3) So, in general, I think we would be well-served by being more numerically humble and accepting that high-dimensional prediction induces substantial uncertainty that classical p-values do not capture well. And while I think we can all agree that the odds that the observed relationship was the result of a sampling error are very low, reporting this as p < 10⁻³⁰⁰ exaggerates that strength of the evidence for that claim, and more importantly adds false weight to the more relevant quantities: the effect size (HR per SD), its robustness across models and populations, and the limited incremental predictive utility of PRSs in the clinical settings.
4) I stand by my claim that this p-value - while primarily a byproduct of large sample size - also represents the use of a model that is not fit for the task. And the fact that a model that would produce an "accurate" p-value (whatever that means) doesn't justify reporting values from one that are both not likely to be correct, and which are of no practical value.