Here is a proof that adding the condition $\nabla I(\mathbf x) \ne 0$ if $I(\mathbf x) = c$ makes the statement true. In fact, a stronger statement will be proved:
Let $M \subseteq \mathbb R^n$ be a nonempty, $k$-dimensional differentiable manifold, and suppose that $F(\mathbf x)$ is tangent to $M$ at every point $\mathbf x \in M$. Then, $M$ is invariant.
The lemma stated in the question follows from the latter due to the level set $I^{-1}(c)$ being an $(n - 1)$-dimensional differentiable manifold (this follows from the implicit function theorem when $I$ has nonzero gradient).
Now to the proof. Fix $\mathbf x_0 \in M$ and let $\varphi: U \longrightarrow V$, where $U$ is a neighborhood of $0$ in $\mathbb R^n$ and $V$ is a neighborhood of $\mathbf x_0$, be a diffeomorphism such that $\varphi(0) = \mathbf x_0$ and $\varphi(E \cap U) = M \cap V$, where $E = \mathrm{Span}\{e_1, \, \dots, \, e_k\}$. Any trajectory $t \mapsto \mathbf x(t)$ whose support lies entirely inside $V$ can be mapped to $U$ via $\varphi^{-1}$, i.e. we can define $\mathbf y(t) = \varphi^{-1}(\mathbf x(t))$. Then, $\mathbf y$ satisfies $$\dot{\mathbf y} = J_{\varphi^{-1}}(\mathbf x)\dot{\mathbf x} = (J_{\varphi}(\mathbf y))^{-1}F(\varphi(\mathbf y)) =: G(\mathbf y).$$ Now, whenever $\mathbf y \in E$ (which means $\varphi(\mathbf y) \in M$), $F(\varphi(\mathbf y))$ belongs to the tangent space to $M$ at $\mathbf x
= \varphi(\textbf y)$, $T_{\mathbf x}M = J_{\varphi}(\mathbf y)(E)$. Thus $F(\mathbf x) = J_{\varphi}(\mathbf y)\mathbf v$ for some $\mathbf v \in E$, and $G(\mathbf y) = (J_{\varphi}(\mathbf y))^{-1}J_{\varphi}(\mathbf y)\mathbf v = \mathbf v$. If we manage to show that the trajctories of $\dot{\mathbf y} = G(\mathbf y)$ starting at any $\tilde{\mathbf y} \in E$ lie (locally) in $E$ (that is, $E$ is invariant for $G$) we are done, since these trajectories are one-to-one with the trajectories of $\dot{\mathbf x} = F(\mathbf x)$ starting at $\tilde{\mathbf x} \in M$ via $\varphi$, and $\varphi$ maps points of $E$ to points of $M$.
To do this, consider any solution of $\dot{\mathbf y} = G(\mathbf y)$ with initial value $\mathbf y(0) = \tilde{\mathbf y}$. If we let $\tilde{\mathbf y} = (\tilde{\mathbf z}, \, 0)$ (where the $0$ is $(n - k)$-dimensional) and, for $\mathbf z \in \mathbb R^k$, $G((\mathbf z, \, 0)) = (H(\mathbf z), \, 0)$ (recall that $G(E \cap U) \subseteq E$), we can consider the system $\begin{cases} \dot{\mathbf z} = H(\mathbf z) \\ \mathbf z(0) = \tilde{\mathbf z} \end{cases}$ in $\mathbb R^k$. Due to the regularity of $H$, local existence and uniqueness of the solution hold. But now, if $\mathbf z(t)$ is the local solution, $\mathbf y(t) = (\mathbf z(t), \, 0)$ is a solution to $\dot{\mathbf y} = F(\mathbf y)$ that lies entirely in $E$, and it is in fact the only one.
Best Answer
You are right in saying that every orbit is an invariant set. However, in general (but not always) when we talk about invariant sets, we talk about sets of nonzero measure. Examples are easy to find. If you take a 2D nonlinear ODE with a sink type fixed point, it will have a basin of attraction that will form an invariant set.