$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
$$
\mbox{In spherical coordinates,}\quad
\delta\pars{\vec{r}}={\delta\pars{r}\delta\pars{\cos\pars{\theta}}\delta\pars{\phi} \over r^{2}}\quad
\mbox{such that}
$$
\begin{align}
\color{#66f}{\large\int_{{\mathbb R}^{3}}\delta\pars{\vec{r}}\,\dd^{3}\vec{r}}
&=\int_{0^{-}}^{\infty}\dd r\,r^{2}\int_{0}^{\pi}\dd\theta\,\sin\pars{\theta}
\int_{0}^{2\pi}\dd\phi\,{\delta\pars{r}\delta\pars{\cos\pars{\theta}}\delta\pars{\phi} \over r^{2}}
\\[3mm]&=\underbrace{\bracks{\int_{0^{-}}^{\infty}\delta\pars{r}\,\dd r}}
_{\ds{=\ 1}}\
\underbrace{\bracks{%
\int_{0}^{\pi}\delta\pars{\cos\pars{\theta}}\sin\pars{\theta}\,\dd\theta}}
_{\ds{=\ 1}}\
\underbrace{\bracks{\int_{0}^{2\pi}\delta\pars{\phi}\,\dd\phi}}_{\ds{=\ 1}}\
\\[3mm]&=\ \color{#66f}{\Large 1}
\end{align}
$$\mbox{Note that}\quad
\int_{{\mathbb R}^{3}}\delta\pars{\vec{r} - \vec{r}_{0}}\,\dd^{3}\vec{r}
=\int_{{\mathbb R}^{3}}\delta\pars{\vec{r}}\,\dd^{3}\vec{r}
$$
Write $\delta_y(x) = \delta(x-y)$. Let $\phi\in \mathcal{D}(\mathbb{R})$ be a test function. Then we have
$$
\langle f\delta_y',\phi\rangle = \langle \delta_y',f\phi\rangle = - \langle \delta_y,(f\phi)'\rangle = -\langle \delta_y,f'\phi + f\phi'\rangle = -f'(y)\phi(y)-f(y)\phi'(y).
$$
Now, clearly $\phi(y)=\langle \delta_y,\phi\rangle$ and $\phi'(y)=-\langle \delta_y',\phi\rangle$, so we obtain
$$
\langle f\delta_y',\phi\rangle = \langle -f'(y)\delta_y+f(y)\delta_y',\phi\rangle,
$$
which means
$$
f\delta_y' = -f'(y)\delta_y + f(y)\delta_y',
$$
or, equivalently,
$$
f(x)\delta'(x-y) = f(y)\delta'(x-y) - f'(y)\delta(x-y).
$$
Note: If you think in terms of integrals, just write
$$
\langle \delta_y,g\rangle = \int_{\mathbb{R}} \delta(x-y)g(x) \; dx
$$
in the above.
Best Answer
The basic reason for this is that $G(\mathbf{r},\mathbf{r_0})$ is locally integrable on the whole $\mathbb{R}^n$, $n\geq 3$, i.e $G\in L^1_{loc}(\mathbb{R}^n)$. To use the terminology of V.S. Vladimirov [1, § 1.6 p. 15], $G$ is a regular generalized function so, by the linearity of and by the theorem on the passage to the limit under the sign of integral, it is represented exactly as $$ \langle G,\varphi \rangle= \int_{\mathbb{R}^n} G(\mathbf{r},\mathbf{r_0})\varphi(\mathbf{r}) d\mathbf{r} \quad \forall \varphi \in C^\infty_0(\mathbb{R}^n) $$ It is therefore a somewhat natural choice to use the structure of this distribution to prove that $$ \langle\nabla^2 G,\varphi \rangle= \langle G,\nabla^2\!\varphi \rangle= \int_{\mathbb{R}^n} G(\mathbf{r},\mathbf{r_0}) \nabla^2\!\varphi(\mathbf{r}) d\mathbf{r}=\varphi(\mathbf{r_0}) \quad \forall \varphi \in C^\infty_0(\mathbb{R}^n) $$
The method of proof of this equality to which user258521 alludes is simply the above formula in disguise. To see this, consider any domain $\Omega\in \mathbb{R}^n$ for which a form of Gauss's theorem holds, and a mollifier $\varphi_\varepsilon \in C^\infty_0(B_\varepsilon)$, i.e. an infinitely smooth function of compact support contained inside the ball $B_\varepsilon$ of radius $\varepsilon>0$ centered in $\mathbf{0}\in \mathbb{R}^n$ such that $\varphi_\varepsilon(\mathbf{r})\to \delta(\mathbf{r})$ in $\mathscr{D}^\prime$ for $\varepsilon\to 0$. Then $$ \varphi_\varepsilon\ast\chi_\Omega (\mathbf{r})= \int_{\mathbb{R}^n} \varphi_\varepsilon(\mathbf{r}-\mathbf{s})\chi_\Omega(\mathbf{s}) d\mathbf{s} \in C^\infty_0(\Omega+B_\varepsilon) \subset C^\infty_0(\mathbb{R}^n) $$ where
By using the above definition for $\langle G,\varphi_\varepsilon\ast\chi_\Omega \rangle$ and letting $\varepsilon\to 0$ $$ \langle G,\varphi_\varepsilon\ast\chi_\Omega \rangle \to \int_\Omega G(\mathbf{r},\mathbf{r_0}) d\mathbf{r}, $$ and consequently applying Gauss's theorem: $$ \langle \nabla^2 G,\varphi_\varepsilon\ast\chi_\Omega \rangle = \langle G,\nabla^2\!\varphi_\varepsilon\ast\chi_\Omega \rangle\to \begin{cases} 1 & \mathbf{r_0}\in\Omega\\ 0 & \mathbf{r_0}\notin\Omega \end{cases} $$
[1] Vladimirov, V. S. (2002), Methods of the theory of generalized functions, Analytical Methods and Special Functions, 6, London–New York: Taylor & Francis, pp. XII+353, ISBN 0-415-27356-0, MR 2012831, Zbl 1078.46029.