The argument from the link that OP included can be mimicked: reducing the problem to complex Nullstellensatz with a bit of analysis. Consider a non-singlular point $\mathbf{p}\in Z$ near which the equation $Q(x_1,\dots,x_n)=0$ can be solved for one of the coordinates, say $x_n$. Hence there exists an open neighborhood $U$ of $\mathbf{p}=(p_1,\dots,p_{n-1},p_n)$ in $\Bbb{C}^n$ and a complex analytic function $\gamma$ with
$$
\forall (x_1,\dots,x_n)\in U: Q(x_1,\dots,x_{n-1},x_n)=0 \Leftrightarrow x_n=\gamma(x_1,\dots,x_{n-1}).
$$
The assumption on $P$ now shows that
$P(x_1,\dots,x_{n-1},\gamma(x_1,\dots,x_{n-1}))=0$ whenever all the arguments are real numbers.
But since $\mathbf{p}\in\Bbb{R}^n$ and $Q$ has real coefficients, all coefficients of the Taylor expansion of $\gamma$ at $(p_1,\dots,p_{n-1})$ are real. Therefore, all Taylor coefficients of
$P(x_1,\dots,x_{n-1},\gamma(x_1,\dots,x_{n-1}))$ at $(p_1,\dots,p_{n-1})$ are real too. But the latter function is identically zero when $(x_1,\dots,x_{n-1})$ comes from a small enough neighborhood of $(p_1,\dots,p_{n-1})$ in $\Bbb{R}^{n-1}$, because then $(x_1,\dots,x_{n-1},\gamma(x_1,\dots,x_{n-1}))\in Z$. Therefore all Taylor coefficients of $P(x_1,\dots,x_{n-1},\gamma(x_1,\dots,x_{n-1}))$ at $(p_1,\dots,p_{n-1})$ are zero. In particular, this quantity is zero as $(x_1,\dots,x_{n-1})$ varies in a small enough neighborhood of $(p_1,\dots,p_{n-1})$ in $\Bbb{C}^{n-1}$. We conclude that there is a perhaps smaller open neighborhood $U'\subseteq U$ of $\mathbf{p}$ in $\Bbb{C}^n$ such that $P$ vanishes on the non-empty open subset
$$
U'\cap\{(x_1,\dots,x_n)\in\Bbb{C}^n\mid Q(x_1,\dots,x_n)=0\}
$$
of the complex zero locus
$Z_{\Bbb{C}}(Q)=\{(x_1,\dots,x_n)\in\Bbb{C}^n\mid Q(x_1,\dots,x_n)=0\}$. Its non-singular locus $Z_{\Bbb{C}}(Q)^*$, a complex submanifold of $\Bbb{C}^n$, is dense in $Z_{\Bbb{C}}(Q)$; and is furthermore connected due to the irreducibility of $Q$ (see Chapter 0 of Griffiths & Harris for these standard facts). Now $P$ restricts to a holomorphic function on the connected complex manifold $Z_{\Bbb{C}}(Q)^*$ which is zero on the non-empty open subset $U'\cap Z_{\Bbb{C}}(Q)^*$. This holomorphic function should thus be identically zero. Consequently, $P$ vanishes on $\overline{Z_{\Bbb{C}}(Q)^*}=Z_{\Bbb{C}}(Q)$. The complex Nullstellensatz now implies that a power of $P$ is a multiple of $Q$ in $\Bbb{C}[x_1,\dots,x_n]$. But $P$ and $Q$ belong to $\Bbb{R}[x_1,\dots,x_n]$ and the latter is irreducible. We conclude that $Q\mid P$ in $\Bbb{R}[x_1,\dots,x_n]$.
Remark) The requirement that the real zero locus of $Q$ contains an $(n-1)$-dimensional submanifold of $\Bbb{R}^n$ is necessary for having $Q\mid P$. Without that, one can merely say that $Q\mid\sum_{i=1}^sR_i^2+P^{2m}$ for suitable polynomials $R_1,\dots,R_s\in\Bbb{R}[x_1,\dots,x_n]$ and positive integer $m$ (this due to Positivstellensatz, see Theorem 5.5 here for instance). This does not imply $Q\mid P$ in general. For instance, take $Q$ to be the Motzkin polynomial $x^4y^2 + x^2y^4 + 1 − 3x^2y^2$, a famous example from the literature on Hilbert's 17th problem. Then $Q\mid x^2y^2(x^2+y^2+1)(x^2+y^2−2)^2 + (x^2−y^2)^2$. But $Q(x,y)=x^4y^2 + x^2y^4 + 1 − 3x^2y^2$ does not divide $P(x,y):=x^2-y^2$. Here, $Z=\{(\pm 1,\pm 1)\}$ is zero-dimensional.