The level sets of a regular functions $f$ which is a first integral of a vector field (i.e. for a first order ODE) are invariant in the sense that $f$ does not vary along the integral curves.
This means that $\nabla f(x)\cdot X(x) = 0 $ for any $x$ in the space domain. Hence if $f$ is not a constant function, this means that the variation of $f$ happens along directions which live in the orthogonal space to the one generated by the vectors of the vector field $X$ associated to the dynamical system $\dot{x}=f(x)$ you defined in the question.
Hence the level sets $\{x\in\mathbb{R}^n : f(x)=c\}=f^{-1}(\{c\})$ are mapped by the flow of the vector field in a subset of the level set itself. In formulas, let $\Phi^t:\mathbb{R}^n\rightarrow \mathbb{R}^n$ be the flow of the system, then we have that $\forall c\in f(\mathbb{R}^n)$ and $\forall t\in\mathbb{R}$ it follows $\phi^t(f^{-1}(\{c\}))\subset f^{-1}(\{c\})$.
The reason behind this invariancy is simply because the gradient $\nabla f$ is orthogonal to the level sets of $f$, hence since the gradient vector is orthogonal again to the vector field $X$ by definition of first integral, then it follows that the level sets are parallel to the vector field $X$, or more precisely to its integral curves. Hence the level sets are invariant with respect to the flow of $X$.