We prove gradient estimates for harmonic functions with respect to a $d$-dimensional unimodal pure-jump Levy process under some mild assumptions on the density of its Levy measure. These assumptions allow for a construction of an unimodal Levy process in $\R^{d+2}$ with the same characteristic exponent as the original process. The relationship between the two processes provides a fruitful source of gradient estimates of transition densities. We also construct another process called a difference process which is very useful in the analysis of differential properties of harmonic functions.