One (but probably not the easiest) way to proceed is this: Just for the heck of it, I'll explain as much as I can in arbitrary dimensions and codimensions. Work with a fixed set of local co-ordinates on the submanifold. Let $f: M^n \rightarrow \mathbb{R}^{N}$. Let $g_{ij}$ be the induced metric given by
$$
\partial_i f \cdot \partial_j f = g_{ij}
$$
Using this it's easy to compute the first variation of $g$. Let $\Gamma^k_{ij}$ be the Christoffel symbols. Verify that
$$
g_{kl}\Gamma^l_{ij} = \partial_kf \cdot \partial^2_{ij}f.
$$
Use this to compute the first variation of $\Gamma^k_{ij}$.
Next, verify that
$$
\partial_{ij}f = \Gamma^k_{ij}\partial_kf + H_{ij},
$$
where $H_{ij}$ is normal to the tangent space of $M$ (as a subspace of $\mathbb{R}^N$. From this you can compute the first variation of $H_{ij}$. The curvature tensor is given by the Gauss equations:
$$
R_{ijkl} = H_{ik}\cdot H_{jl} - H_{il}\cdot H_{jk}.
$$
From this you can compute the first variation of the curvature tensor. If $n = 2$, you can take advantage of the fact that Gauss curvature is equal to a constant times scalar curvature $S = g^{ik}g^{jl}R_{ijkl}$ to calculate its variation.

Now let's assume $N = n+1$. To get the first variation of the second fundamental form, you need a formula for the unit normal. If $n = 2$, then you can use the cross product of $\partial_1f$ and $\partial _2f$ divided by its norm. In higher dimensions, the Gauss map is the unit vector dual to $*df/|*df|_g$, where $*$ is the Hodge star. Using this formula you can get a first variation of the normal vector. The second fundamental form of a hypersurface is now given by
$$
h_{ij} = \nu\cdot H_{ij}.
$$
The mean curvature is easy to get from this. As for the principal curvatures, the first variation of one particular principal curvature is guaranteed to exist only if it is different from all of the others. You can get its variation by implicitly differentiating
$$
0 = \det (H_{ij} - \lambda g_{ij}).
$$