Via a Bismut-Elworthy-Li formula from [KPP23], we derive uniform gradient estimates for transition semigroups associated with stochastic differential equations driven by a large class of cylindrical Lévy processes which includes the important case of cylindrical $α$-stable processes. As the first application, we formulate a Stein's method for quantitative approximation of the invariant measure of these stochastic differential equations in Wasserstein distance. As the second and main application, we study Euler-Maruyama numerical schemes of stochastic differential equations driven by stable Lévy processes with i.i.d. stable components and obtain a uniform-in-time approximation error in Wasserstein distance. Our approximation error has a linear dependence on the stepsize, which is expected to be tight, as can be seen from an explicit calculation for the case of an Ornstein-Uhlenbeck process.