On Wikipedia (https://en.wikipedia.org/wiki/Mean_curvature) it says that mean curvature $H$ in 3D-space can be calculated via $$2H= - \nabla \cdot n$$ but also via $$2H = tr ((II)(I^{-1}))$$. How do these definitions coincide?
On a unit sphere with parametrization $X(\theta,\phi)=\bigl(\sin(\theta)\cos(\phi), \sin(\theta)\sin(\phi, \cos(\theta)\bigr)^T$ for example, I get $$(II)(I^{-1})=\begin{pmatrix}1 &0 \\ 0 &1\end{pmatrix}$$thus $H=1$.Given the normalized normal vector $N=\begin{pmatrix}\sin(\theta)\cos(\phi)\\ \sin(\theta)\sin(\phi)\\\cos(\theta)\end{pmatrix}$, I am not really sure in how to apply the operator $\nabla$ onto $N$.
You made a sign error here. For $\mathbf{n}=\mathbf{e}_r$, you get $X_{\theta\theta}\cdot\mathbf{n}<0$ in $\mathrm{I\!I}$, similarly other components, so $(\mathrm{I\!I})(\mathrm{I}^{-1})=\begin{pmatrix}-1\\&-1\end{pmatrix}$.
What you want to do is to extend $\mathbf{n}$ from $S$ to a neighbourhood of our surface $S$, then calculate its divergence in $\mathbb{R}^3$.
Specifically, every point $\mathbf{p}$ in a sufficiently small neighbourhood of $S$ has a unique closest point $\mathbf{x}$, and we let $\mathbf{n}(\mathbf{p})=\mathbf{n}(\mathbf{x})$, where $\mathbf{n}(\mathbf{x})$ is our outward normal to surface $S$ at $\mathbf{x}\in S$. So in your example, in a small spherical shell we have $\mathbf{n}=\mathbf{e}_r$. Then $$ 2H=-\nabla\cdot\mathbf{n}=-\left.\frac1{r^2}\frac{\partial}{\partial r}(r^2\mathbf{n}\cdot\mathbf{e}_r)\right\vert_{r=1}=-2. $$