This answer is a **heuristic** along the lines of Joro's.

We use $p,q,r$ to denote primes.  Let $S(p,a)$ denote the number of pairs of primes $(q,r)$ with $q,r\leq p$ and $p|(qr+a)$.  We are interested in the case $a=-1$, but in general by the orthogonality relations of the characters we have $$ S(p,a)=\frac{1}{\phi(p)}\sum_{r\leq p}\sum_{q\leq p}\sum_{\chi\ \text{mod}\ p}\overline{\chi}(-1)\chi\left(qr\right).$$   Rearranging this is $$S(p,a)=\frac{1}{p-1}\sum_{\chi\ \text{mod}\ p}\overline{\chi}(a)\left(\sum_{q\leq p}\chi\left(q\right)\right)^{2}.$$

We might hope that this behaves a similar way as when we try to prove the prime number theorem, that is we might hope that the sums are all very small for $\chi$ non principle, and that only the principle character contributes.  With this in mind we expect

$$S(p,a)\approx \frac{1}{p} \text{li}(p)^2.$$

Numerically this is **remarkably** close.  Using the calculation done in Joro'sanswer, letting $a=-1$ and $p=1000003$ we have $$S(p,-1)=6184$$ whereas $$\frac{1}{p} \text{li}(p)^2=6182.307\dots $$

Now all that remains is to understand the sum $$\sum_{q\leq p}\chi\left(q\right)$$ for a character modulo $p$.  However, I believe this is very difficult.