The Lagrangian written in the question is a low-energy effective Lagrangian. By definition such a low-energy effective description only sees the massless degrees of freedom. It is obtained from some UV description by integrating out all the massive degrees of freedom. The Coulomb branch is parametrized by the expectation values of the massless vector multiplet scalars. One cannot see the massive BPS particles in the low energy description (except at special points of the moduli space where some of the generically massive BPS particles become massless. These special points correspond generally to singularities of the moduli space: at these points our low energy description breaks down because we have forgotten some massless states).

For the reason explained in the previous paragraph, the vector multiplets appearing in the Lagrangian are massless. They are also BPS. It is not a contradiction because their charge $Z$ is zero. Indeed the fields of the vector multiplets live in the adjoint representation of the gauge group and as this gauge group is abelian generically on the Coulomb branch they are in fact neutral.

As the Lagrangian is a low energy description, it depends on a choice of vacua. As these vacua are parametrized by the expectation values of the scalars $a$, the parameters of the low energy Lagrangian, and in particular the coupling constant $\tau_{IJ}$, depend on $a$. The one loop argument is correct and give the behaviour of $\tau_{IJ}(a)$ in the region $|a|>>\Lambda$ of the moduli space where the theory is weakly coupled.

Let $G$ be the gauge group of the UV theory, a compact connected Lie group. By definition, at a point of the Coulomb branch one can have non-trivial expectation values for the scalars in the vector multiplets. These scalars are valued in the adjoint representation, i.e. the Lie algebra $\mathfrak{g}$ of $G$ with adjoint action of $G$. Classically, the Coulomb branch is simply the quotient of $\mathfrak{g}$ by the adjoint action of $G$, i.e. the quotient of a Cartan subalgebra $\mathfrak{h}$ of $\mathfrak{g}$ by the Weyl group $W$. For example, for $G=U(n)$, we are simply looking to hermitian matrices up to unitary conjugation, $\mathfrak{h}$ can be taken as the set of real diagonal matrices and $W$ is the symmetric group permuting these $n$ diagonal values. At a point $\phi \in \mathfrak{h}/W$ of the classical Coulomb branch, the gauge group $G$ is classically Higgsed to the stabilizer $G_\phi$ of this element under the adjoint action. Let $T=U(1)^r$ be the maximal torus in $G$ whose Lie algebra is $\mathfrak{h}$. It is obvious that $T$ is included in $G_\phi$: the gauge group cannot be Higgsed to something smaller than $T$ on the Coulomb branch. It is true that classically, for some $\phi$, $G_\phi$ can be bigger than $T$ and non-abelian. The claim is that these elements are "special", "non-generic". More precisely, the following is true: look at the roots of $\mathfrak{g}$, there are finitely many linear forms on $\mathfrak{h}$ and so their vanishing locus is a finite union of hyperplanes in $\mathfrak{h}$. The claim is that $G_\phi=T$ for every $\phi$ in the complement of these hyperplanes. The general proof of this claim needs a bit of structure theory of Lie algebras and Lie groups and I skip it (I can add it later if someone would like to see it). The important example to have in mind is $G=U(n)$. If $\phi$ is a diagonal matrix $(a_1,...,a_n)$ and if the $a_i$'s are all distinct, it is an easy computation to check that $G_\phi=T=U(1)^n$. But if for example $a_i=a_j$ for some $i\neq j$ then the $U(2)$ subgroup rotating the ith and jth coordiantes will be in $G_\phi$. The $n(n-1)/2$ hyperplanes defined by the equations $a_i=a_j$, $i\neq j$, are precisely the root hyperplanes in the Lie algebra of $U(n)$.

The previous description is classical. The key point is that by $N=2$ supersymmetry, no potential can be generated and lift the previous moduli space. So the Coulomb branch will stay globally the same and a statement such that "at a generic point the gauge group is Higgsed to its maximal torus" is expected to remain true. But the precise metric on the moduli space is quantum corrected and so the precise description given by the above classical arguments of what happens at a precise point of the moduli space will not be valid in general. For example in the simplest possible case of the pure $G=SU(2)$ $N=2$ gauge theory, the classical description predicts that the Coulomb branch can be described as a complex plane, that at a point $u\neq 0$, the gauge group is broken to its maximal torus $U(1)$, and that at $u=0$, $SU(2)$ is unbroken. The correct quantum mechanical description, famously derived by Seiberg and Witten, is that it is still true that for a generic $u$ the $SU(2)$ symmetry is broken to $U(1)$, as expected, but also that there are in fact two particular singular points in the moduli space where something special happens and there is in fact no point where the $SU(2)$ symmetry is unbroken.