The claim that the $F_\epsilon$ map any Sobolev space to any other has nothing to do with the properties (i)-(iii), it is just a consequence of $F_\epsilon$ being a smoothing operator. In fact, in the literature the property of mapping any Sobolev space to any other is frequently used to define what it means to be a smoothing operator.
Recall that $W^r=(-\Delta +1)^{-r/2}L^2(M)$, where $-\Delta$ is the positive Laplacian on $M$ with respect to some Riemannian metric (the choice of which does not matter because $M$ is compact). Again because $M$ is compact, we can assume for our purposes that the volume density used to define $L^2(M)$ agrees with the Riemannian volume density of the metric we have chosen to define $\Delta$. For a smoothing operator $A$ with smooth kernel $k$ on $M\times M$, a smooth function $f$ on $M$, and $r,r'\in \mathbb R$ one has at $x\in M$
$$
((-\Delta +1)^{r'/2}A(-\Delta +1)^{r/2}f)(x)=(-\Delta_x +1)^{r'/2}\int_{M}k(x,y)((-\Delta_y +1)^{r/2}f)(y)dvol(y)
$$
$$
=\int_{M}((-\Delta_x +1)^{r'/2}(-\Delta_y +1)^{r/2}k)(x,y)f(y)dvol(y)
$$
because $\Delta_y$ is symmetric in $L^2(M,dvol)$ and we can interchange $\Delta_x$ with integration because $M$ is compact.
Thus, $(-\Delta +1)^{r'/2}A(-\Delta +1)^{r/2}$ is an integral operator with smooth kernel $k_{r,r'}$ on $M\times M$, given by
$$
k_{r,r'}(x,y):=((-\Delta_x +1)^{r'/2}(-\Delta_y +1)^{r/2}k)(x,y).
$$
In particular, it defines a bounded operator on $L^2(M)$, which shows that $A$ is bounded from $W^r$ to $W^{r'}$.
A similar argument applies to $B'\circ A\circ B$ for any two differential operators $B$ and $B'$, with the only difference that if $B$ is not symmetric in $L^2$, then the operator that is applied to $k$ is the formal adjoint of $B$.
For the uniform boundedness of $F_\epsilon$ as an operator family on Sobolev spaces, one can use the following characterization of $W^r$ for integer $r\geq 0$:
$$
W^r=\{u\in \mathscr D'(M): Du\in L^2(M)\;\text{for every differential operator of order }\leq r \}.
$$
Here $\mathscr D'(M)$ is the space of distributions on $M$ and the norm on $W^r$ can be defined (up to equivalence) using a partition of unity and the local norms
$$
\Vert u \Vert_r^2:=\sum_{|I|\leq r}\Vert \partial_I u \Vert_{L^2}^2,
$$
where $u\in C^\infty_c(U)$, $U\subset M$ being a chart domain, and $\partial_I:=\partial_{x_1}^{i_1}\cdots \partial_{x_n}^{i_n}$, $I=(i_1,\ldots,i_n)$, $n=\mathrm{dim}\,M$, are the standard partial derivatives in $U$.
Now, given any $r'\in \mathbb R$ with $r'\geq r$ one has on $U$
$$
\partial_I\circ F_\epsilon \circ(-\Delta +1)^{-r'/2}=
[\partial_I, F_\epsilon] \circ(-\Delta +1)^{-r'/2}+F_\epsilon \circ\partial_I\circ (-\Delta +1)^{-r'/2},
$$
where $F_\epsilon$ and $[\partial_I, F_\epsilon]$ are uniformly bounded on $L^2$ by (i) and (ii), and $(-\Delta +1)^{-r'/2}$ as well as $\partial_I\circ (-\Delta +1)^{-r'/2}$ are bounded on $L^2$ because they are pseudodifferential operators of orders $\leq 0$.
Patching together these local estimates using a partition of unity, one obtains that
$F_\epsilon: W^{-r'}\to W^r$ is uniformly bounded for every integer $r\geq 0$ and any real $r'\geq r$. For non-integer $r\geq 0$, one can use Sobolev interpolation to get the same result.
