First, the divergence in spherical coordinates, expressed in terms of derivatives, would take the form
$$\nabla \cdot \vec A = \frac{1}{r^2} \frac{\partial}{\partial r} [r^2 A^r] + \ldots$$
where $A^r$ is the radial component of the vector field $\vec A$. In this case, that's $1/r^2$, so we naively get 0 for this contribution.
Of course, this formula cannot be valid at the origin (there's a $1/r^2$ in it), so we resort to an alternative definition of divergence, one that is regular at the origin:
$$\nabla \cdot \vec A |_{\vec 0} = \lim_{\delta r \to 0} \frac{1}{4 \pi (\delta r)^3/3} \int_0^{2\pi} \int_0^\pi \vec A \cdot \hat r \, (\delta r)^2 \, \sin \theta \, d\theta \, d\phi$$
This basically follows from the divergence theorem and is a general alternative definition for divergence that eschews derivatives. It also can be done with curl and gradient instead. The limit process, of course, makes this entirely equivalent to a derivative.
Evaluating this definition for $\vec A = \hat r/r^2$ at $r = \delta r$ (as we're integrating over a sphere of this radius, and then taking the limit) gives
$$\left. \nabla \cdot \left( \frac{\hat r}{r^2} \right) \right|_{\vec 0} = \lim_{\delta r \to 0} \frac{3}{\delta r} $$
which obviously diverges (positively, as negative radii have no meaning).
So this is a vector field whose divergence is zero everywhere except the origin, where its divergence...well, diverges. That all certainly sounds like a delta function.
Typically, one uses the divergence theorem directly to verify the stated condition of the delta function: that its integral over any region containing zero is 1. That is, we do
$$4\pi \int_0^r \nabla \cdot \frac{\hat r}{r^2} r^2 \, dr= \int_0^{2\pi} \int_0^\pi \frac{\hat r}{r^2} \cdot \hat r r^2 \sin \theta \, d\theta \, d\phi$$
If $\nabla \cdot \hat r/r^2 = 4 \pi \delta$, then the surface integral on the left should evaluate to $4\pi$ (which it does, obviously). This is not 100% formally sound, but you can always check this by integrating against a test function (a scalar field would do), just as you would in 1d.
Best Answer
For terms like $\frac{\partial}{\partial x} \int u(x,y,z) \mathrm dz$ you could apply the Leibniz integral rule, see https://en.wikipedia.org/wiki/Leibniz_integral_rule.