Surface integral over the surface of a sphere

I have been given an excercise:

$$\int_V\Delta\frac{1}{r}dV \hspace{1cm} \text{with} \hspace{1cm} r=\sqrt{x^2+y^2+z^2}$$

The Volume V is a sphere, centerd at the coordinate system's origin with radius R=1. The differential surface element in spherical coordinates is given by

$$d\textbf{S} = \sin(\vartheta)\textbf{e}_rd\vartheta d\varphi$$

Use "Gauss's theorem" to calculate the Integral.

Now what I did:

$$I = \int_{V}\Delta\frac{1}{r}dV=\int_V \nabla \cdot \nabla\frac{1}{r}dV = \int_{\delta V}\nabla\frac{1}{r}\cdot d\textbf{S} = \int_{\delta V}\nabla\frac{sin(\vartheta)}{r}\cdot \textbf{e}_rd\vartheta d\varphi$$

Now to my question: How do I have to deal with the nabla operator?

Is it

$$I = \int_{\delta V}(\nabla\frac{\sin(\vartheta)}{r})\cdot \textbf{e}_r d\vartheta d\varphi = \int_{\delta V}\frac{-\sin(\vartheta)}{r^2}d\vartheta d\varphi = \frac{-1}{r^2}\int_{\vartheta = 0}^{\pi}\sin(\vartheta)d\vartheta\int_{\varphi = 0}^{2\pi}d\varphi = -4\pi$$

or

$$I = \int_{\delta V}\nabla(\frac{\sin(\vartheta)}{r} \textbf{e}_r)d\vartheta d\varphi = \int_{\delta V}\left(\left(\nabla\frac{\sin(\vartheta)}{r}\right)\cdot \textbf{e}_r + \frac{\sin(\vartheta)}{r}\nabla\cdot\textbf{e}_r\right)d\vartheta d\varphi = \int_{\delta V}\left(-\frac{\sin(\vartheta)}{r^2}+\frac{2\sin(\vartheta)}{r^2}\right)d\vartheta d\varphi = \int_{\delta V}\frac{\sin(\vartheta)}{r^2}d\vartheta d\varphi = \frac{1}{r^2}\int_{\vartheta = 0}^{\pi}\sin(\vartheta)d\vartheta\int_{\varphi = 0}^{2\pi}d\varphi = 4\pi$$

Which of these is correct and why? I myself lean towards the first approach. Thanks in advance.

I don't know about the operations of vector analysis in spherical coordinates. But we don't need them anyway. I write ${\bf x}$ for the points and $|{\bf x}|=:r$ for short.
The function at stake is $$f({\bf x}):={1\over r}\qquad({\bf x}\ne{\bf 0})\ .$$ One computes $${\partial f\over\partial x_i}=-{1\over r^2}{\partial r\over\partial x_i}=-{1\over r^2}\>{x_i\over r}\ ,\tag{1}$$ and continuing in this way one arrives at $$\Delta f({\bf x})\equiv0\qquad({\bf x}\ne{\bf0})\ .\tag{2}$$ Now let $B$ be a ball of radius $R>0$ with center the origin, and put $B\setminus\{{\bf0 }\}=:\dot B$. You are told to compute the quantity $$Q:=\int_{\dot B} \Delta f({\bf x})\>{\rm d}({\bf x})\ .$$ In the light of $(2)$ one has $Q=0$. It is suggested that we compute $Q$ using Gauss' theorem instead. To this end we take from $(1)$ that $$\nabla f({\bf x})=-{1\over r^2}\>{{\bf x}\over r}\ .$$ Let $S$ be the bounding sphere of $B$, oriented outwards. Then Gauss' theorem gives $$Q=\int_B {\rm div}(\nabla f) \>{\rm d}({\bf x})=\int_S \nabla f\cdot {\rm d}\vec S\ .$$ On $S$ we have ${{\bf x}\over R}={\bf n}$ and therefore $$\nabla f\cdot {\rm d}\vec S=-{1\over R^2}{\bf n}\cdot{\bf n}\>{\rm d}S=-{1\over R^2}\>{\rm d}S\ .$$ In this way we obtain $$Q=-{1\over R^2}\int_S{\rm d}S=-{1\over R^2}\>4\pi R^2=-4\pi\ .$$ Where is the mistake? We did apply Gauss' theorem to a situation where it is not applicable. Gauss theorem applies to a $C^1$ vector field defined on a compact body $B$ with piecewise smooth boundary surface $\partial B$. But the domain $\dot B$ of our example is not of this kind. In fact the origin belongs to the boundary of $\dot B$, and the field $\nabla f$ has a heavy singularity there.
Note that $$\Delta\frac{1}{r}=-4\pi\delta(\vec{r})$$ Hence $$I=-4\pi$$