In $\mathbb R^d$, $d \geq 3$, consider the divergence and the non-divergence form operators \begin{equation} \tag{$i$} - \nabla \cdot a \cdot \nabla + b \cdot \nabla, \end{equation} \begin{equation} \tag{$ii$} - a \cdot \nabla^2 + b \cdot \nabla, \end{equation} where $a=I+c \mathsf{f} \otimes \mathsf{f}$, the vector fields $\nabla_i \mathsf{f}$ ($i=1,2,\dots,d$) and $b$ are form-bounded (this includes the sub-critical class $[L^d + L^\infty]^d$ as well as vector fields having critical-order singularities). We characterize quantitative dependence on $c$ and the values of the form-bounds of the $L^q \rightarrow W^{1,qd/(d-2)}$ regularity of the resolvents of the operator realizations of ($i$), ($ii$) in $L^q$, $q \geq 2 \vee ( d-2)$ as (minus) generators of positivity preserving $L^\infty$ contraction $C_0$ semigroups. The latter allows to run an iteration procedure $L^p \rightarrow L^\infty$ that yields associated with ($i$), ($ii$) $L^q$-strong Feller semigroups.