You don't want to do this by brute force computation. Your intuition is correct: The principal directions are the tangent vector to the profile curve (meridian) and any basis vectors for the tangent space to the sphere.
You can verify this by looking at the covariant derivative of the unit normal in these directions. As you move along a profile curve, you are considering a plane curve and the covariant derivative of the normal must be orthogonal to the normal, hence tangent to the curve. Next, by rotational symmetry, as you move along a latitude sphere, you can see that the normal vector to the hypersurface has two components, one in the direction of the normal to the sphere and one pointing along the first axis; by symmetry, the magnitudes in those directions are constant. Thus, when you compute the covariant derivative in the direction of $V$, tangent to the latitude sphere, the first component gives you a multiple of the $V$ and the second component drops out.
Using the usual Meusnier's Theorem analysis that works for a surface in $\Bbb R^3$ [see, for example, pp. 51-2 of my text] you can figure out the principal curvatures. Up to sign, the principal curvature in the direction of the profile curve is the curvature of the curve. More interestingly, by Meusnier's Theorem, the normal curvature in any direction $V$ tangent to the sphere will be the $1/f(t)$ (the curvature of a great circle with tangent vector $V$) multiplied by $\cos\phi$, where $\phi$ is the angle between the normal to the hypersurface and the principal normal of that great circle. The principal normal points toward the center of the sphere, and so with a drawing you can check that (again up to supplements) $\phi$ is the angle between $(-f'(t),1)$ and $(0,1)$ in $\Bbb R^2$. That is, $\cos\phi = 1/\sqrt{1+f'(t)^2}$.
Finally, then, the principal curvatures (up to sign) of your hypersurface are $\dfrac{f''(t)}{(1+f'(t)^2)^{3/2}}$ and $\dfrac 1{f(t)\sqrt{1+f'(t)^2}}$ (with multiplicity $n-1$).
You cannot show that $S_p(Z)$ is zero, you have to show that $S_p(Z)$ is a multiple of $Z$, i.e. $S_p(Z) = \lambda(p,Z) Z$, where $\lambda$ is a function that depends on the direction $Z$ and the point $p$. Once you have this fact, you have to show that $\lambda$ does not depend on the direction $Z$ and that $\lambda$ is constant on the surface. Then you have proven that the surface is umbilical.
Hint: Consider a point $p$ on the surface. Show that, if the curvature of $\gamma$ is not $0$ at $p$, then the normal vector of $\gamma$ is equal to the surface normal.
This first hint should suffice as a starting point. If you need more input, leave a comment.
Addendum 1: further explanation. Your specific question (the case $\kappa(0)=0$) has already been asked before and a hint can be found in this question.
The hint actually shows that we can "work around" the directions in which the normal curvatures are zero.
A point is parabolic by definition if there is only one direction in which $\kappa_n$ is zero. A point is hyperbolic, by definition, if $K< 0$, or equivalently, if there are two directions in which the normal curvature is zero. So in case of a parabolic or hyperbolic point, we already know that in every direction, except in one or two directions, $S_p(Z) = \lambda(p,Z) Z$ holds.
Now take two linearly independent vectors $v$ and $w$, such that the normal sections in the directions $v$ and $w$ and $v+w$ are not zero. By the argument hereabove, we know that $S_p(v)=\lambda(p,v) v$, $S_p(w) = \lambda(p,w) w$ and $S_p(v+w) = \lambda(p,v+w) v+w$. Then by linearity
$$
\lambda(p,v+w)(v+w) = S_p(v+w) = S_p(v) + S_p(w)
= \lambda(p,v) v + \lambda(p,w) w.
$$
We may conclude that $\lambda(p,v)=\lambda(p)$ does not depend on the directions, except maybe in the one or two directions where the normal curvature is zero.
Consider now a direction $Z$ in which $\kappa_n = 0$. Write $Z=av+bw$ with $a,b\in \mathbb{R}$. Then by linearity
$$
S_p(Z) = a S_p(v) + b S_p(w) = a \lambda(p) v + b \lambda(p) w = \lambda(p) Z.
$$
So we conclude that the normal curvature section $\kappa_n$ must be equal in every direction, so the point is umbilical. Hence, in conclusion, we see that $p$ cannot be a parabolic or hyperbolic point.
Addendum 2: asymptotic directions. Consider a point at a surface. If $e_1$ and $e_2$ are the two principal directions with principal curvatures $\kappa_1$ and $\kappa_2$ respectively, then the normal curvature in the direction of $\cos \theta e_1 + \sin\theta e_2$ is
$$
\kappa_n = \kappa_1 \cos^2 \theta + \kappa_2 \sin^2 \theta.
$$
This is sometimes called Euler's formula. If $K = \kappa_1 \kappa_2$ is negative at the point, then $\kappa_1 > 0$ and $\kappa_2 < 0$. Then there are exactly two solutions $\theta$ such that $\kappa_n = 0$, which correspond to two linearly independent directions with zero normal curvature. These directions are called asymptotic directions. Also note that the point is planar (every normal curvature is zero) if and only if $\kappa_1=\kappa_2=0$.
Best Answer
I only know that operator $df(X) \cdot dN(Y)$ is symmetric because the partial derivatives commute. The formal proof of that is just one line in W.P.A. Klingenberg "A Course in Differential Geometry", p.38. Can be this treated as a geometric intuition?
Informally, I would think of $df(X) \cdot dN(Y)$ as being defined completely just by by $df$ (and the dot product in the ambient space, which is $\mathbb{R}^3$ in the question), so the changes in direction $X$ are reflected automatically (via $N=\frac{df}{\left\Vert df\right\Vert }$) in direction $Y$. Maybe, this observation can be made more rigorous.
Edit. I've just noticed an error in the above. Actually, $N=\pm \frac{ds}{\left\Vert ds\right\Vert }$ for a defining function $s$ of the hypersurface.