Probabilistic approach
I have to think about the right reference, but things are very much different in the non-local case.
The solution $u = u^\epsilon$ can be written in probabilistic terms as
$$ u(x) = \mathbb E^x e^{-k \tau_D} ,$$
where $k = \lambda / \epsilon$, $\tau_D$ is the first exit time from the open set $D = \Omega$ (when probability enters, I prefer to use $D$ rather than $\Omega$ for the domain) for the corresponding isotropic stable Lévy process $X_t$, and $\mathbb E^x$ is the expectation corresponding to the process started at $x$.
In the local case, the Brownian motion is very unlikely to exit $D$ quickly. Roughly speaking, for $t$ small, the probability that $\tau_D < t$ is comparable to the probability that $X_t \notin D$, which in turn is roughly $$\exp(-\operatorname{dist}(x,\partial\Omega)^2 / (4 t)).$$ Thus, $u(x)$ can be estimated from below by (again roughly)
$$ \exp(-k t - \operatorname{dist}(x,\partial\Omega)^2 / (4 t)) . $$
Optimizing with respect to $t$ leads to a lower bound of the form
$$ \exp(-\sqrt{k} \operatorname{dist}(x,\partial\Omega) k t) , $$
which is off by a factor $2$ in the exponent compared to Varadhan's result.
Now let us consider the non-local case and the jump-type stable process. In this case the probability that $\tau_D < t$ is quite large: it is of the order $$c_{N,s} t \int_{\mathbb R^N \setminus \Omega} |x - z|^{-N-2s} dz + O(t^2) $$ (this is comparable to $t \operatorname{dist}(x,\partial\Omega)^{-2s}$ given some regularity of the complement of $\Omega$). Therefore, we immediately get a lower bound for $u(x)$ of the form
$$ t e^{-k t} c_{N,s} \int_{\mathbb R^N \setminus \Omega} |x - z|^{-N-2s} dz , $$
and optimization leads to
$$ c_{N,s} k^{-1} \int_{\mathbb R^N \setminus \Omega} |x - z|^{-N-2s} dz , $$
which behaves in a completely different way.
Conjecture
So a natural conjecture is that
$$ \lim_{\epsilon \to 0^+} \frac{u^\epsilon(x)}{\epsilon} = \frac{c_{N,s}}{k} \int_{\mathbb R^N \setminus \Omega} |x - z|^{-N-2s} dz . $$
It should not be very difficult to write down a rigorous argument. In particular, I expect that the non-local case is far simpler than the local one.
My guess is that the above result follows relatively easily from what is known and can be found in literature, but, again, I do not have any particular reference in mind at this moment.
Final remark: my guess is that Varadhan's result, as well as the above conjecture, require an exterior cone condition, or something similar.
PDE approach
Edit: Here is a more PDE-oriented approach.
Suppose that $u$ is the solution of
$$ u + \epsilon (-\Delta)^s u = 0 \qquad \text{in } \Omega $$
and $u = 1$ in the complement of $\Omega$. Let $w = u$ in $\Omega$ and $w = 0$ in the complement of $\Omega$. Then, at least formally, we have
$$ (\epsilon^{-1} I + (-\Delta)^s) w(x) = c_{N,s} \int_{\mathbb R^N \setminus \Omega} \frac{1}{|y - x|^{N + 2 s}} \, dy =: g(x) $$
in $\Omega$, and therefore $w$ is the resolvent for $(-\Delta)^s$ in $\Omega$ applied to $g$:
$$ u(x) = w(x) = \int_\Omega \biggl(\int_0^\infty e^{-t/\epsilon} p_t^\Omega(x, y) dt\biggr) g(y) dy ;$$
here $p_t^\Omega(x,y)$ is the heat kernel for $(-\Delta)^s$ in $\Omega$. (This is a very informal derivation of what is essentially the Ikeda–Watanabe formula in probability.)
Now let us rewrite the expression for $u$ a bit. We have
$$ \frac{u(x)}{\epsilon} = \int_\Omega \biggl(\int_0^\infty e^{-s} p_{\epsilon s}^\Omega(x, y) dt\biggr) g(y) dy .$$
Now let us write $u = u^\epsilon$ and consider $\epsilon \to 0^+$. Intuitively it is clear that the expression in brackets converges to the Dirac delta at $x$ as $\epsilon \to 0^+$. Making this rigorous requires some effort, but since the first part was already informal, let me ignore all details. We find that
$$ \lim_{\epsilon \to 0^+} \frac{u^\epsilon(x)}{\epsilon} = \int_\Omega g(y) \delta_x(dy) = g(x) ,$$
which is precisely the conjectured result.
