Dirac's delta in 3 dimensions: proof of $\nabla^2(\|\boldsymbol{x}-\boldsymbol{x}_0\|^{-1})=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_0)$

You are right, the partial derivatives of distributions on higher-dimensional spaces are defined as you surmised,

$$\frac{\partial T}{\partial x_i} \colon \varphi \mapsto -T\biggl[\frac{\partial \varphi}{\partial x_i}\biggr],$$

more generally

$$D^{\alpha} T \colon \varphi \mapsto (-1)^{\lvert\alpha\rvert} T[D^{\alpha}\varphi]$$

for higher derivatives.

From that you obtain $(\nabla^2 T)[\varphi] = T[\nabla^2\varphi]$, and for the locally integrable function $f(x) = \lVert x-x_0\rVert^{-1}$ we therefore have

$$(\nabla^2 T_f)[\varphi] = T[\nabla^2\varphi] = \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3.\tag{1}$$

Now choose $R > 0$ so large that $\operatorname{supp} \varphi \subset B_R(x_0)$, and choose $0 < \varepsilon < R$. Then

$$\int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = \lim_{\varepsilon \searrow 0} \int_{\varepsilon < \lVert x-x_0\rVert < R} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3$$

since $f$ is locally integrable and $\nabla^2\varphi$ is continuous.

Away from $x_0$, $f$ is smooth, and a slightly tedious computation shows that $\nabla^2 f \equiv 0$ on $\mathbb{R}^3\setminus \{x_0\}$, so

$$\int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = \lim_{\varepsilon \searrow 0} \int_{\varepsilon < \lVert x-x_0\rVert < R} f(x)\cdot \nabla^2\varphi(x) - \varphi(x)\cdot \nabla^2 f(x)\,dx_1\,dx_2\,dx_3,$$

and to that integral you can apply Green's formula to obtain

\begin{align} \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 &= \lim_{\varepsilon \searrow 0}\; \Biggl(\int_{\lVert x-x_0\rVert = R} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS\\ &\qquad\qquad + \int_{\lVert x-x_0\rVert = \varepsilon} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS\Biggr)\\ &= \lim_{\varepsilon \searrow 0} \int_{\lVert x-x_0\rVert = \varepsilon} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS, \end{align}

where $dS$ denotes the surface measure of the sphere, $\frac{\partial}{\partial\nu}$ the directional derivative in direction of the outer normal, and the first integral on the right hand side vanishes since $\varphi$ vanishes in a neighbourhood of the outer sphere. The outer normal on the sphere $\lVert x-x_0\rVert = \varepsilon$ is $-\frac{x-x_0}{\lVert x-x_0\rVert}$,

so we are left with

$$\frac{1}{\varepsilon}\int_{\lVert x-x_0\rVert = \varepsilon} f(x)\langle \nabla\varphi(x),x-x_0\rangle - \varphi(x)\langle \nabla f, x-x_0\rangle\,dS.$$

An easy estimate shows

$$\frac{1}{\varepsilon}\int_{\lVert x-x_0\rVert = \varepsilon} f(x)\langle\nabla\varphi(x),x-x_0\rangle\,dS \leqslant \frac{1}{\varepsilon^2}\int_{\lVert x-x_0\rVert = \varepsilon} \lVert\nabla\varphi\rVert\cdot\varepsilon\,dS \leqslant C\frac{4\pi\varepsilon^2}{\varepsilon} \xrightarrow{\varepsilon \searrow 0} 0.$$

For the other part of the integral, computing $\frac{\partial f}{\partial\nu}$ gives

$$\int_{\lVert x-x_0\rVert = \varepsilon} \varphi(x)\cdot \frac{1}{\varepsilon^2}\,dS,$$

and since $\varphi$ is continuous at $x_0$ we have

$$\lim_{\varepsilon\searrow 0} \int_{\lVert x-x_0\rVert = \varepsilon} \varphi(x)\cdot \frac{1}{\varepsilon^2}\,dS = 4\pi\varphi(x_0).$$

Collecting everything without mucking up the signs, we get

$$(\nabla^2 f)[\varphi] = \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = -4\pi \varphi(x_0) = -4\pi\delta_{x_0}[\varphi].$$

That holds for all test functions $\varphi$, hence $\nabla^2 f = -4\pi\delta_{x_0}$.