Consider the following figure
where $R=\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}=|\mathbf{r}-\mathbf{r}'|$ is the module of the $\mathbf{R}$ vector depends not only on the location of the $P$ point but also on the location $P'$ where the $dV'$ volume is located (fixed once located in the volume $\mathcal{V}$). Obviously if you change $P'=(x',y',z')$ it will also change $\mathbf{R}$. Since the potential \begin{equation} \psi(\mathbf{r})=-G\iiint_{\mathcal{V}} \frac{\rho(x',y',z')dx'dy'dz'}{|\mathbf{r}-\mathbf{r}'|} \end{equation}
we calculate the gradient of the quantity
\begin{equation} \boldsymbol{\nabla}_{(\mathbf r)}\frac{1}{|\mathbf r-\mathbf r'|} \end{equation}
Calculating, respectively, the partial derivatives $\partial_x=\partial/\partial x$, $\partial_y=\partial/\partial y$ and $\partial_z=\partial/\partial z$ compared to the function $1/|\mathbf{r}-\mathbf{r}'|$, we will have
\begin{align*} \frac{\partial}{\partial x}\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}} & = \frac{\partial}{\partial x}\left((x-x')^2+(y-y')^2+(z-z')^2\right)^{-\frac 12}= && \\ &=\left(-\frac 12\right)\Bigl[\ldots\ldots\Bigr]^{-\frac32}\cdot 2\cdot (x-x')= && \\ &=-\frac{x-x'}{R^3} && \\ \end{align*}
Similarly $$\frac{\partial}{\partial y}\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}=-\frac{y-y'}{R^3}$$ $$\frac{\partial}{\partial z}\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}=-\frac{z-z'}{R^3}$$
Hence
$$ \boldsymbol{\nabla}_{(\mathbf r)}\frac{1}{|\mathbf r-\mathbf r'|}=-\frac{\mathbf r-\mathbf r'}{|\mathbf r-\mathbf r'|^3}. $$ from which
\begin{align} \boldsymbol{\nabla}_{(\mathbf{r})}\psi(\mathbf{r}) & = {\mathbf \nabla}_{(\mathbf{r})}\left(-G\iiint_{\mathcal{V}} \frac{\rho(x',y',z')\,dx'dy'dz'}{|\mathbf{r}-\mathbf{r}'|}\right) = && \tag{*}\\ &= -G\iiint_{\mathcal{V}}\left( {\mathbf \nabla}_{(\mathbf{r})}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\right)\rho(x',y',z')\,dx'dy'dz'= && \\ &= -G\iiint_{\mathcal{V}} \left(-\frac{\mathbf{r}-\mathbf{r}'}{|\mathbf{r}-\mathbf{r}'|^3}\right)\rho(x',y',z')\,dx'dy'dz' = && \nonumber\\ &= G\iiint_{\mathcal{V}} \frac{\rho(x',y',z') \, dx'dy'dz'}{R^2}\,\mathbf{\widehat R}=\mathbf{g}(\mathbf{r}) && \nonumber\\ \nonumber \end{align}
I ask with such kindness, considering that I have to prove that $\mathbf{g}(\mathbf{r})=-\boldsymbol{\nabla}\psi(\mathbf r)$ I was not able to find the error of a less missing sign in the various steps of the $(^*)$.
I hope you will appreciate my effort and my question is clear.