By dyadic decomposition, it suffices to obtain bounds on the quantities $$ \int_{|x-a_i| \sim R_i \forall i} \left[ \frac{V(x,a_1,\dots,a_n)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}}\ dx$$ which, when summed over dyadic powers of two $R_1,\dots,R_n > 0$, becomes bounded by $O(1)$, uniformly over $a_1,\dots,a_n$.
By permuting we may assume that $R_1 \leq \dots \leq R_n$. The most important geometric feature of the set $\{a_1,\dots,a_n\}$ turns out to be its diameter; by rescaling we can normalize this diameter to be $1$. The triangle inequality then lets us assume $R_n \gtrsim 1$. In fact the triangle inequality leaves us with just two scenarios: the "far field" scenario $$ R_1 \sim \dots \sim R_n \gtrsim 1 \quad (1)$$ and the "near field" scenario $$ R_1 \leq \dots \leq R_n \sim 1 \quad (2).$$
One could refine the range of possibilities of the near field further by inspecting the finer metric geometry of the $a_1,\dots,a_n$ beyond just using the diameter normalization, but fortunately we do not have to do so for this problem.
As implicitly observed by the OP, the volume $V(x,a_1,\dots,a_n)$ is of size at most $O(R_1 \cdots R_n)$ by multiplying all the lengths emenating from $x$. On the other hand, since $a_1,\dots,a_n$ lie in a diameter one set and $x$ is at distance $O(R_1)$ from $a_1$, we also have the bound $O(R_1)$ by the base-times-height formula. [Note that this is a significantly superior bound in the far field case, which was identified by the OP as the case where the previous estimates were poor.] Finally, the condition $|x-a_1| \sim R_1$ restricts $x$ to a set of measure $O(R_1^n)$. Thus the above integral can be bounded by $$ \lesssim \left[\frac{\min(R_1 \cdots R_n, R_1)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n.$$
In the far field case (1), this simplifies to $O(R_n^{-2})$ which is summable to $O(1)$. [NOTE: in a previous version I instead obtained the estimate $O(\frac{1}{R_1^{\frac{n-3}{n-1}}R_n^{\frac{n+1}{n-1}}})$, which was easily controlled for $n \geq 3$ but (as pointed out in comments) required a different argument to sum for $n=2$ if one did not exploit the additional property $R_1 \sim R_2$ of the far field.] In the near field case (2), we bound the above by $$ \left[\frac{(R_1 \cdots R_n)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n = \prod_{j=1}^n \frac{R_1}{R_j} \leq \frac{R_1}{R_n} \sim R_1$$ which is summable to $O(1)$ since $R_1 \lesssim 1$ (and there are only $O(\log^{n-2} (1/R_1))$ choices of the remaining scales $R_2,\dots,R_n$ for any fixed $R_1$).
EDIT: one can refactor the above proof to avoid dyadic decomposition as follows.
We again normalize $\mathrm{diam}(a_1,\dots,a_n)=1$, and by translation we also normalize $a_1=0$. We now split into the far field $|x| \geq 2$ and the near field $|x| < 2$. In the far field we have $V(x,a_1,\dots,a_n) \lesssim |x|$ by the base-times-height formula and $|x-a_i| \sim |x|$ by the triangle inequality, so the contribution of the far field is $$ \lesssim \int_{|x| \geq 2} \left[ \frac{|x|^2}{|x|^{n(n+1)}}\right]^{\frac{1}{n-1}}\ dx = \int_{|x| \geq 2} \frac{1}{|x|^{n+2}}\ dx = O(1).$$ In the near field we use the Hadamard bound $V(x,a_1,\dots,a_n) \leq \frac{1}{n!} \prod_{i=1}^n |x-a_i|$ as in the OP to bound things by $$ \lesssim \int_{|x| \leq 2} \frac{dx}{\prod_{i=1}^n |x-a_i|}.$$ By the triangle inequality, at least one of the $|x-a_i|$ is comparable to one. Deleting this factor and then applying the AM-GM inequality, one can bound this by $$ \lesssim \int_{|x| \leq 2} \sum_{i=1}^n \frac{dx}{|x-a_i|^{n-1}} = O(1)$$ and the claim follows. (Note how it was important to delete one of the factors in the denominator before applying AM-GM in order to obtain an integrable exponent of $n-1$ instead of the divergent exponent $n$. One can also use Holder's inequality in place of AM-GM if desired.)
This calculation and its refactoring illustrates how dyadic decomposition can be a systematic tool to estimate (up to constants) various unsigned geometric integrals (with every geometric fact about the integrand and the domain of integration reducible to some inequality regarding the various scales and amplitudes involved, which can then be combined and optimized by elementary algebra), but once one gains enough understanding of the dyadic decomposition argument, one can often streamline the proof to erase traces of the dyadic decomposition used, though sometimes at the cost of providing motivation as to where the proof came from.