Subdifferentiable functions satisfy Lusin properties of class $C^1$ or $C^2$

Let $f:\mathbb{R}^n\to\mathbb{R}$ be a function. Assume that for a measurable set $\Omega$ and almost every $x\in\Omega$ there exists a vector $\xi_x\in\mathbb{R}^n$ such that $$\liminf_{h\to 0}\frac{f(x+h)-f(x)-\langle \xi_x, h\rangle}{|h|^2}>-\infty.$$ Then we show that $f$ satisfies a Lusin-type property of order $2$ in $\Omega$, that is to say, for every $\varepsilon>0$ there exists a function $g\in C^2(\mathbb{R}^n)$ such that $\mathcal{L}^{n}\left(\{x\in\Omega : f(x)\neq g(x)\}\right)\leq\varepsilon$. In particular every function which has a nonempty proximal subdifferential almost everywhere also has the Lusin property of class $C^2$. We also obtain a similar result (replacing $C^2$ with $C^1$) for the Fr\'echet subdifferential. Finally we provide some examples showing that this kind of results are no longer true for "Taylor subexpansions" of higher order.