We study gradient estimates of $q$-harmonic functions $u$ of the fractional Schr{ö}dinger operator $Δ^{α/2} + q$, $α\in (0,1]$ in bounded domains $D \subset \R^d$. For nonnegative $u$ we show that if $q$ is H{ö}lder continuous of order $η> 1 - α$ then $\nabla u(x)$ exists for any $x \in D$ and $|\nabla u(x)| \le c u(x)/ (\dist(x,\partial D) \wedge 1)$. The exponent $1 - α$ is critical i.e. when $q$ is only $1 - α$ H{ö}lder continuous $\nabla u(x)$ may not exist. The above gradient estimates are well known for $α\in (1,2]$ under the assumption that $q$ belongs to the Kato class $\calJ^{α- 1}$. The case $α\in (0,1]$ is different. To obtain results for $α\in (0,1]$ we use probabilistic methods. As a corollary, we obtain for $α\in (0,1)$ that a weak solution of $Δ^{α/2}u + q u = 0$ is in fact a strong solution.