WORK IN PROGRESS (This is a community wiki to record progress in this problem).
Easy case. If $\mu$ and $D$ are constants, there is an explicit representation of the solution:
$$u(t, x)=\sqrt{\frac{1}{2\pi Dt}}\int_{-\infty}^\infty \exp\left(\frac{(x-y-\mu t)^2}{Dt}\right)\, u_0(y)\, dy.$$
From this formula one can easily establish all the continuity properties one wants.
Constant $D$ case.
We assume that $D$ is a constant and we ask whether
$$\lim_{D\downarrow 0} u(t, x) = w(t, x),$$
where $\partial_t w = -\partial_x (\mu(t, x)w(t, x))+\frac{D}2\partial^2_x w,$ and $w(0, x)=u(0, x)$. (The convergence is intended in an appropriate sense to be specified).
If $\mu=\mu(x)$ the answer is affirmative: see Hans comment below for a transformation that reduces the problem to the constant $\mu$ case. The abstract framework of semigroup theory also applies; see Pazy, Semigroups of linear operators and applications to PDEs, section 3.4: "The Trotter approximation theorem".
If $\mu=\mu(t)$ the answer is affirmative in the following sense. Taking the spatial Fourier transform termwise, the equation reduces to $\partial_t \hat{u} = \mu(t)(-i\xi \hat{u}(t, \xi))-\frac{D}2\xi^2\hat{u}(t, \xi).$ This is a family of ODEs indexed by $\xi$, and by Gronwall's lemma arguments, as $D\to 0$ its solutions converge uniformly for $t$ in compact intervals. (There is some work to be done here, to show that there is uniformity in $\xi$. This seems true, though).
If $\mu=\mu(t, x)$ is analytic in time, the problem can be reduced to the study of the sequence of PDEs
$$\partial_t u_n=-t^n \partial_x(\mu_n(x)u_n(t, x))+\frac{D}2\partial_x^2 u_n.$$