Let $α(x)$ be a measurable function taking values in $ [α_1,α_2]$ for $0<\A_1\le \A_2<2$, and $κ(x,z)$ be a positive measurable function that is symmetric in $z$ and bounded between two positive constants. Under a uniform Hölder continuous assumptions on $α(x)$ and $x\mapsto κ(x,z)$, we obtain existence, upper and lower bounds, and regularity properties of the heat kernel associated with the following non-local operator of variable order $$ \LL f(x)=\int_{\R^d}\big(f(x+z)-f(x)-\langle\nabla f(x), z\rangle \I_{\{|z|\le 1\}}\big) \frac{κ(x,z)}{|z|^{d+α(x)}}\,dz. $$ In particular, we show that the operator $\LL$ generates a conservative Feller process on $\R^d$ having the strong Feller property, which is usually assumed a priori in the literature to study analytic properties of $\LL$ via probabilistic approaches. Our near-diagonal estimates and lower bound estimates of the heat kernel depend on the local behavior of index function $α(x)$, when $α(x)\equiv \A\in(0,2)$, our results recover some results by Chen and Kumagai (2003) and Chen and Zhang (2016).