In my knot theory notes I have written down the following statement:
Let $L$ be a link with $k$ components, then $P_L(x,y) - 1$ is divisible by $x+y-1$ and $P_L(x, y) + (-1)^k$ is divisible by $x+y+1$, thus HOMFLY polynomial never vanishes.
That HOMFLY never vanishes follows from similar statement for Jones polynomial (see Why does the Jones polynomial never vanish?), but I am still interested in the reference for the divisibility part.
I am not sure which parametrisation this is, probably three-variable one from https://ncatlab.org/nlab/show/HOMFLY-PT+polynomial evaluated in $z = 1$.
For the polynomial $P(L)$ satisfying $x \cdot P(L_+) + y \cdot P(L_-) + P(L_0) = 0$, which is what we get when we set $z=1$ in the three-variable formula, we can show that $P(L)-1$ is divisible by $x+y+1$. This is not the divisibility relation you want, but I expect that
Induct on the number of crossings. When $L$ has no crossings, it consists of $k$ copies of the unknot, and $P(L) = (-x-y)^{k-1}$, so $P(L)-1 = (-x-y)^{k-1} - 1$ has a factor of $-x-y-1$, and is divisible by $x+y+1$ as desired.
For a link $L$ with $n$ crossings, we consider all $2^n$ possible versions of $L$ with positive and negative versions of each crossing, and prove that $x+y+1 \mid P(L)-1$ for all of them at once.
For all of them, by looking at any one of the crossings, we have a skein relation $x \cdot P(L_+) + y \cdot P(L_-) + P(L_0) = 0$, where $L_0$ has $n-1$ crossings; therefore $P(L_0) \equiv 1 \pmod{x+y+1}$ by the induction hypothesis. We have $$ x \cdot (P(L_+) - 1) + y \cdot (P(L_-) - 1) + x + y + P(L_0) = 0 $$ and $x + y + P(L_0) \equiv x + y + 1 \equiv 0 \pmod{x+y+1}$, so $$ x \cdot (P(L_+)-1) \equiv -y \cdot (P(L_-)-1) \pmod{x+y+1}. $$ In particular, if we discover that one of $P(L_+) -1$ or $P(L_-)-1$ is divisible by $x+y+1$, we immediately learn that so is the other (because the factors of $x$ and $-y$ can't make a difference).
But at least one of the $2^n$ links can be simplified by Reidemeister moves to a link with fewer crossings, and so for that link, we know that $P(L)-1$ is divisible by $x+y+1$. This propagates to all the links that differ by one crossing, via the argument above, and then propagates again, until we know it for all $2^n$ links we're looking at.
By induction, it's true for all links.