We study harmonic functions associated to systems of stochastic differential equations of the form $dX_t^i=A_{i1}(X_{t-})dZ_t^1+\cdots+A_{id}(X_{t-})dZ_t^d$, $i\in\{1,\dots,d\}$, where $Z_t^j$ are independent one-dimensional symmetric stable processes with indices $α_j\in(0,2)$, $j\in\{1,\dots,d\}$. In this article we prove Hölder regularity of bounded harmonic functions with respect to solutions to such systems.