The claim holds under the much weaker assumption that the exceptional set satisfies $\mathcal H^{n-1}(E)=0$.
Since the limits of $f$ and $\nabla f$ exist everywhere we get two continuous functions $F\in C(\mathbb R^n)$ and $G\in C(\mathbb R^n)^n$ such that $F$ is differentiable in $\mathbb R^n\setminus E$ with $\nabla F=G$ there. We will show that $F$ is differentiable in each $x\in E$ with Fréchet derivative $G(x)$.
To verify this, we take $\varepsilon>0$ and $\delta>0$ with $\|G(x)-Gy)\|\le \varepsilon/2$ for all $\|y-x\|\le\delta$. The very nice answer to
If $\mathcal{H}^{n-1}(E)=0$ then $\mathbb{R}^n\setminus E$ is connected by @saúl-rm shows that, for $\|y-x\|\le\delta/2$, there is a smooth curve $\gamma:[0,1]\to \mathbb R^n$ of length $\le 2 \|x-y\|$ with $\gamma(0)=x$, $\gamma(1)=y$, and $\gamma(t)\notin E$ for all $t\in (0,1)$. Since the fundamenal theorem of calculus only requires differentiability in the interior of the interval we obtain
$$F(y)-F(x)-\langle G(x),y-x\rangle=\int_0^1\langle \nabla F(\gamma(t))-G(x),\dot\gamma(t)\rangle dt.$$
As the length of $\gamma$ is less than $2\|x-y\|\le\delta$, all $\gamma(t)$ are $\delta$-close to $x$ so that $\|\nabla F(\gamma(t))-G(x)\|=\|G(\gamma(t))-G(x)\|\le \varepsilon/2$. Using Cauchy Schwarz we get
$$\|F(y)-F(x)-\langle G(x),y-x\rangle\|\le \frac{\varepsilon}{2} \int_0^1\|\dot\gamma(t)\|dt \le \varepsilon \|y-x\|$$
which is the desired differentiability of $F$ in $x$.
The case $n=1$ only requires $E$ to be finte, i.e., that $\mathcal H^{n-1}(E)$ is finite. For $n\ge 2$ the assumption $\mathcal H^{n-1}(E)$ finite is certainly not enough for the given proof.