In this paper we consider the following SDE with distributional drift $b$: $$ {\rm d} X_t=σ(X_t){\rm d} B_t+b(X_t){\rm d} t,\ X_0=x\in{\mathbb R}^d, $$ where $σ$ is a bounded continuous and uniformly non-degenerate $d\times d$-matrix-valued function, $B$ is a $d$-dimensional standard Brownian motion. Let $α\in(0,\frac{1}{2}]$, $p\in(\frac{d}{1-α},\infty)$ and $β\in[α,1]$, $q\in(\frac{d}β,\infty)$. Assume $\|({\mathbb I}-Δ)^{-α/2}b\|_p+\|(-Δ)^{β/2}σ\|_q<\infty$. We show the existence and uniqueness of martingale solutions to the above SDE, and obtain sharp two-sided and gradient estimates of the heat kernel associated to the above SDE. Moreover, we study the ergodicity and global regularity of the invariant measures of the associated semigroup under some dissipative assumptions.