A note on what I think to be a misleading trend in some texts, where I have recently stumbled.
As a proper Riemann integral, $$\iiint_{\tau}\nabla\cdot\left(\frac{\mathbf{e}_r}{r^2}\right)d\tau$$ where $\frac{\mathbf{e}_r}{r^2}=\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}$ (it seems that $\mathbf{y}=\mathbf{0}$ in your case), is $0$ if $\mathbf{y}\notin\bar{\tau}$ and does not exist if $\mathbf{y}\in\bar{\tau}$, because the integrand of a Riemann integral, in the usual calculus definitions of it, has to be defined on all the domain.
As the limit $$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)dx_1dx_2dx_3$$ it is $0$ because $\lim_{\varepsilon\to 0}0=0$ and the same holds for the Lebesgue integral $$\int_{\tau}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)d\mu_{\mathbf{x}}$$which is calculated as the preceding limit.
This shows that, under these definitions of the integral and the usual definition of the derivative, the divergence theorem, certainly valid if, instead of $\frac{\mathbf{e}_r}{r^2}$ you had a vector field $\mathbf{F}\in C^1(\mathring{A})$ with $\tau\subset\mathring{A}$ satisfying opportune assumptions, cannot be applied.
Since $\forall\mathbf{x}\in\mathbb{R}^3\setminus\{\mathbf{y}\}\quad\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}=-\nabla\left(\frac{1}{\|\mathbf{x}-\mathbf{y}\|}\right)$ and the divergence of the gradient is the Laplacian $\nabla\cdot\nabla=\nabla^2$, we see that $\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)=-\nabla^2\left(\frac{1}{\|\mathbf{x}-\mathbf{y}\|}\right)$. Then by reading the integral $$\int_{\tau}-\nabla^2\left(\frac{1}{r}\right)\varphi \,d\tau$$ where $\varphi\in C^2(\mathbb{R}^3)$ (typically $\varphi\in C^\infty(\mathbb{R}^3)$) is such that $\forall\mathbf{x}\notin\tau\quad\varphi(\mathbf{x})=0$, in the symbolic way of the Laplacian of the distribution defined by $-\frac{1}{r}$, it can be shown, as it is here, that $-\int_{\tau}\nabla^2\left(\frac{1}{r}\right)\varphi \,d\tau=4\pi\int_{\tau}\delta_{\mathbf{y}}\varphi \,d\tau:=4\pi\varphi(\mathbf{y})$ (where $\delta_{\mathbf{y}}(\mathbf{x}):=\delta(\mathbf{x}-\mathbf{y})$), but that is$$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\frac{-\nabla^2\varphi(\mathbf{x})}{\|\mathbf{x}-\mathbf{y}\|}dx_1dx_2dx_3=\int_{\tau}\frac{-\nabla^2\varphi(\mathbf{x})}{\|\mathbf{x}-\mathbf{y}\|}d\mu_{\mathbf{x}}$$if we use the usual Riemann (one th left) or Lebesgue (on the right) integrals, while $$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)\varphi(\mathbf{x})dx_1dx_2dx_3=\int_{\tau}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)\varphi(\mathbf{x})d\mu_{\mathbf{x}}\equiv 0$$for all $\mathbf{y}$ and all functions $\varphi$.