Gauge theories become constrained Hamiltonian systems when passing from the Lagrangian $L(q,\dot{q},t)$ to the Hamiltonian $H(q,p,t)$ where $p = \frac{\partial L}{\partial \dot{q}}$. Generically, you get a constrained Hamiltonian system whenever the matrix/operator with components
$$ \frac{\partial^2 L}{\partial \dot{q}^i\partial\dot{q}^j}$$
is singular, i.e. non-invertible, i.e $\det(\frac{\partial^2 L}{\partial \dot{q}\partial\dot{q}}) = 0$. As you correctly observe, this is already the case for a massive vector field.
So, let us look at the Lagrangian of a generic vector field:
$$ L(A,\dot{A}) = \int (\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \mu^2 A^\mu A_\mu) \mathrm{d}^3x$$
Regardless of $\mu$, the canonical momenta are
$$ \pi^\mu = F^{\mu 0}$$
so we have the primary constraint1
$$\pi^0 \approx 0\tag{1}$$
always. The canonical Hamiltonian reads
$$ H = \int (\frac{1}{2}\pi^i\pi_i - \frac{1}{4} F^{ij}F_{ij} - A_0 \partial_i \pi^i -\mu^2 A^\mu A_\mu)\mathrm{d}^3 x$$
and, for consistency of the constraint, we incur a secondary constraint
$$ \dot{\pi}^0 = \{\pi^0,H\} = \partial_i\pi^i + \mu^2 A_0 \approx 0 \tag{2}$$
using that the only non-zero Poisson bracket of $\pi_0$ is $\{A_0,\pi^0\} = 1$. The nature of the theory is now very different depending on whether or not $\mu^2 = 0$.
If $\mu^2\neq 0$, then $(2)$ does not, effectively, impose a restriction on $\partial_i \pi^i$. Solving the two constraints means just putting $A_0 = -\frac{1}{\mu^2}\partial_i\pi^i$ and $\pi^0 = 0$, meaning we reduce the phase space dimension by two (effectively forgetting there was a pair of coordinates $(A_0, \pi^0)$) and are in a constraint-free theory. On this reduced phase space, canonical quantization may now proceed as usual, and we do not have gauge degrees of freedom left. In particular, canonical quantization delivers a propagator, since the operator is invertible on the degrees of freedom that are left.
If $\mu^2 = 0$, then $(2)$ is just Gauss' law $\nabla \cdot \vec E = 0$ since $\pi^i = F^{i0} = E^i$. Although it is in principle possible to locally solve this constraint and again pass to a reduced phase space (with lower dimension than before), this is not actually feasible or desirable in practice. Hence, canonical quantization as for unconstrained systems is not possible, and we do not obtain a propagator if we try to calculate one, since there are gauge degrees of freedom left.
1$\approx$ denotes weak equalities which hold on the constraint surface, but are not identically zero throughout the entire phase space.