The algebra in equation (155) of [1] (stated for convex functions of real variable) is elementary and applies to Hilbert spaces as well.

Let $\lambda>0$ and $y$ be such that $\lambda+y>0$. Let $p_0=prox[\lambda f](x)$ and $p_y=prox[(\lambda+y)f](x)$. By definition, the objective function at $y+\lambda$ of the minimizer $p_y$ is smaller than the same objective at $p_0$,
\begin{align}
h(y) &= f_{\lambda+y}(x)-f_\lambda(x) + \frac{y}{2\lambda^2}\|x-p_0\|^2
\\&\le \frac{\|x-p_0\|^2}{2(\lambda+y)} + f(p_0) - f(p_0) - \frac{\|x-p_0\|^2}{2\lambda} + \frac{y}{2\lambda^2}\|x-p_0\|^2
\\&= \|x-p_0\|^2 \frac{y^2}{2\lambda^2(y+\lambda)}.
\end{align}
Because $\lambda\mapsto f_\lambda(x)$ is convex in $\lambda$, $y\mapsto h(y)$ is also convex and $h(0)=0$ gives 
$$h(y)\ge - h(-y)\ge - \|x-p_0\|^2 \frac{y^2}{2\lambda^2(-y+\lambda)}.$$
This gives $h(y)=O(y^2)$ hence the desired derivative.

**Reference**

[1] Chistos Thrampoulidis, Ehsan Abbasi, Babak Hassibi, 
"Precise error analysis of regularized M-estimators in high dimensions", IEEE Transactions on Information Theory 64, No. 8, 5592-5628 (2018), [arxiv:1601.06233v1](https://arxiv.org/abs/1601.06233)**[cs.IT]**, DOI:[10.1109/TIT.2018.2840720](https://doi.org/10.1109/TIT.2018.2840720), [MR3832326](https://mathscinet.ams.org/mathscinet-getitem?mr=3832326), [Zbl 1401.94051](https://zbmath.org/1401.94051).