Circular velocity in axisymmetrici potential If a matter distribution is spherically symmetric, the circular velocity is straightforwardly related to the mass via:

$v_{circ} = \sqrt{GM/r}$

In an axisymmetric potential, on the other hand hand, we have

$v_{circ}(R) = \sqrt{R\frac{d\Phi}{dR}}\neq \sqrt{GM/r}$

with $\Phi$ the gravitational potential.

What I don't quite understand is the following: the potential of the Milky Way is not spherical. How is the measured velocity then reliably related to a mass?
What you do is predict the orbital velocity for a more realistic hypothesised mass distribution. There is an easy analytic solution for a disc. More complicated distributions require a numerical solution.

If you get far enough outside an arbitrary mass distribution, then that solution will asymptotically look similar to $v_{\rm circ} = \sqrt{GM/r}$.

Astronomers get into trouble when explaining the evidence for dark matter from galaxy rotation curves to the layperson. They often use the analogy of the solar system, where the orbital speed goes as $r^{-1/2}$. The trouble is that that a Galaxy is not like the solar system because the mass is distributed and in general cannot be assumed to be concentrated at a point and Newton's shell theorem does not apply.

In practice what you can do is posit a more realistic distribution of matter based on star counts and the gas that can be observed. You can then predict (using a computer model) what the rotation velocity should be at any point if that is all the mass that exists. This model will converge to the simple $r^{-1/2}$ law once you get well outside the bulk of the mass in the model. The model can be compared to the observed velocities of gas, globular clusters, planetary nebulae and other velocity tracers - some of which do indeed exist well outside the bulk of the visible galaxy matter. What you find is that they move much faster than predicted, indicating that considerably more unseen matter must exist inside their galactic orbits. You can then add that (dark) matter in various ways. Non-spherically symmetric and spherically symmetric dark matter profiles (of the same total mass) would produce different orbital velocities at different radii for precisely the reason you identify. By looking at tracers at different radii you can map out the distribution of unseen mass. Once happy with your model you can integrate the density of seen and unseen matter to estimate the total mass.

The bottom line is that the $\sqrt{GM/r}$ is just a "back of the envelope" demonstration of the effect. Real research uses more accurate models or works in a regime far outside the bulk of the visible matter where the "back of the envelope" approximation works.